Метод Ньютона

Обзор

Метод Ньютона - один из самых популярных численных методов, и Бёрден и Фейрес даже называют его самым мощным методом, который используется для решения уравнения \(f(x) = 0\). Этот метод основан на разложении в ряд Тейлора функции \(f(x)\) относительно точки \(x_1\):

\begin{equation} f(x)=f\left(x_{1}\right)+\left(x-x_{1}\right) f^{\prime}\left(x_{1}\right)+\frac{1}{2 !}\left(x-x_{1}\right)^{2} f^{\prime \prime}\left(x_{1}\right)+\cdots \end{equation}

где \(f\), а также его производные первого и второго порядка, \(f'\) и \(f'\) вычисляются в \(x_1\). Если мы возьмем первые два члена разложения в ряд Тейлора, мы получим:

\begin{equation} f(x) \approx f\left(x_{1}\right)+\left(x-x_{1}\right) f^{\prime}\left(x_{1}\right). \end{equation}

Затем мы устанавливаем предыдущее выражение равным нулю (т.е. \(f(x) = 0\)), чтобы найти корень уравнения, которое дает нам:

\begin{equation} f\left(x_{1}\right)+\left(x-x_{1}\right) f^{\prime}\left(x_{1}\right)=0. \end{equation}

Переставляя предыдущее выражение, мы получаем следующее приближение к корню, что дает нам:

\begin{equation} x=x_{2}=x_{1}-\frac{f\left(x_{1}\right)}{f^{\prime}\left(x_{1}\right)} \end{equation}

Таким образом, обобщая предыдущее выражение, мы получаем итерационный метод Ньютона:

\begin{equation} x_{i}=x_{i-1}-\frac{f\left(x_{i-1}\right)}{f^{\prime}\left(x_{i-1}\right)}, i \in \mathbb{N} \end{equation}

где \(x_{i} \rightarrow \bar{x}\) (при \(i \rightarrow \infty\)), а x - приближение к корню функции f (x).

Примечание: Поскольку итерации начинают иметь одинаковые повторяющиеся значения, т.е. как \(x_i = x_{i+1} = \bar{x}\), это указывает на то, что \(f(x)\) сходится к \(\bar{x}\). Таким образом, \(x_i\) - это корень функции \(f(x)\).

Метод

Шаг 1:

Пусть \(\mathbf{x}^{(0)}=\left(x_{1}^{(0)}, x_{2}^{(0)}, \ldots, x_{n}^{(0)}\right)\) - заданный начальный вектор.

Шаг 2:

Вычислите \(J\left(\mathbf{x}^{(0)}\right)\) и \(\mathbf{F}\left(\mathbf{x}^{(0)}\right)\).

Шаг 3:

Теперь нам нужно вычислить вектор \(\mathbf{y}^{(0)}\), где

\begin{equation} \mathbf{y}= \left[ \begin{array}{c} y_{1} \\ y_{2} \\ \vdots \\ y_{n} \end{array} \right] \end{equation}

Чтобы найти \(\mathbf{y}^{(0)}\), мы решаем линейную систему :math:`Jleft(mathbf{x}^{(0)}right) mathbf{y}^{(0)}=-mathbf{F}left(mathbf{x}^{(0)}right)`используя метод исключения Гаусса.

Примечание: преобразовав систему на шаге 3, мы получим \(\mathbf{y}^{(0)}=-J\left(\mathbf{x}^{(0)}\right)^{-1} \mathbf{F}\left(\mathbf{x}^{(0)}\right)\). Значение этого состоит в том, что, поскольку math:mathbf{y}^{(0)}=-Jleft(mathbf{x}^{(0)}right)^{-1} mathbf{F}left(mathbf{x}^{(0)}right), мы можем заменить \(-J\left(\mathbf{x}^{(0)}\right)^{-1} \mathbf{F}\left(\mathbf{x}^{(0)}\right)\) в наша итерационная формула с \(\mathbf{y}^{(0)}\). Этот результат даст, что

\begin{equation} \mathbf{x}^{(k)}=\mathbf{x}^{(k-1)}-J\left(\mathbf{x}^{(k-1)}\right)^{-1} \mathbf{F}\left(\mathbf{x}^{(k-1)}\right)=\mathbf{x}^{(k-1)}-\mathbf{y}^{(k-1)} \end{equation}

Шаг 4:

Как только \(\mathbf{y}^{(0)}\) найден, мы можем приступить к завершению первой итерации, решив для \(\mathbf{x}^{(1)}\). Таким образом, используя результат шага 3, мы имеем

\begin{equation} \mathbf{x}^{(1)}=\mathbf{x}^{(0)}+\mathbf{y}^{(0)}= \left[\begin{array}{c} x_{1}^{(0)} \\ x_{2}^{(0)} \\ \vdots \\ x_{n}^{(0)} \end{array}\right] + \left[\begin{array}{c} y_{1}^{(0)} \\ y_{2}^{(0)} \\ \vdots \\ y_{n}^{(0)} \end{array}\right] \end{equation}

Шаг 5:

После того, как мы вычислили \(\mathbf{x}^{(1)}\), мы повторяем процесс снова, пока \(\mathbf{x}^{(k)}\) не сходится к \(\bar{x}\). Это указывает на то, что мы достигли решения \(\mathbf{F}(\mathbf{x})=\mathbf{0}\), где \(\bar{x}\) - решение системы.

Примечание: когда набор векторов сходится, норма \(\left\|\mathbf{x}^{(k)}-\mathbf{x}^{(k-1)}\right\|=0\). Это означает, что

\begin{equation} \left\|\mathbf{x}^{(k)}-\mathbf{x}^{(k-1)}\right\|=\sqrt{\left(x_{1}^{(k)}-x_{1}^{(k-1)}\right)^{2}+\cdots+\left(x_{n}^{(k)}-x_{n}^{(k-1)}\right)^{2}}=0 \end{equation}

Использование

представьте, что мы хотим решить следующую нелинейную систему уравнений:

\begin{equation}\left\{\begin{array}{l} f(x, y)=x^{2}+y^{2}-5 \\ g(x, y)=y-3 x+5 \end{array}\right.\end{equation}

тогда код будет выглядеть так:

// example_newton.cpp

#include <iostream>
#include "../src/numerary.hpp" // Numerary library

using namespace std;
using namespace numerary;

/* System to solve */
void system(double *x, double *fv, int n) {
    fv[0] = x[0]*x[0] + x[1]*x[1] - 5;
    fv[1] = x[1] - 3*x[0] + 5;
}

/* The main function */
int main() {

    const int maxiter = 100; // Maximum interations for set method (100 by default)
    const double eps = 1.e-7; // eps value for method (1.e-7 by default)
    double *x = new double[2], *fv = new double[2];
    short int is_solved;

    // Initial point
    x[0] = 1; x[1] = 2;

    is_solved = Numerary::nonlinear_systems_of_equations(system, x, fv, 2, "newton", eps, maxiter);

    if (is_solved == 1) {
        cout << "Solved!" << endl;
        cout << "x = " << x[0] << endl;
        cout << "y = " << x[1] << endl;
    } else {
        cout << "Method not allowed!";
    }

    delete[] x;
    delete[] fv;

    return 0;
}