Miejsca Zerowe Funkcji - Metoda Regula Falsi

Mamy daną funkcję f(x) oraz przedział [a,b] poszukiwań pierwiastka. W przedziale tym funkcja musi spełniać następujące warunki:

Funkcja f(x) jest określona - dla każdej wartości argumentu x z przedziału [a,b] potrafimy policzyć wartość funkcji.
Funkcja f(x) jest ciągła - 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.
Funkcja f(x) na krańcach przedziału [a,b] przyjmuje różne znaki. Ponieważ funkcja, zgodnie z poprzednim wymogiem, 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_o w przedziale [a,b], dla którego funkcja przyjmuje wartość pośrednią:

Gdy funkcja f(x) spełnia podane warunki, to w przedziale [a,b] zagwarantowane jest istnienie pierwiastka i możemy go wyszukać algorytmem regula falsi.

obrazek

W języku łacińskim regula falsi oznacza fałszywą prostą. Ideą tej metody jest założenie, iż funkcja w coraz mniejszych przedziałach wokół pierwiastka zaczyna przypominać funkcję liniową. Skoro tak, to przybliżenie pierwiastka otrzymujemy prowadząc linię prostą (sieczną) z punktów krańcowych przedziału. Sieczna przecina oś OX w punkcie x_o, który przyjmujemy za przybliżenie pierwiastka - gdyby funkcja faktycznie była liniowa, otrzymany punkt x_o byłby rzeczywistym pierwiastkiem.

Wzór dla x_o można wyprowadzić na kilka sposobów. My wybierzemy znane twierdzenie Talesa, które mówi, iż jeżeli ramiona kąta przetniemy dwiema prostymi równoległymi, to długości odcinków wyznaczonych przez te proste na jednym ramieniu kąta będą proporcjonalne do długości odpowiednich odcinków wyznaczonych przez te proste na ramieniu drugim. Poniższy rysunek obrazuje tę sytuację:

obrazek

W naszym przypadku postępujemy następująco:

Kąt utworzą odpowiednie odcinki:

pionowo FA = (a,0)-(a,f(a)) powiększony o odcinek FB = (b,0)-(b,-f(b))

poziomo XAB = (a,0)-(b,0)

Prostymi równoległymi będzie cięciwa z punktów krańcowych przedziału, oraz ta sama cięciwa przesunięta pionowo w górę o długość odcinka FB.

Poniższy rysunek obrazuje otrzymaną sytuację:

obrazek

Zgodnie z twierdzeniem Talesa mamy:

Jeśli podstawimy do tego wzoru długości odcinków:

Otrzymamy:

a dalej:

Ostatnie przekształcenie ma na celu otrzymanie wzoru o lepszej "zapamiętywalności". Mnożymy mianownik przez (-1), dzięki czemu staje się on spójny z licznikiem ułamka. Sam ułamek zmienia znak na minus.

Algorytm regula falsi jest bardzo podobny do opisanego w poprzednim rozdziale algorytmu bisekcji. Założenia wstępne dla badanej funkcji w obu algorytmach są identyczne. Różnią się one sposobem wyznaczania punktu x_o. W algorytmie bisekcji punkt ten zawsze wyznaczany był w środku przedziału [a,b]. Z tego powodu algorytm bisekcji jest "nieczuły" na przebieg funkcji - zawsze zbliża się do pierwiastka w ten sam sposób.

W algorytmie regula falsi jest inaczej. Punk x_o wyznaczany jest w zależności od wartości funkcji na krańcach przedziału poszukiwań. W efekcie w każdej iteracji otrzymujemy lepsze przybliżenie do rzeczywistej wartości pierwiastka. Zatem w większości przypadków algorytm regula falsi szybciej osiągnie założoną dokładność pierwiastka od algorytmu bisekcji (chociaż można oczywiście podać kontrprzykłady).

Po wyznaczeniu przybliżonego pierwiastka postępowanie w obu algorytmach jest w zasadzie takie samo. Sprawdzamy, czy wartość modułu różnicy pomiędzy dwoma ostatnimi przybliżeniami pierwiastka jest mniejsza od zadanego minimum. Jeśli tak, obliczenia kończymy zwracając x_o.

