Miejsca Zerowe Funkcji - Metoda siecznych

Mamy daną funkcję f(x), dwa punkty startowe x₁ oraz x₂ i przedział [a,b] poszukiwań pierwiastka, do którego należą punkty x₁, x₂. W przedziale poszukiwań pierwiastka 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 (nie obowiązuje to punktów x₁ i x₂). 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ą:

Dodatkowo w przedziale [a,b] pierwsza pochodna f '(x) jest różna od zera. Nie istnieje zatem minimum lub maksimum lokalne. Ten warunek gwarantuje nam, iż sieczna nie będzie równoległa do osi OX, co uniemożliwiłoby wyznaczenie jej punktu przecięcia z tą osią.

Gdy funkcja f(x) spełnia podane warunki, to w przedziale [a,b] istnieje pierwiastek i możemy go wyszukać metodą siecznych.

Metoda siecznych (ang. secant method) jest ulepszeniem metody regula falsi. W metodzie regula falsi żądamy, aby funkcja f(x) zawsze przyjmowała różne znaki na krańcach przedziału poszukiwań pierwiastka. Różne znaki gwarantują nam istnienie pierwiastka w tym przedziale. Otóż w metodzie siecznych taki wymóg nie jest konieczny (za wyjątkiem pierwszego obiegu). Do wyznaczenia kolejnego przybliżenia pierwiastka funkcji wykorzystujemy dwa poprzednio wyznaczone punkty:

Okazuje się, iż takie podejście do problemu znacznie przyspiesza znajdowanie pierwiastka funkcji. Prześledźmy graficznie tok postępowania:

	Za początkowe punkty x₁ i x₂ przyjmujemy punkty krańcowe przedziału poszukiwań pierwiastka. Ponieważ funkcja zmienia znak w tym przedziale i jest określona oraz ciągła, mamy gwarancję istnienia pierwiastka. Z punktów wykresu funkcji dla x₁ i x₂ prowadzimy sieczną - na rysunku zaznaczoną kolorem czerwonym. Punkt przecięcia się tej siecznej z osią OX daje nam następny punkt x₃.
	Drugą sieczną prowadzimy z punktów wykresu funkcji dla x₂ oraz x₃. Otrzymujemy punkt x₄.
	Trzecia sieczna prowadzona jest z punktów wykresu funkcji dla x₃ i x₄. Otrzymujemy kolejny punkt x₅. Już widać, iż leży on bardzo blisko rzeczywistego pierwiastka funkcji.
	Ostatnią sieczną prowadzimy z punktów wykresu funkcji dla x₄ i x₅. Otrzymany punkt x₆ jest już dostatecznie blisko pierwiastka funkcji. Zwróć uwagę, iż kolejne punkty x_i zbliżają się do siebie.

Jeśli punkty początkowe x₁, x₂ zostaną źle dobrane, to algorytm może nie dać rozwiązania. Dlatego stosujemy w nim dodatkowy licznik obiegów pętli wyznaczającej kolejne przybliżenia pierwiastka funkcji. Jeśli po wykonaniu zadanej liczby obiegów (u nas 64) nie zostanie znalezione rozwiązanie, to algorytm przerywa obliczenia z odpowiednim komunikatem. Obliczenia również będą przerwane w przypadku, gdy algorytm stwierdzi, iż sieczna jest równoległa do osi OX. Poprawnym warunkiem zakończenia algorytmu jest:

do podrozdziału do strony

Specyfikacja problemu

Dane wejściowe

f(x)	–	funkcja rzeczywista, której pierwiastek liczymy. Musi być ciągła i określona w przedziale poszukiwań pierwiastka.
x₁, x₂	–	punkty krańcowe przedziału poszukiwań pierwiastka funkcji f(x). W dalszej fazie obliczeń punkty te przechowują dwie poprzednie wartości x_o - x₁ ostatnią, a x₂ przedostatnią. x₁, x₂ ∈ R

Dane wyjściowe

x_o

–

pierwiastek funkcji f(x). x_o ∈ R

Zmienne pomocnicze i funkcje

