diff options
| -rw-r--r-- | log-fourier_integrator.cpp | 18 |
1 files changed, 13 insertions, 5 deletions
diff --git a/log-fourier_integrator.cpp b/log-fourier_integrator.cpp index f875254..4c22ba2 100644 --- a/log-fourier_integrator.cpp +++ b/log-fourier_integrator.cpp @@ -100,25 +100,33 @@ int main(int argc, char* argv[]) { Real ΔC₀ = 100; unsigned it = 0; while (ΔC > ε) { -// std::vector<Real> dfC(N); + std::vector<Real> dfC(N); std::vector<Real> RddfC(N); for (unsigned n = 0; n < N; n++) { RddfC[n] = Rₜ[n] * ddf(λ, p, s, Cₜ[n]); -// dfC[n] = df(λ, p, s, Cₜ[n]); + dfC[n] = df(λ, p, s, Cₜ[n]); } std::vector<Complex> RddfCt = fft.fourier(RddfC, false); -// std::vector<Complex> dfCt = fft.fourier(dfC, true); + std::vector<Complex> dfCt = fft.fourier(dfC, true); 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); -// Ĉₜ₊₁[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(Ȓₜ₊₁); + + 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<Real> Cₜ₊₁ = fft.inverse(Ĉₜ₊₁); ΔC = 0; |