Реализация метода W4 (расширение Ньютона-Рафсона)

В настоящее время я реализую метод поиска корня Ньютона-Рафсона, который гарантирует сходимость в многомерной среде (не домашняя работа!). В настоящее время он находит корень для x, но не для y. Я также наблюдал странное поведение, когда f1 и f2 приравнивались к одному и тому же числу. Например, после 2000 итераций оба равны ≈ 560,0. Я думаю, что f1 и f2 оба должны приближаться к 0. По крайней мере, так это работает с использованием классического метода Ньютона-Рафсона.

Кто-нибудь может увидеть, что может быть причиной этого? Мне нужна вторая пара глаз.

документ: 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;

    }
}

math linear-algebra numerical-methods newtons-method
person zeno    schedule 13.03.2020    source источник

Ответы (1)


arrow_upward
1
arrow_downward

Кажется, я не знал, как правильно разложить матрицу.

Ниже приведен рабочий пример метода W4 для двумерной системы.

#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 (std::abs(f1) > eps && std::abs(f2) > 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)*J(1,0)/J(1,1),   J(0,1),
                     0.0,                          J(1,1);
        J_U_inv = J_U_inv.inverse().eval();

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


        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;

    }
}
person zeno    schedule 14.03.2020