diff options
| -rw-r--r-- | .gitignore | 1 | ||||
| -rw-r--r-- | Makefile | 15 | ||||
| -rw-r--r-- | fourier.cpp | 153 | ||||
| -rw-r--r-- | fourier.hpp | 47 | ||||
| -rw-r--r-- | fourier_integrator.cpp | 29 | ||||
| -rw-r--r-- | get_energy.cpp | 15 | ||||
| -rw-r--r-- | integrator.cpp | 247 | ||||
| -rw-r--r-- | log-fourier_integrator.cpp | 114 |
8 files changed, 425 insertions, 196 deletions
@@ -6,3 +6,4 @@ vgcore.* walk get_energy integrator +fftw.wisdom @@ -1,6 +1,6 @@ -all: walk correlations integrator fourier_integrator get_energy +all: walk correlations integrator fourier_integrator get_energy log-fourier_integrator -CC := clang++ -std=c++17 -flto -Wno-mathematical-notation-identifier-extension -O3 -march=native -mtune=native -fopenmp +CC := clang++ -std=c++17 -flto -Wno-mathematical-notation-identifier-extension -O3 -march=native -mtune=native -fopenmp -DFFTW_THREADS=1 walk: walk.cpp $(CC) walk.cpp -o walk @@ -8,14 +8,17 @@ 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 -lgsl -o integrator fourier.o: fourier.cpp fourier.hpp $(CC) fourier.cpp -c -o fourier.o fourier_integrator: fourier_integrator.cpp fourier.hpp fourier.o - $(CC) fourier_integrator.cpp fourier.o -lfftw3 -o fourier_integrator + $(CC) fourier_integrator.cpp fourier.o -lfftw3 -lgsl -o fourier_integrator + +log-fourier_integrator: log-fourier_integrator.cpp fourier.hpp fourier.o + $(CC) log-fourier_integrator.cpp fourier.o -lfftw3 -lgsl -o log-fourier_integrator get_energy: get_energy.cpp fourier.hpp fourier.o - $(CC) get_energy.cpp fourier.o -lfftw3 -o get_energy + $(CC) get_energy.cpp fourier.o -lfftw3 -lgsl -o get_energy diff --git a/fourier.cpp b/fourier.cpp index 8943521..cf45ad1 100644 --- a/fourier.cpp +++ b/fourier.cpp @@ -1,4 +1,5 @@ #include "fourier.hpp" +#include <boost/multiprecision/fwd.hpp> inline Real fP(unsigned p, Real q) { return 0.5 * pow(q, p); @@ -24,37 +25,153 @@ Real ddf(Real λ, unsigned p, unsigned s, Real q) { return (1 - λ) * ddfP(p, q) + λ * ddfP(s, q); } -FourierTransform::FourierTransform(unsigned n, Real Δω, Real Δτ, unsigned flags) : a(2 * n), â(n + 1), Δω(Δω), Δτ(Δτ) { - plan_r2c = fftw_plan_dft_r2c_1d(2 * n, a.data(), reinterpret_cast<fftw_complex*>(â.data()), flags); - plan_c2r = fftw_plan_dft_c2r_1d(2 * n, reinterpret_cast<fftw_complex*>(â.data()), a.data(), flags); +FourierTransform::FourierTransform(unsigned n, Real Δω, Real Δτ, unsigned flags) : n(n), Δω(Δω), Δτ(Δτ) { + a = fftw_alloc_real(2 * n); + â = reinterpret_cast<Complex*>(fftw_alloc_complex(n + 1)); +// fftw_init_threads(); +// fftw_plan_with_nthreads(FFTW_THREADS); + fftw_import_wisdom_from_filename("fftw.wisdom"); + plan_r2c = fftw_plan_dft_r2c_1d(2 * n, a, reinterpret_cast<fftw_complex*>(â), flags); + plan_c2r = fftw_plan_dft_c2r_1d(2 * n, reinterpret_cast<fftw_complex*>(â), a, flags); + fftw_export_wisdom_to_filename("fftw.wisdom"); } FourierTransform::~FourierTransform() { fftw_destroy_plan(plan_r2c); fftw_destroy_plan(plan_c2r); + fftw_free(a); + fftw_free(â); fftw_cleanup(); } std::vector<Complex> FourierTransform::fourier(const std::vector<Real>& c) { - a = c; + for (unsigned i = 0; i < 2 * n; i++) { + a[i] = c[i]; + } fftw_execute(plan_r2c); - std::vector<Complex> ĉ(â.size()); - for (unsigned i = 0; i < â.size(); i++) { + std::vector<Complex> ĉ(n + 1); + for (unsigned i = 0; i < n + 1; i++) { ĉ[i] = â[i] * (Δτ * M_PI); } return ĉ; } +std::vector<Complex> FourierTransform::fourier() { + fftw_execute(plan_r2c); + std::vector<Complex> ĉ(n+1); + for (unsigned i = 0; i < n+1; i++) { + ĉ[i] = â[i] * (Δτ * M_PI); + } + 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]; + } fftw_execute(plan_c2r); - std::vector<Real> c(a.size()); - for (unsigned i = 0; i < a.size(); i++) { + std::vector<Real> c(2*n); + for (unsigned i = 0; i < 2*n; i++) { c[i] = a[i] * (Δω / (2 * M_PI)); } return c; } +void FourierTransform::writeToA(unsigned i, Real ai) { + a[i] = ai; +} + +LogarithmicFourierTransform::LogarithmicFourierTransform(unsigned N, Real k, Real Δω, Real Δτ, Real Δs) : N(N), k(k), Δτ(Δτ), Δω(Δω), Δs(Δs) { + τₛ = -0.5 * N; + ωₛ = -0.5 * N; + sₛ = -0.5 * N; + a = reinterpret_cast<Complex*>(fftw_alloc_complex(N)); + â = reinterpret_cast<Complex*>(fftw_alloc_complex(N)); + fftw_import_wisdom_from_filename("fftw.wisdom"); + a_to_â = fftw_plan_dft_1d(N, reinterpret_cast<fftw_complex*>(a), reinterpret_cast<fftw_complex*>(â), 1, 0); + â_to_a = fftw_plan_dft_1d(N, reinterpret_cast<fftw_complex*>(â), reinterpret_cast<fftw_complex*>(a), 1, 0); + fftw_export_wisdom_to_filename("fftw.wisdom"); +} + +LogarithmicFourierTransform::~LogarithmicFourierTransform() { + fftw_destroy_plan(a_to_â); + fftw_destroy_plan(â_to_a); + fftw_free(a); + fftw_free(â); + fftw_cleanup(); +} + +Real LogarithmicFourierTransform::τ(unsigned n) const { + return Δτ * (n + τₛ); +} + +Real LogarithmicFourierTransform::ω(unsigned n) const { + return Δω * (n + ωₛ); +} + +Real LogarithmicFourierTransform::s(unsigned n) const { + return Δs * (n + sₛ); +} + +Real LogarithmicFourierTransform::t(unsigned n) const { + return exp(τ(n)); +} + +Real LogarithmicFourierTransform::ν(unsigned n) const { + return exp(ω(n)); +} + +Complex gamma(Complex z) { + gsl_sf_result logΓ; + gsl_sf_result argΓ; + + gsl_sf_lngamma_complex_e(z.real(), z.imag(), &logΓ, &argΓ); + + return exp(logΓ.val + 1i * argΓ.val); +} + +std::vector<Complex> LogarithmicFourierTransform::fourier(const std::vector<Real>& c, bool symmetric){ + std::vector<Complex> ĉ(N); + std::vector<Real> σs = {1}; + if (symmetric){ + σs.push_back(-1); + } + for (Real σ : σs) { + for (unsigned n = 0; n < N; n++) { + a[n] = c[n] * exp((1 - k) * τ(n)); + } + fftw_execute(a_to_â); + for (unsigned n = 0; n < N; n++) { + â[(n + N / 2) % N] *= pow(-1i * σ, 1i * (s(n)) - k) * gamma(k - 1i * (s(n))); + } + fftw_execute(â_to_a); + for (unsigned n = 0; n < N; n++) { + ĉ[n] += exp(-k * ω(n) * 4 * M_PI) * a[n] / pow(2 * M_PI, 1); + } + } + + return ĉ; +} + std::string fourierFile(std::string prefix, unsigned p, unsigned s, Real λ, Real τ₀, Real y, unsigned log2n, Real τₘₐₓ) { return prefix + "_" + std::to_string(p) + "_" + std::to_string(s) + "_" + std::to_string(λ) + "_" + std::to_string(τ₀) + "_" + std::to_string(y) + "_" + std::to_string(log2n) + "_" + std::to_string(τₘₐₓ) + ".dat"; } @@ -70,23 +187,25 @@ Real energy(const std::vector<Real>& C, const std::vector<Real>& R, unsigned p, } std::tuple<std::vector<Complex>, std::vector<Complex>> RddfCtdfCt(FourierTransform& fft, const std::vector<Real>& C, const std::vector<Real>& R, unsigned p, unsigned s, Real λ) { - std::vector<Real> RddfC(C.size()); for (unsigned i = 0; i < C.size() / 2; i++) { - RddfC[i] = R[i] * ddf(λ, p, s, C[i]); + fft.writeToA(i, R[i] * ddf(λ, p, s, C[i])); + } + for (unsigned i = C.size() / 2; i < C.size(); i++) { + fft.writeToA(i, 0); } - std::vector<Complex> RddfCt = fft.fourier(RddfC); + std::vector<Complex> RddfCt = fft.fourier(); - std::vector<Real> dfC(C.size()); for (unsigned i = 0; i < C.size(); i++) { - dfC[i] = df(λ, p, s, C[i]); + fft.writeToA(i, df(λ, p, s, C[i])); } - std::vector<Complex> dfCt = fft.fourier(dfC); + std::vector<Complex> dfCt = fft.fourier(); return {RddfCt, dfCt}; } -Real estimateZ(FourierTransform& 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 y) { +Real estimateZ(FourierTransform& 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 y) { auto [RddfCt, dfCt] = RddfCtdfCt(fft, C, R, p, s, λ); + Real Γ₀ = 1 + τ₀ / 2; - return ((std::conj(Rt[0]) + pow(y, 2) * (RddfCt[0] * Ct[0] + dfCt[0] * std::conj(Rt[0]))) / Ct[0]).real(); + return ((Γ₀ * std::conj(Rt[0]) + pow(y, 2) * (RddfCt[0] * Ct[0] + dfCt[0] * std::conj(Rt[0]))) / Ct[0]).real(); } diff --git a/fourier.hpp b/fourier.hpp index 3925058..22a57e7 100644 --- a/fourier.hpp +++ b/fourier.hpp @@ -3,6 +3,12 @@ #include <fftw3.h> #include <complex> #include <vector> +#include <tuple> +#include <gsl/gsl_sf_gamma.h> + +#ifndef FFTW_THREADS +#define FFTW_THREADS 1 +#endif using Real = double; using Complex = std::complex<Real>; @@ -14,20 +20,55 @@ Real ddf(Real λ, unsigned p, unsigned s, Real q); class FourierTransform { private: - std::vector<Real> a; - std::vector<Complex> â; + Real* a; + Complex* â; fftw_plan plan_r2c; fftw_plan plan_c2r; + unsigned n; Real Δω; Real Δτ; public: FourierTransform(unsigned n, Real Δω, Real Δτ, unsigned flags = 0); ~FourierTransform(); std::vector<Complex> fourier(const std::vector<Real>& c); + 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); }; +class LogarithmicFourierTransform { +private: + Complex* a; + Complex* â; + fftw_plan a_to_â; + fftw_plan â_to_a; + unsigned N; + Real k; + Real Δτ; + Real Δω; + Real Δs; + Real τₛ; + Real ωₛ; + Real sₛ; +public: + LogarithmicFourierTransform(unsigned N, Real k, Real Δω, Real Δτ, Real Δs); + ~LogarithmicFourierTransform(); + Real τ(unsigned n) const; + Real ω(unsigned n) const; + Real t(unsigned n) const; + Real ν(unsigned n) const; + Real s(unsigned n) const; + std::vector<Complex> fourier(const std::vector<Real>& c, bool symmetric); + 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); +}; + +Complex gamma(Complex z); + std::string fourierFile(std::string prefix, unsigned p, unsigned s, Real λ, Real τ₀, Real y, unsigned log2n, Real τₘₐₓ); Real energy(const std::vector<Real>& C, const std::vector<Real>& R, unsigned p, unsigned s, Real λ, Real y, Real Δτ); std::tuple<std::vector<Complex>, std::vector<Complex>> RddfCtdfCt(FourierTransform& fft, const std::vector<Real>& C, const std::vector<Real>& R, unsigned p, unsigned s, Real λ); -Real estimateZ(FourierTransform& 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 y); +Real estimateZ(FourierTransform& 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 y); diff --git a/fourier_integrator.cpp b/fourier_integrator.cpp index d83fb9d..4f4eff5 100644 --- a/fourier_integrator.cpp +++ b/fourier_integrator.cpp @@ -17,7 +17,7 @@ int main(int argc, char* argv[]) { unsigned maxIterations = 1000; Real ε = 1e-14; - Real γ = 0; + Real γ = 1; bool loadData = false; @@ -65,11 +65,11 @@ int main(int argc, char* argv[]) { unsigned n = pow(2, log2n); - Real Δτ = τₘₐₓ / M_PI / n; - Real Δω = M_PI / τₘₐₓ; + Real Δτ = (1 + τ₀ / 2) * τₘₐₓ / M_PI / n; + Real Δω = M_PI / ((1 + τ₀ / 2) * τₘₐₓ); - Real z = (-1+sqrt(1+2*τ₀)) / (2 * τ₀); - Real Γ₀ = 1; + Real z = 0.5; + Real Γ₀ = 1 + τ₀ / 2; std::vector<Real> C(2 * n); std::vector<Real> R(2 * n); @@ -84,7 +84,11 @@ int main(int argc, char* argv[]) { // start from the exact solution for τ₀ = 0 for (unsigned i = 0; i < n; i++) { Real τ = i * Δτ * M_PI; - C[i] = Γ₀ / 2 * (exp(-z * τ) - z * τ₀ * exp(-τ / τ₀)) / (z - pow(z, 3) * pow(τ₀, 2)); + if (τ₀ > 0) { + C[i] = Γ₀ / 2 * (exp(-z * τ) - z * τ₀ * exp(-τ / τ₀)) / (z - pow(z, 3) * pow(τ₀, 2)); + } else { + C[i] = Γ₀ / 2 * exp(-z * τ) / z; + } if (i > 0) { C[2 * n - i] = C[i]; } @@ -94,16 +98,19 @@ int main(int argc, char* argv[]) { Rt = fft.fourier(R); } else { std::ifstream cfile(fourierFile("C", p, s, λ, τ₀, y, log2n, τₘₐₓ), std::ios::binary); - cfile.read((char*)(C.data()), C.size() * sizeof(Real)); + cfile.read((char*)(C.data()), (C.size() / 2) * sizeof(Real)); cfile.close(); + for (unsigned i = 1; i < n; i++) { + C[2 * n - i] = C[i]; + } std::ifstream rfile(fourierFile("R", p, s, λ, τ₀, y, log2n, τₘₐₓ), std::ios::binary); - rfile.read((char*)(R.data()), R.size() * sizeof(Real)); + rfile.read((char*)(R.data()), (R.size() / 2) * sizeof(Real)); rfile.close(); Ct = fft.fourier(C); Rt = fft.fourier(R); - z = estimateZ(fft, C, Ct, R, Rt, p, s, λ, y); + z = estimateZ(fft, C, Ct, R, Rt, p, s, λ, τ₀, y); } while (y += Δy, y <= yₘₐₓ) { @@ -165,11 +172,11 @@ int main(int argc, char* argv[]) { std::cerr << "y " << y << " " << e << " " << z << std::endl; std::ofstream outfile(fourierFile("C", p, s, λ, τ₀, y, log2n, τₘₐₓ), std::ios::out | std::ios::binary); - outfile.write((const char*)(C.data()), C.size() * sizeof(Real)); + outfile.write((const char*)(C.data()), (C.size() / 2) * sizeof(Real)); outfile.close(); std::ofstream outfileR(fourierFile("R", p, s, λ, τ₀, y, log2n, τₘₐₓ), std::ios::out | std::ios::binary); - outfileR.write((const char*)(R.data()), R.size() * sizeof(Real)); + outfileR.write((const char*)(R.data()), (R.size() / 2) * sizeof(Real)); outfileR.close(); } diff --git a/get_energy.cpp b/get_energy.cpp index b6e835f..c33c04e 100644 --- a/get_energy.cpp +++ b/get_energy.cpp @@ -50,12 +50,12 @@ int main(int argc, char* argv[]) { unsigned n = pow(2, log2n); - Real Δτ = τₘₐₓ / M_PI / n; - Real Δω = M_PI / τₘₐₓ; + Real Δτ = (1 + τ₀ / 2) * τₘₐₓ / M_PI / n; + Real Δω = M_PI / ((1 + τ₀ / 2) * τₘₐₓ); Real y = y₀; - FourierTransform fft(n, Δω, Δτ, FFTW_ESTIMATE); + FourierTransform fft(n, Δω, Δτ); std::vector<Real> C(2 * n); std::vector<Real> R(2 * n); @@ -63,11 +63,14 @@ int main(int argc, char* argv[]) { while (y += Δy, y <= yₘₐₓ) { std::ifstream cfile(fourierFile("C", p, s, λ, τ₀, y, log2n, τₘₐₓ), std::ios::binary); if (cfile.is_open()) { - cfile.read((char*)(C.data()), C.size() * sizeof(Real)); + cfile.read((char*)(C.data()), (C.size() / 2) * sizeof(Real)); cfile.close(); + for (unsigned i = 1; i < n; i++) { + C[2 * n - i] = C[i]; + } std::ifstream rfile(fourierFile("R", p, s, λ, τ₀, y, log2n, τₘₐₓ), std::ios::binary); - rfile.read((char*)(R.data()), R.size() * sizeof(Real)); + rfile.read((char*)(R.data()), (R.size() / 2) * sizeof(Real)); rfile.close(); Real e = energy(C, R, p, s, λ, y, Δτ); @@ -75,7 +78,7 @@ int main(int argc, char* argv[]) { std::vector<Complex> Ct = fft.fourier(C); std::vector<Complex> Rt = fft.fourier(R); - Real z = estimateZ(fft, C, Ct, R, Rt, p, s, λ, y); + Real z = estimateZ(fft, C, Ct, R, Rt, p, s, λ, τ₀, y); std::cout << y << " " << e << " " << Ct[0].real() << " " << z << std::endl; } diff --git a/integrator.cpp b/integrator.cpp index 8b13a9b..7e5c512 100644 --- a/integrator.cpp +++ b/integrator.cpp @@ -1,124 +1,51 @@ +#include "fourier.hpp" +#include <fstream> +#include <fftw3.h> #include <getopt.h> -#include <vector> -#include <cmath> #include <iostream> -using Real = double; - -unsigned p = 2; - -Real f(Real q) { - return 0.5 * pow(q, p); -} - -Real df(Real q) { - return 0.5 * p * pow(q, p - 1); -} - -Real ddf(Real q) { - return 0.5 * p * (p - 1) * pow(q, p - 2); -} - -Real integrate(const std::vector<Real>& C, signed τ = std::numeric_limits<unsigned>::max()) { +Real energy(const std::vector<Real>& C, std::vector<Real>& R, Real λ, unsigned p, unsigned s, Real Δτ) { Real I = 0; - if (τ > C.size() - 1) { - τ = C.size() - 1; - } -#pragma omp parallel for reduction(+:I) - for (unsigned σ = 0; σ < τ; σ++) { - unsigned τ_σ = τ - σ; - Real Cτ_σ = (C[τ_σ] + C[τ_σ - 1]) / 2; - Real dCσ = C[σ + 1] - C[σ]; - - I += df(Cτ_σ) * dCσ; - } - return I; -} - -Real integratePast(const std::vector<Real>& C, signed τ) { - Real I = 0; -#pragma omp parallel for reduction(+:I) - for (signed σ = -C.size() + τ + 3; σ < τ - 2; σ++) { - signed τ_σ = τ - σ; - - Real Cτ_σ = (C[abs(τ_σ)] + C[abs(τ_σ) - 1]) / 2; - Real Cσ = (C[abs(σ) + 1] + C[abs(σ)]) / 2; - Real dddC; - if (τ_σ != 0) { - dddC = (τ_σ / abs(τ_σ)) * (C[abs(τ_σ)+2] - 2 * C[abs(τ_σ)+1] + 2 * C[abs(τ_σ)-1] - C[abs(τ_σ)-2]) / 2; - } else { - dddC = 0; - } - - I += dddC * ddf(Cτ_σ) * Cσ; - } -#pragma omp parallel for reduction(+:I) - for (signed σ = -C.size() + τ + 3; σ < -1; σ++) { - signed τ_σ = τ - σ; - - Real Cτ_σ = (C[abs(τ_σ)] + C[abs(τ_σ) - 1]) / 2; - Real dddC; - if (σ != 0) { - dddC = -(σ / abs(σ)) * (C[abs(σ)+2] - 2 * C[abs(σ)+1] + 2 * C[abs(σ)-1] - C[abs(σ)-2]) / 2; - } else { - dddC = 0; - } - - I += dddC * df(Cτ_σ); - } - return I; -} - -Real integrateDelay(const std::vector<Real>& C, unsigned τ, Real Δτ, Real τ₀) { - Real I = 0; -#pragma omp parallel for reduction(+:I) - for (signed σ = 2; σ < C.size() - τ - 2; σ++) { - unsigned τ_σ = τ + σ; - Real dC = -(C[σ+1] - C[σ-1]) / 2; - Real dddC = -(C[σ+2] - 2 * C[σ+1] + 2 * C[σ-1] - C[σ-2]) / 2; - - I += (dC - pow(τ₀ / Δτ, 2) * dddC) * exp(-(τ_σ * Δτ / τ₀)); - } - return I / τ₀; -} - -Real energy(const std::vector<Real>& C, Real Δτ, Real τ₀) { - Real I = 0; - for (unsigned σ = 0; σ < C.size() - 1; σ++) { - Real Cσ = (C[σ] + C[σ + 1]) / 2; - Real dC = (C[σ + 1] - C[σ]) / Δτ; - - Real dddC = 0; - if (σ > 1 && σ < C.size() - 2 && C.size() > 3) { - dddC = (C[σ+1] - 3 * C[σ] + 3 * C[σ-1] - C[σ-2]) / pow(Δτ, 3); - } - I += Δτ * df(Cσ) * (dC - pow(τ₀, 2) * dddC); + for (unsigned σ = 0; σ < C.size(); σ++) { + I += Δτ * df(λ, p, s, C[σ]) * R[σ]; } return I; } int main(int argc, char* argv[]) { - Real Δτ = 1e-3; + unsigned p = 3; + unsigned s = 4; + Real λ = 0.5; Real τₘₐₓ = 1e3; Real τ₀ = 0; - Real y = 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:y: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 'y': - y = atof(optarg); + case '0': + β₀ = atof(optarg); + break; + case 'b': + βₘₐₓ = atof(optarg); + break; + case 'd': + Δβ = atof(optarg); break; case 'I': iterations = (unsigned)atof(optarg); @@ -128,71 +55,85 @@ int main(int argc, char* argv[]) { } } - Real z = 0.4794707565634420155347; - Real Γ₀ = 1; + unsigned N = pow(2, log2n); - Real τ = 0; - std::vector<Real> C; - C.reserve(τₘₐₓ / Δτ + 1); + Real Δτ = (1 + τ₀ / 2) * τₘₐₓ / M_PI / N; + Real Δω = M_PI / ((1 + τ₀ / 2) * τₘₐₓ); - C.push_back(1); + FourierTransform fft(N, Δω, Δτ, FFTW_ESTIMATE); -// while (std::cout << τ << " " << C.back() << std::endl, τ < τₘₐₓ) { - while (τ < τₘₐₓ) { - τ += Δτ; - Real dC = -(z - 2 * pow(y, 2)) * C.back() - 2 / Γ₀ * pow(y, 2) * integrate(C); - C.push_back(C.back() + Δτ * dC); + Real Γ₀ = 1; + Real μ = 1; + if (τ₀ > 0) { + μ = (sqrt(1+4*Γ₀*τ₀) - 1) / (2 * τ₀); } - - for (unsigned it = 0; it < iterations; it++) { - - τ = 0; - std::vector<Real> C2; - C2.reserve(τₘₐₓ / Δτ + 1); - C2.push_back(1); - while (τ < τₘₐₓ) { - τ += Δτ; - Real dC = -(z - 2 * pow(y, 2)) * C2.back() + integrateDelay(C, C2.size() - 1, Δτ, τ₀) - 2 / Γ₀ * pow(y, 2) * (integrate(C2) - pow(τ₀ / Δτ, 2) * integratePast(C, C2.size()-1)); - C2.push_back(C2.back() + Δτ * dC); + Real τ = 0; + 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(-μ * τ); + } - Real error = 0; - - for (unsigned i = 0; i < std::min(C.size(), C2.size()); i++) { - error += pow(C[i] - C2[i], 2); + 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 * Δτ; + } + + 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; } - std::cerr << "Iteration " << it << ": " << sqrt(error / C.size()) << " " << z << std::endl; - - C = C2; - } - /* - Real zNew = (2.0 * ((C[2] - 2 * C[1] + C[0]) / pow(Δτ, 2) - pow(τ₀, 2) * (C[4] - 4 * C[3] + 6 * C[2] - 4 * C[1] + C[0]) / pow(Δτ, 4))); - Real zNew = (2.0 * ((C[2] - 2 * C[1] + C[0]) / pow(Δτ, 2) - pow(τ₀, 2) * (C[4] - 4 * C[3] + 6 * C[2] - 4 * C[1] + C[0]) / pow(Δτ, 4))); -// Real zNew = (2.0 * ((83 * C[6] - 245 * C[5] + 101 * C[4] + 254 * C[3] - 31 * C[2] - 377 * C[1] + 215 * C[0]) / (132 * pow(Δτ, 2)) - pow(τ₀, 2) * (3 * C[6] - 7 * C[5] + C[4] + 6 * C[3] + C[2] - 7 * C[1] + 3 * C[0]) / (11 * pow(Δτ, 4)))); - z = z / zNew; - τ = 0; - C.clear(); - C.reserve(τₘₐₓ / Δτ + 1); - - C.push_back(1); - -// while (std::cout << τ << " " << C.back() << std::endl, τ < τₘₐₓ) { - while (τ < τₘₐₓ) { - τ += Δτ; - Real dC = -z * C.back() - 2 / Γ₀ * pow(y, 2) * integrate(C); - C.push_back(C.back() + Δτ * dC); - } - */ + 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 * y / Γ₀ * energy(C, Δτ, τ₀) << std::endl; - return 0; } diff --git a/log-fourier_integrator.cpp b/log-fourier_integrator.cpp new file mode 100644 index 0000000..247270b --- /dev/null +++ b/log-fourier_integrator.cpp @@ -0,0 +1,114 @@ +#include "fourier.hpp" +#include <getopt.h> +#include <iostream> +#include <fstream> + +int main(int argc, char* argv[]) { + unsigned p = 2; + unsigned s = 2; + Real λ = 0.5; + Real τ₀ = 0; + Real β₀ = 0; + Real yₘₐₓ = 0.5; + Real Δy = 0.05; + + unsigned log2n = 8; + Real τₘₐₓ = 20; + + unsigned maxIterations = 1000; + Real ε = 1e-14; + Real γ = 1; + + int opt; + + while ((opt = getopt(argc, argv, "p:s:2:T:t:0:y:d:I:g:l")) != -1) { + switch (opt) { + case 'p': + p = atoi(optarg); + break; + case 's': + s = atoi(optarg); + break; + case '2': + log2n = atoi(optarg); + break; + case 'T': + τₘₐₓ = atof(optarg); + break; + case 't': + τ₀ = atof(optarg); + break; + case '0': + β₀ = atof(optarg); + break; + case 'y': + yₘₐₓ = atof(optarg); + break; + case 'd': + Δy = atof(optarg); + break; + case 'I': + maxIterations = (unsigned)atof(optarg); + break; + case 'g': + γ = atof(optarg); + break; + default: + exit(1); + } + } + + unsigned N = pow(2, log2n); + + Real Δτ = 1e-2; + Real Δω = 1e-2; + Real Δs = 1e-2; + Real k = 0.1; + + Real Γ₀ = 1.0; + Real μ = Γ₀; + if (τ₀ > 0) { + μ = (sqrt(1+4*Γ₀*τ₀)-1)/(2*τ₀); + } + + std::vector<Real> C(N); + std::vector<Real> R(N); + + LogarithmicFourierTransform fft(N, k, Δω, Δτ, Δs); + + for (unsigned n = 0; n < N; n++) { + if (τ₀ > 0) { + C[n] = Γ₀ * (exp(-μ * fft.t(n)) - μ * τ₀ * exp(-fft.t(n) / τ₀)) / (μ - pow(μ, 3) * pow(τ₀, 2)); + } else { + C[n] = Γ₀ * exp(-μ * fft.t(n)) / μ; + } + R[n] = exp(-μ * fft.t(n)); + } + + std::vector<Complex> Ct = fft.fourier(C, true); + std::vector<Complex> Rt = fft.fourier(R, false); + + /* + for (unsigned n = 0; n < N; n++) { + std::cout << fft.t(n) << " " << C[n] << std::endl; + } + */ + + for (unsigned n = 0; n < N; n++) { + std::cout << fft.ν(n) << " " << Ct[n].real() << std::endl; + } + + /* + // start from the exact solution for τ₀ = 0 + for (unsigned i = 0; i < N + 1; i++) { + Real ω = i * Δω; + Ct[i] = 2 * Γ₀ / (pow(μ, 2) + pow(ω, 2)) / (1 + pow(τ₀ * ω, 2)); + Rt[i] = 1.0 / (μ + 1i * ω); + Γt[i] = Γ₀ / (1 + pow(τ₀ * ω, 2)); + } + C = fft.inverse(Ct); + R = fft.inverse(Rt); + */ + + return 0; +} |