Testing alternate iterate

1 files changed, 7 insertions, 28 deletions

diff --git a/log-fourier_integrator.cpp b/log-fourier_integrator.cpp
index b21f275..879580d 100644
--- a/log-fourier_integrator.cpp
+++ b/log-fourier_integrator.cpp
@@ -16,7 +16,7 @@ int main(int argc, char* argv[]) {
   Real k = 0.1;
 
   /* Iteration parameters */
-  Real ε = 1e-13;
+  Real ε = 1e-14;
   Real γ = 1;
   Real βₘₐₓ = 0.7;
   Real Δβ = 0.01;
@@ -63,11 +63,6 @@ int main(int argc, char* argv[]) {
 
   Real Γ₀ = 1.0 + τ₀;
   Real μ = 1.0;
-/*  Real μ = Γ₀;
-  if (τ₀ > 0) {
-    μ = (sqrt(1+4*Γ₀*τ₀)-1)/(2*τ₀);
-  }
-*/
 
   std::vector<Real> Cₜ₋₁(N);
   std::vector<Real> Rₜ₋₁(N);
@@ -92,12 +87,9 @@ int main(int argc, char* argv[]) {
   std::vector<Complex> Ĉₜ = Ĉₜ₋₁;
   std::vector<Complex> Ȓₜ = Ȓₜ₋₁;
 
-  Real fac = 1;
   Real β = 0;
   while (β < βₘₐₓ) {
     Real ΔC = 100;
-    Real ΔC₀ = 100;
-    unsigned it = 0;
     while (ΔC > ε) {
       std::vector<Real> dfC(N);
       std::vector<Real> RddfC(N);
@@ -108,22 +100,25 @@ int main(int argc, char* argv[]) {
       std::vector<Complex> RddfCt = fft.fourier(RddfC, false);
       std::vector<Complex> dfCt = fft.fourier(dfC, true);
 
-      std::vector<Complex> Ȓₜ₊₁(N);
       std::vector<Complex> Ĉₜ₊₁(N);
-
+      std::vector<Complex> Ȓₜ₊₁(N);
       for (unsigned n = 0; n < N; n++) {
         Ȓₜ₊₁[n] = (1.0 + pow(β, 2) * RddfCt[n] * Ȓₜ[n]) / (μ + 1i * fft.ν(n));
+        Ĉₜ₊₁[n] = (2 * Γ₀ * std::conj(Ȓₜ[n]) / (1 + pow(τ₀ * fft.ν(n), 2)) + pow(β, 2) * (RddfCt[n] * Ĉₜ[n] + dfCt[n] * std::conj(Ȓₜ[n]))) / (μ + 1i * fft.ν(n));
       }
 
       std::vector<Real> Rₜ₊₁ = fft.inverse(Ȓₜ₊₁);
 
+      /*
       for (unsigned n = 0; n < N; n++) {
         RddfC[n] = Rₜ₊₁[n] * ddf(λ, p, s, Cₜ[n]);
       }
       RddfCt = fft.fourier(RddfC, false);
+
       for (unsigned n = 0; n < N; n++) {
         Ĉₜ₊₁[n] = (2 * Γ₀ * std::conj(Ȓₜ₊₁[n]) / (1 + pow(τ₀ * fft.ν(n), 2)) + pow(β, 2) * (RddfCt[n] * Ĉₜ[n] + dfCt[n] * std::conj(Ȓₜ₊₁[n]))) / (μ + 1i * fft.ν(n));
       }
+      */
       std::vector<Real> Cₜ₊₁ = fft.inverse(Ĉₜ₊₁);
 
       μ *= pow(tanh(Cₜ₊₁[0]-1)+1, 0.05);
@@ -142,22 +137,6 @@ int main(int argc, char* argv[]) {
         Ȓₜ[i] += γ * (Ȓₜ₊₁[i] - Ȓₜ[i]);
       }
 
-      /*
-      if (ΔC < ΔC₀) {
-        ΔC₀ = ΔC;
-        it = 0;
-        γ = std::min(1.001 * γ, 1.0);
-      } else {
-        it++;
-      }
-
-      if (it > 100) {
-        γ = std::max(0.5 * γ, 1e-3);
-        it = 0;
-        ΔC₀ = ΔC;
-      }
-      */
-
       std::cerr << "\x1b[2K" << "\r";
       std::cerr << β << " " << μ << " " << ΔC << " " << γ << " " << Cₜ[0];
     }
@@ -178,7 +157,7 @@ int main(int argc, char* argv[]) {
     }
     E *= β;
 
-      std::cerr << "\x1b[2K" << "\r";
+    std::cerr << "\x1b[2K" << "\r";
     std::cerr << β << " " << μ << " " << Ĉₜ[0].real() << " " << E << " " << γ << std::endl;
     β += Δβ;
     Cₜ₋₁ = Cₜ;