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

©2020 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
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:

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,
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: 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
©2020 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.