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_0\), \(x_1\) i \(x_2\), 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_0\), \(x_1\) i \(x_2\) 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 wybieramy 3 punkty \(x_0\), \(x_1\) i \(x_2\):

Znajdujemy parabolę, która przechodzi przez 3 wybrane punkty wykresu \((x_0,f(x_0))\), \((x_1,f(x_1))\) i \((x_2,f(x_2))\):

Parabola ma wzór:

Musimy znaleźć wartości współczynników \(a_2\), \(a_1\) i \(a_0\). Tworzymy układ 3 równań:

Z równania 3) wyliczamy współczynnik \(a_0\):

Wyliczony współczynnik \(a_0\) wstawiamy do równań 1) i 2):

Z równania 2) wyliczamy współczynnik \(a_1\):

I wstawiamy go do równania 1) w celu wyliczenia współczynnika \(a_2\):

Aby uprościć rachunki, dokonamy podstawień:

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

Za \(x_3\) przyjmujemy ten z pierwiastków \(x_{31}\), \(x_{32}\) paraboli, który leży bliżej \(x_2\). Następnie sprawdzamy, czy wartość funkcji w punkcie \(x_3\) 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_3 - x_2|}\) jest dostatecznie mała. Jeśli nie, to za nowe trzy punkty przyjmujemy \(x_1\), \(x_2\) oraz wyliczony punkt \(x_3\) 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

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)2026 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.10;
double x1   = -1.05;
double x2   = -1.00;
// Maksymalna liczba obiegów
int max_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ą
  int n = max_n;
  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) << max_n - 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
max_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
n = max_n
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"{(max_n - n):4}")

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