“ W4方法的实现”（牛顿-拉普森扩展）

Question

我目前正在实施一种Newton-raphson寻根方法，该方法可确保在多维环境（而不是家庭作业！）中收敛。当前，它找到x的根，而不是y。我还观察到一个奇怪的行为，其中f1和f2等于相同的数字。例如，在2000次迭代之后，两者均为≈560.0。我认为f1和f2都需要接近0。至少，这是使用经典的newton-raphson方法的工作方式。

任何人都可以看到导致这种情况的原因吗？我需要另一双眼睛。

论文：https://arxiv.org/pdf/1809.04495.pdf和附录：https://arxiv.org/pdf/1809.04358.pdf（第D.2节->包括所附数学）

注：U，L是雅可比行列的上和下三角矩阵（偏导数矩阵）。

我当前的实现如下所示（使用了Eigen，但是很清楚它的作用）。目前有些奇怪

#include "../../Eigen/Eigen/Core"
#include "../../Eigen/Eigen/LU"
#include <iostream>

int main(){
    double eps = 1e-4;
    Eigen::Vector2d p(0.0, 0.0);

    double x = 0.1;
    double y = 1.0;
    double f1 = 1e9;
    double f2 = 1e9;

    unsigned int count = 0;

    while (count < 2000 && f1 > eps){
        std::cout << "count : " << count << std::endl;

        f1 = x*x - 10*x + y*y - 10*y + 34;
        f2 = x*x - 22*x + y*y - 10*y + 130;

        std::cout << "f1: " << f1 << ", f2: " << f2 << std::endl;

        double A = 2*x - 10;
        double B = 2*y - 10;
        double C = 2*x - 22;
        double D = 2*y - 10;

        Eigen::Matrix2d J;
        J << A, B, C, D;

        Eigen::Matrix2d J_U_inv;
        J_U_inv << J(0,0), J(0,1), 0.0, J(1,1);
        J_U_inv = J_U_inv.inverse();

        Eigen::Matrix2d J_L_inv;
        J_L_inv << J(0,0), 0.0, J(1,0), J(1,1);
        J_L_inv = J_L_inv.inverse();

        Eigen::Vector2d f3(f1, f2);
        Eigen::Vector2d T(x, y);

        if (count == 0){
            p = -0.5 * J_U_inv * f3;
        }

        Eigen::Vector2d E = T + 0.5 * J_L_inv * p;

        p = -0.5 * J_U_inv * f3;

        x = E(0);
        y = E(1);

        std::cout << "x, y: " << x << ", " << y << std::endl;

        ++count;

    }
}

我目前正在实施一种Newton-raphson寻根方法，该方法可确保在多维环境（而不是家庭作业！）中收敛。当前，它找到x的根，而不是y。我也是...

Answer 1

似乎我不知道进行矩阵分解的正确方法。