diff options
| -rw-r--r-- | Makefile | 4 | ||||
| -rw-r--r-- | fourier.cpp | 19 | ||||
| -rw-r--r-- | fourier.hpp | 1 | ||||
| -rw-r--r-- | integrator.cpp | 96 |
4 files changed, 92 insertions, 28 deletions
@@ -8,8 +8,8 @@ walk: walk.cpp correlations: correlations.cpp $(CC) correlations.cpp -o correlations -integrator: integrator.cpp - $(CC) integrator.cpp -o integrator +integrator: integrator.cpp fourier.o fourier.hpp + $(CC) integrator.cpp fourier.o -lfftw3 -lfftw3_omp -o integrator fourier.o: fourier.cpp fourier.hpp $(CC) fourier.cpp -c -o fourier.o diff --git a/fourier.cpp b/fourier.cpp index 07f8fdc..bc4b633 100644 --- a/fourier.cpp +++ b/fourier.cpp @@ -64,6 +64,25 @@ std::vector<Complex> FourierTransform::fourier() { return ĉ; } +std::vector<Real> FourierTransform::convolve(const std::vector<Real>& Γh, const std::vector<Real>& R) { + a[0] = 0; + for (unsigned i = 1; i < n; i++) { + a[i] = R[i]; + a[2 * n - i] = -R[i]; + } + fftw_execute(plan_r2c); + for (unsigned i = 1; i < n + 1; i++) { + â[i] *= Γh[i] * (Δτ * M_PI); + } + fftw_execute(plan_c2r); + std::vector<Real> dC(n); + for (unsigned i = 0; i < n; i++) { + dC[i] = a[i] * (Δω / (2 * M_PI)); + } + + return dC; +} + std::vector<Real> FourierTransform::inverse(const std::vector<Complex>& ĉ) { for (unsigned i = 0; i < n + 1; i++) { â[i] = ĉ[i]; diff --git a/fourier.hpp b/fourier.hpp index 9451f69..791953b 100644 --- a/fourier.hpp +++ b/fourier.hpp @@ -33,6 +33,7 @@ public: std::vector<Complex> fourier(); std::vector<Real> inverse(const std::vector<Complex>& ĉ); void writeToA(unsigned i, Real ai); + std::vector<Real> convolve(const std::vector<Real>& Γh, const std::vector<Real>& R); }; std::string fourierFile(std::string prefix, unsigned p, unsigned s, Real λ, Real τ₀, Real y, unsigned log2n, Real τₘₐₓ); diff --git a/integrator.cpp b/integrator.cpp index faba06d..7e5c512 100644 --- a/integrator.cpp +++ b/integrator.cpp @@ -1,4 +1,6 @@ #include "fourier.hpp" +#include <fstream> +#include <fftw3.h> #include <getopt.h> #include <iostream> @@ -14,27 +16,36 @@ int main(int argc, char* argv[]) { unsigned p = 3; unsigned s = 4; Real λ = 0.5; - Real Δτ = 1e-3; Real τₘₐₓ = 1e3; Real τ₀ = 0; - Real β = 0.5; + Real β₀ = 0; + Real βₘₐₓ = 1; + Real Δβ = 1e-2; unsigned iterations = 10; + unsigned log2n = 8; + Real ε = 1e-14; int opt; - while ((opt = getopt(argc, argv, "d:T:t:b:I:")) != -1) { + while ((opt = getopt(argc, argv, "T:2:t:0:b:d:I:")) != -1) { switch (opt) { - case 'd': - Δτ = atof(optarg); - break; case 'T': τₘₐₓ = atof(optarg); break; + case '2': + log2n = atof(optarg); + break; case 't': τ₀ = atof(optarg); break; + case '0': + β₀ = atof(optarg); + break; case 'b': - β = atof(optarg); + βₘₐₓ = atof(optarg); + break; + case 'd': + Δβ = atof(optarg); break; case 'I': iterations = (unsigned)atof(optarg); @@ -44,6 +55,13 @@ int main(int argc, char* argv[]) { } } + unsigned N = pow(2, log2n); + + Real Δτ = (1 + τ₀ / 2) * τₘₐₓ / M_PI / N; + Real Δω = M_PI / ((1 + τ₀ / 2) * τₘₐₓ); + + FourierTransform fft(N, Δω, Δτ, FFTW_ESTIMATE); + Real Γ₀ = 1; Real μ = 1; if (τ₀ > 0) { @@ -51,45 +69,71 @@ int main(int argc, char* argv[]) { } Real τ = 0; - unsigned N = τₘₐₓ / Δτ + 1; std::vector<Real> C(N); std::vector<Real> R(N); std::vector<Real> Γ(N); + std::vector<Real> Γh(N+1); + + Γh[0] = Γ₀; for (unsigned i = 0; i < N; i++) { Real τ = i * Δτ; + Real ω = (i + 1) * Δω * M_PI; if (τ₀ > 0) { C[i] = (Γ₀ / μ) * (exp(-μ * τ) - μ * τ₀ * exp(-τ / τ₀)) / (1 - pow(μ * τ₀, 2)); Γ[i] = (Γ₀ / τ₀) * exp(-τ / τ₀); } else { C[i] = (Γ₀ / μ) * exp(-μ * τ); } + Γh[i+1] = Γ₀ / (1 + pow(ω * τ₀, 2)); R[i] = exp(-μ * τ); } - for (unsigned it = 0; it < iterations; it++) { - /* First step: integrate R from C */ - std::vector<Real> R₊(N); - R₊[0] = 1; - for (unsigned i = 1; i < N; i++) { - Real I = 0; - for (unsigned j = 0; j <= i; j++) { - I += R[i - j] * ddf(λ, p, s, C[i - j]) * R[j] * Δτ; + for (Real β = β₀; β < βₘₐₓ; β += Δβ) { + Real Rerr = 100; + while (sqrt(Rerr / N) > ε) { + /* First step: integrate R from C */ + std::vector<Real> R₊(N); + R₊[0] = 1; + for (unsigned i = 1; i < N; i++) { + Real I = 0; + for (unsigned j = 0; j <= i; j++) { + I += R[i - j] * ddf(λ, p, s, C[i - j]) * R[j] * Δτ; + } + Real dR = -μ * R₊[i - 1] + pow(β, 2) * I; + R₊[i] = R₊[i - 1] + dR * Δτ; } - Real dR = -μ * R[i] + pow(β, 2) * I; - R₊[i] = R₊[i - 1] + dR * Δτ; + + Rerr = 0; + for (unsigned i = 0; i < N; i++) { + Rerr += pow(R[i] - R₊[i], 2); + } + + R = R₊; + + /* Second step: integrate C from R */ + std::vector<Real> dC = fft.convolve(Γh, R); + Real Cₜ₊₁ = 0; + for (unsigned i = 0; i < N; i++) { + Real Cₜ = Cₜ₊₁ + dC[N - i - 1] * Δτ; + C[N - i - 1] = Cₜ; + Cₜ₊₁ = Cₜ; + } + + /* Third step: adjust μ */ + μ *= C[0]; + + std::cerr << β << " " << sqrt(Rerr / N) << std::endl; } - /* Second step: integrate C from R */ - } + std::ofstream outfile(fourierFile("Ci", p, s, λ, τ₀, β, log2n, τₘₐₓ), std::ios::out | std::ios::binary); + outfile.write((const char*)(C.data()), N * sizeof(Real)); + outfile.close(); - τ = 0; - for (Real Ci : C) { - std::cout << τ << " " << Ci << std::endl; - τ += Δτ; + std::ofstream outfileR(fourierFile("Ri", p, s, λ, τ₀, β, log2n, τₘₐₓ), std::ios::out | std::ios::binary); + outfileR.write((const char*)(R.data()), N * sizeof(Real)); + outfileR.close(); } - std::cerr << - 2 * β / Γ₀ * energy(C, R, λ, p, s, Δτ) << std::endl; - return 0; } |