Obliczamy wartość funkcji w punkcie x_o i sprawdzamy, czy jest ona dostatecznie bliska zeru. Jeśli tak, zwracamy x_o i kończymy. Jeśli nie, za nowy przedział przyjmujemy tą część przedziału [a,b] rozdzielonego przez x_o, w której funkcja zmienia znak. Całą procedurę powtarzamy, aż do osiągnięcia pożądanego wyniku.

Na początek: podrozdziału strony

Specyfikacja problemu

Dane wejściowe

f(x)	–	funkcja, której pierwiastek liczymy. Musi spełniać warunki podane na początku tego rozdziału.
a,b	–	punkty krańcowe przedziału poszukiwań pierwiastka funkcji f(x). a,b ∈ R

Dane wyjściowe

x_o

–

pierwiastek funkcji f(x)

Zmienne pomocnicze i funkcje

f_a , f_b , f_o	– wartości funkcji odpowiednio w punktach a, b, x_o. f_a , f_b , f_o ∈ R
x₁	– przechowuje poprzednią wartość przybliżenia x_o
ε_o	– określa dokładność porównania z zerem.	ε_o = 0.0000000001
ε_x	– określa dokładność wyznaczania pierwiastka x_o.	ε_x = 0.0000000001

Lista kroków

	K01:	Czytaj a i b
	K02:	f_a ← f (a); f_b ← f (b); x₁ ← a; x_o ← b
	K03:	Jeśli f_a · f_b > 0, to pisz "Funkcja nie spełnia założeń" i zakończ
	K04:	Dopóki \| x₁ - x_o \| > ε_x: wykonuj kroki K05 ... K08
	K05:	x₁ ← x_o
	K06:
	K07:	Jeśli \| f_o \| < ε_o, to idź do kroku K09
	K08:	Jeśli f_a · f_o < 0, to b ← x_o; f_b ← f_o, inaczej a ← x_o; f_a ← f_o
	K09:	Pisz x_o
	K10:	Zakończ

Schemat blokowy

obrazek

Algorytm regula falsi jest bardzo podobny do algorytmu bisekcji. Zmiany zaznaczono na schemacie blokowym zielonym kolorem elementów.

Odczytujemy zakres poszukiwań pierwiastka danej funkcji. W zmiennych f_a i f_b umieszczamy wartość funkcji na krańcach tego zakresu. Inicjujemy x₁ oraz x_o, tak aby algorytm nie zatrzymał się przy pierwszym obiegu pętli.

Pierwszy test ma na celu sprawdzenie warunku różnych znaków wartości funkcji na krańcach przedziału [a,b]. Różne znaki gwarantują nam istnienie pierwiastka w danym zakresie. Jeśli funkcja na krańcach ma ten sam znak, wypisujemy odpowiedni komunikat i kończymy algorytm.

Rozpoczynamy pętlę wyliczania kolejnych przybliżeń pierwiastka funkcji. Pętla wykonywana jest do momentu, gdy dwa ostatnio obliczone pierwiastki różnią się o ε_x.

Na początku pętli do x₁ wprowadzamy poprzednio wyliczony pierwiastek i obliczamy nowy. Następne obliczamy wartość funkcji w x_o umieszczając ją w zmiennej f_o. Sprawdzamy, czy wartość f_o jest dostatecznie bliska zeru (wpada w otoczenie zera o promieniu ε_o), Jeśli tak, to przerywamy pętlę, co spowoduje wypisanie x_o i zakończenie algorytmu.

Za nowy przedział [a,b] przyjmujemy tę część przedziału podzielonego przez x_o, w której funkcja zmienia znak. Testujemy iloczyn f_a przez f_o. Jeśli jest ujemny, to funkcja zmienia znak przyjmuje w [a,x_o]. Za nowy koniec b przyjmujemy x_o i kontynuujemy pętlę. W przeciwnym razie zmiana znaku występuje w [x_o,b]. Za nowy początek a przyjmujemy x_o i kontynuujemy pętlę. W obu przypadkach przepisujemy również f_o odpowiednio do f_a lub f_b w celu uniknięcia ponownych obliczeń wartości funkcji w nowych punktach krańcowych przedziału.

Na początek: podrozdziału 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 regula falsi
//---------------------------------------------
// (C)2006 mgr J.Wałaszek       I LO w Tarnowie

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

using namespace std;

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

