Rewrote the original integrator to use a new approach.HEADmaster

4 files changed, 92 insertions, 28 deletions

diff --git a/Makefile b/Makefile
index 547e336..688a99c 100644
--- a/Makefile
+++ b/Makefile
@@ -8,8 +8,8 @@ walk: walk.cpp
 correlations: correlations.cpp
 	$(CC) correlations.cpp -o correlations
 
-integrator: integrator.cpp
-	$(CC) integrator.cpp -o integrator
+integrator: integrator.cpp fourier.o fourier.hpp
+	$(CC) integrator.cpp fourier.o -lfftw3 -lfftw3_omp -o integrator
 
 fourier.o: fourier.cpp fourier.hpp
 	$(CC) fourier.cpp -c -o fourier.o
diff --git a/fourier.cpp b/fourier.cpp
index 07f8fdc..bc4b633 100644
--- a/fourier.cpp
+++ b/fourier.cpp
@@ -64,6 +64,25 @@ std::vector<Complex> FourierTransform::fourier() {
   return ĉ;
 }
 
+std::vector<Real> FourierTransform::convolve(const std::vector<Real>& Γh, const std::vector<Real>& R) {
+  a[0] = 0;
+  for (unsigned i = 1; i < n; i++) {
+    a[i] = R[i];
+    a[2 * n - i] = -R[i];
+  }
+  fftw_execute(plan_r2c);
+  for (unsigned i = 1; i < n + 1; i++) {
+    â[i] *= Γh[i] * (Δτ * M_PI);
+  }
+  fftw_execute(plan_c2r);
+  std::vector<Real> dC(n);
+  for (unsigned i = 0; i < n; i++) {
+    dC[i] = a[i] * (Δω / (2 * M_PI));
+  }
+
+  return dC;
+}
+
 std::vector<Real> FourierTransform::inverse(const std::vector<Complex>& ĉ) {
   for (unsigned i = 0; i < n + 1; i++) {
     â[i] = ĉ[i];
diff --git a/fourier.hpp b/fourier.hpp
index 9451f69..791953b 100644
--- a/fourier.hpp
+++ b/fourier.hpp
@@ -33,6 +33,7 @@ public:
   std::vector<Complex> fourier();
   std::vector<Real> inverse(const std::vector<Complex>& ĉ);
   void writeToA(unsigned i, Real ai);
+  std::vector<Real> convolve(const std::vector<Real>& Γh, const std::vector<Real>& R);
 };
 
 std::string fourierFile(std::string prefix, unsigned p, unsigned s, Real λ, Real τ₀, Real y, unsigned log2n, Real τₘₐₓ);
diff --git a/integrator.cpp b/integrator.cpp
index faba06d..7e5c512 100644
--- a/integrator.cpp
+++ b/integrator.cpp
@@ -1,4 +1,6 @@
 #include "fourier.hpp"
+#include <fstream>
+#include <fftw3.h>
 #include <getopt.h>
 #include <iostream>
 
@@ -14,27 +16,36 @@ int main(int argc, char* argv[]) {
   unsigned p = 3;
   unsigned s = 4;
   Real λ = 0.5;
-  Real Δτ = 1e-3;
   Real τₘₐₓ = 1e3;
   Real τ₀ = 0;
-  Real β = 0.5;
+  Real β₀ = 0;
+  Real βₘₐₓ = 1;
+  Real Δβ = 1e-2;
   unsigned iterations = 10;
+  unsigned log2n = 8;
+  Real ε = 1e-14;
 
   int opt;
 
-  while ((opt = getopt(argc, argv, "d:T:t:b:I:")) != -1) {
+  while ((opt = getopt(argc, argv, "T:2:t:0:b:d:I:")) != -1) {
     switch (opt) {
-    case 'd':
-      Δτ = atof(optarg);
-      break;
     case 'T':
       τₘₐₓ = atof(optarg);
       break;
+    case '2':
+      log2n = atof(optarg);
+      break;
     case 't':
       τ₀ = atof(optarg);
       break;
+    case '0':
+      β₀ = atof(optarg);
+      break;
     case 'b':
-      β = atof(optarg);
+      βₘₐₓ = atof(optarg);
+      break;
+    case 'd':
+      Δβ = atof(optarg);
       break;
     case 'I':
       iterations = (unsigned)atof(optarg);
@@ -44,6 +55,13 @@ int main(int argc, char* argv[]) {
     }
   }
 
+  unsigned N = pow(2, log2n);
+
+  Real Δτ = (1 + τ₀ / 2) * τₘₐₓ / M_PI / N;
+  Real Δω = M_PI / ((1 + τ₀ / 2) * τₘₐₓ);
+
+  FourierTransform fft(N, Δω, Δτ, FFTW_ESTIMATE);
+
   Real Γ₀ = 1;
   Real μ = 1;
   if (τ₀ > 0) {
@@ -51,45 +69,71 @@ int main(int argc, char* argv[]) {
   }
 
   Real τ = 0;
-  unsigned N = τₘₐₓ / Δτ + 1;
   std::vector<Real> C(N);
   std::vector<Real> R(N);
   std::vector<Real> Γ(N);
+  std::vector<Real> Γh(N+1);
+
+  Γh[0] = Γ₀;
 
   for (unsigned i = 0; i < N; i++) {
     Real τ = i * Δτ;
+    Real ω = (i + 1) * Δω * M_PI;
     if (τ₀ > 0) {
       C[i] = (Γ₀ / μ) * (exp(-μ * τ) - μ * τ₀ * exp(-τ / τ₀)) / (1 - pow(μ * τ₀, 2));
       Γ[i] = (Γ₀ / τ₀) * exp(-τ / τ₀);
     } else {
       C[i] = (Γ₀ / μ) * exp(-μ * τ);
     }
+    Γh[i+1] = Γ₀ / (1 + pow(ω * τ₀, 2));
     R[i] = exp(-μ * τ);
   }
 
-  for (unsigned it = 0; it < iterations; it++) {
-    /* First step: integrate R from C */
-    std::vector<Real> R₊(N);
-    R₊[0] = 1;
-    for (unsigned i = 1; i < N; i++) {
-      Real I = 0;
-      for (unsigned j = 0; j <= i; j++) {
-        I += R[i - j] * ddf(λ, p, s, C[i - j]) * R[j] * Δτ;
+  for (Real β = β₀; β < βₘₐₓ; β += Δβ) {
+    Real Rerr = 100;
+    while (sqrt(Rerr / N) > ε) {
+      /* First step: integrate R from C */
+      std::vector<Real> R₊(N);
+      R₊[0] = 1;
+      for (unsigned i = 1; i < N; i++) {
+        Real I = 0;
+        for (unsigned j = 0; j <= i; j++) {
+          I += R[i - j] * ddf(λ, p, s, C[i - j]) * R[j] * Δτ;
+        }
+        Real dR = -μ * R₊[i - 1] + pow(β, 2) * I;
+        R₊[i] = R₊[i - 1] + dR * Δτ;
       }
-      Real dR = -μ * R[i] + pow(β, 2) * I;
-      R₊[i] = R₊[i - 1] + dR * Δτ;
+
+      Rerr = 0;
+      for (unsigned i = 0; i < N; i++) {
+        Rerr += pow(R[i] - R₊[i], 2);
+      }
+
+      R = R₊;
+
+      /* Second step: integrate C from R */
+      std::vector<Real> dC = fft.convolve(Γh, R);
+      Real Cₜ₊₁ = 0;
+      for (unsigned i = 0; i < N; i++) {
+        Real Cₜ = Cₜ₊₁ + dC[N - i - 1] * Δτ;
+        C[N - i - 1] = Cₜ;
+        Cₜ₊₁ = Cₜ;
+      }
+
+      /* Third step: adjust μ */
+      μ *= C[0];
+
+      std::cerr << β << " " << sqrt(Rerr / N) << std::endl;
     }
 
-    /* Second step: integrate C from R */
-  }
+    std::ofstream outfile(fourierFile("Ci", p, s, λ, τ₀, β, log2n, τₘₐₓ), std::ios::out | std::ios::binary);
+    outfile.write((const char*)(C.data()), N * sizeof(Real));
+    outfile.close();
 
-  τ = 0;
-  for (Real Ci : C) {
-    std::cout << τ << " " << Ci << std::endl;
-    τ += Δτ;
+    std::ofstream outfileR(fourierFile("Ri", p, s, λ, τ₀, β, log2n, τₘₐₓ), std::ios::out | std::ios::binary);
+    outfileR.write((const char*)(R.data()), N * sizeof(Real));
+    outfileR.close();
   }
 
-  std::cerr << - 2 * β / Γ₀ * energy(C, R, λ, p, s, Δτ) << std::endl;
-
   return 0;
 }