diff options
| -rw-r--r-- | langevin.cpp | 28 | ||||
| -rw-r--r-- | stokes.hpp | 50 |
2 files changed, 56 insertions, 22 deletions
diff --git a/langevin.cpp b/langevin.cpp index e1464ec..e466860 100644 --- a/langevin.cpp +++ b/langevin.cpp @@ -182,30 +182,26 @@ int main(int argc, char* argv[]) { Rope<Complex> stokes(10, z1, z2); - std::cout << stokes.cost(Jp) << std::endl; + std::cout << stokes.cost(Jp) << " " << stokes.length() << " " << stokes.error(Jp) << std::endl; - stokes.relax(Jp, 10000, 1e-4); + stokes.relax(Jp, 1e5, 1e-1); - std::cout << stokes.cost(Jp) << std::endl; + std::cout << stokes.cost(Jp) << " " << stokes.length() << " " << stokes.error(Jp) << std::endl; - Rope<Complex> stokes2 = stokes.interpolate(); + stokes = stokes.interpolate(); + stokes.relax(Jp, 1e5, 1e-1); - stokes2.relax(Jp, 10000, 1e-4); + std::cout << stokes.cost(Jp) << " " << stokes.length() << " " << stokes.error(Jp) << std::endl; - std::cout << stokes2.cost(Jp) << std::endl; + stokes = stokes.interpolate(); + stokes.relax(Jp, 1e5, 1e-1); - stokes2 = stokes2.interpolate(); + std::cout << stokes.cost(Jp) << " " << stokes.length() << " " << stokes.error(Jp) << std::endl; - stokes2.relax(Jp, 10000, 1e-4); - - std::cout << stokes2.cost(Jp) << std::endl; - - stokes2 = stokes2.interpolate(); - - stokes2.relax(Jp, 10000, 1e-4); - - std::cout << stokes2.cost(Jp) << std::endl; + stokes = stokes.interpolate(); + stokes.relax(Jp, 1e5, 1e-1); + std::cout << stokes.cost(Jp) << " " << stokes.length() << " " << stokes.error(Jp) << std::endl; return 0; } @@ -59,6 +59,36 @@ class Rope { return z.size() - 2; } + double length() const { + double l = 0; + + for (unsigned i = 0; i < z.size() - 1; i++) { + l += (z[i + 1] - z[i]).norm(); + } + + return l; + } + + template <int p> + double error(const Tensor<Scalar, p>& J) const { + Scalar H0, HN; + std::tie(H0, std::ignore, std::ignore) = hamGradHess(J, z.front()); + std::tie(HN, std::ignore, std::ignore) = hamGradHess(J, z.back()); + + double ImH = imag((H0 + HN) / 2.0); + + double err = 0; + + for (unsigned i = 1; i < z.size() - 1; i++) { + Scalar Hi; + std::tie(Hi, std::ignore, std::ignore) = hamGradHess(J, z[i]); + + err += pow(imag(Hi) - ImH, 2); + } + + return sqrt(err); + } + Rope(unsigned N, const Vector<Scalar>& z1, const Vector<Scalar>& z2) : z(N + 2) { for (unsigned i = 0; i < N + 2; i++) { z[i] = normalize(z1 + (z2 - z1) * ((double)i / (N + 1.0))); @@ -86,6 +116,16 @@ class Rope { δz[i] = variation(z[i], dz(i), ddz(i), dH, ddH); } + // We return a δz with average norm of one. + double mag = 0; + for (unsigned i = 1; i < z.size() - 1; i++) { + mag += δz[i].norm(); + } + + for (unsigned i = 1; i < z.size() - 1; i++) { + δz[i] /= mag / n(); + } + return δz; } @@ -95,10 +135,12 @@ class Rope { Rope rNew = *this; double δ = δ0; + // We rescale the step size by the distance between points. + double Δl = length() / (n() + 1); while (rNew.cost(J) >= cost(J)) { for (unsigned i = 1; i < z.size() - 1; i++) { - rNew.z[i] = normalize(z[i] - δ * Δz[i].conjugate()); + rNew.z[i] = normalize(z[i] - (δ * Δl) * Δz[i].conjugate()); } δ /= 2; @@ -114,11 +156,7 @@ class Rope { } void spread() { - double l = 0; - - for (unsigned i = 0; i < z.size() - 1; i++) { - l += (z[i + 1] - z[i]).norm(); - } + double l = length(); double a = 0; unsigned pos = 0; |