// 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,x1,fa,fb,f0;

  cout << setprecision(8)     // 8 cyfr po przecinku
       << fixed;              // format stałoprzecinkowy

  cout << "Obliczanie pierwiastka funkcji - metoda regula falsi\n"
          "f(x) = x^3*(x+sin(x^2-1)-1)-1\n"
          "----------------------------------------------------\n"
          "(C)2006 mgr Jerzy Walaszek           I LO w Tarnowie\n\n"
          "Podaj zakres poszukiwan pierwiastka:\n\n";
  cout << "a = "; cin >> a;
  cout << "b = "; cin >> b;
  cout << "\n----------------------------------------------------\n\n"
          "WYNIK:\n\n";
  fa = f(a); fb = f(b); x1 = a; x0 = b;
  if(fa * fb > 0) cout << "Funkcja nie spelnia zalozen\n";
  else
  {
    while(fabs(x1 - x0) > EPSX)
    {
      x1 = x0;
      x0 = a - fa * (b - a) / (fb - fa);
      f0 = f(x0);
      if(fabs(f0) < EPS0) break;
      if(fa * f0 < 0)
      {
        b = x0; fb = f0;
      }
      else
      {
        a = x0; fa = f0;
      }
    }
    cout << "x0 = " << setw(15) << x0 << endl;
  }
  cout << "\n------------------------------------------------\n\n";
  system("pause");
  return 0;
}

Pascal

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

program mzf1;

uses math;

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

// 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
  Result := x * x * x * (x + sin(x * x - 1) - 1) - 1;
end;

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

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

begin
  writeln('Obliczanie pierwiastka funkcji - metoda regula falsi');
  writeln('f(x) = x^3*(x+sin(x^2-1)-1)-1');
  writeln('----------------------------------------------------');
  writeln('(C)2006 mgr Jerzy Walaszek           I LO w Tarnowie');
  writeln;
  writeln('Podaj zakres poszukiwan pierwiastka:');
  writeln;
  write('a = '); readln(a);
  write('b = '); readln(b);
  writeln;
  writeln('----------------------------------------------------');
  writeln('WYNIK:');
  writeln;
  fa := f(a); fb := f(b); x1 := a;  x0 := b; 
  if fa * fb > 0 then writeln('Funkcja nie spelnia zalozen')
  else
  begin
    while abs(x1 - x0) > EPSX do
    begin
      x1 := x0;
      x0 := a - fa * (b - a) / (fb - fa);
      f0 := f(x0); 
      if abs(f0) < EPS0 then break;
      if fa * f0 < 0 then
      begin
        b := x0; fb := f0;
      end
      else
      begin
        a := x0; fa := f0;
      end;
    end;
    writeln('x0 = ',x0:15:8);
  end;
  writeln;
  writeln('----------------------------------------------------');
  writeln('Koniec. Nacisnij klawisz Enter...');
  readln;
end.

Basic

' Program znajduje miejsce zerowe funkcji f(x)
' za pomocą algorytmu regula falsi
'---------------------------------------------
' (C)2006 mgr J.Wałaszek       I LO w Tarnowie

Declare Function f(x As Double) As Double

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

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

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

Print "Obliczanie pierwiastka funkcji - metoda regula falsi"
Print "f(x) = x^3*(x+sin(x^2-1)-1)-1"
Print "----------------------------------------------------"
Print "(C)2006 mgr Jerzy Walaszek           I LO w Tarnowie"
Print
Print "Podaj zakres poszukiwan pierwiastka:"
Print
Input "a = ", a
Input "b = ", b
Print
Print "----------------------------------------------------"
Print
Print "WYNIK:"
Print
fa = f(a) : fb = f(b) : x1 = a : x0 = b
If fa * fb > 0 Then
   Print "Funkcja nie spelnia zalozen"
Else
   While Abs(x1 - x0) > EPSX
      x1 = x0
      x0 = a - fa * (b - a) / (fb - fa)
      f0 = f(x0)
      If Abs(f0) < EPS0 Then Exit While
      If fa * f0 < 0 Then
         b = x0 : fb = f0
      Else
         a = x0 : fa = f0
      End If
   Wend
   
   Print Using "x0 = ######.########"; x0
End If

Print
Print "------------------------------------------------"
Print
Print "Koniec. Nacisnij klawisz 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
  Return x * x * x * (x + Sin(x * x - 1) - 1) - 1
