From 717558730554dbb470fbda0700f14fc43595063d Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Fri, 18 Apr 2025 20:04:58 -0300
Subject: Got the log-fourier method working, but integrator is not working
 correctly

---
 log-fourier_integrator.cpp | 110 +++++++++++++++++++++++++++------------------
 1 file changed, 66 insertions(+), 44 deletions(-)

(limited to 'log-fourier_integrator.cpp')

diff --git a/log-fourier_integrator.cpp b/log-fourier_integrator.cpp
index 247270b..177db46 100644
--- a/log-fourier_integrator.cpp
+++ b/log-fourier_integrator.cpp
@@ -1,27 +1,25 @@
 #include "fourier.hpp"
 #include <getopt.h>
 #include <iostream>
-#include <fstream>
 
 int main(int argc, char* argv[]) {
   unsigned p = 2;
   unsigned s = 2;
   Real λ = 0.5;
   Real τ₀ = 0;
-  Real β₀ = 0;
-  Real yₘₐₓ = 0.5;
-  Real Δy = 0.05;
 
   unsigned log2n = 8;
-  Real τₘₐₓ = 20;
+  Real Δτ = 0.1;
+  Real k = 0.1;
 
-  unsigned maxIterations = 1000;
-  Real ε = 1e-14;
+  Real ε = 1e-16;
   Real γ = 1;
+  Real βₘₐₓ = 0.7;
+  Real Δβ = 0.01;
 
   int opt;
 
-  while ((opt = getopt(argc, argv, "p:s:2:T:t:0:y:d:I:g:l")) != -1) {
+  while ((opt = getopt(argc, argv, "p:s:2:T:t:b:d:g:k:D:")) != -1) {
     switch (opt) {
     case 'p':
       p = atoi(optarg);
@@ -32,27 +30,24 @@ int main(int argc, char* argv[]) {
     case '2':
       log2n = atoi(optarg);
       break;
-    case 'T':
-      τₘₐₓ = atof(optarg);
-      break;
     case 't':
       τ₀ = atof(optarg);
       break;
-    case '0':
-      β₀ = atof(optarg);
-      break;
-    case 'y':
-      yₘₐₓ = atof(optarg);
+    case 'b':
+      βₘₐₓ = atof(optarg);
       break;
     case 'd':
-      Δy = atof(optarg);
-      break;
-    case 'I':
-      maxIterations = (unsigned)atof(optarg);
+      Δβ = atof(optarg);
       break;
     case 'g':
       γ = atof(optarg);
       break;
+    case 'k':
+      k = atof(optarg);
+      break;
+    case 'D':
+      Δτ = atof(optarg);
+      break;
     default:
       exit(1);
     }
@@ -60,10 +55,7 @@ int main(int argc, char* argv[]) {
 
   unsigned N = pow(2, log2n);
 
-  Real Δτ = 1e-2;
-  Real Δω = 1e-2;
-  Real Δs = 1e-2;
-  Real k = 0.1;
+  LogarithmicFourierTransform fft(N, k, Δτ, 4);
 
   Real Γ₀ = 1.0;
   Real μ = Γ₀;
@@ -73,9 +65,10 @@ int main(int argc, char* argv[]) {
 
   std::vector<Real> C(N);
   std::vector<Real> R(N);
+  std::vector<Complex> Ct(N);
+  std::vector<Complex> Rt(N);
 
-  LogarithmicFourierTransform fft(N, k, Δω, Δτ, Δs);
-
+  // start from the exact solution for β = 0
   for (unsigned n = 0; n < N; n++) {
     if (τ₀ > 0) {
       C[n] = Γ₀ * (exp(-μ * fft.t(n)) - μ * τ₀ * exp(-fft.t(n) / τ₀)) / (μ - pow(μ, 3) * pow(τ₀, 2));
@@ -83,32 +76,61 @@ int main(int argc, char* argv[]) {
       C[n] = Γ₀ * exp(-μ * fft.t(n)) / μ;
     }
     R[n] = exp(-μ * fft.t(n));
+
+    Ct[n] = 2 * Γ₀ / (pow(μ, 2) + pow(fft.ν(n), 2)) / (1 + pow(τ₀ * fft.ν(n), 2));
+    Rt[n] = 1.0 / (μ + 1i * fft.ν(n));
   }
 
-  std::vector<Complex> Ct = fft.fourier(C, true);
-  std::vector<Complex> Rt = fft.fourier(R, false);
+  Real β = 0;
+  while (β < βₘₐₓ) {
+  Real ΔC = 100;
+  while (ΔC > ε) {
+    std::vector<Real> RddfC(N);
+    std::vector<Real> dfC(N);
+    for (unsigned n = 0; n < N; n++) {
+      RddfC[n] = R[n] * ddf(λ, p, s, C[n]);
+      dfC[n] = df(λ, p, s, C[n]);
+    }
+    std::vector<Complex> RddfCt = fft.fourier(RddfC, false);
+    std::vector<Complex> dfCt = fft.fourier(dfC, true);
 
-  /*
-  for (unsigned n = 0; n < N; n++) {
-    std::cout << fft.t(n) << " " << C[n] << std::endl;
+    std::vector<Complex> Rtnew(N);
+    std::vector<Complex> Ctnew(N);
+
+    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));
+    }
+
+    std::vector<Real> Cnew = fft.inverse(Ctnew);
+    std::vector<Real> Rnew = fft.inverse(Rtnew);
+
+    ΔC = 0;
+    for (unsigned i = 0; i < N; i++) {
+      ΔC += std::norm(Ct[i] - Ctnew[i]);
+      ΔC += std::norm(Rt[i] - Rtnew[i]);
+    }
+    ΔC = sqrt(ΔC) / (2*N);
+
+    for (unsigned i = 0; i < N; i++) {
+      C[i] += γ * (Cnew[i] - C[i]);
+      R[i] += γ * (Rnew[i] - R[i]);
+      Ct[i] += γ * (Ctnew[i] - Ct[i]);
+      Rt[i] += γ * (Rtnew[i] - Rt[i]);
+    }
+
+    μ *= C[0];
+
+//    std::cerr << ΔC << std::endl;
   }
-  */
 
-  for (unsigned n = 0; n < N; n++) {
-    std::cout << fft.ν(n) << " " << Ct[n].real() << std::endl;
+  std::cerr << β << " " << μ << " " << Ct[0].real() << std::endl;
+   β += Δβ;
   }
 
-  /*
-  // start from the exact solution for τ₀ = 0
-  for (unsigned i = 0; i < N + 1; i++) {
-    Real ω = i * Δω;
-    Ct[i] = 2 * Γ₀ / (pow(μ, 2) + pow(ω, 2)) / (1 + pow(τ₀ * ω, 2));
-    Rt[i] = 1.0 / (μ + 1i * ω);
-    Γt[i] = Γ₀ / (1 + pow(τ₀ * ω, 2));
+  for (unsigned i = 0; i < N; i++) {
+    std::cout << fft.t(i) << " " << Rt[i].imag() << std::endl;
   }
-  C = fft.inverse(Ct);
-  R = fft.inverse(Rt);
-  */
 
   return 0;
 }
-- 
cgit v1.2.3-70-g09d2