Miejsca Zerowe Funkcji - Metoda Połowienia

W naszym serwisie jest nowszy artykuł o obliczaniu pierwiastków funkcji: "Metody numeryczne".

Metoda połowienia

Mamy daną funkcję f(x) oraz przedział ⟨a;b⟩, w którym będziemy poszukiwali miejsca zerowego (czyli pierwiastka funkcji f(x)). Aby można było zastosować algorytm połowienia (zwany również algorytmem bisekcji), w przedziale ⟨a;b⟩ muszą być spełnione poniższe warunki:

Funkcja f(x) jest określona - uczniowie często nie rozumieją tego pojęcia. Określoność funkcji oznacza, iż dla każdej wartości argumentu x z przedziału ⟨a;b⟩ potrafimy policzyć wartość funkcji. W trakcie pracy algorytm połowienia oblicza wartości funkcji dla argumentów należących do tego przedziału ⟨a;b⟩. Jeśli przypadkowo trafi na punkt, dla którego wartość funkcji jest nieokreślona, obliczenia się załamią. W praktyce konsekwencje mogą być tragiczne, np. zmiana lotu samolotu z poziomego na pionowy w kierunku ziemi...

Dla przykładu rozważmy prostą funkcję:

Ile wynosi wartość tej funkcji dla x = 1? Musimy dzielić przez 0, a jak wiadomo jest to zadanie niewykonalne. W punkcie x = 1 tak podana funkcja ma nieokreśloną wartość. Co gorsza punkt ten wypada akurat w środku przedziału ⟨a;b⟩ i algorytm bisekcji trafi na niego przy pierwszym obiegu.
Funkcja f(x) jest ciągła. Ciągłość funkcji oznacza z kolei, iż jej wartości nie "wykonują" nagłych skoków. Funkcja przebiega przez wszystkie wartości pośrednie - nie istnieją zatem przerwy w kolejnych wartościach funkcji. Dla przykładu rozważmy taką oto funkcję:

Funkcja w przedziale ⟨-2;1⟩ posiada następujący wykres:
Nieciągłość występuje w punkcie x = 0, czyli w miejscu zmiany przepisu funkcji. Zwróć uwagę, iż funkcja jest określona w tym punkcie - nieciągłość i nieokreśloność to dwie różne sprawy - nie myl ich!!!
Funkcja f(x) na krańcach przedziału ⟨a;b⟩ przyjmuje różne znaki. Ponieważ funkcja, zgodnie z poprzednim podpunktem, jest ciągła, to przyjmuje w przedziale ⟨a;b⟩, wszystkie wartości pośrednie pomiędzy f(a) i f(b). Wartości te mają różne znaki (czyli leżą po różnych stronach osi OX), zatem musi być taki punkt x₀ w przedziale ⟨a;b⟩, dla którego funkcja przyjmuje wartość pośrednią:

Gdy funkcja f(x) spełnia powyższe trzy warunki, to w przedziale ⟨a;b⟩ zagwarantowane jest istnienie pierwiastka i możemy go wyszukać algorytmem połowienia (bisekcji). Zasada jest następująca:
Wyznaczamy punkt x₀ jako środek przedziału ⟨a;b⟩ zgodnie ze wzorem:

Obliczamy wartość funkcji w punkcie x₀. Jeśli długość przedziału jest mniejsza od założonej dokładności wyliczeń pierwiastka, to wartość x₀ jest poszukiwanym przybliżeniem. Kończymy algorytm. W przeciwnym razie sprawdzamy, czy f(x₀) znajduje się dostatecznie blisko 0:

Jeśli nierówność jest spełniona, to x₀ jest poszukiwaną wartością pierwiastka. Zwracamy wynik i kończymy algorytm. W przeciwnym razie za nowy przedział poszukiwań pierwiastka przyjmujemy tą połówkę ⟨a;x₀⟩ lub ⟨x₀;b⟩,, w której funkcja zmienia znak na krańcach. Algorytm powtarzamy od początku dotąd, aż znajdziemy pierwiastek lub przedział ⟨a;b⟩ osiągnie założoną długość (może to być również epsilon). Wtedy kończymy zwracając ostatnio wyliczone x₀.

do podrozdziału do strony

Specyfikacja problemu

Dane wejściowe

f(x)	–	funkcja, której pierwiastek liczymy. Musi być ciągła i określona w przedziale poszukiwań pierwiastka.
a,b	–	punkty krańcowe przedziału poszukiwań pierwiastka funkcji f(x); a,b ∈ R

Dane wyjściowe

x₀

–

pierwiastek funkcji f(x)

Zmienne pomocnicze i funkcje

