diff options
| -rw-r--r-- | dynamics.hpp | 8 | ||||
| -rw-r--r-- | langevin.cpp | 7 |
2 files changed, 9 insertions, 6 deletions
diff --git a/dynamics.hpp b/dynamics.hpp index 639803b..21afd30 100644 --- a/dynamics.hpp +++ b/dynamics.hpp @@ -55,13 +55,13 @@ Vector<Real> findMinimum(const Tensor<Real, p>& J, const Vector<Real>& z0, Real Matrix<Real> M = ddH - (dH * z.transpose() + z.dot(dH) * Matrix<Real>::Identity(z.size(), z.size()) + (ddH * z) * z.transpose()) / (Real)z.size() + 2.0 * z * z.transpose(); while (g.norm() / z.size() > ε && λ < 1e8) { - Vector<Real> dz = (M + λ * Matrix<Real>::Identity(z.size(), z.size())).partialPivLu().solve(g); + Vector<Real> dz = (M + λ * (Matrix<Real>)abs(M.diagonal().array()).matrix().asDiagonal()).partialPivLu().solve(g); dz -= z.dot(dz) * z / (Real)z.size(); Vector<Real> zNew = normalize(z - dz); auto [HNew, dHNew, ddHNew] = hamGradHess(J, zNew); - if (HNew <= H * 1.01) { + if (HNew * 1.0001 <= H) { z = zNew; H = HNew; dH = dHNew; @@ -70,12 +70,10 @@ Vector<Real> findMinimum(const Tensor<Real, p>& J, const Vector<Real>& z0, Real g = dH - z.dot(dH) * z / (Real)z.size(); M = ddH - (dH * z.transpose() + z.dot(dH) * Matrix<Real>::Identity(z.size(), z.size()) + (ddH * z) * z.transpose()) / (Real)z.size() + 2.0 * z * z.transpose(); - λ /= 1.001; + λ /= 2; } else { λ *= 1.5; } - - std::cout << "error : " << H << " " << g.norm() / z.size() << " " << λ << std::endl; } return z; diff --git a/langevin.cpp b/langevin.cpp index a775774..1181173 100644 --- a/langevin.cpp +++ b/langevin.cpp @@ -146,8 +146,13 @@ int main(int argc, char* argv[]) { Vector<Real> zMin = randomMinimum(ReJ, Red, r, ε); auto [Hr, dHr, ddHr] = hamGradHess(ReJ, zMin); - Eigen::EigenSolver<Matrix<Real>> eigenS(ddHr - ((ddHr * zMin) * zMin.transpose()) / (Real)zMin.size()); + Eigen::EigenSolver<Matrix<Real>> eigenS(ddHr - (dHr * zMin.transpose() + zMin.dot(dHr) * Matrix<Real>::Identity(zMin.size(), zMin.size()) + (ddHr * zMin) * zMin.transpose()) / (Real)zMin.size() + 2.0 * zMin * zMin.transpose()); std::cout << eigenS.eigenvalues().transpose() << std::endl; + for (unsigned i = 0; i < N; i++) { + Vector<Real> zNew = normalize(zMin + 0.01 * eigenS.eigenvectors().col(i).real()); + std::cout << getHamiltonian(ReJ, zNew) - Hr << " " << real(eigenS.eigenvectors().col(i).dot(zMin)) << std::endl; + } + std::cout << std::endl; getchar(); complex_normal_distribution<Real> d(0, 1, 0); |