Pierwiastki funkcji

Metoda bisekcji (ang. bisection method) , zwana również metodą połowienia lub wyszukiwaniem binarnym pozwala stosunkowo szybko znaleźć pierwiastek dowolnej funkcji w zadanym przedziale poszukiwań <a;b>. Aby można było zastosować metodę bisekcji, funkcja musi spełniać kilka warunków początkowych:

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

Metoda bisekcji stosowana jest do przybliżonego rozwiązywania równania f(x) = 0 w przedziale <a;b>. Funkcja f(x) musi być określona i ciągła w tym przedziale oraz musi posiadać różne znaki na krańcach przedziału <a;b>, co gwarantuje istnienie przynajmniej jednego pierwiastka w tym przedziale. Rozwiązanie znajdowane jest za pomocą kolejnych przybliżeń. Z tego powodu należy określić dokładność, z którą chcemy otrzymać pierwiastek funkcji oraz dokładność wyznaczania samej funkcji.

W każdym przybliżeniu algorytm wyznacza środek x₀ przedziału <a,b> jako średnią arytmetyczną krańców. Następnie sprawdzane jest, czy odległość tego środka od krańców przedziału jest mniejsza od założonej dokładności wyliczania pierwiastka. Jeśli tak, to algorytm kończy pracę z wynikiem w x₀. Jeśli nie, to wyznaczana jest wartość funkcji w punkcie x₀ i sprawdza się, czy odległość tej wartości od 0 jest mniejsza od założonej dokładności wyznaczania funkcji. Jeśli tak, to algorytm kończy pracę z wynikiem w x₀.

W przeciwnym razie punkt x₀ dzieli przedział na dwie równe połowy: <a;x₀> i <x₀;b>. Algorytm za nowy przedział przyjmuję tę połówkę, w której funkcja zmienia znak na krańcach i kontynuuje wyznaczanie pierwiastka funkcji.

Z opisu wynika, iż w każdym obiegu szerokość przedziału maleje dwukrotnie. Dzięki temu pierwiastek jest wyznaczany coraz dokładniej. Po dziesięciu obiegach szerokość przedziału maleje 2¹⁰ = 1024 razy. Po 20 obiegach szerokość maleje 2²⁰ = 1048576 razy. Jest to postęp wykładniczy, zatem bardzo szybki.

Algorytm bisekcji

Dane wejściowe:

Dane wyjściowe

Zmienne pomocnicze

Lista kroków

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 bisekcji
// (C)2019 mgr Jerzy Wałaszek
// Metody numeryczne 0032
//--------------------------------------

#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ść x
double EPSX = 1e-14;
// 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;

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

  // Licznik obiegów pętli
  int i = 0;

  // 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ą bisekcji\n"
  "-------------------------------------"
  "-----------------------\n\n";

  // Obliczamy i zapamiętujemy
  // wartości funkcji
  // na krańcach przedziału
  fa = f(a);
  fb = f(b);

  // Sprawdzamy, czy na krańcach
  // przedziału 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)
    {
      // Zwiększamy licznik
      // obiegów pętli
      i++;

      // Wyznaczamy środek przedziału
      x0 = (a + b) / 2;

      // Sprawdzamy, czy szerokość
      // przedziału jest mniejsza
      // od zadanej dokładności
      if(fabs(a - x0) < EPSX)
        break; // Jeśli tak, to kończymy

      // 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ł przyjmujemy
      // tę z połówek <a;x0>, <x0;b>,
      // w której funkcja ma różne
      // znaki na krańcach
      if(fa * fx < 0)
        b = x0;
      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 x EPSX = "
    << setw(20) << EPSX << endl
    << "Dokładność dla y EPSY = "
    << setw(20) << EPSY << endl
    << "Liczba obiegów      i = "
    << setw(4) << i;

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

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

import math

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

# 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ść x
EPSX = 1e-14
# Dokładność y
EPSY = 1e-14
# Początek przedziału
# poszukiwań pierwiastka
a = -1.1
# Koniec przedziału
# poszukiwań pierwiastka
b = -1.0

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

# Licznik obiegów pętli
i = 0

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

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

# Obliczamy i zapamiętujemy
# wartości funkcji
# na krańcach przedziału
fa = f(a)
fb = f(b)

# Sprawdzamy, czy na krańcach
# przedziału 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:

        # Zwiększamy licznik
        # obiegów pętli
        i += 1

        # Wyznaczamy środek przedziału
        x0 = (a + b) / 2

        # Sprawdzamy, czy szerokość
        # przedziału jest mniejsza
        # od zadanej dokładności
        if math.fabs(a - x0) < EPSX:
            break

        # 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

        # Za nowy przedział przyjmujemy
        # tę z połówek <a,x0>, <x0,b>,
        # w której funkcja ma różne
        # znaki na krańcach
        if fa * fx < 0:
            b = x0
        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 x EPSX = {EPSX:20.15f}")
    print(
  f"Dokładność dla y EPSY = {EPSY:20.15f}")
    print(
  f"Liczba obiegów      i = {i:4}")

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