f_a , f_b, f₀	–	wartości funkcji odpowiednio w punktach a, b, x₀; f_a, f_b, f₀ ∈ R
ε₀	–	określa dokładność porównania z zerem; ε₀ = 0.0000000001
ε_x	–	określa dokładność wyznaczania pierwiastka x₀; ε_x = 0.0000000001

Lista kroków

	K01:	Czytaj a i b
	K02:	f_a ← f(a) ; f_b ← f(b)
	K03:	Jeśli f_a · f_b > 0, to pisz "Brak pierwiastka" i zakończ
	K04:	Dopóki \| a - b \| > ε_x: wykonuj kroki K05...K07
	K05:
	K06:	Jeśli \| f₀ \| < ε₀, to idź do kroku K08
	K07:	Jeśli f_a · f₀ < 0, to b ← x₀ inaczej a ← x₀; f_a ← f₀
	K08:	Pisz x₀i zakończ

Schemat blokowy

Wykonanie algorytmu rozpoczynamy od odczytu zakresu poszukiwań pierwiastka danej funkcji. W zmiennych f_a i f_b umieszczamy wartość funkcji na krańcach tego zakresu

Pierwszy test ma na celu sprawdzenie warunku różnych znaków wartości funkcji na krańcach zakresu poszukiwań pierwiastka. Różne znaki gwarantują nam istnienie pierwiastka w danym przedziale poszukiwań. Zatem jeśli funkcja na krańcach nie przybiera różnych znaków, wypisujemy odpowiedni komunikat i kończymy algorytm.

Rozpoczynamy pętlę wyliczania kolejnych przybliżeń pierwiastka funkcji. Pętla wykonuje się do momentu, aż przedział poszukiwań pierwiastka osiągnie długość ε_x.

Wewnątrz pętli wyznaczamy punkt x₀ leżący w środku przedziału ⟨a;b⟩ oraz obliczamy wartość funkcji w punkcie x₀ i umieszczamy ją w zmiennej f₀. Jeśli wartość f₀ jest dostatecznie bliska zeru (wpada w otoczenie zera o promieniu ε₀), przerywamy pętlę, wypisujemy x₀ i kończymy algorytm.

W przeciwnym razie za nowy przedział ⟨a;b⟩ przyjmujemy połówkę wyznaczoną przez x₀, w której funkcja zmienia znak. Testujemy iloczyn f_a przez f₀. Jeśli iloczyn jest ujemny, to różne znaki funkcja przyjmuje w połówce ⟨a;b⟩. Zatem za nowy koniec b przyjmujemy x₀ i kontynuujemy pętlę. W przeciwnym razie zmiana znaku występuje w drugiej połówce przedziału - ⟨a;b⟩. Za nowy początek a przyjmujemy x₀. Dodatkowo za f_a przyjmujemy f₀ - zwolni nas to od konieczności ponownego wyliczania wartości funkcji w nowym punkcie a. Po tych czynnościach kontynuujemy pętlę wyznaczania kolejnych przybliżeń pierwiastka x₀.

do podrozdziału do strony

Programy wyznaczają miejsce zerowe funkcji:

Pierwiastków należy poszukiwać w przedziałach ⟨-1;0⟩ i ⟨1;2⟩.

C++

// Program znajduje miejsce
// zerowe funkcji f(x) za
// pomocą algorytmu
// połowienia - bisekcji
// ------------------------
//  (C)2006 mgr J.Wałaszek

#include <iostream>
#include <iomanip>
#include <cmath>
#include <cstdlib>

using namespace std;

// dokładność porównania z zerem
const double EPS0 = 0.0000000001;
// dokładność wyznaczenia pierwiastka
const double EPSX = 0.0000000001;

// Funkcja, której miejsce
// zerowe obliczamy
// f(x) = x^3*(x+sin(x^2-1)-1)-1
// <-1,0> i <1,2>
double f(double x)
{
  return x * x * x *
        (x + sin(x * x - 1) - 1)
        - 1;
}

//---------------
// Program główny
//---------------

int main()
{

  double a,b,x0,fa,fb,f0;

  cout << setprecision(5)
       << fixed
       << "Obliczanie pierwiastka funkcji\n"
       << "       metoda bisekcji\n"
       << "f(x) = x^3*(x+sin(x^2-1)-1)-1\n"
       << "------------------------------\n"
       << "  (C)2006 mgr Jerzy Walaszek\n\n"
       << "Podaj przedział poszukiwania:\n\n";
  cout << "a = "; cin >> a;
  cout << "b = "; cin >> b;
  cout << "\nWYNIK:\n\n";
  fa = f(a); fb = f(b);
  if(fa * fb > 0)
    cout << "Brak pierwiastka\n";
  else
  {
    while(fabs(a - b) > EPSX)
    {
      x0 = (a + b) / 2; f0 = f(x0);
      if(fabs(f0) < EPS0) break;
      if(fa * f0 < 0) b = x0;
      else
      {
        a = x0; fa = f0;
      }
    }
    cout << "x0 = "
         << setw(15) << x0 << endl;
  }
  cout << endl;
  system("pause");
  return 0;
}

