Pierwiastki funkcji

Metoda siecznych (ang. secant method) jest ulepszeniem metody regula falsi. Warunki początkowe są zatem podobne:

Na osi OX wybieramy pewien przedział <a;b> o następujących własnościach:

Funkcja jest 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 jest 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 ma 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>:

Dodatkowo w przedziale <a;b> pierwsza pochodna funkcji jest różna od zera. Warunek ten gwarantuje nam, iż w tym przedziale funkcja nie posiada minimum/maksimum lokalnego, a co za tym idzie jej sieczna nie będzie równoległa do osi OX, ponieważ algorytm bazuje na znajdowaniu punktów przecięcia siecznych z osią OX.

W przedziale <a;b> spełniającym warunki początkowe metody wybieramy dwa punkty początkowe x₁ i x₂, które leżą w tym przedziale. Będą to dwa pierwsze przybliżenia pierwiastka:

Z punktów wykresu o współrzędnych (x₁,f(x₁)) i (x₂,f(x₂)) prowadzimy sieczną i wyznaczamy punkt x₃, który jest punktem przecięcia tej siecznej z osią OX:

Do wyznaczenia punktu x₃ stosujemy ten sam wzór, co w metodzie regula falsi:

Sprawdzamy, czy w wyznaczonym punkcie wartość funkcji jest dostatecznie bliska zeru. Jeśli tak, to kończymy z wynikiem x₃. Jeśli nie, to za nowe przybliżenia przyjmujemy dwa ostatnie x₂ i x₃ i kontynuujemy obliczenia wg tej samej metody:

Wzory możemy uogólnić:

Metoda siecznych jest zwykle szybko zbieżna do pierwiastka funkcji. Jednak po wybraniu złych punktów początkowych może się zdarzyć, iż nie będzie ona zbieżna. Dlatego należy zastosować licznik kolejnych przybliżeń, a po przekroczeniu zadanej ich liczby algorytm powinien zatrzymać się z błędem.

Istnieją dwa warunki poprawnego zakończenia obliczeń:

1. Gdy wartość funkcji w wyznaczonym punkcie x_n jest dostatecznie bliska zeru.
2. Gdy odległość dwóch ostatnio wyliczonych przybliżeń pierwiastka jest mniejsza od założonej dokładności wyznaczania pierwiastka.

Algorytm metody siecznych

Dane wejściowe:

Dane wyjściowe

Zmienne pomocnicze

Lista kroków

K01:	f1 ← f(x1) f2 ← f(x2)	Obliczamy wartości funkcji w podanych punktach startowych
K02:	Jeśli n = 0, to zakończ z błędem inaczej n ← n - 1	Sprawdzamy licznik obiegów
K03:		Wyznaczamy punkt przecięcia siecznej z osią OX
K04:	fn  ← f(xn)	Obliczamy wartość funkcji w wyznaczonym punkcie przecięcia
K05:	Jeśli (| fn | < εy) lub (| xn - x2 | < εx), to zakończ	Sprawdzamy warunek zakończenia
K06:	x1 ← x2, f1 ← f2 x2 ←xn, f2 ← fn	Uaktualniamy dwa ostatnie przybliżenia
K07:	Idź do kroku K02	Kontynuujemy obliczenia w następnym obiegu pętli.

Poniższy program jest przykładową realizacją algorytmu siecznych. 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,0>. Jako przybliżenia początkowe wybieramy krańce tego przedziału.

// Pierwiastek funkcji: metoda siecznych
// (C)2019 mgr Jerzy Wałaszek
// Metody numeryczne 0034
//--------------------------------------

#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ść wyznaczania pierwiastka
double epsx = 1e-14;
// Dokładność wyznaczania zera
double epsy = 1e-14;
// Punkty startowe
double x1   = -1.1;
double x2   = -1.0;
// Maksymalna liczba obiegów
int n       = 64;

// Program główny
//---------------
int main()
{
  // Zmienne
  double f1,f2,fn,xn;
  // Zmienna informująca o poprawnym
  // zakończeniu metody
  bool result = false;

  SetConsoleOutputCP(CP_UTF8); 
  SetConsoleCP(CP_UTF8);

  cout << setprecision(15)
       << fixed;
  cout <<
  "Obliczanie przybliżonego pierwiastka "
  "funkcji metodą siecznych\n"
  "-------------------------------------"
  "------------------------\n\n";

  // Obliczamy i zapamiętujemy
  // wartości funkcji na krańcach
  // przedziału <a;b>
  f1 = f(x1);
  f2 = f(x2);

  while(--n)
  {
    // Obliczamy punkt przecięcia
    // siecznej z osią OX
    xn = (f1 * x2 - f2 * x1) / (f1 - f2);

    // Obliczamy wartość funkcji
    // w punkcie przecięcia
    fn = f(xn);

    // sprawdzamy warunki zakończenia
    if((fabs(fn)< epsy) ||
       (fabs(xn - x2) < epsx))
    {
      result = true;
      break;
    }

    // Uaktualniamy przybliżenia
    x1 = x2; f1 = f2;
    x2 = xn; f2 = fn;
  }

  if(!result)
    cout <<
    "Przekroczono liczbę obiegów !!!\n\n";

  cout << "Pierwiastek        xn = "
       << setw(20) << xn << endl
       << "Wartość funkcji f(xn) = "
       << setw(20) << fn << endl
       << "Dokładność dla x epsx = "
       << setw(20) << epsx << 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ą siecznych
-------------------------------------------------------------

Pierwiastek        xn =   -1.077713513691340
Wartość funkcji f(xn) =    0.000000000000000
Dokładność dla x epsx =    0.000000000000010
Dokładność dla y epsy =    0.000000000000010
Liczba obiegów        =    5

Porównaj otrzymany wynik z dwiema poprzednimi metodami:

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

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

# Pierwiastek funkcji: metoda siecznych
# (c)2026 mgr Jerzy Wałaszek
# Metody numeryczne 0034
#--------------------------------------

import math

# 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ść wyznaczania pierwiastka
epsx = 1e-14
# Dokładność wyznaczania zera
epsy = 1e-14
# Punkty startowe
x1 = -1.1
x2 = -1.0
# Maksymalna liczba obiegów
n = 64

# Program główny
#---------------
   
# Zmienna informująca o poprawnym
# zakończeniu metody
result = False
    
print("Obliczanie przybliżonego "
      "pierwiastka funkcji "
      "metodą siecznych\n"
      "-------------------------"
      "--------------------"
      "----------------\n")
    
# Obliczamy i zapamiętujemy
# wartości funkcji na krańcach
f1 = f(x1)
f2 = f(x2)
    
while n > 1:
    n -= 1
        
    # Obliczamy punkt przecięcia
    # siecznej z osią OX
    xn = ((f1 * x2 - f2 * x1)
         / (f1 - f2))
        
    # Obliczamy wartość funkcji
    # w punkcie przecięcia
    fn = f(xn)
        
    # Sprawdzamy warunki zakończenia
    if ((abs(fn) < epsy) or
       (abs(xn - x2) < epsx)):
        result = True
        break
            
    # Uaktualniamy przybliżenia
    x1 = x2
    f1 = f2
    x2 = xn
    f2 = fn
        
if not result:
    print(
    "Przekroczono liczbę obiegów !!!\n")
        
print(
f"Pierwiastek        xn = {xn:20.15f}")
print(
f"Wartość funkcji f(xn) = {fn: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        = {64 - n:4}")
print("\n")
input("Naciśnij Enter...")