Pierwiastki funkcji

Metoda Newtona (ang. Newton's method) polega na kolejnych przybliżeniach pierwiastka funkcji przez wyznaczanie przecięć stycznej do wykresu funkcji z osią OX. Z tego powodu, podobnie jak w metodzie siecznych, warunki początkowe są następujące:

Na osi OX wybieramy pewien przedział [a,b] o następujących własnościach:

Funkcja jest określona w przedziale [a,b]

Funkcja jest ciągła w przedziale [a,b]

Na krańcach przedziału [a,b] funkcja ma różne znaki

Dodatkowo w przedziale [a,b] pierwsza pochodna funkcji jest różna od zera. Warunek ten gwarantuje nam, iż w tym przedziale funkcja nie posiada minimum/maksimum lokalnego, a co za tym idzie jej styczna nie będzie równoległa do osi OX, ponieważ algorytm bazuje na znajdowaniu punktów przecięcia stycznych z osią OX.

Wyjdźmy od wzoru na punkty przecięcia siecznych z osią OX, który wyprowadziliśmy w metodzie siecznych:

Przekształćmy ten wzór następująco:

We wzorze otrzymaliśmy odwrotność ilorazu różnicowego funkcji f(). Policzmy granicę tego ułamka, gdy x_{n - 2} zmierza do x_{n - 1}:

W granicy otrzymamy odwrotność pierwszej pochodnej funkcji w punkcie poprzednio wyznaczonym. Otrzymujemy zatem wzór na kolejne przybliżenie pierwiastka funkcji:

Wzór ten nosi nazwę wzoru Newtona. Jest on bardzo prosty, jednakże wymaga znajomości pierwszej pochodnej funkcji.

Gdy punkt x_{n - 2} zbliża się do punktu x_{n - 1}, to sieczna staje się styczną do wykresu funkcji w punkcie x_{n - 1}. W metodzie Newtona potrzebujemy tylko jednego punktu startowego w przedziale [a,b]. Zasada działania tej metody jest następująca:

W przedziale [a,b] spełniającym warunki początkowe wybieramy dowolny punkt startowy x₀:

Za pomocą wzoru Newtona wyznaczamy punkt przecięcia stycznej w punkcie (x₀,f(x₀)) z osią OX:

Sprawdzamy teraz, czy wartość funkcji w wyliczonym punkcie x₁ jest dostatecznie bliska zeru. Jeśli tak, to kończymy. Jeśli nie, to całą procedurę powtarzamy dla x₁, wyznaczając kolejny punkt x₂:

Sprawdzamy, czy funkcja w punkcie x₂ jest dostatecznie bliska zeru. Jeśli tak, to kończymy. Jeśli nie, wyliczamy następne punkty x₃, x₄, ... , x_n, aż warunek będzie spełniony lub zostanie przekroczony limit aproksymacji.

Kolejne punkty x_i leżą coraz bliżej pierwiastka funkcji. Metoda Newtona jest zwykle szybko zbieżna. Odległości pomiędzy kolejnymi dwoma punktami x_{i - 1} i x_i maleje. Gdy spadnie poniżej dokładności wyznaczania pierwiastka, obliczenia można przerwać. Wynika z tego, iż w metodzie Newtona obliczenia kończymy w trzech przypadkach:

1. Funkcja ma wartość dostatecznie bliską 0 w wyznaczonym punkcie x_i: | f(x_i) | < ε_y.
2. Odległość dwóch kolejnych aproksymacji pierwiastka spadnie poniżej założonej dokładności obliczeń: | x_{i - 1} - x_i | < ε_x.
3. Przekroczono zadaną liczbę aproksymacji bez osiągnięcia wyniku. Jest to sytuacja błędna, która może powstać przy nieprawidłowym doborze punktu startowego.

Wadą metody Newtona jest konieczność wyliczania pierwszej pochodnej w każdym obiegu aproksymacji pierwiastka.

Algorytm metody Newtona

Dane wejściowe:

εx	–	dokładność przybliżenia pierwiastka
εy	–	dokładność przybliżenia do zera
x0	–	początkowe przybliżenia pierwiastka
f()	–	funkcja, której pierwiastek będzie obliczany
f'()	–	pochodna funkcji
n	–	maksymalna liczba aproksymacji

Dane wyjściowe

x0

–

przybliżony pierwiastek funkcji

Zmienne pomocnicze

x1	–	poprzednia wartość x0
f0	–	wartość funkcji w punkcie x0
f1	–	wartość pierwszej pochodnej w punkcie x0

Lista kroków

K01:	n ← n - 1	Zmniejszamy licznik obiegów
K02:	Jeśli n = 0, zakończ z błędem	Sprawdzamy trzeci warunek zakończenia
K03:	f0 ← f(x0)	Obliczamy wartość funkcji w punkcie x0
K04:	Jeśli | f0 | < εy, to zakończ	Sprawdzamy pierwszy warunek zakończenia
K05:	f1 ← f'(x0)	Obliczamy wartość pierwszej pochodnej
K06:	x1 ← x0	Zapamiętujemy x0
K07:		Wyznaczamy kolejne przybliżenie pierwiastka
K08:	Jeśli | x1 - x0 | < εx, zakończ	Sprawdzamy drugi warunek zakończenia
K09:	Idź do kroku K01	Kontynuujemy obliczenia

Poniższy program jest przykładową realizacją algorytmu Newtona. Program oblicza pierwiastek funkcji:

Wzór jej pochodnej jest następujący (nie wyprowadzamy tego wzoru, metody obliczania pochodnych powinieneś poznać na matematyce. Dobrą alternatywą jest opanowanie MathCAD'a):

Wykres tej funkcji w przedziale [-3,1] jest następujący:

Z wykresu wynika, iż pierwiastek znajduje się w przedziale [-1,1 , -1]. Jako przybliżenia początkowe wybieramy x₀ = -1.

// Pierwiastek funkcji - metoda Newtona
// (C)2019 mgr Jerzy Wałaszek
// Metody numeryczne
//------------------------------------------

#include <iostream>
#include <cmath>
#include <iomanip>

using namespace std;

// Tutaj definiujemy funkcję, której pierwiastek jest wyliczany
//-------------------------------------------------------------

double f(double x)
{
  return sin(x*x-x+1/3.0)+0.5*x;
}

// Tutaj definiujemy pochodną funkcji
//-----------------------------------

double df(double x)
{
  return cos(x*x-x+1/3.0)*(2*x-1)+0.5;
}

// Tutaj definiujemy parametry początkowe

double epsx = 1e-14; // Dokładność wyznaczania pierwiastka
double epsy = 1e-14; // Dokładność wyznaczania zera
double x0   = -1;    // Punkt startowy
int n       = 64;    // Maksymalna liczba obiegów

// Program główny
//---------------

int main()
{
    // Zmienne

    double f0,f1,x1;
    bool result = false;

    setlocale(LC_ALL,"");

    cout << setprecision(15) << fixed;
    cout << "Obliczanie przybliżonego pierwiastka funkcji metodą Newtona" << endl
         << "-----------------------------------------------------------" << endl << endl;
    while(--n)
    {
       // Obliczamy wartość funkcji w punkcie x0

       f0 = f(x0);

       // Sprawdzamy, czy funkcja jest dostatecznie bliska zeru

       if(fabs(f0) < epsy)
       {
           result = true;
           break;
       }

       // Obliczamy wartość pierwszej pochodnej funkcji

       f1 = df(x0);

       // Zapamiętujemy bieżące przybliżenie

       x1 = x0;

       // Obliczamy kolejne przybliżenie

       x0 -= f0/f1;

       // Sprawdzamy, czy odległość pomiędzy dwoma ostatnimi przybliżeniami
       // jest mniejsza od założonej dokładności

       if(fabs(x1 - x0) < epsx)
       {
           result = true;
           break;
       }

       // Kontynuujemy obliczenia
    }

    if(!result) cout << "Zakończono z błędem!" << endl << endl;

    cout << "Pierwiastek        x0 = " << setw(18) << x0 << endl
         << "Wartość funkcji f(x0) = " << setw(18) << f0 << endl
         << "Dokładność dla x epsx = " << setw(18) << epsx << endl
         << "Dokładność dla y epsy = " << setw(18) << epsy << endl
         << "Liczba obiegów      n = " << 64 - n;

    cout << endl << endl;

    return 0;
}

Porównaj otrzymany wynik z czterema poprzednimi metodami:

Z powyższego porównania wynika, że metoda Newtona przegrywa z metodą Mullera. Jednak dokładna analiza numeryczna pokazuje, że w przypadku ogólnym metoda Newtona jest najszybsza z opisanych wcześniej metod. Jej podstawową wadą jest konieczność obliczania wartości pierwszej pochodnej funkcji w każdym obiegu aproksymacyjnym.