From ae23cdf557611d213b6aa03597ac4f8dd139f98b Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Sat, 19 Apr 2025 14:33:13 -0300 Subject: Changed interative scheme and many variable names --- log-fourier_integrator.cpp | 148 +++++++++++++++++++++++---------------------- 1 file changed, 75 insertions(+), 73 deletions(-) (limited to 'log-fourier_integrator.cpp') diff --git a/log-fourier_integrator.cpp b/log-fourier_integrator.cpp index fd8b220..bf83da5 100644 --- a/log-fourier_integrator.cpp +++ b/log-fourier_integrator.cpp @@ -64,98 +64,100 @@ int main(int argc, char* argv[]) { μ = (sqrt(1+4*Γ₀*τ₀)-1)/(2*τ₀); } - std::vector C(N); - std::vector R(N); - std::vector Ct(N); - std::vector Rt(N); + std::vector Cₜ(N); + std::vector Rₜ(N); + std::vector Ĉₜ(N); + std::vector Ȓₜ(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)); - Ct[n] = 2 * Γ₀ / (pow(μ, 2) + pow(fft.ν(n), 2)) / (1 + pow(τ₀ * fft.ν(n), 2)); - Rt[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)); } Real β = 0; while (β < βₘₐₓ) { - Real ΔC = 100; - while (ΔC > ε) { - std::vector RddfC(N); - std::vector 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 RddfCt = fft.fourier(RddfC, false); - std::vector dfCt = fft.fourier(dfC, true); - - std::vector Rtnew(N); - std::vector 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 * Γ₀ * Rtnew[n].imag() / (1 + pow(τ₀ * fft.ν(n), 2)) / fft.ν(n); - } - std::vector Rnew = fft.inverse(Rtnew); - for (unsigned n = 0; n < N; n++) { - RddfC[n] = Rnew[n] * ddf(λ, p, s, C[n]); - } - RddfCt = fft.fourier(RddfC, false); - - for (unsigned n = 0; n < N; n++) { - Ctnew[n] = (2 * Γ₀ * std::conj(Rtnew[n]) / (1 + pow(τ₀ * fft.ν(n), 2)) + pow(β, 2) * (RddfCt[n] * Ct[n] + dfCt[n] * std::conj(Rtnew[n]))) / (μ + 1i * fft.ν(n)); - } - + Real ΔC = 100; + while (ΔC > ε) { + std::vector RddfC(N); +// std::vector 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 RddfCt = fft.fourier(RddfC, false); +// std::vector dfCt = fft.fourier(dfC, true); + + std::vector Ȓₜ₊₁(N); + std::vector Ĉₜ₊₁(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 Rₜ₊₁ = fft.inverse(Ȓₜ₊₁); + /* + for (unsigned n = 0; n < N; n++) { + RddfC[n] = Rₜ₊₁[n] * ddf(λ, p, s, Cₜ[n]); + } + RddfCt = fft.fourier(RddfC, false); + + for (unsigned n = 0; n < N; 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 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); + + 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]; - std::vector Cnew = fft.inverse(Ctnew); - - ΔC = 0; - for (unsigned i = 0; i < N; i++) { - ΔC += std::norm(Ct[i] - Ctnew[i]); - ΔC += std::norm(Rt[i] - Rtnew[i]); +// std::cerr << ΔC << std::endl; } - Δ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]); + /* Integrate the energy using Simpson's rule */ + Real E = 0; + for (unsigned n = 0; n < N/2-1; n++) { + Real h₂ₙ = fft.t(2*n+1) - fft.t(2*n); + Real h₂ₙ₊₁ = fft.t(2*n+2) - fft.t(2*n+1); + Real f₂ₙ = Rₜ[2*n] * df(λ, p, s, Cₜ[2*n]); + Real f₂ₙ₊₁ = Rₜ[2*n+1] * df(λ, p, s, Cₜ[2*n+1]); + Real f₂ₙ₊₂ = Rₜ[2*n+2] * df(λ, p, s, Cₜ[2*n+2]); + E += (h₂ₙ + h₂ₙ₊₁) / 6 * ( + (2 - h₂ₙ₊₁ / h₂ₙ) * f₂ₙ + + pow(h₂ₙ + h₂ₙ₊₁, 2) / (h₂ₙ * h₂ₙ₊₁) * f₂ₙ₊₁ + + (2 - h₂ₙ / h₂ₙ₊₁) * f₂ₙ₊₂ + ); } + E *= β; - μ *= C[0]; - -// std::cerr << ΔC << std::endl; - } - - /* Integrate the energy using Simpson's rule */ - Real E = 0; - for (unsigned i = 0; i < N/2-1; i++) { - Real h₂ᵢ = fft.t(2*i+1) - fft.t(2*i); - Real h₂ᵢ₊₁ = fft.t(2*i+2) - fft.t(2*i+1); - Real f₂ᵢ = R[2*i] * df(λ, p, s, C[2*i]); - Real f₂ᵢ₊₁ = R[2*i+1] * df(λ, p, s, C[2*i+1]); - Real f₂ᵢ₊₂ = R[2*i+2] * df(λ, p, s, C[2*i+2]); - E += (h₂ᵢ + h₂ᵢ₊₁) / 6 * ( - (2 - h₂ᵢ₊₁ / h₂ᵢ) * f₂ᵢ - + pow(h₂ᵢ + h₂ᵢ₊₁, 2) / (h₂ᵢ * h₂ᵢ₊₁) * f₂ᵢ₊₁ - + (2 - h₂ᵢ / h₂ᵢ₊₁) * f₂ᵢ₊₂ - ); - } - E *= β; - - std::cerr << β << " " << μ << " " << Ct[0].real() << " " << E << std::endl; - β += Δβ; + std::cerr << β << " " << μ << " " << Ĉₜ[0].real() << " " << E << std::endl; + β += Δβ; } for (unsigned i = 0; i < N; i++) { - std::cout << fft.t(i) << " " << C[i] << std::endl; + std::cout << fft.t(i) << " " << Cₜ[i] << std::endl; } return 0; -- cgit v1.2.3-70-g09d2