Fixed bug and got integrator working

1 files changed, 21 insertions, 4 deletions

diff --git a/log-fourier_integrator.cpp b/log-fourier_integrator.cpp
index 177db46..eda4945 100644
--- a/log-fourier_integrator.cpp
+++ b/log-fourier_integrator.cpp
@@ -12,7 +12,7 @@ int main(int argc, char* argv[]) {
   Real Δτ = 0.1;
   Real k = 0.1;
 
-  Real ε = 1e-16;
+  Real ε = 1e-13;
   Real γ = 1;
   Real βₘₐₓ = 0.7;
   Real Δβ = 0.01;
@@ -99,7 +99,8 @@ int main(int argc, char* argv[]) {
 
     for (unsigned n = 0; n < N; n++) {
       Rtnew[n] = (1.0 + pow(β, 2) * RddfCt[n] * Rt[n]) / (μ + 1i * fft.ν(n));
-      Ctnew[n] = (2 * Γ₀ * std::conj(Rt[n]) / (1 + pow(τ₀ * fft.ν(n), 2)) + pow(β, 2) * (RddfCt[n] * Ct[n] + dfCt[n] * std::conj(Rt[n]))) / (μ + 1i * fft.ν(n));
+      Ctnew[n] = (2 * Γ₀ * std::conj(Rtnew[n]) / (1 + pow(τ₀ * fft.ν(n), 2)) + pow(β, 2) * (RddfCt[n] * Ct[n] + dfCt[n] * std::conj(Rtnew[n]))) / (μ + 1i * fft.ν(n));
+//      Ctnew[n] = - 2 * Γ₀ * Rtnew[n].imag() / (1 + pow(τ₀ * fft.ν(n), 2)) / fft.ν(n);
     }
 
     std::vector<Real> Cnew = fft.inverse(Ctnew);
@@ -124,12 +125,28 @@ int main(int argc, char* argv[]) {
 //    std::cerr << ΔC << std::endl;
   }
 
-  std::cerr << β << " " << μ << " " << Ct[0].real() << std::endl;
+  /* Integrate the energy using Simpson's rule */
+  Real E = 0;
+  for (unsigned i = 0; i < N/2-1; i++) {
+    Real h₂ᵢ   = fft.t(2*i+1) - fft.t(2*i);
+    Real h₂ᵢ₊₁ = fft.t(2*i+2) - fft.t(2*i+1);
+    Real f₂ᵢ   = R[2*i]   * df(λ, p, s, C[2*i]);
+    Real f₂ᵢ₊₁ = R[2*i+1] * df(λ, p, s, C[2*i+1]);
+    Real f₂ᵢ₊₂ = R[2*i+2] * df(λ, p, s, C[2*i+2]);
+    E += (h₂ᵢ + h₂ᵢ₊₁) / 6 * (
+          (2 - h₂ᵢ₊₁ / h₂ᵢ) * f₂ᵢ
+          + pow(h₂ᵢ + h₂ᵢ₊₁, 2) / (h₂ᵢ * h₂ᵢ₊₁) * f₂ᵢ₊₁
+          + (2 - h₂ᵢ / h₂ᵢ₊₁) * f₂ᵢ₊₂
+        );
+  }
+  E *= β;
+
+  std::cerr << β << " " << μ << " " << Ct[0].real() << " " << E << std::endl;
    β += Δβ;
   }
 
   for (unsigned i = 0; i < N; i++) {
-    std::cout << fft.t(i) << " " << Rt[i].imag() << std::endl;
+    std::cout << fft.t(i) << " " << C[i] << std::endl;
   }
 
   return 0;