diff options
| -rw-r--r-- | integrator.cpp | 48 |
1 files changed, 36 insertions, 12 deletions
diff --git a/integrator.cpp b/integrator.cpp index 349e2ff..a94f37e 100644 --- a/integrator.cpp +++ b/integrator.cpp @@ -19,11 +19,14 @@ Real ddf(Real q) { return 0.5 * p * (p - 1) * pow(q, p - 2); } -Real integrate(const std::vector<Real>& C) { +Real integrate(const std::vector<Real>& C, signed τ = std::numeric_limits<unsigned>::max()) { Real I = 0; + if (τ > C.size() - 1) { + τ = C.size() - 1; + } #pragma omp parallel for reduction(+:I) - for (unsigned σ = 0; σ < C.size() - 1; σ++) { - unsigned τ_σ = C.size() - 1 - σ; + for (unsigned σ = 0; σ < τ; σ++) { + unsigned τ_σ = τ - σ; Real Cτ_σ = (C[τ_σ] + C[τ_σ - 1]) / 2; Real dCσ = C[σ + 1] - C[σ]; @@ -99,10 +102,11 @@ int main(int argc, char* argv[]) { Real τₘₐₓ = 1e3; Real τ₀ = 0; Real y = 0.5; + unsigned iterations = 10; int opt; - while ((opt = getopt(argc, argv, "d:T:t:y:")) != -1) { + while ((opt = getopt(argc, argv, "d:T:t:y:I:")) != -1) { switch (opt) { case 'd': Δτ = atof(optarg); @@ -116,15 +120,16 @@ int main(int argc, char* argv[]) { case 'y': y = atof(optarg); break; + case 'I': + iterations = (unsigned)atof(optarg); + break; default: exit(1); } } -// Real z = (sqrt(1 + 2 * τ₀) - 1) / (2 * τ₀); -// Real z = 0.4833773593561142778087; - Real z = 0.4794707565634420155347 - 2 * pow(y, 2); - Real Γ₀ = 1.0; + Real z = 0.5; + Real Γ₀ = 1; Real τ = 0; std::vector<Real> C; @@ -135,19 +140,20 @@ int main(int argc, char* argv[]) { // while (std::cout << τ << " " << C.back() << std::endl, τ < τₘₐₓ) { while (τ < τₘₐₓ) { τ += Δτ; - Real dC = -z * C.back() - 2 / Γ₀ * pow(y, 2) * integrate(C); + Real dC = -(z - 2 * pow(y, 2)) * C.back() - 2 / Γ₀ * pow(y, 2) * integrate(C); C.push_back(C.back() + Δτ * dC); } - for (unsigned it = 0; it < 20; it++) { + for (unsigned it = 0; it < iterations; 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)); + 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); } @@ -157,10 +163,28 @@ int main(int argc, char* argv[]) { error += pow(C[i] - C2[i], 2); } - std::cerr << "Iteration " << it << ": " << sqrt(error / C.size()) << 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); + } + */ τ = 0; for (Real Ci : C) { |