Pierwiastki funkcji

Algorytm regula falsi ( łac. fałszywa prosta) (ang. false position method ) posiada takie same warunki początkowe jak opisany wcześniej algorytm bisekcji:

Funkcja musi być określona w przedziale [a,b]

Funkcja musi być ciągła w przedziale [a,b]

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

Nazwa "fałszywej prostej" bierze się stąd, iż w coraz węższych przedziałach wykres funkcji zaczyna przypominać wykres prostej:

Naszą funkcją jest:

W przedziale [-3, 1] funkcja posiada wykres:

Zmniejszamy przedział do [-1,5 , -1]. Wykres jest następujący:

Zwróć uwagę, iż w tym przedziale wykres jest bardzo zbliżony do wykresu prostej. Ponieważ nie jest to dokładna prosta, stąd nazwa "fałszywa prosta". W metodzie traktujemy funkcję tak, jakby była funkcją, której wykres jest prostą. Prosta ta przecina oś OX w punkcie, który jest przybliżeniem pierwiastka. Wyliczamy punkt przecięcia fałszywej prostej x₀ z osią OX, obliczamy wartość funkcji w tym punkcie. Jeśli jest ona dostatecznie bliska zeru, to kończymy z wynikiem x₀. Jeśli nie trafiliśmy, to punkt x₀ dzieli przedział poszukiwań pierwiastka na dwie części ( tutaj nie będą to zwykle połówki jak w metodzie bisekcji ). Za nowy przedział przyjmujemy tę część, w której funkcja ma różne znaki na krańcach.

Dlaczego metoda regula falsi jest zwykle szybsza od metody bisekcji? Otóż wykorzystywany jest tutaj przebieg funkcji. W metodzie bisekcji x₀ było wyznaczane zawsze w środku przedziału poszukiwań pierwiastka. Tutaj punkt ten jest wyliczany jako punkt przecięcia fałszywej prostej z osią OX. Gdy przedział maleje, prosta ta zaczyna coraz bardziej przypominać faktyczny przebieg funkcji, a zatem x₀ jest wyznaczany bardziej trafnie niż zawsze w środku przedziału. Dlatego algorytm wymaga mniej kroków.

Podstawowym problemem jest tutaj znalezienie wzoru na x₀. Wbrew pozorom jest to łatwe pod warunkiem, iż znasz geometrię.

Fałszywą prostą prowadzimy zawsze, tak aby przechodziła przez punkty krańcowe wykresu: a,f(a)) i (b,f(b). Popatrz na rysunek powyżej. Powstają tutaj dwa trójkąty podobne:

p – podstawa trójkąta niebieskiego w – wysokość trójkąta niebieskiego p' – podstawa trójkąta żółtego w' – wysokość trójkąta żółtego

Długości boków trójkątów są następujące:

Ponieważ trójkąty są podobne, to odpowiadające sobie boki spełniają proporcje:

Podstawiamy i wyliczamy x₀:

Mamy wszystkie elementy.

Algorytm regula falsi

Dane wejściowe:

εy	–	dokładność obliczania funkcji
a	–	początek przedziału poszukiwań pierwiastka
b	–	koniec przedziału poszukiwań pierwiastka
f()	–	funkcja, której pierwiastek będzie obliczany
n	–	maksymalna liczba aproksymacji

Dane wyjściowe

x0

–

przybliżony pierwiastek funkcji

Zmienne pomocnicze

fa	–	wartość funkcji w punkcie a
fb	–	wartość funkcji w punkcie b
fx	–	wartość funkcji w punkcie x0

Lista kroków

K01:	fa ← f ( a ) fb ← f ( b )	Obliczamy wartości funkcji na krańcach przedziału [a,b].
K02:	Jeśli fa · fb > 0, to Zakończ z błędem	Funkcja nie ma różnych znaków na krańcach przedziału [a,b].
K03:		Wyznaczamy x0.
K04:	Jeśli n = 0, to zakończ inaczej n ← n - 1	Przekroczono dozwoloną liczbę iteracji. Wynik w x0.
K05:	fx ← f ( x0 )	Obliczamy wartość funkcji w środku przedziału
K06:	Jeśli | fx | < εy, to zakończ	Sprawdzamy, czy funkcja ma wartość 0 w środku przedziału. Wynik w x0.
K07:	Jeśli fx · fa < 0, to b ← x0, fb ← fx Inaczej a ← x0, fa ← fx	Za nowy przedział [a,b] przyjmujemy tę połówkę [a,x0], [x0,b], w której funkcja zmienia znak
K08:	Idź do kroku K03	Kontynuujemy obliczenia

Zwróć uwagę, że algorytm regula falsi jest praktycznie taki sam jak algorytm bisekcji. Różnice występują w krokach 3, 4 i 7.

Poniższy program jest przykładową realizacją algorytmu bisekcji. 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]. Taki więc przedział został przyjęty w programie.

// Pierwiastek funkcji - metoda regula falsi
// (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 epsy = 1e-14; // Dokładność y
double a    = -1.1;  // Początek przedziału poszukiwań pierwiastka
double b    = -1.0;  // Koniec przedziału poszukiwań pierwiastka
int n       = 64;    // Maksymalna liczba obiegów

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

int main()
{
    // Zmienne

    double fa,fb,fx,x0;

    // Licznik obiegów pętli

    // Zmienna informująca o znalezieniu wyniku

    bool result = false;

    setlocale(LC_ALL,"");

    cout << setprecision(15) << fixed;
    cout << "Obliczanie przybliżonego pierwiastka funkcji metodą regula falsi" << endl
         << "----------------------------------------------------------------" << endl << endl;

    // Obliczamy i zapamiętujemy wartości funkcji na krańcach przedziału [a,b]

    fa = f(a);
    fb = f(b);

    //Sprawdzamy, czy na krańcach przedziału [a,b] wartości funkcji mają różne znaki

    if(fa * fb > 0) cout << "BŁĄD!!! Funkcja nie ma różnych znaków na krańcach przedziału";
    else
    {
        result = true;

        // W pętli wyznaczamy kolejne przybliżenia pierwiastka

        while(true)
        {
            // Wyznaczamy x0

            x0 = (fa * b - fb * a) / (fa - fb);

            // Sprawdzamy, wykonano maksymalną liczbę iteracji.

            if(!n) break; // Jeśli tak, to kończymy
            else n--;

            // Obliczamy i zapamiętujemy wartość funkcji w punkcie x0

            fx = f(x0);

            // Sprawdzamy, czy wartość funkcji jest dostatecznie bliska zeru

            if(fabs(fx) < epsy) break; // Jeśli tak, to kończymy

            // Za nowy przedział [a,b] przyjmujemy tą z części [a,x0], [x0,b],
            // w której funkcja ma różne znaki na krańcach

            if(fa * fx < 0)
            {
                b = x0;
                fb = fx;
            }
            else
            {
                a = x0;
                fa = fx;
            }
        }
    }

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

    cout << endl << endl;

    return 0;
}

Zwróć uwagę na liczbę obiegów oraz wartość funkcji dla wyliczonego pierwiastka. Porównaj je z wartościami otrzymanymi w metodzie bisekcji:

Metoda regula falsi ze względu na swoją prostotę jest dobrym wyborem przy obliczaniu pierwiastków funkcji.