Pierwiastki funkcji

Metoda Mullera ( ang. Muller's method ) jest ulepszeniem metody siecznych. Różnica polega głównie na tym, iż funkcja jest aproksymowana nie prostą, lecz parabolą. Dzięki temu osiąga się lepsze typowanie położenia pierwiastka funkcji, zatem pierwiastek ten zostaje znaleziony w mniejszej liczbie kroków. Dodatkową zaletą metody Mullera jest to, iż po pewnej modyfikacji pozwala ona znajdować pierwiastki funkcji zespolonych.

Danymi początkowymi są trzy wybrane współrzędne na osi OX x₀, x₁ i x₃, które powinny leżeć w pobliżu pierwiastka funkcji i które są wstępnymi przybliżeniami pierwiastka funkcji. Nie ma potrzeby, aby wartości funkcji dla tych współrzędnych miały różne znaki, jednakże powinny one leżeć  wewnątrz przedziału, na krańcach którego funkcja ma różne znaki, co gwarantuje istnienie w nim pierwiastka rzeczywistego. Zatem warunki początkowe można określić jako:

1. Istnieje przedział [a,b], na krańcach którego funkcja ma różne znaki.
2. Funkcja jest określona w przedziale [a,b].
3. Funkcja jest ciągła w przedziale [a,b].
4. Przybliżenia początkowe x₀, x₁ i x₂ powinny należeć do przedziału [a,b].

Ideą metody Mullera jest przybliżanie funkcji za pomocą paraboli, która przebiega przez wybrane trzy punkty wykresu funkcji.

W przedziale poszukiwań pierwiastka [a,b] wybieramy 3 punkty x₀, x₁ i x₂:

Znajdujemy parabolę, która przechodzi przez 3 wybrane punkty wykresu (x₀,f(x₀)), (x₁,f(x₁)) i (x₂,f(x₂)):

Parabola ma wzór:

Musimy znaleźć wartości współczynników wielomianu a₂, a₁ i a₀. Tworzymy układ 3 równań:

Z równania 3 wyliczamy współczynnik a₀:

Wyliczony współczynnik a₀ wstawiamy do równań 1 i 2:

Z równania 2 wyliczamy współczynnik a₁:

I wstawiamy go do równania 1 w celu wyliczenia współczynnika a₂:

Aby uprościć rachunki, dokonamy podstawień:

Mając współczynniki, obliczamy pierwiastki paraboli:

Za x₃ przyjmujemy ten z pierwiastków paraboli, który leży bliżej x₂. Następnie sprawdzamy, czy wartość funkcji w punkcie x₃ jest dostatecznie bliska 0. Jeśli tak, to kończymy. Kończymy również, jeśli różnica pomiędzy dwoma ostatnimi przybliżeniami pierwiastka (x₃ - x₂) jest dostatecznie mała. Jeśli nie, to za nowe trzy punkty przyjmujemy x₁, x₂ oraz wyliczony punkt x₃ i obliczenia powtarzamy.

Ponieważ może się zdarzyć, iż przy złym doborze punktów początkowych pierwiastek nie może być znaleziony metodą Mullera, to należy ograniczyć liczbę iteracji wewnątrz algorytmu. Również obliczenia przerywamy, jeśli parabola nie przecina osi OX.

Algorytm metody Mullera

Dane wejściowe:

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

Dane wyjściowe

x3

–

przybliżony pierwiastek funkcji

Zmienne pomocnicze

f0,f1,f2	–	wartości funkcji w punktach x0,x1,x2
dx		różnice x-ów
dy	–	różnice y-ów
h	–	iloraz różnicowy
a	–	współczynniki paraboli
x31,x32	–	pierwiastki paraboli
d	–	wyróżnik delta paraboli

Lista kroków

K01:	f0 ← f ( x0 ) f1 ← f ( x1 ) f2 ← f ( x2 )	Obliczamy wartości funkcji w podanych punktach startowych
K02:	Jeśli n = 0, to zakończ z błędem inaczej n ← n - 1	Sprawdzamy licznik obiegów
K03:	dx01 ← x0 - x1 dx02 ← x0 - x2 dx12 ← x1 - x2 dy02 ← f0 - f2 dy12 ← f1 - f2	Podstawienia
K04:		Współczynniki paraboli
K05:	Jeśli d < 0, to zakończ z błędem	Wyróżnik delta paraboli
K06:		Pierwiastki paraboli
K07:	Jeśli | x31 - x2 | < | x32 - x2 |, to x3 ← x31 inaczej x3 ← x32	Wybieramy pierwiastek bliższy x2
K08:	f3 ← f ( x3 )	Obliczamy wartość funkcji w punkcie x3
K09:	Jeśli ( | f3 | < εy ) lub ( | x3 - x2 | <εx ), to zakończ	Sprawdzamy warunki zakończenia obliczeń
K10:	x0 ← x1; f0 ← f1 x1 ← x2; f1 ← f2 x2 ← x3; f2 ← f3	Zapamiętujemy ostatnie 3 przybliżenia pierwiastka
K11:	Idź do kroku K02	Powtarzamy obliczenia

Poniższy program jest przykładową realizacją algorytmu Mullera. 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 krańce i środek tego przedziału.

// Pierwiastek funkcji - metoda Mullera
// (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 x0   = -1.1;  // Punkty startowe
double x1   = -1.05;
double x2   = -1.0;
int n       = 64;    // Maksymalna liczba obiegów

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

int main()
{
    // Zmienne

    double x31,x32,x3,f0,f1,f2,f3;
    double dx01,dx02,dx12,dy02,dy12,h12;
    double a0,a1,a2,d;
    bool result = false;

    setlocale(LC_ALL,"");

    cout << setprecision(15) << fixed;
    cout << "Obliczanie przybliżonego pierwiastka funkcji metodą Mullera" << endl
         << "-----------------------------------------------------------" << endl << endl;
    // Obliczamy wartości funkcji w punktach startowych

    f0 = f(x0);
    f1 = f(x1);
    f2 = f(x2);

    // Rozpoczynamy pętlę obliczeniową

    while(true)
    {
        if(!--n)
        {
           cout << "Przekroczono dozwoloną liczbę obiegów !!!" << endl;
           break;
        }

        // Podstawienia

        dx01 = x0 - x1;
        dx02 = x0 - x2;
        dx12 = x1 - x2;
        dy02 = f0 - f2;
        dy12 = f1 - f2;
        h12  = dy12 / dx12;

        // Obliczanie współczynników paraboli

        a2 = (dy02 / dx02 - h12) / dx01;
        a1 = h12 - a2 * (x1 + x2);
        a0 = f2 - x2 * (a2 * x2 + a1);

        // Pierwiastki paraboli

        d = a1 * a1 - 4 * a2 * a0;

        if (d < 0)
        {
            cout << "Źle dobrane punkty startowe, brak pierwiastka rzeczywistego !!!" << endl;
            break;
        }
        d = sqrt(d);
        x31 = (-a1 - d) / 2 / a2;
        x32 = (-a1 + d) / 2 / a2;

        // Wybieramy pierwiastek bliższy x2

        if(fabs(x31 - x2) < fabs(x32 - x2)) x3 = x31; else x3 = x32;

        // Obliczamy wartość funkcji w x3

        f3 = f(x3);

        // Sprawdzamy, czy osiągnięto wymaganą dokładność

        if((fabs(f3) < epsy) || (fabs(x3 - x2) < epsx))
        {
            result = true;
            break;
        }

        // Za nowe przybliżenia wybieramy trzy ostatnie

        x0 = x1; f0 = f1;
        x1 = x2; f1 = f2;
        x2 = x3; f2 = f3;
    }

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

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