diff options
| -rw-r--r-- | log-fourier.cpp | 2 | ||||
| -rw-r--r-- | log-fourier_integrator.cpp | 14 |
2 files changed, 10 insertions, 6 deletions
diff --git a/log-fourier.cpp b/log-fourier.cpp index f7d6f4b..47d5c5c 100644 --- a/log-fourier.cpp +++ b/log-fourier.cpp @@ -173,7 +173,7 @@ std::tuple<std::vector<Complex>, std::vector<Complex>> RddfCtdfCt(LogarithmicFou Real estimateZ(LogarithmicFourierTransform& fft, const std::vector<Real>& C, const std::vector<Complex>& Ct, const std::vector<Real>& R, const std::vector<Complex>& Rt, unsigned p, unsigned s, Real λ, Real τ₀, Real β) { auto [RddfCt, dfCt] = RddfCtdfCt(fft, C, R, p, s, λ); - Real Γ₀ = 1 + τ₀; + Real Γ₀ = 1.0; return ((2 * Γ₀ * std::conj(Rt[0]) + pow(β, 2) * (RddfCt[0] * Ct[0] + dfCt[0] * std::conj(Rt[0]))) / Ct[0]).real(); } diff --git a/log-fourier_integrator.cpp b/log-fourier_integrator.cpp index cf9819a..412b4d9 100644 --- a/log-fourier_integrator.cpp +++ b/log-fourier_integrator.cpp @@ -71,8 +71,11 @@ int main(int argc, char* argv[]) { LogarithmicFourierTransform fft(N, k, Δτ, 4); - Real Γ₀ = 1.0 + τ₀; - Real μₜ₋₁ = 1.0; + Real Γ₀ = 1.0; + Real μₜ₋₁ = Γ₀; + if (τ₀ > 0) { + μₜ₋₁ = (sqrt(1+4*Γ₀*τ₀)-1)/(2*τ₀); + } std::vector<Real> Cₜ₋₁(N); std::vector<Real> Rₜ₋₁(N); @@ -82,10 +85,10 @@ int main(int argc, char* argv[]) { if (!loadData) { /* Start from the exact solution for β = 0 */ for (unsigned n = 0; n < N; n++) { - if (τ₀ != 1) { + if (τ₀ > 0) { Cₜ₋₁[n] = Γ₀ * (exp(-μₜ₋₁ * fft.t(n)) - μₜ₋₁ * τ₀ * exp(-fft.t(n) / τ₀)) / (μₜ₋₁ - pow(μₜ₋₁, 3) * pow(τ₀, 2)); } else { - Cₜ₋₁[n] = Γ₀ * exp(-fft.t(n)) * (1 + fft.t(n)); + Cₜ₋₁[n] = Γ₀ * exp(-μₜ₋₁ * fft.t(n)) / μₜ₋₁; } Rₜ₋₁[n] = exp(-μₜ₋₁ * fft.t(n)); @@ -139,12 +142,13 @@ int main(int argc, char* argv[]) { } if (std::isnan(Cₜ[0])) { + Δβ /= 2; + β -= Δβ; Cₜ = Cₜ₋₁; Rₜ = Rₜ₋₁; Ĉₜ = Ĉₜ₋₁; Ȓₜ = Ȓₜ₋₁; μₜ = μₜ₋₁; - γ /= 2; } else { Real E = energy(fft, Cₜ, Rₜ, p, s, λ, β); |