Serwis Edukacyjny
w I-LO w Tarnowie
obrazek

Materiały dla uczniów liceum

  Wyjście       Spis treści       Wstecz       Dalej  

obrazek

Autor artykułu: mgr Jerzy Wałaszek

©2019 mgr Jerzy Wałaszek
I LO w Tarnowie

obrazek

Pierwiastki funkcji

Metoda Steffensena

SPIS TREŚCI
Podrozdziały

Warunki początkowe

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]

Określoność funkcji oznacza, że dla każdego argumentu x z przedziału [a,b] istnieje wartość funkcji. Warunek ten jest konieczny, ponieważ algorytm bisekcji wybiera punkty w przedziale [a,b] i oblicza dla nich wartość funkcji. Jeśli trafi na punkt nieokreśloności, w którym nie można policzyć wartości funkcji, to cała metoda załamie się.

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

Ciągłość funkcji gwarantuje, iż jej wartości nie wykonują nagłych skoków i dla dowolnych dwóch wartości funkcji w tym przedziale znajdziemy wszystkie wartości pośrednie.

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

Ten warunek wraz z poprzednim gwarantuje, że w przedziale [a,b] istnieje taki argument x0, dla którego funkcja ma wartość 0, która to wartość jest wartością pośrednią pomiędzy wartościami funkcji na krańcach przedziału [a,b]:

obrazek

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.

Na początek:  podrozdziału   strony 

Algorytm i program

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 x0. 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 - x0| < ε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 x0 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
x 1 przybliżenie pierwiastka z poprzedniego kroku

Lista kroków

K01: n ← n - 1 Zmniejszamy licznik
K02: Jeśli n = 0, 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:

obrazek

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

obrazek

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

Przykładowy program w języku C++
// 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;
}
Wynik
Obliczanie przybliżonego pierwiastka funkcji metodą Steffensena
---------------------------------------------------------------

Pierwiastek         x = -1.077713513691340
Wartość funkcji  f(x) = -0.000000000000001
Dokładność dla x epsx =  0.000000000000010
Dokładność dla y epsy =  0.000000000000010
Liczba obiegów      n = 7

Porównaj otrzymany wynik z pięcioma poprzednimi metodami:

Obliczanie przybliżonego pierwiastka funkcji metodą bisekcji
------------------------------------------------------------

Pierwiastek        x0 = -1.077713513691339
Wartość funkcji f(x0) =  0.000000000000002
Dokładność dla x epsx =  0.000000000000010
Dokładność dla y epsy =  0.000000000000010
Liczba obiegów      i = 44
 
Obliczanie przybliżonego pierwiastka funkcji metodą regula falsi
----------------------------------------------------------------

Pierwiastek x0        = -1.077713513691339
Wartość funkcji f(x0) =  0.000000000000001
Dokładność dla y epsy =  0.000000000000010
Liczba obiegów        = 9
 
Obliczanie przybliżonego pierwiastka funkcji metodą siecznych
-------------------------------------------------------------

Pierwiastek        xn = -1.077713513691340
Wartość funkcji f(xn) =  0.000000000000000
Dokładność dla x epsx =  0.000000000000010
Dokładność dla y epsy =  0.000000000000010
Liczba obiegów        = 5
 
Obliczanie przybliżonego pierwiastka funkcji metodą Mullera
-----------------------------------------------------------

Pierwiastek        x3 = -1.077713513691340
Wartość funkcji f(x3) = 0.000000000000000
Dokładność dla x epsx = 0.000000000000010
Dokładność dla y epsy = 0.000000000000010
Liczba obiegów      n = 4
 
Obliczanie przybliżonego pierwiastka funkcji metodą Newtona
-----------------------------------------------------------

Pierwiastek        x0 = -1.077713513691340
Wartość funkcji f(x0) = -0.000000000000001
Dokładność dla x epsx =  0.000000000001000
Dokładność dla y epsy =  0.000000000000010
Liczba obiegów      n =  5
Na początek:  podrozdziału   strony 

Ekstrapolacja Aitkena

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 xr:

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 xr:

obrazek

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