Pascal

// Program znajduje miejsce
// zerowe funkcji f(x) za
// pomocą algorytmu
// połowienia - bisekcji
// ------------------------
//  (C)2006 mgr J.Wałaszek

program mzf1;

uses math;

const
  // dokładność porównania z zerem
  EPS0 = 0.0000000001;
  // dokładność wyznaczenia pierwiastka
  EPSX = 0.0000000001;

// Funkcja, której miejsce
// zerowe obliczamy
// f(x) = x^3*(x+sin(x^2-1)-1)-1
// <-1,0> i <1,2>
// -----------------------------

function f(x : real) : double;
begin
  f := x * x * x * (x +
       sin(x * x - 1) - 1) - 1;
end;

//---------------
// Program główny
//---------------

var
  a,b,x0,fa,fb,f0 : double;

begin
  writeln('Obliczanie pierwiastka funkcji');
  writeln('       metoda bisekcji');
  writeln('f(x) = x^3*(x+sin(x^2-1)-1)-1');
  writeln('------------------------------');
  writeln('  (C)2006 mgr Jerzy Walaszek');
  writeln;
  writeln('Podaj przedzial poszukiwania:');
  writeln;
  write('a = '); readln(a);
  write('b = '); readln(b);
  writeln;
  writeln('WYNIK:');
  writeln;
  fa := f(a); fb := f(b);
  if fa * fb > 0 then
    writeln('Brak pierwiastka')
  else
  begin
    while abs(a - b) > EPSX do
    begin
      x0 := (a + b) / 2; f0 := f(x0);
      if abs(f0) < EPS0 then break;
      if fa * f0 < 0 then b := x0
      else
      begin
        a := x0; fa := f0;
      end;
    end;
    writeln('x0 = ',x0:15:8);
  end;
  writeln;
  writeln('Nacisnij Enter...');
  readln;
end.

Basic

' Program znajduje miejsce
'  zerowe funkcji f(x) za
'     pomocą algorytmu
'  połowienia  - bisekcji
' ------------------------
'  (C)2006 mgr J.Wałaszek

Declare Function f(x As Double) As Double

' dokładność porównania z zerem
Const EPS0 As Double = 0.0000000001
' dokładność wyznaczenia pierwiastka
Const EPSX As Double = 0.0000000001

'---------------
' Program główny
'---------------

Dim As double a, b, x0, fa, fb, f0

Print "Obliczanie pierwiastka funkcji"
Print "       metoda bisekcji"
Print "f(x) = x^3*(x+sin(x^2-1)-1)-1"
Print "------------------------------"
Print "  (C)2006 mgr Jerzy Walaszek"
Print
Print "Podaj przedzial poszukiwania:"
Print
Input "a = ", a
Input "b = ", b
Print
Print "WYNIK:"
Print
fa = f(a) : fb = f(b)
If fa * fb > 0 Then
  Print "Brak pierwiastka"
Else
  While Abs(a - b) > EPSX
    x0 = (a + b) / 2 : f0 = f(x0)
    If Abs(f0) < EPS0 Then Exit While
    If fa * f0 < 0 Then
      b = x0
    Else
      a = x0 : fa = f0
    End If
  Wend
  Print Using "x0 = ######.########"; x0
End If
Print
Print "Nacisnij Enter..."
Sleep
End

' Funkcja, której miejsce
' zerowe obliczamy
' f(x) = x^3*(x+sin(x^2-1)-1)-1
' <-1,0> i <1,2>
'------------------------------
Function f(x As Double) As Double
  f = x*x*x*(x+Sin(x*x-1)-1)- 1
End Function

Python (dodatek)

# Program znajduje miejsce
#  zerowe funkcji f(x) za
#     pomocą algorytmu
#  połowienia  - bisekcji
# ------------------------
#  (C)2026 mgr J.Wałaszek

import math

# Funkcja, której miejsce
# zerowe obliczamy
# f(x) = x^3*(x+sin(x^2-1)-1)-1
# <-1,0> i <1,2>
#------------------------------
def f(x):
    return x * x * x * \
           (x+math.sin(x*x-1)-1)- 1

# Dokładność porównania z zerem
EPS0 = 0.0000000001
# Dokładność wyznaczenia pierwiastka
EPSX = 0.0000000001

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

