Many changes to the original integrator

1 files changed, 77 insertions, 25 deletions

diff --git a/fourier_integrator.cpp b/fourier_integrator.cpp
index 4f4eff5..9f8db97 100644
--- a/fourier_integrator.cpp
+++ b/fourier_integrator.cpp
@@ -8,9 +8,9 @@ int main(int argc, char* argv[]) {
   unsigned s = 2;
   Real λ = 0.5;
   Real τ₀ = 0;
-  Real y₀ = 0;
-  Real yₘₐₓ = 0.5;
-  Real Δy = 0.05;
+  Real β₀ = 0;
+  Real βₘₐₓ = 0.5;
+  Real Δβ = 0.05;
 
   unsigned log2n = 8;
   Real τₘₐₓ = 20;
@@ -41,13 +41,13 @@ int main(int argc, char* argv[]) {
       τ₀ = atof(optarg);
       break;
     case '0':
-      y₀ = atof(optarg);
+      β₀ = atof(optarg);
       break;
     case 'y':
-      yₘₐₓ = atof(optarg);
+      βₘₐₓ = atof(optarg);
       break;
     case 'd':
-      Δy = atof(optarg);
+      Δβ = atof(optarg);
       break;
     case 'I':
       maxIterations = (unsigned)atof(optarg);
@@ -68,8 +68,11 @@ int main(int argc, char* argv[]) {
   Real Δτ = (1 + τ₀ / 2) * τₘₐₓ / M_PI / n;
   Real Δω = M_PI / ((1 + τ₀ / 2) * τₘₐₓ);
 
-  Real z = 0.5;
-  Real Γ₀ = 1 + τ₀ / 2;
+  Real Γ₀ = 1.0;
+  Real μ = Γ₀;
+  if (τ₀ > 0) {
+    μ = (sqrt(1+4*Γ₀*τ₀)-1)/(2*τ₀);
+  }
 
   std::vector<Real> C(2 * n);
   std::vector<Real> R(2 * n);
@@ -78,42 +81,52 @@ int main(int argc, char* argv[]) {
   std::vector<Complex> Ct;
   std::vector<Complex> Rt;
 
-  Real y = y₀;
+  Real β = β₀;
 
   if (!loadData) {
     // start from the exact solution for τ₀ = 0
     for (unsigned i = 0; i < n; i++) {
       Real τ = i * Δτ * M_PI;
       if (τ₀ > 0) {
-        C[i] = Γ₀ / 2 * (exp(-z * τ) - z * τ₀ * exp(-τ / τ₀)) / (z - pow(z, 3) * pow(τ₀, 2));
+        C[i] = Γ₀ * (exp(-μ * τ) - μ * τ₀ * exp(-τ / τ₀)) / (μ - pow(μ, 3) * pow(τ₀, 2));
       } else {
-        C[i] = Γ₀ / 2 * exp(-z * τ) / z;
+        C[i] = Γ₀ * exp(-μ * τ) / μ;
       }
       if (i > 0) {
         C[2 * n - i] = C[i];
       }
-      R[i] = exp(-z * τ);
+      R[i] = exp(-μ * τ);
     }
     Ct = fft.fourier(C);
     Rt = fft.fourier(R);
   } else {
-    std::ifstream cfile(fourierFile("C", p, s, λ, τ₀, y, log2n, τₘₐₓ), std::ios::binary);
+    std::ifstream cfile(fourierFile("C", p, s, λ, τ₀, β, log2n, τₘₐₓ), std::ios::binary);
     cfile.read((char*)(C.data()), (C.size() / 2) * sizeof(Real));
     cfile.close();
     for (unsigned i = 1; i < n; i++) {
       C[2 * n - i] = C[i];
     }
-    std::ifstream rfile(fourierFile("R", p, s, λ, τ₀, y, log2n, τₘₐₓ), std::ios::binary);
+    std::ifstream rfile(fourierFile("R", p, s, λ, τ₀, β, log2n, τₘₐₓ), std::ios::binary);
     rfile.read((char*)(R.data()), (R.size() / 2) * sizeof(Real));
     rfile.close();
 
     Ct = fft.fourier(C);
     Rt = fft.fourier(R);
 
-    z = estimateZ(fft, C, Ct, R, Rt, p, s, λ, τ₀, y);
+    μ = estimateZ(fft, C, Ct, R, Rt, p, s, λ, τ₀, β);
   }
 
-  while (y += Δy, y <= yₘₐₓ) {
+  std::vector<Real> Cb = C;
+  std::vector<Real> Rb = R;
+  std::vector<Complex> Ctb = Ct;
+  std::vector<Complex> Rtb = Rt;
+
+  Real fac = 1;
+  while (β += Δβ, β <= βₘₐₓ) {
+    Real Δμ = 1e-2;
+    Real μ₁ = 0;
+    Real μ₂ = 0;
+    while (true) {
     Real ΔC = 1;
     Real ΔCprev = 1000;
     unsigned it = 0;
@@ -123,12 +136,12 @@ int main(int argc, char* argv[]) {
 
       for (unsigned i = 0; i < Rt.size(); i++) {
         Real ω = i * Δω;
-        Rt[i] = (1.0 + pow(y, 2) * RddfCt[i] * Rt[i]) / (z + 1i * ω);
+        Rt[i] = (1.0 + pow(β, 2) * RddfCt[i] * Rt[i]) / (μ + 1i * ω);
       }
 
       for (unsigned i = 0; i < Ct.size(); i++) {
         Real ω = i * Δω;
-        Ct[i] = (Γ₀ * std::conj(Rt[i]) / (1 + pow(τ₀ * ω, 2)) + pow(y, 2) * (RddfCt[i] * Ct[i] + dfCt[i] * std::conj(Rt[i]))) / (z + 1i * ω);
+        Ct[i] = (2 * Γ₀ * std::conj(Rt[i]) / (1 + pow(τ₀ * ω, 2)) + pow(β, 2) * (RddfCt[i] * Ct[i] + dfCt[i] * std::conj(Rt[i]))) / (μ + 1i * ω);
       }
 
       std::vector<Real> Cnew = fft.inverse(Ct);
@@ -151,8 +164,7 @@ int main(int argc, char* argv[]) {
         R[i] += γ * (Rnew[i] - R[i]);
       }
 
-      z *= Cnew[0];
-
+      /*
       if (it % maxIterations == 0) {
         if (ΔC < ΔCprev) {
           γ = std::min(1.1 * γ, 1.0);
@@ -162,20 +174,60 @@ int main(int argc, char* argv[]) {
 
         ΔCprev = ΔC;
       }
+      */
 
-      std::cerr << it << " " << p << " " << s << " " << τ₀ << " " << y << " " << sqrt(2 * ΔC / C.size()) << " " << γ << std::endl;
+      std::cerr << μ << " " << p << " " << s << " " << τ₀ << " " << β << " " << sqrt(2 * ΔC / C.size()) << " " << γ << " " << C[0];
+      std::cerr << "\r";
 
     }
+      if (std::isnan(C[0])) {
+        C = Cb;
+        R = Rb;
+        Ct = Ctb;
+        Rt = Rtb;
+        μ /= sqrt(sqrt(fac*std::tanh(Cb[0]-1)+1));
+        fac /= 2;
+        μ₁ = 0;
+        μ₂ = 0;
+      } else {
+        Cb = C;
+        Rb = R;
+        Ctb = Ct;
+        Rtb = Rt;
+      if (pow(C[0] - 1, 2) < ε) {
+        break;
+      }
+      if (μ₁ == 0 || μ₂ == 0) {
+        if (C[0] > 1 && μ₁ == 0) {
+          /* We found a lower bound */
+          μ₁ = μ;
+        }
+        if (C[0] < 1 && μ₂ == 0) {
+          /* We found an upper bound */
+          μ₂ = μ;
+        }
+        μ *= sqrt(sqrt(fac*std::tanh(C[0]-1)+1));
+      } else {
+        /* Once the bounds are set, we can use bisection */
+        if (C[0] > 1) {
+          μ₁ = μ;
+        } else {
+          μ₂ = μ;
+        }
+        μ = (μ₁ + μ₂) / 2;
+      }
+      }
+    }
 
-    Real e = energy(C, R, p, s, λ, y, Δτ);
+    Real e = energy(C, R, p, s, λ, β, Δτ);
 
-    std::cerr << "y " << y << " " << e << " " << z << std::endl;
+    std::cerr << "y " << β << " " << e << " " << μ << std::endl;
 
-    std::ofstream outfile(fourierFile("C", p, s, λ, τ₀, y, log2n, τₘₐₓ), std::ios::out | std::ios::binary);
+    std::ofstream outfile(fourierFile("C", p, s, λ, τ₀, β, log2n, τₘₐₓ), std::ios::out | std::ios::binary);
     outfile.write((const char*)(C.data()), (C.size() / 2) * sizeof(Real));
     outfile.close();
 
-    std::ofstream outfileR(fourierFile("R", p, s, λ, τ₀, y, log2n, τₘₐₓ), std::ios::out | std::ios::binary);
+    std::ofstream outfileR(fourierFile("R", p, s, λ, τ₀, β, log2n, τₘₐₓ), std::ios::out | std::ios::binary);
     outfileR.write((const char*)(R.data()), (R.size() / 2) * sizeof(Real));
     outfileR.close();
   }