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ł , na krańcach którego funkcja ma różne znaki.
2. Funkcja jest określona w przedziale.
3. Funkcja jest ciągła w przedziale.
4. Przybliżenia początkowe x₀, x₁ i x₂ powinny należeć do przedziału.

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

W przedziale poszukiwań pierwiastka 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:

Dane wyjściowe

Zmienne pomocnicze

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 <windows.h>
#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

// Dokładność wyznaczania pierwiastka
double epsx = 1e-14;
// Dokładność wyznaczania zera
double epsy = 1e-14;
// Punkty startowe
double x0   = -1.1;
double x1   = -1.05;
double x2   = -1.0;
// Maksymalna liczba obiegów
int n       = 64;

// 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;

  SetConsoleOutputCP(CP_UTF8); 
  SetConsoleCP(CP_UTF8);

  cout << setprecision(15)
       << fixed;
  cout <<
  "Obliczanie przybliżonego pierwiastka "
  "funkcji metodą Mullera\n"
  "-------------------------------------"
  "----------------------\n\n";

  // 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 limit "
              "obiegów !!!\n";
      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 !!!\n";
      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!\n\n";

  cout
  << "Pierwiastek        x3 = "
  << setw(20) << x3 << endl
  << "Wartość funkcji f(x3) = "
  << setw(20) << f3 << endl
  << "Dokładność dla x epsx = "
  << setw(20) << epsx << endl
  << "Dokładność dla y epsy = "
  << setw(20) << epsy << endl
  << "Liczba obiegów      n = "
  << setw(4) << 64 - n;

  cout << endl << endl;
  system("pause");
  return 0;
}

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

Porównaj otrzymany wynik z trzema 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.077713513691340
Wartość funkcji f(x0) =    0.000000000000000
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

# Pierwiastek funkcji: metoda Mullera
# (C)2026 mgr Jerzy Wałaszek
# Metody numeryczne
# ------------------------------------

import math

# Tutaj definiujemy funkcję,
# której pierwiastek jest wyliczany
#----------------------------------
def f(x):
    return (math.sin(x * x - x + 1 / 3.0)
            + 0.5 * x)

# Parametry początkowe
epsx = 1e-14
epsy = 1e-14
x0 = -1.1
x1 = -1.05
x2 = -1.0
n = 64

# Program główny
result = False

print(
    "Obliczanie przybliżonego "
    "pierwiastka\n"
    "funkcji metodą Mullera\n"
    "-------------------------"
    "----------------------\n"
)

# Wartości startowe
f0 = f(x0)
f1 = f(x1)
f2 = f(x2)

# Pętla obliczeniowa
while True:
    n -= 1
    if not n:
        print("Przekroczono limit "
              "obiegów !!!")
        break

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

    # Współczynniki 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:
        print("Źle dobrane punkty "
              "startowe, "
              "brak pierwiastka "
              "rzeczywistego !!!")
        break

    d = math.sqrt(d)
    x31 = (-a1 - d) / 2 / a2
    x32 = (-a1 + d) / 2 / a2

    # Wybór pierwiastka bliższego x2
    if (math.fabs(x31 - x2) <
        math.fabs(x32 - x2)):
        x3 = x31
    else:
        x3 = x32

    # Wartość funkcji w x3
    f3 = f(x3)

    # Kryterium stopu
    if (math.fabs(f3) < epsy or 
        math.fabs(x3 - x2) < epsx):
        result = True
        break

    # Przypisanie nowych przybliżeń
    x0, f0 = x1, f1
    x1, f1 = x2, f2
    x2, f2 = x3, f3

if not result:
    print("Zakończono z błędem!\n")

print(f"Pierwiastek        x3 = "
      f"{x3:20.15f}\n"
      f"Wartość funkcji f(x3) = "
      f"{f3:20.15f}\n"
      f"Dokładność dla x epsx = "
      f"{epsx:20.15f}\n"
      f"Dokładność dla y epsy = "
      f"{epsy:20.15f}\n"
      f"Liczba obiegów      n = "
      f"{64 - n:4}")

print("\n")
input("Naciśnij Enter...")