diff options
| -rw-r--r-- | log-fourier_integrator.cpp | 9 |
1 files changed, 7 insertions, 2 deletions
diff --git a/log-fourier_integrator.cpp b/log-fourier_integrator.cpp index 2f80f4c..8a2fe5c 100644 --- a/log-fourier_integrator.cpp +++ b/log-fourier_integrator.cpp @@ -94,6 +94,7 @@ int main(int argc, char* argv[]) { std::vector<Real> Rₜ₋₁(N); std::vector<Complex> Ĉₜ₋₁(N); std::vector<Complex> Ȓₜ₋₁(N); + std::vector<Real> Γ(N); if (!loadData) { /* Start from the exact solution for β = 0 */ @@ -111,6 +112,7 @@ int main(int argc, char* argv[]) { Ĉₜ₋₁[n] = 2 * Γ₀ / (pow(μ₀, 2) + pow(fft.ν(n), 2)) / (1 + pow(τ₀ * fft.ν(n), 2)); Ȓₜ₋₁[n] = (Real)1.0 / (μ₀ + II * fft.ν(n)); + Γ[n] = Γ₀ / (1 + std::pow(τ₀ * fft.ν(n), 2)); } } else { logFourierLoad(Cₜ₋₁, Rₜ₋₁, Ĉₜ₋₁, Ȓₜ₋₁, p, s, λ, τ₀, β₀, log2n, Δτ, logShift); @@ -138,7 +140,7 @@ int main(int argc, char* argv[]) { while (std::abs(C₀ - 1) > ε) { for (unsigned n = 0; n < N; n++) { - Ĉₜ₊₁[n] = ((2 * Γ₀ * std::conj(Ȓₜ[n]) / (1 + std::pow(τ₀ * fft.ν(n), 2)) + std::pow(β, 2) * (RddfCt[n] * Ĉₜ[n] + dfCt[n] * std::conj(Ȓₜ[n]))) / (μₜ + II * fft.ν(n))).real(); + Ĉₜ₊₁[n] = ((2 * Γ[n] * std::conj(Ȓₜ[n]) + std::pow(β, 2) * (RddfCt[n] * Ĉₜ[n] + dfCt[n] * std::conj(Ȓₜ[n]))) / (μₜ + II * fft.ν(n))).real(); } C₀ = C0(fft, Ĉₜ₊₁); if (C₀ > 1) { @@ -155,7 +157,7 @@ int main(int argc, char* argv[]) { ΔCₜ = 0; for (unsigned n = 0; n < N; n++) { - ΔCₜ += std::norm((2 * Γ₀ * std::conj(Ȓₜ[n]) / (1 + std::pow(τ₀ * fft.ν(n), 2)) + std::pow(β, 2) * (RddfCt[n] * Ĉₜ[n] + dfCt[n] * std::conj(Ȓₜ[n]))) - Ĉₜ[n] * (μₜ + II * fft.ν(n))); + ΔCₜ += std::norm((2 * Γ[n] * std::conj(Ȓₜ[n]) + std::pow(β, 2) * (RddfCt[n] * Ĉₜ[n] + dfCt[n] * std::conj(Ȓₜ[n]))) - Ĉₜ[n] * (μₜ + II * fft.ν(n))); ΔCₜ += std::norm(((Real)1.0 + std::pow(β, 2) * RddfCt[n] * Ȓₜ[n]) - Ȓₜ[n] * (μₜ + II * fft.ν(n))); } ΔCₜ = sqrt(ΔCₜ) / (2*N); @@ -179,6 +181,9 @@ int main(int argc, char* argv[]) { } if (std::isnan(Cₜ[0])) { + β -= Δβ; + Δβ *= 0.1; + β += Δβ; Cₜ = Cₜ₋₁; Rₜ = Rₜ₋₁; Ĉₜ = Ĉₜ₋₁; |