End Function

JavaScript

<html>
  <head>
  </head>
  <body>
    <div align="center">
      <form style="BORDER-RIGHT: #ff9933 1px outset; PADDING-RIGHT: 4px;
                   BORDER-TOP: #ff9933 1px outset; PADDING-LEFT: 4px;
                   PADDING-BOTTOM: 1px; BORDER-LEFT: #ff9933 1px outset;
                   PADDING-TOP: 1px; BORDER-BOTTOM: #ff9933 1px outset;
                   BACKGROUND-COLOR: #ffcc66" name="frmbincode">
        <h3 style="TEXT-ALIGN: center">
          Obliczanie pierwiastka funkcji metodą regula falsi
        </h3>
        <p style="TEXT-ALIGN: center">
          <i>f(x) = x<sup>3</sup>(x + sin(x<sup>2</sup> - 1) - 1) - 1</i>
        </p>
        <p style="TEXT-ALIGN: center">
          (C)2006 mgr Jerzy Wałaszek I LO w Tarnowie
        </p>
        <hr>
        <p style="text-align: center">
          Wpisz do pól edycyjnych zakres przedziału poszukiwań pierwiastka
        </p>
        <div align="center">
          <table border="0" id="table144" cellpadding="8"
                 style="border-collapse: collapse">
            <tr>
              <td>
                a = <input type="text" name="wsp_a" size="20" value="1"
                           style="text-align: right">
              </td>
              <td>
                b = <input type="text" name="wsp_b" size="20" value="2"
                           style="text-align: right">
              </td>
              <td>
                <input type="button" value="Szukaj pierwiastka" name="B1"
                       onclick="main()">
              </td>
            </tr>
          </table>
        </div>
        <div id="out" align=center>...</div>
      </form> 

<script language=javascript>

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

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

// 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)
{
  return x * x * x * (x + Math.sin(x * x - 1) - 1) - 1;
}

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

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

  a = parseFloat(document.frmbincode.wsp_a.value);
  b = parseFloat(document.frmbincode.wsp_b.value);
  if(isNaN(a) || isNaN(b))
    t = "<font color=red><b>Złe krańce zakresu poszukiwań pierwiastka!</b></font>";
  else
  {
    t  = "x<sub>o</sub> = ";
    fa = f(a); fb = f(b); x1 = a; x0 = b;
    if(fa * fb > 0)
      t = "<font color=red><b>Funkcja nie spelnia zalozen</b></font>";
    else
    {
      while(Math.abs(x1 - x0) > EPSX)
      {
        x1 = x0;
        x0 = a - fa * (b - a) / (fb - fa);
        f0 = f(x0);
        if(Math.abs(f0) < EPS0) break;
        if(fa * f0 < 0)
        {
          b = x0; fb = f0;
        }
        else
        {
          a = x0; fa = f0;
        }
      }
      t += x0;
    }
  }
  document.getElementById("out").innerHTML = t;
}

</script>
    </div>
  </body>
</html>

Wynik:

Obliczanie pierwiastka funkcji - metoda regula falsi
f(x) = x^3*(x+sin(x^2-1)-1)-1
----------------------------------------------------
(C)2006 mgr Jerzy Wałaszek I LO w Tarnowie

Podaj zakres poszukiwań pierwiastka:

a = 1
b = 2

----------------------------------------------------

WYNIK:

x0 = 1,18983299

----------------------------------------------------
Koniec. Naciśnij klawisz Enter...

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

Na początek: podrozdziału strony

Materiały tylko do użytku dydaktycznego. Ich kopiowanie i powielanie jest dozwolone
pod warunkiem podania źródła oraz niepobierania za to pieniędzy.

Pytania proszę przesyłać na adres email: i-lo@eduinf.waw.pl

Serwis wykorzystuje pliki cookies. Jeśli nie chcesz ich otrzymywać, zablokuj je w swojej przeglądarce.

Informacje dodatkowe.

Metoda Regula Falsi

Algorytm

Opis Algorytmu

Specyfikacja problemu

Dane wejściowe

Dane wyjściowe

Zmienne pomocnicze i funkcje

Lista kroków

Schemat blokowy

Programy

Obliczanie pierwiastka funkcji metodą regula falsi