diff options
| -rw-r--r-- | Makefile | 2 | ||||
| -rw-r--r-- | integrator.cpp | 95 |
2 files changed, 86 insertions, 11 deletions
@@ -1,6 +1,6 @@ all: walk correlations integrator -CC := clang++ -std=c++17 -Wno-mathematical-notation-identifier-extension -O3 -march=native -mtune=native -fopenmp +CC := clang++ -std=c++17 -Wno-mathematical-notation-identifier-extension -O3 -march=native -mtune=native -fopenmp -g walk: walk.cpp $(CC) walk.cpp -o walk diff --git a/integrator.cpp b/integrator.cpp index e184de6..349e2ff 100644 --- a/integrator.cpp +++ b/integrator.cpp @@ -19,23 +19,66 @@ Real ddf(Real q) { return 0.5 * p * (p - 1) * pow(q, p - 2); } -Real integrate(const std::vector<Real>& C, Real Δτ, Real τ₀) { +Real integrate(const std::vector<Real>& C) { Real I = 0; #pragma omp parallel for reduction(+:I) for (unsigned σ = 0; σ < C.size() - 1; σ++) { unsigned τ_σ = C.size() - 1 - σ; Real Cτ_σ = (C[τ_σ] + C[τ_σ - 1]) / 2; - Real dCσ = (C[σ + 1] - C[σ]) / Δτ; + Real dCσ = C[σ + 1] - C[σ]; - Real dddC = 0; - if (σ > 3 && σ < C.size() && C.size() > 3) { - dddC += (C[τ_σ] - 3 * C[τ_σ+1] + 3 * C[τ_σ+2] - C[τ_σ+3]) / pow(Δτ, 3); + 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 += Δτ * df(Cτ_σ) * (dCσ - pow(τ₀, 2) * dddC); + + 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; σ++) { @@ -78,8 +121,10 @@ int main(int argc, char* argv[]) { } } - Real z = 0.5; - Real Γ₀ = 0.5 * (2 * τ₀) / (sqrt(1 + 2 * τ₀) - 1); +// Real z = (sqrt(1 + 2 * τ₀) - 1) / (2 * τ₀); +// Real z = 0.4833773593561142778087; + Real z = 0.4794707565634420155347 - 2 * pow(y, 2); + Real Γ₀ = 1.0; Real τ = 0; std::vector<Real> C; @@ -87,12 +132,42 @@ int main(int argc, char* argv[]) { C.push_back(1); - while (std::cout << τ << " " << C.back() << std::endl, τ < τₘₐₓ) { +// while (std::cout << τ << " " << C.back() << std::endl, τ < τₘₐₓ) { + while (τ < τₘₐₓ) { τ += Δτ; - Real dC = -z * C.back() - 2 / Γ₀ * pow(y, 2) * integrate(C, Δτ, τ₀); + Real dC = -z * C.back() - 2 / Γ₀ * pow(y, 2) * integrate(C); C.push_back(C.back() + Δτ * dC); } + + for (unsigned it = 0; it < 20; it++) { + τ = 0; + std::vector<Real> C2; + C2.reserve(τₘₐₓ / Δτ + 1); + C2.push_back(1); + while (τ < τₘₐₓ) { + τ += Δτ; + Real dC = -z * 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); + } + + std::cerr << "Iteration " << it << ": " << sqrt(error / C.size()) << std::endl; + + C = C2; + } + + τ = 0; + for (Real Ci : C) { + std::cout << τ << " " << Ci << std::endl; + τ += Δτ; + } + std::cerr << - 2 * y / Γ₀ * energy(C, Δτ, τ₀) << std::endl; return 0; |