Jest szybciej zbieżny do xr 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 x0. Następnie za pomocą metody Steffensena wyliczamy dwa kolejne punkty x1 i x2. Punkty te są konieczne do wyznaczenia x3 ekstrapolacją Aitkena. Punkt x3 wyliczamy ze wzoru:

Sprawdzamy, czy funkcja w punkcie x3 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 x0 przyjmujemy wyliczony punkt x3 i powtarzamy obliczenia. Przy złym wyborze x0 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, 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: obrazek Obliczamy różnicę dzieloną
K07: obrazek 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: obrazek Obliczamy różnicę dzieloną
K12: obrazek 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:

obrazek

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

obrazek

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

Przykładowy program w języku C++
// 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;
}
Wynik
Obliczanie przybliżonego pierwiastka funkcji metodą Steffensena/Aitkena
-----------------------------------------------------------------------

Pierwiastek         x = -1.077713513691340
Wartość funkcji  f(x) =  0.000000000000000
Dokładność dla x epsx =  0.000000000000010
Dokładność dla y epsy =  0.000000000000010
Liczba obiegów      n = 4

Porównaj otrzymany wynik z sześcioma poprzednimi metodami:

Obliczanie przybliżonego pierwiastka funkcji metodą bisekcji
------------------------------------------------------------

Pierwiastek        x0 = -1.077713513691339
Wartość funkcji f(x0) =  0.000000000000002
Dokładność dla x epsx =  0.000000000000010
Dokładność dla y epsy =  0.000000000000010
Liczba obiegów      i = 44
 
Obliczanie przybliżonego pierwiastka funkcji metodą regula falsi
----------------------------------------------------------------

Pierwiastek x0        = -1.077713513691339
Wartość funkcji f(x0) =  0.000000000000001
Dokładność dla y epsy =  0.000000000000010
Liczba obiegów        = 9
 
Obliczanie przybliżonego pierwiastka funkcji metodą siecznych
-------------------------------------------------------------

Pierwiastek        xn = -1.077713513691340
Wartość funkcji f(xn) =  0.000000000000000
Dokładność dla x epsx =  0.000000000000010
Dokładność dla y epsy =  0.000000000000010
Liczba obiegów        = 5
 
Obliczanie przybliżonego pierwiastka funkcji metodą Mullera
-----------------------------------------------------------

Pierwiastek        x3 = -1.077713513691340
Wartość funkcji f(x3) = 0.000000000000000
Dokładność dla x epsx = 0.000000000000010
Dokładność dla y epsy = 0.000000000000010
Liczba obiegów      n = 4
 
Obliczanie przybliżonego pierwiastka funkcji metodą Newtona
-----------------------------------------------------------

Pierwiastek        x0 = -1.077713513691340
Wartość funkcji f(x0) = -0.000000000000001
Dokładność dla x epsx =  0.000000000001000
Dokładność dla y epsy =  0.000000000000010
Liczba obiegów      n =  5
 
Obliczanie przybliżonego pierwiastka funkcji metodą Steffensena
---------------------------------------------------------------

Pierwiastek         x = -1.077713513691340
Wartość funkcji  f(x) = -0.000000000000001
Dokładność dla x epsx =  0.000000000000010
Dokładność dla y epsy =  0.000000000000010
Liczba obiegów      n = 7

Każdy obieg aproksymacyjny Aitkena wymaga dwukrotnego wyliczania pierwiastka metodą Steffensona (dla x1 i x2). 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.

Na początek:  podrozdziału   strony 

Zespół Przedmiotowy
Chemii-Fizyki-Informatyki

w I Liceum Ogólnokształcącym
im. Kazimierza Brodzińskiego
w Tarnowie
ul. Piłsudskiego 4
©2019 mgr Jerzy Wałaszek

Materiały tylko do użytku dydaktycznego. Ich kopiowanie i powielanie jest dozwolone
pod warunkiem podania źródła oraz niepobierania za to pieniędzy.

Pytania proszę przesyłać na adres email: i-lo@eduinf.waw.pl

Serwis wykorzystuje pliki cookies. Jeśli nie chcesz ich otrzymywać, zablokuj je w swojej przeglądarce.
Informacje dodatkowe.