From d93dc5dba5d9e00298038c3dac3b51aa554d2ca8 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Sat, 10 May 2025 13:46:24 -0300 Subject: New time scaling --- log-fourier.cpp | 29 +++++++++++++++-------------- log-fourier.hpp | 7 ++++--- log-fourier_integrator.cpp | 39 ++++++++++++++++++++------------------- log_get_energy.cpp | 12 +++++++----- 4 files changed, 46 insertions(+), 41 deletions(-) diff --git a/log-fourier.cpp b/log-fourier.cpp index e60439a..d5ac17d 100644 --- a/log-fourier.cpp +++ b/log-fourier.cpp @@ -4,9 +4,9 @@ #include #include -LogarithmicFourierTransform::LogarithmicFourierTransform(unsigned N, Real k, Real Δτ, unsigned pad, Real shift) : N(N), pad(pad), k(k), Δτ(Δτ) { - τₛ = -shift * N; - ωₛ = -(1-shift) * N; +LogarithmicFourierTransform::LogarithmicFourierTransform(unsigned N, Real k, Real Δτ, unsigned pad, Real shift) : N(N), pad(pad), k(k), Δτ(Δτ), shift(shift) { + τₛ = -0.5 * N; + ωₛ = -0.5 * N; sₛ = -0.5 * pad * N; a = reinterpret_cast(FFTW_ALLOC_COMPLEX(pad*N)); â = reinterpret_cast(FFTW_ALLOC_COMPLEX(pad*N)); @@ -37,11 +37,11 @@ Real LogarithmicFourierTransform::s(unsigned n) const { } Real LogarithmicFourierTransform::t(unsigned n) const { - return exp(τ(n)); + return std::exp(τ(n)) / shift; } Real LogarithmicFourierTransform::ν(unsigned n) const { - return exp(ω(n)); + return std::exp(ω(n)) * shift; } Complex Γ(Complex z) { @@ -50,7 +50,7 @@ Complex Γ(Complex z) { gsl_sf_lngamma_complex_e((double)z.real(), (double)z.imag(), &logΓ, &argΓ); - return exp((Real)logΓ.val + II * (Real)argΓ.val); + return std::exp((Real)logΓ.val + II * (Real)argΓ.val); } std::vector LogarithmicFourierTransform::fourier(const std::vector& c, bool symmetric) { @@ -63,16 +63,16 @@ std::vector LogarithmicFourierTransform::fourier(const std::vector= (pad - 1) * N) { - a[n] = c[pad*N-n-1] * exp((1 - k) * τ(pad*N-n-1)); + a[n] = c[pad*N-n-1] * std::exp((1 - k) * τ(pad*N-n-1)); } else { a[n] = 0; } } FFTW_EXECUTE(a_to_â); for (unsigned n = 0; n < pad*N; n++) { - â[(pad*N / 2 + n) % (pad*N)] *= std::exp(II*(0.5 * N + τₛ) * s(n) / Δτ) * std::pow(II * σ, II * s(n) - k) * Γ(k - II * s(n)); + â[(pad*N / 2 + n) % (pad*N)] *= std::pow(II * σ, II * s(n) - k) * Γ(k - II * s(n)); } FFTW_EXECUTE(â_to_a); for (unsigned n = 0; n < N; n++) { @@ -82,6 +82,8 @@ std::vector LogarithmicFourierTransform::fourier(const std::vector LogarithmicFourierTransform::inverse(const std::vector= (pad - 1) * N) { - a[n] = (ĉ[pad*N-n-1].real() + II * σ * ĉ[pad*N-n-1].imag()) * std::exp((1 - k) * ω(pad*N-n-1)); + a[n] = (ĉ[pad*N-n-1].real() - II * σ * ĉ[pad*N-n-1].imag()) * std::exp((1 - k) * ω(pad*N-n-1)); } else { a[n] = 0; } } FFTW_EXECUTE(a_to_â); for (unsigned n = 0; n < pad*N; n++) { - â[(pad*N / 2 + n) % (pad*N)] *= std::exp(-II*(0.5 * N + τₛ) * s(n) / Δτ) * std::pow(-II * σ, II * s(n) - k) * Γ(k - II * s(n)); + â[(pad*N / 2 + n) % (pad*N)] *= std::pow(-II * σ, II * s(n) - k) * Γ(k - II * s(n)); } FFTW_EXECUTE(â_to_a); for (unsigned n = 0; n < N; n++) { @@ -112,6 +114,7 @@ std::vector LogarithmicFourierTransform::inverse(const std::vector& C, con Real energy(const LogarithmicFourierTransform& fft, std::vector& C, const std::vector& R, unsigned p, unsigned s, Real λ, Real β) { unsigned n₀ = 0; - /* for (unsigned n = 0; n < C.size(); n++) { if (C[n] > 1 || R[n] > 1) n₀ = n % 2 == 0 ? n / 2 : (n + 1) / 2; } - */ Real E = fft.t(2*n₀) * df(λ, p, s, 1); for (unsigned n = n₀; n < C.size()/2-1; n++) { Real R₂ₙ = R[2*n]; @@ -203,7 +204,7 @@ Real energy(const LogarithmicFourierTransform& fft, std::vector& C, const Real C₂ₙ₊₁ = C[2*n+1]; Real C₂ₙ₊₂ = C[2*n+2]; - //if (C₂ₙ₊₂ < 0 || R₂ₙ₊₂ < 0) break; + if (C₂ₙ₊₂ < 0 || R₂ₙ₊₂ < 0) break; Real h₂ₙ = fft.t(2*n+1) - fft.t(2*n); Real h₂ₙ₊₁ = fft.t(2*n+2) - fft.t(2*n+1); diff --git a/log-fourier.hpp b/log-fourier.hpp index b1e4bd1..b5bb4c0 100644 --- a/log-fourier.hpp +++ b/log-fourier.hpp @@ -20,6 +20,7 @@ private: Real τₛ; Real ωₛ; Real sₛ; + Real shift; public: LogarithmicFourierTransform(unsigned N, Real k, Real Δτ, unsigned pad = 4, Real shift = 0.5); ~LogarithmicFourierTransform(); @@ -32,11 +33,11 @@ public: std::vector inverse(const std::vector& ĉ); }; -std::string logFourierFile(std::string prefix, unsigned p, unsigned s, Real λ, Real τ₀, Real β, unsigned log2n, Real Δτ, Real k); +std::string logFourierFile(std::string prefix, unsigned p, unsigned s, Real λ, Real τ₀, Real β, unsigned log2n, Real Δτ, Real shift); -void logFourierSave(const std::vector& C, const std::vector& R, const std::vector& Ct, const std::vector& Rt, unsigned p, unsigned s, Real λ, Real τ₀, Real β, unsigned log2n, Real Δτ, Real k); +void logFourierSave(const std::vector& C, const std::vector& R, const std::vector& Ct, const std::vector& Rt, unsigned p, unsigned s, Real λ, Real τ₀, Real β, unsigned log2n, Real Δτ, Real shift); -bool logFourierLoad(std::vector& C, std::vector& R, std::vector& Ct, std::vector& Rt, unsigned p, unsigned s, Real λ, Real τ₀, Real β, unsigned log2n, Real Δτ, Real k); +bool logFourierLoad(std::vector& C, std::vector& R, std::vector& Ct, std::vector& Rt, unsigned p, unsigned s, Real λ, Real τ₀, Real β, unsigned log2n, Real Δτ, Real shift); std::tuple, std::vector> RddfCtdfCt(LogarithmicFourierTransform& fft, const std::vector& C, const std::vector& R, unsigned p, unsigned s, Real λ); diff --git a/log-fourier_integrator.cpp b/log-fourier_integrator.cpp index 7b24a55..9db21f0 100644 --- a/log-fourier_integrator.cpp +++ b/log-fourier_integrator.cpp @@ -13,7 +13,7 @@ int main(int argc, char* argv[]) { unsigned log2n = 8; Real Δτ = 0.1; Real k = -0.01; - Real shift = 0.5; + Real logShift = 0; /* Iteration parameters */ Real ε = 1e-15; @@ -54,12 +54,12 @@ int main(int argc, char* argv[]) { case 'k': k = atof(optarg); break; + case 'h': + logShift = atof(optarg); + break; case 'D': Δτ = atof(optarg); break; - case 'h': - shift = atof(optarg); - break; case 'e': ε = atof(optarg); break; @@ -85,13 +85,14 @@ int main(int argc, char* argv[]) { unsigned N = pow(2, log2n); - LogarithmicFourierTransform fft(N, k, Δτ, pad, shift); - Real Γ₀ = 1; - Real μₜ₋₁ = Γ₀; - if (τ₀ > 0) { - μₜ₋₁ = (sqrt(1+4*Γ₀*τ₀)-1)/(2*τ₀); - } + Real μ₀ = τ₀ > 0 ? (sqrt(1+4*Γ₀*τ₀)-1)/(2*τ₀) : Γ₀; + + LogarithmicFourierTransform fft(N, k, Δτ, pad, μ₀ * pow(10, logShift)); + + std::cerr << "Starting, μ₀ = " << μ₀ << ", range " << fft.t(0) << " " << fft.t(N-1) << std::endl; + + Real μₜ₋₁ = μ₀; std::vector Cₜ₋₁(N); std::vector Rₜ₋₁(N); @@ -103,20 +104,20 @@ int main(int argc, char* argv[]) { for (unsigned n = 0; n < N; n++) { if (τ₀ > 0) { if (τ₀ == 2) { - Cₜ₋₁[n] = Γ₀ * exp(-fft.t(n) / 2) * (1 + fft.t(n) / 2); + Cₜ₋₁[n] = Γ₀ * std::exp(-fft.t(n) / 2) * (1 + fft.t(n) / 2); } else { - Cₜ₋₁[n] = Γ₀ * (exp(-μₜ₋₁ * fft.t(n)) - μₜ₋₁ * τ₀ * exp(-fft.t(n) / τ₀)) / (μₜ₋₁ - pow(μₜ₋₁, 3) * pow(τ₀, 2)); + Cₜ₋₁[n] = Γ₀ * (std::exp(-μ₀ * fft.t(n)) - μ₀ * τ₀ * std::exp(-fft.t(n) / τ₀)) / (μ₀ - pow(μ₀, 3) * pow(τ₀, 2)); } } else { - Cₜ₋₁[n] = Γ₀ * exp(-μₜ₋₁ * fft.t(n)) / μₜ₋₁; + Cₜ₋₁[n] = Γ₀ * std::exp(-μ₀ * fft.t(n)) / μ₀; } - Rₜ₋₁[n] = exp(-μₜ₋₁ * fft.t(n)); + Rₜ₋₁[n] = std::exp(-μ₀ * fft.t(n)); - Ĉₜ₋₁[n] = 2 * Γ₀ / (pow(μₜ₋₁, 2) + pow(fft.ν(n), 2)) / (1 + pow(τ₀ * fft.ν(n), 2)); - Ȓₜ₋₁[n] = (Real)1.0 / (μₜ₋₁ + II * fft.ν(n)); + Ĉₜ₋₁[n] = 2 * Γ₀ / (pow(μ₀, 2) + pow(fft.ν(n), 2)) / (1 + pow(τ₀ * fft.ν(n), 2)); + Ȓₜ₋₁[n] = (Real)1.0 / (μ₀ + II * fft.ν(n)); } } else { - logFourierLoad(Cₜ₋₁, Rₜ₋₁, Ĉₜ₋₁, Ȓₜ₋₁, p, s, λ, τ₀, β₀, log2n, Δτ, shift); + logFourierLoad(Cₜ₋₁, Rₜ₋₁, Ĉₜ₋₁, Ȓₜ₋₁, p, s, λ, τ₀, β₀, log2n, Δτ, logShift); μₜ₋₁ = estimateZ(fft, Cₜ₋₁, Ĉₜ₋₁, Rₜ₋₁, Ȓₜ₋₁, p, s, λ, τ₀, β₀); } @@ -139,7 +140,7 @@ int main(int argc, char* argv[]) { std::vector Ȓₜ₊₁(N); for (unsigned n = 0; n < N; n++) { Ȓₜ₊₁[n] = ((Real)1.0 + std::pow(β, 2) * RddfCt[n] * Ȓₜ[n]) / (μₜ + II * fft.ν(n)); - Ĉₜ₊₁[n] = (2 * Γ₀ * std::conj(Ȓₜ[n]) / (1 + std::pow(τ₀ * fft.ν(n), 2)) + std::pow(β, 2) * (RddfCt[n] * Ĉₜ[n] + dfCt[n] * std::conj(Ȓₜ[n]))) / (μₜ + II * fft.ν(n)); + Ĉₜ₊₁[n] = ((2 * Γ₀ * std::conj(Ȓₜ[n]) / (1 + std::pow(τ₀ * fft.ν(n), 2)) + std::pow(β, 2) * (RddfCt[n] * Ĉₜ[n] + dfCt[n] * std::conj(Ȓₜ[n]))) / (μₜ + II * fft.ν(n))).real(); } std::vector Rₜ₊₁ = fft.inverse(Ȓₜ₊₁); std::vector Cₜ₊₁ = fft.inverse(Ĉₜ₊₁); @@ -190,7 +191,7 @@ int main(int argc, char* argv[]) { std::cerr << "\x1b[2K" << "\r"; std::cerr << β << " " << μₜ << " " << Ĉₜ[0].real() << " " << E << " " << γ << std::endl; - logFourierSave(Cₜ, Rₜ, Ĉₜ, Ȓₜ, p, s, λ, τ₀, β, log2n, Δτ, shift); + logFourierSave(Cₜ, Rₜ, Ĉₜ, Ȓₜ, p, s, λ, τ₀, β, log2n, Δτ, logShift); β += Δβ; Cₜ₋₁ = Cₜ; diff --git a/log_get_energy.cpp b/log_get_energy.cpp index b01034c..a183861 100644 --- a/log_get_energy.cpp +++ b/log_get_energy.cpp @@ -16,7 +16,7 @@ int main(int argc, char* argv[]) { Real Δτ = 0.1; Real k = 0.1; unsigned pad = 2; - Real shift = 0.5; + Real logShift = 0; /* Iteration parameters */ Real β₀ = 0; @@ -25,7 +25,7 @@ int main(int argc, char* argv[]) { int opt; - while ((opt = getopt(argc, argv, "p:s:2:T:t:b:d:k:D:0:h:")) != -1) { + while ((opt = getopt(argc, argv, "p:s:2:T:t:b:d:k:h:D:0:")) != -1) { switch (opt) { case 'p': p = atoi(optarg); @@ -49,7 +49,7 @@ int main(int argc, char* argv[]) { k = atof(optarg); break; case 'h': - shift = atof(optarg); + logShift = atof(optarg); break; case 'D': Δτ = atof(optarg); @@ -63,8 +63,10 @@ int main(int argc, char* argv[]) { } unsigned N = pow(2, log2n); + Real Γ₀ = 1; + Real μ₀ = τ₀ > 0 ? (sqrt(1+4*Γ₀*τ₀)-1)/(2*τ₀) : Γ₀; - LogarithmicFourierTransform fft(N, k, Δτ, pad, shift); + LogarithmicFourierTransform fft(N, k, Δτ, pad, μ₀ * pow(10, logShift)); std::vector C(N); std::vector R(N); @@ -76,7 +78,7 @@ int main(int argc, char* argv[]) { std::cout << std::setprecision(16); while (β += Δβ, β <= βₘₐₓ) { - if (logFourierLoad(C, R, Ct, Rt, p, s, λ, τ₀, β, log2n, Δτ, shift)) { + if (logFourierLoad(C, R, Ct, Rt, p, s, λ, τ₀, β, log2n, Δτ, logShift)) { Real e = energy(fft, C, R, p, s, λ, β); -- cgit v1.2.3-70-g09d2