diff options
| -rw-r--r-- | log-fourier_integrator.cpp | 50 |
1 files changed, 38 insertions, 12 deletions
diff --git a/log-fourier_integrator.cpp b/log-fourier_integrator.cpp index 7b689f1..ca378bb 100644 --- a/log-fourier_integrator.cpp +++ b/log-fourier_integrator.cpp @@ -67,24 +67,29 @@ int main(int argc, char* argv[]) { μ = (sqrt(1+4*Γ₀*τ₀)-1)/(2*τ₀); } - std::vector<Real> Cₜ(N); - std::vector<Real> Rₜ(N); - std::vector<Complex> Ĉₜ(N); - std::vector<Complex> Ȓₜ(N); + std::vector<Real> Cₜ₋₁(N); + std::vector<Real> Rₜ₋₁(N); + std::vector<Complex> Ĉₜ₋₁(N); + std::vector<Complex> Ȓₜ₋₁(N); /* 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)); + Cₜ₋₁[n] = Γ₀ * (exp(-μ * fft.t(n)) - μ * τ₀ * exp(-fft.t(n) / τ₀)) / (μ - pow(μ, 3) * pow(τ₀, 2)); } else { - Cₜ[n] = Γ₀ * exp(-μ * fft.t(n)) / μ; + Cₜ₋₁[n] = Γ₀ * exp(-μ * fft.t(n)) / μ; } - Rₜ[n] = exp(-μ * fft.t(n)); + Rₜ₋₁[n] = exp(-μ * fft.t(n)); - Ĉₜ[n] = 2 * Γ₀ / (pow(μ, 2) + pow(fft.ν(n), 2)) / (1 + pow(τ₀ * fft.ν(n), 2)); - Ȓₜ[n] = 1.0 / (μ + 1i * fft.ν(n)); + Ĉₜ₋₁[n] = 2 * Γ₀ / (pow(μ, 2) + pow(fft.ν(n), 2)) / (1 + pow(τ₀ * fft.ν(n), 2)); + Ȓₜ₋₁[n] = 1.0 / (μ + 1i * fft.ν(n)); } + std::vector<Real> Cₜ = Cₜ₋₁; + std::vector<Real> Rₜ = Rₜ₋₁; + std::vector<Complex> Ĉₜ = Ĉₜ₋₁; + std::vector<Complex> Ȓₜ = Ȓₜ₋₁; + Real β = 0; while (β < βₘₐₓ) { Real μ₁ = 0; @@ -92,6 +97,7 @@ int main(int argc, char* argv[]) { while (true) { Real ΔC = 100; Real ΔC₀ = 100; + unsigned it = 0; while (ΔC > ε) { std::vector<Real> dfC(N); std::vector<Real> RddfC(N); @@ -107,8 +113,8 @@ int main(int argc, char* argv[]) { for (unsigned n = 0; n < N; n++) { Ȓₜ₊₁[n] = (1.0 + pow(β, 2) * RddfCt[n] * Ȓₜ[n]) / (μ + 1i * fft.ν(n)); -// Ĉₜ₊₁[n] = - 2 * Γ₀ * Ȓₜ₊₁[n].imag() / (1 + pow(τ₀ * fft.ν(n), 2)) / fft.ν(n); - Ĉₜ₊₁[n] = (2 * Γ₀ * std::conj(Ȓₜ[n]) / (1 + pow(τ₀ * fft.ν(n), 2)) + pow(β, 2) * (RddfCt[n] * Ĉₜ[n] + dfCt[n] * std::conj(Ȓₜ[n]))) / (μ + 1i * fft.ν(n)); + Ĉₜ₊₁[n] = - 2 * Γ₀ * Ȓₜ₊₁[n].imag() / (1 + pow(τ₀ * fft.ν(n), 2)) / fft.ν(n); +// Ĉₜ₊₁[n] = (2 * Γ₀ * std::conj(Ȓₜ[n]) / (1 + pow(τ₀ * fft.ν(n), 2)) + pow(β, 2) * (RddfCt[n] * Ĉₜ[n] + dfCt[n] * std::conj(Ȓₜ[n]))) / (μ + 1i * fft.ν(n)); } std::vector<Real> Rₜ₊₁ = fft.inverse(Ȓₜ₊₁); @@ -130,14 +136,29 @@ int main(int argc, char* argv[]) { if (ΔC < ΔC₀) { ΔC₀ = ΔC; - γ = std::min(1.01 * γ, 1.0); + it = 0; + γ = std::min(1.001 * γ, 1.0); } else { + it++; + } + + if (it > 100) { γ = std::max(0.5 * γ, 1e-3); + it = 0; + ΔC₀ = ΔC; } std::cerr << β << " " << μ << " " << ΔC << " " << γ; std::cerr << "\r"; } + if (std::isnan(Cₜ[0])) { + Cₜ = Cₜ₋₁; + Rₜ = Rₜ₋₁; + Ĉₜ = Ĉₜ₋₁; + Ȓₜ = Ȓₜ₋₁; + Δβ /= 2; + β -= Δβ; + } else { if (pow(Cₜ[0] - 1, 2) < ε) { break; } @@ -160,6 +181,7 @@ int main(int argc, char* argv[]) { } μ = (μ₁ + μ₂) / 2; } + } } /* Integrate the energy using Simpson's rule */ @@ -180,6 +202,10 @@ int main(int argc, char* argv[]) { std::cerr << β << " " << μ << " " << Ĉₜ[0].real() << " " << E << " " << γ << std::endl; β += Δβ; + Cₜ₋₁ = Cₜ; + Rₜ₋₁ = Rₜ; + Ĉₜ₋₁ = Ĉₜ; + Ȓₜ₋₁ = Ȓₜ; } for (unsigned i = 0; i < N; i++) { |