Move Δy when there is a failure

1 files changed, 12 insertions, 6 deletions

diff --git a/fourier_integrator.cpp b/fourier_integrator.cpp
index 2ee304d..8588354 100644
--- a/fourier_integrator.cpp
+++ b/fourier_integrator.cpp
@@ -71,7 +71,7 @@ int main(int argc, char* argv[]) {
   unsigned log2n = 8;
   Real τₘₐₓ = 20 / M_PI;
 
-  unsigned maxIterations = 1000;
+  unsigned maxIterations = 100;
   Real ε = 1e-14;
 
   int opt;
@@ -132,7 +132,9 @@ int main(int argc, char* argv[]) {
   Real y = 0;
 
   while (y += Δy, y <= yₘₐₓ) {
-    for (unsigned it = 0; it < maxIterations; it++) {
+    unsigned it = 0;
+    while (true) {
+      it++;
       std::vector<Real> RddfC(C.size());
       for (unsigned i = 0; i < C.size(); i++) {
         RddfC[i] = R[i] * ddf(p, C[i]);
@@ -172,11 +174,11 @@ int main(int argc, char* argv[]) {
 
       z *= Cnew[0];
 
-  Real energy = 0;
+      Real energy = 0;
 
-  for (unsigned i = 0; i < n; i++) {
-    energy += y * R[i] * df(p, C[i]) * M_PI * Δτ;
-  }
+      for (unsigned i = 0; i < n; i++) {
+        energy += y * R[i] * df(p, C[i]) * M_PI * Δτ;
+      }
 
       C = Cnew;
       R = Rnew;
@@ -184,6 +186,10 @@ int main(int argc, char* argv[]) {
       std::cerr << "y " << y << " iteration " << it << ": " << sqrt(ΔC / C.size()) << " " << pow(y, 2) / z << " " << energy << std::endl;
 
       if (sqrt(ΔC / C.size()) < ε) break;
+      if (it > maxIterations) {
+        y -= Δy;
+        Δy /= 2;
+      }
     }
   }