f_o, f₁ , f₂	– wartości funkcji odpowiednio w punktach x_o, x₁, x₂. f_o, f₁ , f₂ ∈ R
i	– licznik obiegów pętli. Obiegi są zliczane wstecz od wartości i = 64. i ∈ C
ε_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 x₁ i x₂
	K02:	f₁ ← f(x₁); f₂ ← f(x₂); i ← 64
	K03:	Dopóki (i > 0) ∧ (\| x₁ - x₂ \| > ε_x): wykonuj kroki K04...K08
	K04:	Jeśli \| f₁ - f₂ \| < ε_o, to pisz "Złe punkty startowe" i zakończ
	K05:
	K06:	Jeśli \| f_o \| < ε_o, to idź do kroku K10
	K07:	x₂ ← x₁; f₂ ← f₁; x₁ ← x_o; f₁ ← f_o
	K08:	i ← i - 1
	K09:	Jeśli i = 0, to pisz "Przekroczony limit obiegów" i zakończ
	K10:	Pisz x_oi zakończ

Schemat blokowy

obrazek

Wykonanie algorytmu siecznych rozpoczynamy od odczytu dwóch punktów x₁ oraz x₂. Punkty te posłużą do wyznaczenia pierwszej siecznej. Nie muszą one obejmować przedziału, w którym funkcja zmienia znak, jednakże jest to wskazane.

Dla obu punktów obliczamy wartości funkcji umieszczając je w zmiennych odpowiednio f₁ i f₂. Ponieważ algorytm siecznych może nie dać rozwiązania, stosujemy proste zabezpieczenie w postaci zliczania wykonanych obiegów pętli. Zmienna i pełni rolę licznika tych obiegów. Obiegi zliczane są wstecz. Na początku w zmiennej i określamy maksymalną ilość obiegów, po których algorytm zaprzestanie wykonywać obliczenia pierwiastka.

Rozpoczynamy pętlę obliczającą kolejne przybliżenia pierwiastka funkcji. Pętla kończy się w jednym z trzech przypadków:

Licznik i osiąga zero - wypisujemy odpowiedni komunikat i kończymy algorytm - pierwiastek nie został znaleziony, a x_o zawiera zwykle jakieś przypadkowe wartości. Jest to sygnał, iż początkowe punkty x₁ i x₂ były zbyt odległe od pierwiastka funkcji.
Ostatnie dwa pierwiastki przybliżone x₁ i x₂ różnią się od siebie o ε_x lub mniej - w takim przypadku x_o jest przybliżonym pierwiastkiem. Wypisujemy x_o i kończymy algorytm.
Moduł różnicy wartości funkcji f₁ i f₂ dla dwóch ostatnich pierwiastków x₁ i x₂ jest równy zero. Ponieważ różnica ta występuje w mianowniku ułamka we wzorze obliczającym nowy pierwiastek x_o, musimy przerwać obliczenia. Wypisujemy odpowiedni komunikat i kończymy algorytm.

Na początku pętli wyznaczamy nowe przybliżenie pierwiastka x_o oraz wartość funkcji w tym punkcie f_o. Do wyznaczenia x_o stosujemy wzór siecznych, który wykorzystuje dwa ostatnio wyliczone przybliżenia pierwiastka x₁ i x₂.

Jeśli f_o jest dostatecznie bliskie zeru, to x_o jest przybliżonym pierwiastkiem, wypisujemy x_o i kończymy algorytm.

W przeciwnym razie przesuwamy x₁ i x_o odpowiednio do x₂ i x₁ wraz z wartościami funkcji w tych punktach i kontynuujemy pętlę. Przesunięcie wartości funkcji f₁ i f_o zwalnia nas od konieczności ich ponownego wyliczania.

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 siecznych
//---------------------------------------------
// (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 x0,x1,x2,f0,f1,f2;
  int    i;

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

  cout << "Obliczanie pierwiastka funkcji - metoda siecznych\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 dwa wstepne punkty x1 i x2:\n\n";
  cout << "x1 = "; cin >> x1;
  cout << "x2 = "; cin >> x2;
  cout << "\n-------------------------------------------------\n\n"
          "WYNIK:\n\n";
  f1 = f(x1); f2 = f(x2); i = 64;
  while(i && (fabs(x1 - x2) > EPSX))
  {
    if(fabs(f1 - f2) < EPS0)
    {
      cout << "Zle punkty startowe\n";
      i = 0;
      break;
    }
    x0 = x1 - f1 * (x1 - x2) / (f1 - f2);
    f0 = f(x0);
    if(fabs(f0) < EPS0) break;
    x2 = x1; f2 = f1;
    x1 = x0; f1 = f0;
    if(!(--i)) cout << "Przekroczony limit obiegow\n";
  }
  if(i) 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 siecznych
