diff options
| -rw-r--r-- | langevin.cpp | 9 | ||||
| -rw-r--r-- | stokes.hpp | 42 |
2 files changed, 35 insertions, 16 deletions
diff --git a/langevin.cpp b/langevin.cpp index e466860..c3f8436 100644 --- a/langevin.cpp +++ b/langevin.cpp @@ -167,7 +167,7 @@ int main(int argc, char* argv[]) { prevdist = abs(imag(H1 - H2)); εJ = 1e4 * prevdist; - if (abs(imag(H1 - H2)) < 1e-10 && dist > 1e-2) { + if (abs(imag(H1 - H2)) < 1e-13 && dist > 1e-2) { std::cerr << "Found distinct critical points with sufficiently similar Im H." << std::endl; break; } @@ -180,7 +180,7 @@ int main(int argc, char* argv[]) { } } - Rope<Complex> stokes(10, z1, z2); + Rope<Complex> stokes(10, z1, z2, Jp); std::cout << stokes.cost(Jp) << " " << stokes.length() << " " << stokes.error(Jp) << std::endl; @@ -203,5 +203,10 @@ int main(int argc, char* argv[]) { std::cout << stokes.cost(Jp) << " " << stokes.length() << " " << stokes.error(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; } @@ -22,9 +22,9 @@ Vector<Scalar> variation(const Vector<Scalar>& z, const Vector<Scalar>& z´, con Matrix<Scalar>::Identity(ddH.rows(), ddH.cols()) - z.conjugate() * z.transpose() / z² ); - Vector<Scalar> dLdz = - Reż·z´ * ( + Vector<Scalar> dLdz = - ( dżc * z´ + dż * z´.conjugate() - (dż * ż.conjugate() + dżc * ż) * Reż·z´ / ż² - ) / (ż² * z´²); + ) / sqrt(ż² * z´²) / 2; Vector<Scalar> ż´ = -(ddH * z´).conjugate() + ((ddH * z´).dot(z) / z²) * z.conjugate() + ( dH.dot(z) * z´.conjugate() + dH.dot(z´) * z.conjugate() - ( @@ -35,17 +35,16 @@ Vector<Scalar> variation(const Vector<Scalar>& z, const Vector<Scalar>& z´, con double dReż·z´ = real(ż.dot(z´´) + ż´.dot(z´)); Vector<Scalar> ddtdLdz´ = - ( - Reż·z´ * ( + ( ż´.conjugate() - ( Reż·z´ * z´´.conjugate() + dReż·z´ * z´.conjugate() - (Reż·z´ / z´²) * (z´´.dot(z´) + z´.dot(z´´)) * z´.conjugate() ) / z´² ) - + dReż·z´ * (ż.conjugate() - Reż·z´ / z´² * z´.conjugate()) - - Reż·z´ * ( + - 0.5 * ( (ż.dot(ż´) + ż´.dot(ż)) / ż² + (z´´.dot(z´) + z´.dot(z´´)) / z´² ) * (ż.conjugate() - Reż·z´ / z´² * z´.conjugate()) - ) / (ż² * z´²); + ) / sqrt(ż² * z´²) / 2; return dLdz - ddtdLdz´; } @@ -55,6 +54,27 @@ class Rope { public: std::vector<Vector<Scalar>> z; + Rope(unsigned N) : z(N + 2) {} + + template <int p> + Rope(unsigned N, const Vector<Scalar>& z1, const Vector<Scalar>& z2, const Tensor<Scalar, p>& J) : z(N + 2) { + Scalar H1, H2; + std::tie(H1, std::ignore, std::ignore) = hamGradHess(J, z1); + std::tie(H2, std::ignore, std::ignore) = hamGradHess(J, z2); + + if (real(H1) > real(H2)) { + z.front() = z1; + z.back() = z2; + } else { + z.front() = z2; + z.back() = z1; + } + + for (unsigned i = 0; i < N + 2; i++) { + z[i] = normalize(z.front() + (z.back() - z.front()) * ((double)i / (N + 1.0))); + } + } + unsigned n() const { return z.size() - 2; } @@ -89,12 +109,6 @@ class Rope { 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))); - } - } - Vector<Scalar> dz(unsigned i) const { return z[i + 1] - z[i - 1]; } @@ -208,14 +222,14 @@ class Rope { Vector<Scalar> zD = zDot(z[i], dH); - c += 1.0 - pow(real(zD.dot(dz(i))), 2) / zD.squaredNorm() / dz(i).squaredNorm(); + c += 1.0 - real(zD.dot(dz(i))) / zD.norm() / dz(i).norm(); } return c; } Rope<Scalar> interpolate() const { - Rope<Scalar> r(2 * n() + 1, z[0], z[z.size() - 1]); + Rope<Scalar> r(2 * n() + 1); for (unsigned i = 0; i < z.size(); i++) { r.z[2 * i] = z[i]; |