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>

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 musi być 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 musi mieć różne znaki

Ten warunek wraz z poprzednim gwarantuje, że w przedziale <a;b> istnieje taki argument x₀, 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>:

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:

Dane wyjściowe

Zmienne pomocnicze

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

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

// Program główny
//---------------
int main()
{
  // Zmienne
  double fa,fb,fx,x0;

  // Zmienna informująca
  // o znalezieniu wyniku
  bool result = false;

  SetConsoleOutputCP(CP_UTF8); 
  SetConsoleCP(CP_UTF8);

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

  // 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, czy 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(20) << x0 << endl
    << "Wartość funkcji f(x0) = "
    << setw(20) << f(x0) << endl
    << "Dokładność dla y epsy = "
    << setw(20) << epsy << endl
    << "Liczba obiegów        = "
    << setw(4) << 64 - n;

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

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

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

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

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

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

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

# Tutaj definiujemy
# parametry początkowe

# Dokładność y
epsy = 1e-14
# Początek przedziału
# poszukiwań pierwiastka
a = -1.1
# Koniec przedziału
# poszukiwań pierwiastka
b = -1.0
# Maksymalna liczba obiegów
n = 64

# Program główny
#---------------

# Zmienna informująca
# o znalezieniu wyniku
result = False

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

# 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:
    print(
      "BŁĄD!!! Funkcja nie ma różnych "
      "znaków na krańcach przedziału")
else:
    result = True
    x0 = 0.0

    # W pętli wyznaczamy kolejne
    # przybliżenia pierwiastka
    while True:
        # Wyznaczamy x0
        x0 = (fa * b - fb * a) / (fa - fb)

        # Sprawdzamy, czy wykonano
        # maksymalną liczbę iteracji.
        if not n:
            break  # Jeśli tak, to kończymy
        else:
            n -= 1

        # Obliczamy i zapamiętujemy
        # wartość funkcji w punkcie x0
        fx = f(x0)

        # Sprawdzamy, czy wartość funkcji
        # jest dostatecznie bliska zeru
        if math.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:
    print(f"Pierwiastek        x0 = {x0:20.15f}")
    print(f"Wartość funkcji f(x0) = {f(x0):20.15f}")
    print(f"Dokładność dla y epsy = {epsy:20.15f}")
    print(f"Liczba obiegów        = {64 - n:4}")

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