//---------------------------------------------
// (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
  x0,x1,x2,f0,f1,f2 : double;
  i                 : integer;

begin
  writeln('Obliczanie pierwiastka funkcji - metoda siecznych');
  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 dwa wstepne punkty x1 i x2:');
  writeln;
  write('x1 = '); readln(x1);
  write('x2 = '); readln(x2);
  writeln;
  writeln('-------------------------------------------------');
  writeln('WYNIK:');
  writeln;
  f1 := f(x1); f2 := f(x2); i := 64;
  while (i > 0) and (abs(x1 - x2) > EPSX) do
  begin
    if abs(f1 - f2) < EPS0 then
    begin
      writeln('Zle punkty startowe');
      i := 0;
      break;
    end;
    x0 := x1 - f1 * (x1 - x2) / (f1 - f2);
    f0 := f(x0);
    if abs(f0) < EPS0 then break;
    x2 := x1; f2 := f1;
    x1 := x0; f1 := f0;
    dec(i);
    if i = 0 then writeln('Przekroczony limit obiegow');
  end;
  if i > 0 then writeln('x0 = ',x0:15:8);
  writeln;
  writeln('-------------------------------------------------');
  writeln('Koniec. Nacisnij klawisz Enter...');
  readln;
end.

Basic

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

Declare Function f(ByVal 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 x0, x1, x2, f0, f1, f2
Dim i As Integer

print "Obliczanie pierwiastka funkcji - metoda siecznych"
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 dwa wstepne punkty x1 i x2:"
Print
Input "x1 = ", x1
Input "x2 = ", x2
Print
print "-------------------------------------------------"
Print
print "WYNIK:"
Print

f1 = f(x1) : f2 = f(x2) : i = 64

While (i > 0) And (Abs(x1 - x2) > EPSX)
   If Abs(f1 - f2) < EPS0 Then
      Print "Złe punkty startowe"
      i = 0
      Exit While
   End If
   x0 = x1 - f1 * (x1 - x2) / (f1 - f2)
   f0 = f(x0)
   If Abs(f0) < EPS0 Then Exit While
   x2 = x1 : f2 = f1
   x1 = x0 : f1 = f0
   i -= 1
   If i = 0 Then Print "Przekroczony limit obiegów"
Wend

If i > 0 Then Print Using "x0 = ######.########"; x0

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(ByVal 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ą siecznych
        </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 dwa punkty dla pierwszej siecznej
        </p>
        <div align="center">
          <table border="0" cellpadding="8" style="border-collapse: collapse">
            <tr>
              <td>
                x<sub>1</sub> = <input type="text" name="wsp_x1" size="20"
                                       value="1" style="text-align: right">
              </td>
              <td>
                x<sub>2</sub> = <input type="text" name="wsp_x2" 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 siecznych
//---------------------------------------------
// (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 x0,x1,x2,f0,f1,f2,i,t;

  x1 = parseFloat(document.frmbincode.wsp_x1.value);
  x2 = parseFloat(document.frmbincode.wsp_x2.value);
  if(isNaN(x1) || isNaN(x2))
    t = "<font color=red><b>Nieprawidłowe dane!</b></font>";
  else 
  {
    f1 = f(x1); f2 = f(x2); i = 64;
    while(i && (Math.abs(x1 - x2) > EPSX))
    {
      if(Math.abs(f1 - f2) < EPS0)
      {
        t = "<font color=red><b>Złe punkty startowe!</b></font>";
        i = 0;
        break;
      }
      x0 = x1 - f1 * (x1 - x2) / (f1 - f2);
      f0 = f(x0);
      if(Math.abs(f0) < EPS0) break;
      x2 = x1; f2 = f1;
      x1 = x0; f1 = f0;
      if(!(--i))
        t = "<font color=red><b>Przekroczony limit obiegów!</b></font>";
    }
    if(i) t = "x0 = " + x0;
  }
  document.getElementById("out").innerHTML = t;
}

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

Wynik:

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

Podaj dwa wstępne punkty x1 i x2:

x1 = 1
x2 = 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].

do podrozdziału do 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 siecznych

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ą siecznych