More work on the integrator.

1 files changed, 50 insertions, 41 deletions

diff --git a/fourier_integrator.cpp b/fourier_integrator.cpp
index 40a958e..2ee304d 100644
--- a/fourier_integrator.cpp
+++ b/fourier_integrator.cpp
@@ -65,16 +65,18 @@ public:
 int main(int argc, char* argv[]) {
   unsigned p = 2;
   Real τ₀ = 0;
-  Real y = 0.5;
+  Real yₘₐₓ = 0.5;
+  Real Δy = 0.05;
 
   unsigned log2n = 8;
   Real τₘₐₓ = 20 / M_PI;
 
-  unsigned iterations = 10;
+  unsigned maxIterations = 1000;
+  Real ε = 1e-14;
 
   int opt;
 
-  while ((opt = getopt(argc, argv, "p:2:T:t:y:I:")) != -1) {
+  while ((opt = getopt(argc, argv, "p:2:T:t:y:d:I:")) != -1) {
     switch (opt) {
     case 'p':
       p = atoi(optarg);
@@ -89,10 +91,13 @@ int main(int argc, char* argv[]) {
       τ₀ = atof(optarg);
       break;
     case 'y':
-      y = atof(optarg);
+      yₘₐₓ = atof(optarg);
+      break;
+    case 'd':
+      Δy = atof(optarg);
       break;
     case 'I':
-      iterations = (unsigned)atof(optarg);
+      maxIterations = (unsigned)atof(optarg);
       break;
     default:
       exit(1);
@@ -124,51 +129,48 @@ int main(int argc, char* argv[]) {
   std::vector<Complex> Ct = fft.fourier(C);
   std::vector<Complex> Rt = fft.fourier(R);
 
-  for (unsigned it = 0; it < iterations; it++) {
-    std::vector<Real> RddfC(C.size());
-    for (unsigned i = 0; i < C.size(); i++) {
-      RddfC[i] = R[i] * ddf(p, C[i]);
-    }
-    std::vector<Complex> RddfCt = fft.fourier(RddfC);
-
-    std::vector<Real> dfC(C.size());
-    for (unsigned i = 0; i < C.size(); i++) {
-      dfC[i] = df(p, C[i]);
-    }
-    std::vector<Complex> dfCt = fft.fourier(dfC);
+  Real y = 0;
 
-    std::vector<Complex> Ctnew(Ct.size());
-    std::vector<Complex> Rtnew(Rt.size());
+  while (y += Δy, y <= yₘₐₓ) {
+    for (unsigned it = 0; it < maxIterations; it++) {
+      std::vector<Real> RddfC(C.size());
+      for (unsigned i = 0; i < C.size(); i++) {
+        RddfC[i] = R[i] * ddf(p, C[i]);
+      }
+      std::vector<Complex> RddfCt = fft.fourier(RddfC);
 
-    for (unsigned i = 0; i < Rt.size(); i++) {
-      Real ω = i * Δω;
-      Rtnew[i] = (1.0 + pow(y, 2) * RddfCt[i] * Rt[i]) / (z + 1i * ω);
-    }
+      std::vector<Real> dfC(C.size());
+      for (unsigned i = 0; i < C.size(); i++) {
+        dfC[i] = df(p, C[i]);
+      }
+      std::vector<Complex> dfCt = fft.fourier(dfC);
 
-    Rt = Rtnew;
+      std::vector<Complex> Ctnew(Ct.size());
+      std::vector<Complex> Rtnew(Rt.size());
 
-    for (unsigned i = 0; i < Rt.size(); i++) {
-      Real ω = i * Δω;
-      Ctnew[i] = (Γ₀ * std::conj(Rt[i]) / (1 + pow(τ₀ * ω, 2)) + pow(y, 2) * (RddfCt[i] * Ct[i] + dfCt[i] * std::conj(Rt[i]))) / (z + 1i * ω);
-    }
+      for (unsigned i = 0; i < Rt.size(); i++) {
+        Real ω = i * Δω;
+        Rtnew[i] = (1.0 + pow(y, 2) * RddfCt[i] * Rt[i]) / (z + 1i * ω);
+      }
 
-    Ct = Ctnew;
+      Rt = Rtnew;
 
-    std::vector<Real> Cnew = fft.inverse(Ctnew);
-    std::vector<Real> Rnew = fft.inverse(Rtnew);
+      for (unsigned i = 0; i < Rt.size(); i++) {
+        Real ω = i * Δω;
+        Ctnew[i] = (Γ₀ * std::conj(Rt[i]) / (1 + pow(τ₀ * ω, 2)) + pow(y, 2) * (RddfCt[i] * Ct[i] + dfCt[i] * std::conj(Rt[i]))) / (z + 1i * ω);
+      }
 
-    Real ΔC = 0;
-    for (unsigned i = 0; i < Cnew.size(); i++) {
-      ΔC += pow(Cnew[i] - C[i], 2);
-    }
+      Ct = Ctnew;
 
-    z *= Cnew[0];
+      std::vector<Real> Cnew = fft.inverse(Ctnew);
+      std::vector<Real> Rnew = fft.inverse(Rtnew);
 
-    std::cerr << "Iteration " << it << ": " << sqrt(ΔC / C.size()) << " " << z << std::endl;
+      Real ΔC = 0;
+      for (unsigned i = 0; i < Cnew.size(); i++) {
+        ΔC += pow(Cnew[i] - C[i], 2);
+      }
 
-    C = Cnew;
-    R = Rnew;
-  }
+      z *= Cnew[0];
 
   Real energy = 0;
 
@@ -176,7 +178,14 @@ int main(int argc, char* argv[]) {
     energy += y * R[i] * df(p, C[i]) * M_PI * Δτ;
   }
 
-  std::cerr << energy << std::endl;
+      C = Cnew;
+      R = Rnew;
+
+      std::cerr << "y " << y << " iteration " << it << ": " << sqrt(ΔC / C.size()) << " " << pow(y, 2) / z << " " << energy << std::endl;
+
+      if (sqrt(ΔC / C.size()) < ε) break;
+    }
+  }
 
   for (unsigned i = 0; i < C.size(); i++) {
     std::cout << i * Δτ * M_PI << " " << C[i] << std::endl;