summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--fourier.cpp58
-rw-r--r--fourier.hpp8
-rw-r--r--log-fourier_integrator.cpp110
3 files changed, 111 insertions, 65 deletions
diff --git a/fourier.cpp b/fourier.cpp
index cf45ad1..bb1c92f 100644
--- a/fourier.cpp
+++ b/fourier.cpp
@@ -1,5 +1,5 @@
#include "fourier.hpp"
-#include <boost/multiprecision/fwd.hpp>
+#include <fftw3.h>
inline Real fP(unsigned p, Real q) {
return 0.5 * pow(q, p);
@@ -100,15 +100,15 @@ 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) {
+LogarithmicFourierTransform::LogarithmicFourierTransform(unsigned N, Real k, Real Δτ, unsigned pad) : N(N), pad(pad), k(k), Δτ(Δτ) {
τₛ = -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));
+ sₛ = -0.5 * pad * N;
+ a = reinterpret_cast<Complex*>(fftw_alloc_complex(pad*N));
+ â = reinterpret_cast<Complex*>(fftw_alloc_complex(pad*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);
+ a_to_â = fftw_plan_dft_1d(pad*N, reinterpret_cast<fftw_complex*>(a), reinterpret_cast<fftw_complex*>(â), FFTW_BACKWARD, 0);
+ â_to_a = fftw_plan_dft_1d(pad*N, reinterpret_cast<fftw_complex*>(â), reinterpret_cast<fftw_complex*>(a), FFTW_BACKWARD, 0);
fftw_export_wisdom_to_filename("fftw.wisdom");
}
@@ -125,11 +125,11 @@ Real LogarithmicFourierTransform::τ(unsigned n) const {
}
Real LogarithmicFourierTransform::ω(unsigned n) const {
- return Δω * (n + ωₛ);
+ return Δτ * (n + ωₛ);
}
Real LogarithmicFourierTransform::s(unsigned n) const {
- return Δs * (n + sₛ);
+ return (n + sₛ) * 2*M_PI / (pad * N * Δτ);
}
Real LogarithmicFourierTransform::t(unsigned n) const {
@@ -149,29 +149,57 @@ Complex gamma(Complex z) {
return exp(logΓ.val + 1i * argΓ.val);
}
-std::vector<Complex> LogarithmicFourierTransform::fourier(const std::vector<Real>& c, bool symmetric){
+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));
+ for (unsigned n = 0; n < pad*N; n++) {
+ if (n < N) {
+ a[n] = c[n] * exp((1 - k) * τ(n));
+ } else {
+ a[n] = 0;
+ }
}
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)));
+ for (unsigned n = 0; n < pad*N; n++) {
+ â[(pad*N / 2 + n) % (pad*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);
+ ĉ[n] += exp(-k * ω(n)) * a[(pad - 1)*N+n].real() / (Real)(pad*N);
}
}
return ĉ;
}
+std::vector<Real> LogarithmicFourierTransform::inverse(const std::vector<Complex>& ĉ) {
+ std::vector<Real> c(N);
+ std::vector<Real> σs = {1, -1};
+ for (Real σ : σs) {
+ for (unsigned n = 0; n < pad * N; n++) {
+ if (n < N) {
+ a[n] = ĉ[n] * exp((1 - k) * ω(n));
+ } else {
+ a[n] = 0;
+ }
+ }
+ fftw_execute(a_to_â);
+ for (unsigned n = 0; n < pad*N; n++) {
+ â[(pad*N / 2 + n) % (pad*N)] *= pow(-1i * σ, 1i * s(n) - k) * gamma(k - 1i * s(n));
+ }
+ fftw_execute(â_to_a);
+ for (unsigned n = 0; n < N; n++) {
+ c[n] += exp(-k * τ(n)) * a[(pad - 1)*N+n].real() / (Real)(pad*N) / (2 * M_PI);
+ }
+ }
+
+ return c;
+}
+
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";
}
diff --git a/fourier.hpp b/fourier.hpp
index 22a57e7..1708667 100644
--- a/fourier.hpp
+++ b/fourier.hpp
@@ -44,15 +44,14 @@ private:
fftw_plan a_to_â;
fftw_plan â_to_a;
unsigned N;
+ unsigned pad;
Real k;
Real Δτ;
- Real Δω;
- Real Δs;
Real τₛ;
Real ωₛ;
Real sₛ;
public:
- LogarithmicFourierTransform(unsigned N, Real k, Real Δω, Real Δτ, Real Δs);
+ LogarithmicFourierTransform(unsigned N, Real k, Real Δτ, unsigned pad = 4);
~LogarithmicFourierTransform();
Real τ(unsigned n) const;
Real ω(unsigned n) const;
@@ -60,10 +59,7 @@ public:
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);
diff --git a/log-fourier_integrator.cpp b/log-fourier_integrator.cpp
index 247270b..177db46 100644
--- a/log-fourier_integrator.cpp
+++ b/log-fourier_integrator.cpp
@@ -1,27 +1,25 @@
#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;
+ Real Δτ = 0.1;
+ Real k = 0.1;
- unsigned maxIterations = 1000;
- Real ε = 1e-14;
+ Real ε = 1e-16;
Real γ = 1;
+ Real βₘₐₓ = 0.7;
+ Real Δβ = 0.01;
int opt;
- while ((opt = getopt(argc, argv, "p:s:2:T:t:0:y:d:I:g:l")) != -1) {
+ while ((opt = getopt(argc, argv, "p:s:2:T:t:b:d:g:k:D:")) != -1) {
switch (opt) {
case 'p':
p = atoi(optarg);
@@ -32,27 +30,24 @@ int main(int argc, char* argv[]) {
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);
+ case 'b':
+ βₘₐₓ = atof(optarg);
break;
case 'd':
- Δy = atof(optarg);
- break;
- case 'I':
- maxIterations = (unsigned)atof(optarg);
+ Δβ = atof(optarg);
break;
case 'g':
γ = atof(optarg);
break;
+ case 'k':
+ k = atof(optarg);
+ break;
+ case 'D':
+ Δτ = atof(optarg);
+ break;
default:
exit(1);
}
@@ -60,10 +55,7 @@ int main(int argc, char* argv[]) {
unsigned N = pow(2, log2n);
- Real Δτ = 1e-2;
- Real Δω = 1e-2;
- Real Δs = 1e-2;
- Real k = 0.1;
+ LogarithmicFourierTransform fft(N, k, Δτ, 4);
Real Γ₀ = 1.0;
Real μ = Γ₀;
@@ -73,9 +65,10 @@ int main(int argc, char* argv[]) {
std::vector<Real> C(N);
std::vector<Real> R(N);
+ std::vector<Complex> Ct(N);
+ std::vector<Complex> Rt(N);
- LogarithmicFourierTransform fft(N, k, Δω, Δτ, Δs);
-
+ // start from the exact solution for β = 0
for (unsigned n = 0; n < N; n++) {
if (τ₀ > 0) {
C[n] = Γ₀ * (exp(-μ * fft.t(n)) - μ * τ₀ * exp(-fft.t(n) / τ₀)) / (μ - pow(μ, 3) * pow(τ₀, 2));
@@ -83,32 +76,61 @@ int main(int argc, char* argv[]) {
C[n] = Γ₀ * exp(-μ * fft.t(n)) / μ;
}
R[n] = exp(-μ * fft.t(n));
+
+ Ct[n] = 2 * Γ₀ / (pow(μ, 2) + pow(fft.ν(n), 2)) / (1 + pow(τ₀ * fft.ν(n), 2));
+ Rt[n] = 1.0 / (μ + 1i * fft.ν(n));
}
- std::vector<Complex> Ct = fft.fourier(C, true);
- std::vector<Complex> Rt = fft.fourier(R, false);
+ Real β = 0;
+ while (β < βₘₐₓ) {
+ Real ΔC = 100;
+ while (ΔC > ε) {
+ std::vector<Real> RddfC(N);
+ std::vector<Real> dfC(N);
+ for (unsigned n = 0; n < N; n++) {
+ RddfC[n] = R[n] * ddf(λ, 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);
- /*
- for (unsigned n = 0; n < N; n++) {
- std::cout << fft.t(n) << " " << C[n] << std::endl;
+ std::vector<Complex> Rtnew(N);
+ std::vector<Complex> Ctnew(N);
+
+ 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));
+ }
+
+ std::vector<Real> Cnew = fft.inverse(Ctnew);
+ std::vector<Real> Rnew = fft.inverse(Rtnew);
+
+ ΔC = 0;
+ for (unsigned i = 0; i < N; i++) {
+ ΔC += std::norm(Ct[i] - Ctnew[i]);
+ ΔC += std::norm(Rt[i] - Rtnew[i]);
+ }
+ ΔC = sqrt(ΔC) / (2*N);
+
+ for (unsigned i = 0; i < N; i++) {
+ C[i] += γ * (Cnew[i] - C[i]);
+ R[i] += γ * (Rnew[i] - R[i]);
+ Ct[i] += γ * (Ctnew[i] - Ct[i]);
+ Rt[i] += γ * (Rtnew[i] - Rt[i]);
+ }
+
+ μ *= C[0];
+
+// std::cerr << ΔC << std::endl;
}
- */
- for (unsigned n = 0; n < N; n++) {
- std::cout << fft.ν(n) << " " << Ct[n].real() << std::endl;
+ std::cerr << β << " " << μ << " " << Ct[0].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));
+ for (unsigned i = 0; i < N; i++) {
+ std::cout << fft.t(i) << " " << Rt[i].imag() << std::endl;
}
- C = fft.inverse(Ct);
- R = fft.inverse(Rt);
- */
return 0;
}