diff options
| -rw-r--r-- | log-fourier_integrator.cpp | 93 |
1 files changed, 49 insertions, 44 deletions
diff --git a/log-fourier_integrator.cpp b/log-fourier_integrator.cpp index 7569245..94987f3 100644 --- a/log-fourier_integrator.cpp +++ b/log-fourier_integrator.cpp @@ -87,59 +87,64 @@ int main(int argc, char* argv[]) { Real β = 0; while (β < βₘₐₓ) { - Real ΔC = 100; - Real ΔC₀ = 100; - unsigned it = 0; - while (ΔC > ε) { - std::vector<Real> RddfC(N); - for (unsigned n = 0; n < N; n++) { - RddfC[n] = Rₜ[n] * ddf(λ, p, s, Cₜ[n]); - } - std::vector<Complex> RddfCt = fft.fourier(RddfC, false); + while (true) { + Real ΔC = 100; + Real ΔC₀ = 100; + unsigned it = 0; + while (ΔC > ε) { + std::vector<Real> RddfC(N); + for (unsigned n = 0; n < N; n++) { + RddfC[n] = Rₜ[n] * ddf(λ, p, s, Cₜ[n]); + } + std::vector<Complex> RddfCt = fft.fourier(RddfC, false); - std::vector<Complex> Ȓₜ₊₁(N); - std::vector<Complex> Ĉₜ₊₁(N); + std::vector<Complex> Ȓₜ₊₁(N); + std::vector<Complex> Ĉₜ₊₁(N); - 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); - } + 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); + } - std::vector<Real> Rₜ₊₁ = fft.inverse(Ȓₜ₊₁); - std::vector<Real> Cₜ₊₁ = fft.inverse(Ĉₜ₊₁); + std::vector<Real> Rₜ₊₁ = fft.inverse(Ȓₜ₊₁); + std::vector<Real> Cₜ₊₁ = fft.inverse(Ĉₜ₊₁); - ΔC = 0; - for (unsigned i = 0; i < N; i++) { - ΔC += std::norm(Ĉₜ[i] - Ĉₜ₊₁[i]); - ΔC += std::norm(Ȓₜ[i] - Ȓₜ₊₁[i]); - } - ΔC = sqrt(ΔC) / (2*N); + ΔC = 0; + for (unsigned i = 0; i < N; i++) { + ΔC += std::norm(Ĉₜ[i] - Ĉₜ₊₁[i]); + ΔC += std::norm(Ȓₜ[i] - Ȓₜ₊₁[i]); + } + ΔC = sqrt(ΔC) / (2*N); - for (unsigned i = 0; i < N; i++) { - Cₜ[i] += γ * (Cₜ₊₁[i] - Cₜ[i]); - Rₜ[i] += γ * (Rₜ₊₁[i] - Rₜ[i]); - Ĉₜ[i] += γ * (Ĉₜ₊₁[i] - Ĉₜ[i]); - Ȓₜ[i] += γ * (Ȓₜ₊₁[i] - Ȓₜ[i]); - } + for (unsigned i = 0; i < N; i++) { + Cₜ[i] += γ * (Cₜ₊₁[i] - Cₜ[i]); + Rₜ[i] += γ * (Rₜ₊₁[i] - Rₜ[i]); + Ĉₜ[i] += γ * (Ĉₜ₊₁[i] - Ĉₜ[i]); + Ȓₜ[i] += γ * (Ȓₜ₊₁[i] - Ȓₜ[i]); + } - μ = γ * μ * Cₜ[0] + (1-γ) * μ; + if (ΔC < ΔC₀) { + ΔC₀ = ΔC; + it = 0; + γ = std::min(1.001 * γ, 1.0); + } else { + it++; + } - if (ΔC < ΔC₀) { - ΔC₀ = ΔC; - it = 0; - γ = std::min(1.001 * γ, 1.0); - } else { - it++; - } + if (it > 100) { + γ = std::max(0.5 * γ, 1e-3); + it = 0; + ΔC₀ = ΔC; + } - if (it > 100) { - γ = std::max(0.5 * γ, 1e-3); - it = 0; - ΔC₀ = ΔC; + std::cerr << β << " " << μ << " " << ΔC << " " << γ; + std::cerr << "\r"; + } + if (std::abs(Cₜ[0] - 1) < ε) { + break; + } else { + μ *= Cₜ[0]; } - - std::cerr << β << " " << μ << " " << ΔC << " " << γ; - std::cerr << "\r"; } /* Integrate the energy using Simpson's rule */ |