summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--.gitignore9
-rw-r--r--Makefile10
-rw-r--r--fourier.cpp84
-rw-r--r--fourier.hpp14
-rw-r--r--fourier_integrator.cpp54
-rw-r--r--get_energy.cpp23
-rw-r--r--integrator.cpp247
7 files changed, 233 insertions, 208 deletions
diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..13e1cac
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,9 @@
+*.dat
+*.o
+correlations
+fourier_integrator
+vgcore.*
+walk
+get_energy
+integrator
+fftw.wisdom
diff --git a/Makefile b/Makefile
index d4c15bf..688a99c 100644
--- a/Makefile
+++ b/Makefile
@@ -1,6 +1,6 @@
all: walk correlations integrator fourier_integrator get_energy
-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,14 @@ 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
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 -lfftw3_omp -o 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 -lfftw3_omp -o get_energy
diff --git a/fourier.cpp b/fourier.cpp
index 7808989..bc4b633 100644
--- a/fourier.cpp
+++ b/fourier.cpp
@@ -24,37 +24,81 @@ 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> ĉ(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> ĉ(â.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<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;
+}
+
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,17 +114,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(); i++) {
- RddfC[i] = R[i] * ddf(λ, p, s, C[i]);
+ for (unsigned i = 0; i < C.size() / 2; 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 τ₀, 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();
+}
diff --git a/fourier.hpp b/fourier.hpp
index 86335f7..791953b 100644
--- a/fourier.hpp
+++ b/fourier.hpp
@@ -3,6 +3,11 @@
#include <fftw3.h>
#include <complex>
#include <vector>
+#include <tuple>
+
+#ifndef FFTW_THREADS
+#define FFTW_THREADS 1
+#endif
using Real = double;
using Complex = std::complex<Real>;
@@ -14,19 +19,24 @@ 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);
};
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 τ₀, Real y);
diff --git a/fourier_integrator.cpp b/fourier_integrator.cpp
index 86c3b24..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,24 +98,26 @@ 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);
- auto [RddfCt, dfCt] = RddfCtdfCt(fft, C, R, p, s, λ);
-
- z = ((Γ₀ * std::conj(Rt[0]) + pow(y, 2) * (RddfCt[0] * Ct[0] + dfCt[0] * std::conj(Rt[0]))) / Ct[0]).real();
+ z = estimateZ(fft, C, Ct, R, Rt, p, s, λ, τ₀, y);
}
while (y += Δy, y <= yₘₐₓ) {
- Real ΔC = 1;;
+ Real ΔC = 1;
+ Real ΔCprev = 1000;
unsigned it = 0;
- while (sqrt(ΔC / C.size()) > ε) {
+ while (sqrt(2 * ΔC / C.size()) > ε) {
it++;
auto [RddfCt, dfCt] = RddfCtdfCt(fft, C, R, p, s, λ);
@@ -132,7 +138,7 @@ int main(int argc, char* argv[]) {
}
ΔC = 0;
- for (unsigned i = 0; i < Cnew.size(); i++) {
+ for (unsigned i = 0; i < Cnew.size() / 2; i++) {
ΔC += pow(Cnew[i] - C[i], 2);
ΔC += pow(Rnew[i] - R[i], 2);
}
@@ -141,18 +147,24 @@ int main(int argc, char* argv[]) {
C[i] += γ * (Cnew[i] - C[i]);
}
- for (unsigned i = 0; i < Rnew.size(); i++) {
+ for (unsigned i = 0; i < Rnew.size() / 2; i++) {
R[i] += γ * (Rnew[i] - R[i]);
}
z *= Cnew[0];
- std::cerr << it << " " << p << " " << s << " " << τ₀ << " " << sqrt(ΔC / C.size()) << std::endl;
+ if (it % maxIterations == 0) {
+ if (ΔC < ΔCprev) {
+ γ = std::min(1.1 * γ, 1.0);
+ } else {
+ γ /= 2;
+ }
- if (it > maxIterations) {
- it = 0;
- γ /= 2;
+ ΔCprev = ΔC;
}
+
+ std::cerr << it << " " << p << " " << s << " " << τ₀ << " " << y << " " << sqrt(2 * ΔC / C.size()) << " " << γ << std::endl;
+
}
Real e = energy(C, R, p, s, λ, y, Δτ);
@@ -160,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 a388e0c..c33c04e 100644
--- a/get_energy.cpp
+++ b/get_energy.cpp
@@ -50,24 +50,27 @@ 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, Δω, Δτ);
- while (y += Δy, y <= yₘₐₓ) {
- std::vector<Real> C(2 * n);
- std::vector<Real> R(2 * n);
+ std::vector<Real> C(2 * n);
+ std::vector<Real> R(2 * n);
+ 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,9 +78,7 @@ int main(int argc, char* argv[]) {
std::vector<Complex> Ct = fft.fourier(C);
std::vector<Complex> Rt = fft.fourier(R);
- auto [RddfCt, dfCt] = RddfCtdfCt(fft, C, R, p, s, λ);
-
- Real z = ((std::conj(Rt[0]) + pow(y, 2) * (RddfCt[0] * Ct[0] + dfCt[0] * std::conj(Rt[0]))) / Ct[0]).real();
+ 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;
}