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

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

1 голос
/ 14 марта 2020

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

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

бумага: https://arxiv.org/pdf/1809.04495.pdf и приложение: https://arxiv.org/pdf/1809.04358.pdf (раздел D.2 -> включает прилагаемую математику)

image

equation

Примечание: U, L - верхняя и нижняя три angular матрицы якобиана (матрица частичные производные).

Моя текущая реализация выглядит следующим образом (используется 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;

    }
}

1 Ответ

1 голос
/ 14 марта 2020

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

Ниже приведен рабочий пример метода 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;

    }
}
...