print("Obliczanie pierwiastka funkcji")
print("       metoda bisekcji")
print("f(x) = x^3*(x+sin(x^2-1)-1)-1")
print("------------------------------")
print("  (C)2026 mgr Jerzy Wałaszek")
print()
print("Podaj przedział poszukiwania:")
print()
a = float(input("a = "))
b = float(input("b = "))
print()
print("WYNIK:")
print()
fa, fb = f(a), f(b)
if fa * fb > 0:
    print("Brak pierwiastka")
else:
    while abs(a - b) > EPSX:
        x0 = (a + b) / 2
        f0 = f(x0)
        if abs(f0) < EPS0: break
        if fa * f0 < 0:
            b = x0
        else:
            a, fa = x0, f0
    print("x0 = %15.8f" % x0)
print()
input("Naciśnij Enter...")

Wynik:

Obliczanie pierwiastka funkcji
       metoda bisekcji        
f(x) = x^3*(x+sin(x^2-1)-1)-1 
------------------------------
  (C)2026 mgr Jerzy Wałaszek  

Podaj przedział poszukiwania: 

a = 1
b = 2

WYNIK:

x0 =      1.18983299

Naciśnij Enter...

JavaScript

<html>
  <head>
  </head>
  <body>
<div style="overflow-x: auto;"
     align="center">
  <table
  border="0"
  cellpadding="4"
  style="border-collapse:
         collapse">
    <tr>
      <td nowrap>
        <form
        name="frmbincode"
        style="text-align: center;
               background-color:
               #E7E7DA;">
          <b><big>
          Obliczanie<br>
          pierwiastka funkcji<br>
          metodą Bisekcji
          </big></b><br><br>
          f(x) = x<sup>3</sup>(x +
          sin(x<sup>2</sup>-1)-1)-1<br>
          <br>
          &nbsp;&nbsp;
          (C)2026
          mgr Jerzy Wałaszek
          I LO w Tarnowie
          &nbsp;&nbsp;
          <hr>
          Wpisz do pól edycyjnych<br>
          krańce przedziału<br>
          poszukiwań pierwiastka<br>
          <br>
          a = <input
              type="text"
              name="inp_a"
              size="16"
              value="1"
              style="text-align:
              right">
          <br>
          b = <input
              type="text"
              name="inp_b"
              size="16"
              value="2"
              style="text-align:
              right">
          <hr>
          <input
          type="button"
          value="Szukaj"
          name="B1"
          onclick="main()">
          <hr>
          <b>Wynik:</b>
          <div id="out">.</div>
        </form>
      </td>
    </tr>
  </table>
</div>

<script language=javascript>

// Program znajduje
// miejsce zerowe
// funkcji f(x)
// za pomocą algorytmu
// połowienia - bisekcji
//----------------------
// (C)2006 mgr J.Wałaszek
// I LO w Tarnowie

// dokładność porównania
// z zerem
var EPS0 = 0.0000000001
// dokładność wyznaczenia
// pierwiastka
var EPSX = 0.0000000001

// Funkcja, której miejsce
// zerowe obliczamy.
// f(x) = x^3*(x+sin(x^2-1)-1)-1
// Pierwiastki w <-1,0> i <1,2>
function f(x)
{
  return x * x * x * (x +
         Math.sin(x * x - 1) - 1) - 1;
}

//---------------
// Program główny
//---------------

function main()
{
  var a,b,x0,fa,fb,f0,t

  a = parseFloat(document
      .frmbincode.inp_a.value)
  b = parseFloat(document
      .frmbincode.inp_b.value)
  if(isNaN(a) || isNaN(b))
    t = "<font color=red><b>" +
        "Złe krańce " +
        "przedziału " + 
        "poszukiwań " +
        "pierwiastka!" +
        "</b></font>"
  else
  {
    t  = "x<sub>0</sub> = "
    fa = f(a); fb = f(b)
    if(fa * fb > 0)
      t = "<font color=red><b>" +
          "Brak pierwiastka" +
          "</b></font>"
    else
    {
      while(Math.abs(a-b)>EPSX)
      {
        x0 = (a + b) / 2
        f0 = f(x0)
        if(Math.abs(f0) < EPS0)
          break
        if(fa * f0 < 0)
          b = x0
        else
        {
          a = x0
          fa = f0
        }
      }
      t += x0
    }
  }
  t = t
  document.getElementById("out")
  .innerHTML = t
}
</script>
  </body>
</html>

Tutaj możesz przetestować działanie prezentowanego skryptu. Przykładowa funkcja posiada pierwiastki w przedziałach ⟨-1;0⟩ oraz ⟨1;2⟩.

do podrozdziału do strony

Metoda połowienia

Algorytm

Opis algorytmu

Specyfikacja problemu

Dane wejściowe

Dane wyjściowe

Zmienne pomocnicze i funkcje

Lista kroków

Schemat blokowy

Programy