Pierwiastki funkcji

Metoda Steffensena (ang. Steffensen's method) jest bardzo podobna do metody Newtona, jednak nie wymaga liczenia pochodnej funkcji. Warunki początkowe są takie same jak dla metody Newtona:

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.

Metoda Steffensena wykorzystuje do obliczania następnego przybliżenia pierwiastka funkcji podobny wzór jak metoda Newtona. Różnica jest taka, iż pochodna jest zastąpiona funkcją g(x). Zasada działania tej metody jest następująca:

Wybieramy przedział [a,b] spełniający warunki początkowe. W przedziale [a,b] wybieramy punkt startowy x₀. Następne przybliżenie wyliczamy wg wzoru:

We wzorze występuje funkcja g(x), której wartość wyliczamy ze wzoru:

Przekształćmy ten wzór:

Podstawmy h = f(x):

Funkcja g(x) jest średnim nachyleniem funkcji f(x) pomiędzy punktami wykresu (x,f(x)) i (x+h,f(x+h)). Pierwszy punkt x należy do ciągu aproksymacji pierwiastka. Drugi punkt x+h jest punktem pomocniczym. Funkcja ta zwana jest również różnicą dzieloną pierwszego rzędu funkcji f(x) (ang. first-order divided difference). Przyrost h nie jest stały i maleje wraz z kolejnymi przybliżeniami x do pierwiastka funkcji f(x).

Obliczenia przerywamy, podobnie jak w metodzie Newtona, w jednym z trzech przypadków:

1. Funkcja ma wartość dostatecznie bliską 0 w wyznaczonym punkcie x: | f(x) | < ε_y.
2. Odległość dwóch ostatnich aproksymacji pierwiastka spadnie poniżej założonej dokładności obliczeń: |x - x₀| < ε_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.

Jeśli żaden z tych przypadków nie wystąpił, to wyliczone x przepisujemy do x₀ i kontynuujemy wyliczanie następnego przybliżenia pierwiastka funkcji.

Metoda Steffensena wyznacza pierwiastek funkcji równie szybko jak metoda Newtona, jednakże wymaga ona dwukrotnego wyliczania wartości funkcji f(x) w każdym obiegu. Zaletą jednak jest to, iż nie ma potrzeby wyznaczania pochodnej funkcji f(x). Jeśli obliczanie funkcji f(x) jest czasochłonne, to w praktyce może się okazać, że metody teoretycznie wolniejsze (np. metoda siecznych) dadzą rozwiązanie szybciej, ponieważ są one mniej złożone obliczeniowo.

Algorytm metody Steffensena

Dane wejściowe:

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

Dane wyjściowe

x

–

przybliżony pierwiastek funkcji

Zmienne pomocnicze

h	–	wartość funkcji w punkcie x
g	–	wartość różnicy dzielonej
x1	–	przybliżenie pierwiastka z poprzedniego kroku

Lista kroków

K01:	n ← n - 1	Zmniejszamy licznik
K02:	Jeśli n = 0, to zakończ z błędem	Sprawdzamy trzeci warunek zakończenia
K03:	h ← f ( x )	Obliczamy wartość funkcji w punkcie x
K04:	Jeśli | h | < εy, to zakończ	Pierwszy warunek zakończenia
K05:		Obliczamy różnicę dzieloną
K06:	x1 ← x	Zapamiętujemy bieżące przybliżenie pierwiastka
K07:		Wyliczamy kolejne przybliżenie pierwiastka
K08:	Jeśli | x1 - x | < εx, to zakończ	Drugi warunek zakończenia
K09:	Idź do kroku K01	Kontynuujemy obliczenia z nowym x

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

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 Steffensena
// (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 parametry początkowe

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

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

int main()
{
    // Zmienne

    double x1,h,g;
    bool result = false;

    setlocale(LC_ALL,"");

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

       h = f(x);

       // Sprawdzamy, czy funkcja jest dostatecznie bliska zeru

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

       // Obliczamy różnicę dzieloną

       g = (f(x+h) - h) / h;

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

       x1 = x;

       // Obliczamy kolejne przybliżenie

       x -= h/g;

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

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

       // Kontynuujemy obliczenia
    }

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

    cout << "Pierwiastek         x = " << setw(18) << x << endl
         << "Wartość funkcji  f(x) = " << setw(18) << h << 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 pięcioma poprzednimi metodami:

Działanie metody Steffensena można przyspieszyć przez zastosowanie tzw. ekstrapolacji Aitkena.

Załóżmy, że mamy ciąg kolejnych przybliżeń pierwiastka funkcji:

Ciąg ten jest zbieżny do pierwiastka funkcji x_r:

Jeśli ciąg ten jest liniowo zbieżny do pierwiastka, to:

gdzie C to pewna wartość stała.

Dla odpowiednio dużego k możemy zapisać:

Wyliczmy przybliżoną wartość pierwiastka x_r:

Okazuje się, że ciąg o wyrazach:

Jest szybciej zbieżny do x_r od ciągu pierwotnego. Jest to tzw. ekstrapolacja Aitkena. Możemy ją połączyć z metodą Steffensena w sposób następujący:

Wybieramy punkt startowy x₀. Następnie za pomocą metody Steffensena wyliczamy dwa kolejne punkty x₁ i x₂. Punkty te są konieczne do wyznaczenia x₃ ekstrapolacją Aitkena. Punkt x₃ wyliczamy ze wzoru:

Sprawdzamy, czy funkcja w punkcie x₃ jest dostatecznie bliska 0. Jeśli tak, kończymy.

Sprawdzamy, czy różnica ostatnich dwóch aproksymacji jest dostatecznie mała. Jeśli tak, kończymy.

W przeciwnym razie za x₀ przyjmujemy wyliczony punkt x₃ i powtarzamy obliczenia. Przy złym wyborze x₀ algorytm może nie być zbieżny do rozwiązania. Zatem należy ograniczyć liczbę dozwolonych aproksymacji tak, jak w poprzednich metodach.

Algorytm jest następujący:

Algorytm ekstrapolacji Aitkena

Dane wejściowe:

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

Dane wyjściowe

x

–

przybliżony pierwiastek funkcji

Zmienne pomocnicze

h	–	wartość funkcji w punkcie x
g	–	wartość różnicy dzielonej
x1,x2	–	przybliżenia pierwiastka z poprzedniego kroku

Lista kroków

K01:	n ← n - 1	Zmniejszamy licznik
K02:	Jeśli n = 0, to zakończ z błędem	Sprawdzamy warunek zakończenia
K03:	h ← f ( x )	Wyliczamy wartość funkcji w punkcie x
K04:	Jeśli | h | < εy, to zakończ	Sprawdzamy warunek zakończenia
K05:	x0 ← x	Zapamiętujemy x
K06:		Obliczamy różnicę dzieloną
K07:		Wyliczamy kolejne przybliżenie pierwiastka
K08:	x1 ← x	Zapamiętujemy x
K09:	h ← f ( x )	Wyliczamy wartość funkcji w punkcie x
K10:	Jeśli | h | < εy, to zakończ	Sprawdzamy warunek zakończenia
K11:		Obliczamy różnicę dzieloną
K12:		Wyliczamy kolejne przybliżenie pierwiastka
K13:	x2 ← x	Zapamiętujemy x
K14:		Obliczamy kolejne x wg ekstrapolacji Aitkena
K15:	Jeśli | x - x2 | < εx, to zakończ	Sprawdzamy warunek zakończenia
K16:	Idź do kroku K01	Kontynuujemy z x ekstrapolowanym

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

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 Steffensena/Aitkena
// (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 parametry początkowe

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

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

int main()
{
    // Zmienne

    double x0,x1,x2,h,g;
    bool result = false;

    setlocale(LC_ALL,"");

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

       h = f(x);

       // Sprawdzamy, czy funkcja jest dostatecznie bliska zeru

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

       // Obliczamy różnicę dzieloną

       g = (f(x+h) - h) / h;

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

       x0 = x;

       // Obliczamy kolejne przybliżenie

       x -= h/g;

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

       x1 = x;

       // Obliczamy wartość funkcji w punkcie x

       h = f(x);

       // Sprawdzamy, czy funkcja jest dostatecznie bliska zeru

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

       // Obliczamy różnicę dzieloną

       g = (f(x+h) - h) / h;

       // Obliczamy kolejne przybliżenie

       x -= h/g;

       // Zapamietujemy x

       x2 = x;

       // Kolejne przybliżenie wyliczamy ekstrapolacją Aitkena

       x = x0 - (x1 - x0) * (x1 - x0) / (x2 - 2 * x1 + x0);

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

       if(fabs(x - x2) < epsx)
       {
           result = true;
           break;
       }

       // Kontynuujemy obliczenia z nowym punktem x
    }

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

    cout << "Pierwiastek         x = " << setw(18) << x << endl
         << "Wartość funkcji  f(x) = " << setw(18) << f(x) << 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 sześcioma poprzednimi metodami:

Każdy obieg aproksymacyjny Aitkena wymaga dwukrotnego wyliczania pierwiastka metodą Steffensona (dla x₁ i x₂). W efekcie koszt obliczeniowy obu metod jest podobny, pomimo iż liczba obiegów jest mniejsza. Z podanego porównania wynika, iż metoda siecznych daje całkiem przyzwoite wyniki w małej liczbie obiegów, dodatkowo jest bardzo prosta obliczeniowo.