Cleaned up the code a bit and rescaled time for the persistant case.

1 files changed, 26 insertions, 39 deletions

diff --git a/integrator.cpp b/integrator.cpp
index 7d0ad1d..e184de6 100644
--- a/integrator.cpp
+++ b/integrator.cpp
@@ -1,21 +1,10 @@
 #include <getopt.h>
 #include <vector>
-#include <array>
 #include <cmath>
 #include <iostream>
 
 using Real = double;
 
-class Point : public std::array<Real, 2> {
-public:
-  Real τ() const {
-    return front();
-  }
-  Real C() const {
-    return back();
-  }
-};
-
 unsigned p = 2;
 
 Real f(Real q) {
@@ -30,38 +19,34 @@ Real ddf(Real q) {
   return 0.5 * p * (p - 1) * pow(q, p - 2);
 }
 
-Real integrate(const std::vector<Point>& Cₜ, double τ₀) {
+Real integrate(const std::vector<Real>& C, Real Δτ, Real τ₀) {
   Real I = 0;
 #pragma omp parallel for reduction(+:I)
-  for (unsigned i = 0; i < Cₜ.size() - 1; i++) {
-    unsigned j = Cₜ.size() - 1 - i;
-    Real Δτ = Cₜ[i + 1].τ() - Cₜ[i].τ();
-    Real C = (Cₜ[j].C() + Cₜ[j - 1].C()) / 2;
-    Real dC = (Cₜ[i + 1].C() - Cₜ[i].C()) / Δτ;
+  for (unsigned σ = 0; σ < C.size() - 1; σ++) {
+    unsigned τ_σ = C.size() - 1 - σ;
+    Real Cτ_σ = (C[τ_σ] + C[τ_σ - 1]) / 2;
+    Real dCσ = (C[σ + 1] - C[σ]) / Δτ;
 
     Real dddC = 0;
-    if (i > 4 && i < Cₜ.size() && Cₜ.size() > 4) {
-      dddC += (Cₜ[j].C() - 3 * Cₜ[j+1].C() + 3 * Cₜ[j+2].C() - Cₜ[j+3].C()) / pow(Δτ, 3);
-      dddC += (Cₜ[j+1].C() - 3 * Cₜ[j+2].C() + 3 * Cₜ[j+3].C() - Cₜ[j+4].C()) / pow(Δτ, 3);
-      dddC /= 2;
+    if (σ > 3 && σ < C.size() && C.size() > 3) {
+      dddC += (C[τ_σ] - 3 * C[τ_σ+1] + 3 * C[τ_σ+2] - C[τ_σ+3]) / pow(Δτ, 3);
     }
-    I += Δτ * df(C) * (dC - pow(τ₀, 2) * dddC);
+    I += Δτ * df(Cτ_σ) * (dCσ - pow(τ₀, 2) * dddC);
   }
   return I;
 }
 
-Real energy(const std::vector<Point>& Ct, Real τ₀) {
+Real energy(const std::vector<Real>& C, Real Δτ, Real τ₀) {
   Real I = 0;
-  for (unsigned i = 0; i < Ct.size() - 1; i++) {
-    Real Δτ = Ct[i + 1].τ() - Ct[i].τ();
-    Real C = (Ct[i].C() + Ct[i + 1].C()) / 2;
-    Real dC = (Ct[i + 1].C() - Ct[i].C()) / Δτ;
+  for (unsigned σ = 0; σ < C.size() - 1; σ++) {
+    Real Cσ = (C[σ] + C[σ + 1]) / 2;
+    Real dC = (C[σ + 1] - C[σ]) / Δτ;
 
     Real dddC = 0;
-    if (i > 1 && i < Ct.size() - 2 && Ct.size() > 3) {
-      dddC = (Ct[i+1].C() - 3 * Ct[i].C() + 3 * Ct[i-1].C() - Ct[i-2].C()) / pow(Δτ, 3);
+    if (σ > 1 && σ < C.size() - 2 && C.size() > 3) {
+      dddC = (C[σ+1] - 3 * C[σ] + 3 * C[σ-1] - C[σ-2]) / pow(Δτ, 3);
     }
-    I += Δτ * df(C) * (dC - pow(τ₀, 2) * dddC);
+    I += Δτ * df(Cσ) * (dC - pow(τ₀, 2) * dddC);
   }
   return I;
 }
@@ -93,20 +78,22 @@ int main(int argc, char* argv[]) {
     }
   }
 
-  Real z = (sqrt(1 + 2 * τ₀) - 1) / (2 * τ₀);
+  Real z = 0.5;
+  Real Γ₀ = 0.5 * (2 * τ₀) / (sqrt(1 + 2 * τ₀) - 1);
 
-  std::vector<Point> Cₜ;
-  Cₜ.reserve(τₘₐₓ / Δτ + 1);
+  Real τ = 0;
+  std::vector<Real> C;
+  C.reserve(τₘₐₓ / Δτ + 1);
 
-  Cₜ.push_back({0, 1});
+  C.push_back(1);
 
-  while (Cₜ.back().τ() < τₘₐₓ) {
-    Real dC = -z * Cₜ.back().C() - 2 * pow(y, 2) * integrate(Cₜ, τ₀);
-    Cₜ.push_back({Cₜ.back().τ() + Δτ, Cₜ.back().C() + Δτ * dC});
-    std::cout << Cₜ.back().τ() << " " << Cₜ.back().C() << std::endl;
+  while (std::cout << τ << " " << C.back() << std::endl, τ < τₘₐₓ) {
+    τ += Δτ;
+    Real dC = -z * C.back() - 2 / Γ₀ * pow(y, 2) * integrate(C, Δτ, τ₀);
+    C.push_back(C.back() + Δτ * dC);
   }
 
-  std::cerr << - 2 * y * energy(Cₜ, τ₀) << std::endl;
+  std::cerr << - 2 * y / Γ₀ * energy(C, Δτ, τ₀) << std::endl;
 
   return 0;
 }