Began trying to find a convergent iteration for determining z

1 files changed, 36 insertions, 12 deletions

diff --git a/integrator.cpp b/integrator.cpp
index 349e2ff..a94f37e 100644
--- a/integrator.cpp
+++ b/integrator.cpp
@@ -19,11 +19,14 @@ Real ddf(Real q) {
   return 0.5 * p * (p - 1) * pow(q, p - 2);
 }
 
-Real integrate(const std::vector<Real>& C) {
+Real integrate(const std::vector<Real>& C, signed τ = std::numeric_limits<unsigned>::max()) {
   Real I = 0;
+  if (τ > C.size() - 1) {
+    τ = C.size() - 1;
+  }
 #pragma omp parallel for reduction(+:I)
-  for (unsigned σ = 0; σ < C.size() - 1; σ++) {
-    unsigned τ_σ = C.size() - 1 - σ;
+  for (unsigned σ = 0; σ < τ; σ++) {
+    unsigned τ_σ = τ - σ;
     Real Cτ_σ = (C[τ_σ] + C[τ_σ - 1]) / 2;
     Real dCσ = C[σ + 1] - C[σ];
 
@@ -99,10 +102,11 @@ int main(int argc, char* argv[]) {
   Real τₘₐₓ = 1e3;
   Real τ₀ = 0;
   Real y = 0.5;
+  unsigned iterations = 10;
 
   int opt;
 
-  while ((opt = getopt(argc, argv, "d:T:t:y:")) != -1) {
+  while ((opt = getopt(argc, argv, "d:T:t:y:I:")) != -1) {
     switch (opt) {
     case 'd':
       Δτ = atof(optarg);
@@ -116,15 +120,16 @@ int main(int argc, char* argv[]) {
     case 'y':
       y = atof(optarg);
       break;
+    case 'I':
+      iterations = (unsigned)atof(optarg);
+      break;
     default:
       exit(1);
     }
   }
 
-//  Real z = (sqrt(1 + 2 * τ₀) - 1) / (2 * τ₀);
-//  Real z = 0.4833773593561142778087;
-  Real z = 0.4794707565634420155347 - 2 * pow(y, 2);
-  Real Γ₀ = 1.0;
+  Real z = 0.5;
+  Real Γ₀ = 1;
 
   Real τ = 0;
   std::vector<Real> C;
@@ -135,19 +140,20 @@ int main(int argc, char* argv[]) {
 //  while (std::cout << τ << " " << C.back() << std::endl, τ < τₘₐₓ) {
   while (τ < τₘₐₓ) {
     τ += Δτ;
-    Real dC = -z * C.back() - 2 / Γ₀ * pow(y, 2) * integrate(C);
+    Real dC = -(z - 2 * pow(y, 2)) * C.back() - 2 / Γ₀ * pow(y, 2) * integrate(C);
     C.push_back(C.back() + Δτ * dC);
   }
 
 
-  for (unsigned it = 0; it < 20; it++) {
+  for (unsigned it = 0; it < iterations; it++) {
+
     τ = 0;
     std::vector<Real> C2;
     C2.reserve(τₘₐₓ / Δτ + 1);
     C2.push_back(1);
     while (τ < τₘₐₓ) {
       τ += Δτ;
-      Real dC = -z * C2.back() + integrateDelay(C, C2.size() - 1, Δτ, τ₀) - 2 / Γ₀ * pow(y, 2) * (integrate(C2) - pow(τ₀ / Δτ, 2) * integratePast(C, C2.size()-1));
+      Real dC = -(z - 2 * pow(y, 2)) * C2.back() + integrateDelay(C, C2.size() - 1, Δτ, τ₀) - 2 / Γ₀ * pow(y, 2) * (integrate(C2) - pow(τ₀ / Δτ, 2) * integratePast(C, C2.size()-1));
       C2.push_back(C2.back() + Δτ * dC);
     }
 
@@ -157,10 +163,28 @@ int main(int argc, char* argv[]) {
       error += pow(C[i] - C2[i], 2);
     }
 
-    std::cerr << "Iteration " << it << ": " << sqrt(error / C.size()) << std::endl;
+    std::cerr << "Iteration " << it << ": " << sqrt(error / C.size()) << " " << z << std::endl;
 
     C = C2;
   }
+  /*
+    Real zNew = (2.0 * ((C[2] - 2 * C[1] + C[0]) / pow(Δτ, 2) - pow(τ₀, 2) * (C[4] - 4 * C[3] + 6 * C[2] - 4 * C[1] + C[0]) / pow(Δτ, 4)));
+    Real zNew = (2.0 * ((C[2] - 2 * C[1] + C[0]) / pow(Δτ, 2) - pow(τ₀, 2) * (C[4] - 4 * C[3] + 6 * C[2] - 4 * C[1] + C[0]) / pow(Δτ, 4)));
+//    Real zNew = (2.0 * ((83 * C[6] - 245 * C[5] + 101 * C[4] + 254 * C[3] - 31 * C[2] - 377 * C[1] + 215 * C[0]) / (132 * pow(Δτ, 2)) - pow(τ₀, 2) * (3 * C[6] - 7 * C[5] + C[4] + 6 * C[3] + C[2] - 7 * C[1] + 3 * C[0]) / (11 * pow(Δτ, 4))));
+    z = z / zNew;
+  τ = 0;
+  C.clear();
+  C.reserve(τₘₐₓ / Δτ + 1);
+
+  C.push_back(1);
+
+//  while (std::cout << τ << " " << C.back() << std::endl, τ < τₘₐₓ) {
+  while (τ < τₘₐₓ) {
+    τ += Δτ;
+    Real dC = -z * C.back() - 2 / Γ₀ * pow(y, 2) * integrate(C);
+    C.push_back(C.back() + Δτ * dC);
+  }
+  */
 
   τ = 0;
   for (Real Ci : C) {