summaryrefslogtreecommitdiff
path: root/integrator.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'integrator.cpp')
-rw-r--r--integrator.cpp185
1 files changed, 41 insertions, 144 deletions
diff --git a/integrator.cpp b/integrator.cpp
index 8b13a9b..faba06d 100644
--- a/integrator.cpp
+++ b/integrator.cpp
@@ -1,112 +1,28 @@
+#include "fourier.hpp"
#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[]) {
+ unsigned p = 3;
+ unsigned s = 4;
+ Real λ = 0.5;
Real Δτ = 1e-3;
Real τₘₐₓ = 1e3;
Real τ₀ = 0;
- Real y = 0.5;
+ Real β = 0.5;
unsigned iterations = 10;
int opt;
- while ((opt = getopt(argc, argv, "d:T:t:y:I:")) != -1) {
+ while ((opt = getopt(argc, argv, "d:T:t:b:I:")) != -1) {
switch (opt) {
case 'd':
Δτ = atof(optarg);
@@ -117,8 +33,8 @@ int main(int argc, char* argv[]) {
case 't':
τ₀ = atof(optarg);
break;
- case 'y':
- y = atof(optarg);
+ case 'b':
+ β = atof(optarg);
break;
case 'I':
iterations = (unsigned)atof(optarg);
@@ -128,63 +44,44 @@ int main(int argc, char* argv[]) {
}
}
- Real z = 0.4794707565634420155347;
Real Γ₀ = 1;
+ Real μ = 1;
+ if (τ₀ > 0) {
+ μ = (sqrt(1+4*Γ₀*τ₀) - 1) / (2 * τ₀);
+ }
Real τ = 0;
- std::vector<Real> C;
- C.reserve(τₘₐₓ / Δτ + 1);
-
- C.push_back(1);
-
-// 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);
+ unsigned N = τₘₐₓ / Δτ + 1;
+ std::vector<Real> C(N);
+ std::vector<Real> R(N);
+ std::vector<Real> Γ(N);
+
+ for (unsigned i = 0; i < N; i++) {
+ Real τ = i * Δτ;
+ if (τ₀ > 0) {
+ C[i] = (Γ₀ / μ) * (exp(-μ * τ) - μ * τ₀ * exp(-τ / τ₀)) / (1 - pow(μ * τ₀, 2));
+ Γ[i] = (Γ₀ / τ₀) * exp(-τ / τ₀);
+ } else {
+ C[i] = (Γ₀ / μ) * exp(-μ * τ);
+ }
+ R[i] = exp(-μ * τ);
}
-
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 error = 0;
-
- for (unsigned i = 0; i < std::min(C.size(), C2.size()); i++) {
- error += pow(C[i] - C2[i], 2);
+ /* 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] + pow(β, 2) * I;
+ R₊[i] = R₊[i - 1] + dR * Δτ;
}
- 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);
+ /* Second step: integrate C from R */
}
- */
τ = 0;
for (Real Ci : C) {
@@ -192,7 +89,7 @@ int main(int argc, char* argv[]) {
τ += Δτ;
}
- std::cerr << - 2 * y / Γ₀ * energy(C, Δτ, τ₀) << std::endl;
+ std::cerr << - 2 * β / Γ₀ * energy(C, R, λ, p, s, Δτ) << std::endl;
return 0;
}