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_0\), 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_0\) z osią OX, obliczamy wartość funkcji w tym punkcie. Jeśli jest ona dostatecznie bliska zeru, to kończymy z wynikiem \(x_0\). Jeśli nie trafiliśmy, to punkt \(x_0\) 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_0\) 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_0\) jest wyznaczany bardziej trafnie niż zawsze w środku przedziału. Dlatego algorytm wymaga mniej kroków.

Podstawowym problemem jest tutaj znalezienie wzoru na \(x_0\). 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_0\):

Mamy wszystkie elementy.

Algorytm regula falsi

Dane wejściowe:

Dane wyjściowe

Zmienne pomocnicze

Lista kroków

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
#-----------------------------------------

from math import sin

# Tutaj definiujemy funkcję,
# której pierwiastek jest wyliczany
#----------------------------------
def f(x):
    return (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 abs(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}\n"
  f"Wartość funkcji f(x0) = {f(x0):20.15f}\n"
  f"Dokładność dla y EPSY = {EPSY:20.15f}\n"
  f"Liczba obiegów        = {64 - n:4}\n")

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