diff options
| -rw-r--r-- | fourier.cpp | 2 | ||||
| -rw-r--r-- | log-fourier_integrator.cpp | 25 |
2 files changed, 22 insertions, 5 deletions
diff --git a/fourier.cpp b/fourier.cpp index bb1c92f..15cad52 100644 --- a/fourier.cpp +++ b/fourier.cpp @@ -169,7 +169,7 @@ std::vector<Complex> LogarithmicFourierTransform::fourier(const std::vector<Real } fftw_execute(â_to_a); for (unsigned n = 0; n < N; n++) { - ĉ[n] += exp(-k * ω(n)) * a[(pad - 1)*N+n].real() / (Real)(pad*N); + ĉ[n] += exp(-k * ω(n)) * a[(pad - 1)*N+n] / (Real)(pad*N); } } diff --git a/log-fourier_integrator.cpp b/log-fourier_integrator.cpp index 177db46..eda4945 100644 --- a/log-fourier_integrator.cpp +++ b/log-fourier_integrator.cpp @@ -12,7 +12,7 @@ int main(int argc, char* argv[]) { Real Δτ = 0.1; Real k = 0.1; - Real ε = 1e-16; + Real ε = 1e-13; Real γ = 1; Real βₘₐₓ = 0.7; Real Δβ = 0.01; @@ -99,7 +99,8 @@ int main(int argc, char* argv[]) { for (unsigned n = 0; n < N; n++) { Rtnew[n] = (1.0 + pow(β, 2) * RddfCt[n] * Rt[n]) / (μ + 1i * fft.ν(n)); - Ctnew[n] = (2 * Γ₀ * std::conj(Rt[n]) / (1 + pow(τ₀ * fft.ν(n), 2)) + pow(β, 2) * (RddfCt[n] * Ct[n] + dfCt[n] * std::conj(Rt[n]))) / (μ + 1i * fft.ν(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)); +// Ctnew[n] = - 2 * Γ₀ * Rtnew[n].imag() / (1 + pow(τ₀ * fft.ν(n), 2)) / fft.ν(n); } std::vector<Real> Cnew = fft.inverse(Ctnew); @@ -124,12 +125,28 @@ int main(int argc, char* argv[]) { // std::cerr << ΔC << std::endl; } - std::cerr << β << " " << μ << " " << Ct[0].real() << 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; β += Δβ; } for (unsigned i = 0; i < N; i++) { - std::cout << fft.t(i) << " " << Rt[i].imag() << std::endl; + std::cout << fft.t(i) << " " << C[i] << std::endl; } return 0; |