Fixed several bugs.

1 files changed, 24 insertions, 10 deletions

diff --git a/integrator.cpp b/integrator.cpp
index bb22a36..ea5c477 100644
--- a/integrator.cpp
+++ b/integrator.cpp
@@ -16,29 +16,43 @@ public:
   }
 };
 
+unsigned p = 3;
+
 Real f(Real q) {
-  return 0.5 * pow(q, 2);
+  return 0.5 * pow(q, p);
 }
 
 Real df(Real q) {
-  return q;
+  return 0.5 * p * pow(q, p - 1);
 }
 
 Real ddf(Real q) {
-  return 1;
+  return 0.5 * p * (p - 1) * pow(q, p - 2);
 }
 
 Real integrate(const std::vector<Point>& Cₜ) {
   Real I = 0;
   for (unsigned i = 0; i < Cₜ.size() - 1; i++) {
+    unsigned j = Cₜ.size() - 1 - i;
     Real Δτ = Cₜ[i + 1].τ() - Cₜ[i].τ();
-    Real C = (Cₜ[i + 1].C() + Cₜ[i].C()) / 2;
+    Real C = (Cₜ[j].C() + Cₜ[j - 1].C()) / 2;
     Real dC = (Cₜ[i + 1].C() - Cₜ[i].C()) / Δτ;
     I += Δτ * df(C) * dC;
   }
   return I;
 }
 
+Real energy(const std::vector<Point>& Ct) {
+  Real I = 0;
+  for (unsigned i = 0; i < Ct.size() - 1; i++) {
+    Real Δτ = Ct[i + 1].τ() - Ct[i].τ();
+    Real C = (Ct[i].C() + Ct[i + 1].C()) / 2;
+    Real dC = (Ct[i + 1].C() - Ct[i].C()) / Δτ;
+    I += Δτ * df(C) * dC;
+  }
+  return I;
+}
+
 int main(int argc, char* argv[]) {
   Real Δτ = 1e-3;
   Real τₘₐₓ = 1e3;
@@ -66,20 +80,20 @@ int main(int argc, char* argv[]) {
     }
   }
 
+  Real z = 0.5 + 0.0 * pow(y, 2);
+
   std::vector<Point> Cₜ;
-  Cₜ.reserve(τₘₐₓ / Δτ);
+  Cₜ.reserve(τₘₐₓ / Δτ + 1);
 
   Cₜ.push_back({0, 1});
-  Cₜ.push_back({Δτ, 1 - Δτ});
 
   while (Cₜ.back().τ() < τₘₐₓ) {
-    Real dC = -Cₜ.back().C() - 2 * pow(y, 2) * integrate(Cₜ);
+    Real dC = -z * Cₜ.back().C() - 2 * pow(y, 2) * integrate(Cₜ);
     Cₜ.push_back({Cₜ.back().τ() + Δτ, Cₜ.back().C() + Δτ * dC});
+    std::cout << Cₜ.back().τ() << " " << Cₜ.back().C() << std::endl;
   }
 
-  for (const Point& p : Cₜ) {
-    std::cout << p.τ() << " " << p.C() << std::endl;
-  }
+  std::cerr << - 2 * y * energy(Cₜ) << std::endl;
 
   return 0;
 }