Metody numeryczne

Metody numeryczne są to sposoby rozwiązywania złożonych problemów matematycznych za pomocą narzędzi obliczeniowych udostępnianych przez popularne języki programowania.
Metody numeryczne są jedną z tych dziedzin matematyki stosowanej, których zastosowanie w praktyce jest powszechne. Wykorzystywane są wówczas, gdy badany problem nie ma w ogóle rozwiązania analitycznego (danego wzorami), lub korzystanie z takich rozwiązań jest uciążliwe ze względu na ich złożoność lub z innych powodów (np. stosowanie eliminacji Gaussa zamiast wyliczania rozwiązań układu równań metodą wyznaczników jest stosowana dlatego, że jest lepiej uwarunkowana numerycznie, a nie dlatego, że brak jest wzoru). Otrzymywane tą drogą wyniki są na ogół przybliżone, jednak dokładność obliczeń może być z góry określona i dobiera się ją zależnie od potrzeb.
Inaczej mówiąc: metody numeryczne polegają na uzyskiwaniu wyniku poprzez sekwencję kolejnych przybliżeń. W efekcie otrzymany wynik będzie cechował się prawie zawsze pewnym błędem, chociaż ten błąd może być dowolnie mały.
Metody numeryczne znalazły zastosowanie wszędzie tam, gdzie dojście do wyniku innymi sposobami jest niemożliwe lub bardzo trudne.
Takim przykładem z życia szkolnego może być np. obliczenie powierzchni jakiejś nieregularnej figury lub pola pod krzywą.
Z metodami numerycznymi związane są takie pojęcia jak: konwergencja i dywergencja oraz dokładność.
Konwergencja określa zdolność danej metody numerycznej do "zmierzania" w kierunku wyniku. Użytkownik jest zainteresowany, aby zastosowana metoda cechowała się maksymalną konwergencją. Niestety czasami dla pewnych specyficznych danych zastosowana metoda numeryczna może zachowywać się co najmniej dziwnie - zamiast prowadzić do wyniku, daje rezultaty coraz bardziej od wyniku odległe. Tę niepożądaną cechę nazywa się dywergencją. Czasem dywergencja dla pewnego małego zakresu danych jest ceną, jaką płacimy za korzyści wynikające z zastosowanej metody.
Dokładność - ważna cecha metod numerycznych, jest zwykle określana przed rozpoczęciem obliczeń. Im dokładniejszy wynik chce się otrzymać, tym więcej czasu trzeba poświęcić na obliczenia.
Dla niektórych metod wynik może być uzyskany z dowolnie dużą dokładnością, lecz niektóre metody dają wynik obarczony pewnym stałym błędem, zatem zwiększanie precyzji obliczeń nic tu nie pomoże. Czemu stosuje się takie metody - ponieważ są bardzo szybkie i proste w implementacji.
W metodach numerycznych ważny jest kompromis pomiędzy czasem obliczeń, ich dokładnością oraz uniwersalnością metody.
Uwaga, w większości zagadnień numerycznych operuje się na liczbach rzeczywistych.
Programista, który zamierza napisać program rozwiązujący zagadnienia numeryczne musi wybrać odpowiadający mu typ rzeczywisty. Najlepiej gdyby to był typ dający największą dokładność obliczeń (np. long double w języku C++). Niestety, konsekwencją takiego wyboru jest wydłużenie czasu obliczeń i zwiększenie wymaganej przez program pamięci operacyjnej. Gdy problem numeryczny jest bardzo złożony (np. w trakcie obliczeń trzeba jeszcze przydzielić dodatkową pamięć), użycie liczb o dużej precyzji może doprowadzić do szybkiego wyczerpania zasobów maszyny obliczeniowej i w konsekwencji uniemożliwić wykonanie obliczeń.

Szukanie miejsc zerowych funkcji metodami numerycznymi

Aby poznać i zrealizować wybrane metody numeryczne wybrałam jako przykłady szukanie miejsca zerowego konkretnej funkcji z określoną dokładnością. Zastosowane implementacje mają poważne zalety - są proste w zrozumieniu, nie wymagają wyższej matematyki, a co najważniejsze dobrze ilustrują problem.
Funkcja f(x) = x³ - 3x - 20, ma jedno miejsce zerowe, którego położenie można określić, jednak jest ono liczbą niewymierną. Szukanie miejsca zerowego tej funkcji z dokładnością 0,0001 posłuży do zaprezentowania, zilustrowania i porównania kilku wybranych metod numerycznych.

Wyznaczanie miejsc zerowych funkcji jednej zmiennej (f(x) = x³ - 3x - 20)

Metoda bisekcji

Wyznaczanie miejsc zerowych funkcji jednej zmiennej (f(x) = x³ - 3x - 20)

Znaną metodę połowienia do poszukiwania elementu w uporządkowanym ciągu można również użyć do znajdowania miejsca zerowego funkcji w zadanym przedziale.
Dana jest funkcja f(x) oraz końce przedziału domkniętego [a,b].
W tym przedziale chcemy znaleźć punkt c, który jest miejscem zerowym tej funkcji, czyli f(c)=0
Tak postawiony problem ma rozwiązanie, jeśli funkcja f(x) jest ciągła w przedziale [a,b] oraz spełniony jest warunek, że wartości funkcji na końcach przedziału mają przeciwne znaki, tj. f(a) * f(b) <= 0.
Z tego wynika, że wykres funkcji f(x) musi w tym przedziale przynajmniej raz przeciąć oś Ox - ten punkt przecięcia jest właśnie szukanym miejsem zerowym funkcji. Jest to założenie oparte na twierdzeniu Bolzano-Cauchy`ego: Jeżeli funkcja jest ciągła w przedziale domkniętym i na jego końcach przyjmuje wartości różnych znaków, to między punktami przedziału znajduje się co najmniej jeden pierwiastek równania.
Metoda ta zawsze jest konwergentna - zbieżna, co gwarantuje odnalezienie pierwiastka równania, czasem jednak kosztem liczby iteracji.

Ogólny opis metody bisekcji

Użytkownik podaje wartości końców przedziału, w którym spodziewa się znaleźć rozwiązanie. Jego program sprawdza, czy wartości funkcji dla tych krańców są przeciwnego znaku. Jeśli tak tzn, że w tym przedziale na pewno funkcja przecina oś x. Program wg metody bisekcji, dzieli ten przedział na pół i sprawdza, czy środek tego przedziału jest szukanym miejscem zerowym (z założoną dokładnością). Jeśli tak to kończy swoje działanie i zwraca go (ten środek) jako wynik, jeśli nie to wybiera (w zależności od znaku funkcji w połowie przedziału) nowy, pomniejszony o połowę przedział i kontynuuje obliczenia.

Specyfikacja

Dane:

• funkcja f(x) ciągła w przedziale domkniętym [a,b] i spełniająca f(a) * f(b) < 0;
• epsilon - przyjęta dokładność obliczeń (epsilon>0)

Wyniki:

• c przybliżenie zera funkcji f w przedziale [a,b], spełniające warunki:
  | f ( c ) | <= epsilon (założone epsilon powinno być znacząco większe od używanej dokładności liczb rzeczywistych)

Algorytm:

1. Pobierz wartości końców przedziału [a,b]. Sprawdź czy f(a) * f(b) < 0
2. Jeśli f(a) * f(b) < 0, idź do polecenia 3, wpp wróć do polecenia 1.
3. Oblicz c = (a + b) / 2
4. Policz f(c).
5. Jeśli |f(c)| <= epsilon lub |a - b| <= epsilon zakończ obliczenia - c jest szukanym miejscem zerowym funkcji, wpp sprawdź znak f(c).
   Jeśli sign f(c) = sign f(a) podstaw a = c, wpp b = c. Wróć do polecenia 3
   Uwaga, Dzieki drugiemu warunkowi mamy pewność, że uniknie się zapętlenia. 

Metoda siecznych

Wyznaczanie miejsc zerowych funkcji jednej zmiennej (y = x³ - 3x - 20)

Metoda siecznych (interpolacji liniowej) polega na przyjęciu, że funkcja na dostatecznie małym odcinku w przybliżeniu zmienia się w sposób liniowy. Możemy wtedy na odcinku krzywą y=f(x) zastąpić sieczną. Za przybliżoną wartość pierwiastka przyjmujemy punkt przecięcia siecznej z osią odciętych OX.

Ogólny opis metody siecznych (i metody reguła falsi)

Metoda siecznych wymaga od użytkownika podania dwóch przybliżeń x_o i x₁, dla których wartości funkcji mogą, ale nie muszą leżeć po przeciwnych stronach osi x (jak w metodzie bisekcji i w metodzie reguła falsi). Jeśli jednak te pierwsze, podane przez użytkownika przybliżenia, będą jednocześnie końcami przedziału, w którym funkcja ma miejsce zerowe, to konwergencja metody będzie wówczas wysoka. Przez te podane dwa punkty (użytkownik podaje wartości x_o i x₁, program oblicza wartości funkcji), prowadzona jest prosta, miejsce przecięcia tej prostej z osią x jest przybliżonym wynikiem szukanego miejsca zerowego, o ile bezwzględna wartość funkcji w tym punkcie jest mniejsza od założonej dokładności.
Uwaga, niektórzy autorzy utożsamiają metodę siecznych z metodą reguła falsi, gdyż obliczenia w obu wypadkach przeprowadza się w ten sam sposób i według tego samego wzoru, jedyna różnica polega właśnie na warunku jaki spełnia wybrany do szukania miejsca zerowego przedział: w metodzie siecznych te dwa przybliżenia podane przez użytkownika są przedziałem w którym może znajdować się pierwiastek równania (lub w pobliżu tego przedziału), zaś w metodzie reguła falsi musi leżeć w nim ten pierwiastek, czyli należy sprawdzić, że f(x₀) * f(x₁) < 0. Porównaj [2] str. 173 i 176

Specyfikacja

Dane:

• funkcja f(x) ciągła w przedziale domkniętym [x₀,x₁] i spełniająca f(x₀) * f(x₁) < 0 (ten drugi warunek oznacza, że tak naprawdę przyjmujemy założenia do metody reguła falsi, dla metody siecznej ważne jest, aby wartość f(x₀) <> f(x₁), by to były różne punkty startowe);
• epsilon - przyjęta dokładność obliczeń (epsilon>0)

Wyniki:

• x₂ przybliżenie miejsca zerowego funkcji f, spełniające warunki: |f(x₂)| <= epsilon; wartość funkcji f powinna być liczona ze znacząco większą dokładnością liczb rzeczywistych niż epsilon.

Uwaga 1 na podstawie [3] str. 125 - "...niestety zdarzają się przypadki, gdy może ona (m. siecznych) nie być zbieżna, np. gdy początkowe przybliżenia nie leżą dostatecznie blisko pierwiastka. W metodzie tej istotne znaczenie ma maksymalna graniczna dokładność (wynikająca z przyjętej arytmetyki). Gdy bowiem różnica x_n+1 - x_n jest tego samego rzędu co oszacowanie błędu, jakim jest obarczona, następne przybliżenie może już być całkowicie błędne. Dlatego też za dodatkowe kryterium przerwania iteracji należy przyjmować wartości |f(x_n)|, tak aby tworzyły one ciąg malejący (w końcowej fazie obliczeń). ... jeśli różnica między kolejnymi przybliżeniami zamiast maleć zaczyna szybko wzrastać. W takim przypadku należy przeprowadzić powtórną lokalizację pierwiastka.."
Uwaga 2 na podstawie [3] str. 126 - "Przy stosowaniu metody siecznych jest niezbędne, aby pierwsza iteracja zaczynała się z punktów, w których funkcja ma różne znaki, w przeciwnym razie możemy wykryć nieistniejący pierwiastek, co jest szczególnie niebezpieczne przy obliczeniach prowadzonych na maszynach cyfrowych"
Algorytm:

1. Pobierz wartości dwóch pierwszych przybliżeń [x₀,x₁]. Sprawdź czy f(x₀) * f(₁) < 0, jeśli nie pobierz inne przybliżenia. Sprawdzanie, czy podawane przez użytkownika dwa pierwsze przybliżenia są jednocześnie końcami przedziału, w którym znajduje się miejsce zerowe funkcji, następuje na początku algorytmu.
2. Oblicz^* x₂ = x₁ - f (x₁) * (x₁ - x₀) / (f(x₁) - f(x₀)) .
3. Policz f(x₂).
4. Jeśli |f(x₂)| <= epsilon zakończ obliczenia - x₂ jest szukanym miejscem zerowym funkcji, wpp podstaw jako nowe x₀ = x₁ i nowe x₁ = x₂. Wróć do polecenia 2 - to w metodzie siecznych; w metodzie reguła falsi należy podobnie jak w metodzie bisekcji wybrać przedział, w którym znajduje się szukany pierwiastek funkcji. 
5. ^* Wyprowadzenie wzoru na x₂

• x₀ i x₁ - dwa kolejne przybliżenia miejsca zerowego funkcji, zarazem końce przedziału, w którym to miejsce się znajduje,
• przez te dwa punkty prowadzimy sieczną o wzorze w(x) = a * x + b,
• wartości funkcji i siecznej są w tych punktach te same, czyli
  w(x₀) = f(x₀) oraz
  w(x₁) = f(x₁), zatem
  f(x₀) = a * x₀ + b oraz
  f(x₁) = a * x₁ + b
• mając te dwa równania można wyliczyć a i b (odjąć stronami i obliczyć a, następnie z jednego z równań policzyć b) - reszta jest dana,
  
  
• x₂ jest to miejsce przecięcia siecznej z osią x, czyli
  0 = a * x₂ + b (a i b wcześniej policzone), stąd można wyliczyć x₂,
• x₂ = x₁ - f (x₁) * (x₁ - x₀) / (f(x₁) - f(x₀))  

• Metoda stycznych (Newtona)

  Wyznaczanie miejsc zerowych funkcji jednej zmiennej (y = x³ - 3x - 20)

  Metoda stycznych (Newtona) opiera się na założeniu, że ponieważ pochodna w punkcie jest granicą ilorazu różnicowego^* przy epsilon dążącym do 0, to wartość ilorazu różnicowego liczona dla dostatecznie małych wartości epsilon jest bliska dokładnej wartości pochodnej.
  ^*Przybliżone wyznaczanie pochodnej na podstawie definicji ilorazu różnicowego funkcji:
  
  
  
  

  Ogólny opis metody stycznych (Newtona)

  Metoda stycznych wymaga od użytkownika podania pierwszego przybliżenia x₀ (leżącego dostatecznie blisko pierwiastka funkcji), dla którego obliczana jest wartość funkcji. Jeśli bezwzględna wartość funkcji (co do modułu) w tym punkcie jest mniejsza od założonej dokładności, to jest ona szukanym przybliżonym miejscem zerowym, jeśli nie to w tym punkcie prowadzona jest styczna do funkcji (obliczana jest pochodna funkcji - musi być różna od zera) i miejsce przecięcia tej stycznej z osią x przyjmowane jest za kolejne przybliżenie pierwiastka funkcji. Powtarza się wszystkie obliczenia, aż do uzyskania założonej dokładności.
  

  Specyfikacja

  Dane:
  
• funkcja f(x) ciągła w otoczeniu punktu x₀,
• x₀ - początkowe przybliżenie pierwiastka, spełniające warunek, iż pochodna funkcji dla tego punktu jest różna od zera,
• epsilon - przyjęta dokładność obliczeń (epsilon>0)

Wyniki:

• x₁ przybliżenie miejsca zerowego funkcji f, spełniające warunek: |f(x₁)| <= epsilon

Algorytm:

1. Pobierz wartości pierwszego przybliżenia x₀. Oblicz pochodną f'(x₀). Sprawdź czy f'(x₀) <> 0, jeśli nie podaj inne przybliżenia. Sprawdzanie, czy pochodna jest różna od zera jest przygotowaniem do następnego kroku przy założeniu, że użytkownik od razu nie poda rozwiązania z założoną dokładnością. Idź do polecenia 3.
2. Oblicz pochodną f'(x₀). Sprawdź czy f'(x₀) <> 0, jeśli nie zakończ obliczenia.
3. Oblicz^* x₁ = x₀ - f(x₀) / f'(x₀).
4. Policz f(x₁).
5. Jeśli |f(x₁)| <= epsilon zakończ obliczenia - x₁ jest szukanym miejscem zerowym funkcji, wpp podstaw jako nowe x₀ = x₁. Wróć do polecenia 2.

Uwaga - kiedy ta metoda gwarantuje zakończenie obliczeń?
Gdy w przedziale, w którym położony jest pierwiastek, funkcja f(x) ma w końcach różne znaki oraz druga pochodna ma stały znak. Uzasadniając, można skorzystać z twierdzenia podanego w [3] na str. 130.
Twierdzenie. Jeżeli mamy przedział < a; b > taki, że:
(i) f(a) i f(b) mają przeciwne znaki,
(ii) f"(x) jest ciągła i nie zmienia znaku na < a; b >,
(iii) styczne do krzywej y=f(x) poprowadzone w punktach o odciętych a i b przecinają oś Ox wewnątrz przedziału < a; b >,
wówczas równanie f(x)=0 ma dokładnie jeden pierwiastek alfa w przedziale < a; b > i metoda stycznych (Newtona) jest zbieżna do alfa dla dowolnego punktu startowego x₀ należącego do < a; b >.
^* Wyprowadzenie wzoru na x₁

• styczna do funkcji f(x) w punkcie x₀ przechodzi przez punkt [x₀, f(x₀) ] i przecina oś x w punkcie o współrzędnych [x₁, 0],
• styczna do funkcji ma w punkcie styczności taki sam kąt nachylenia wykresu jak funkcja,
• pochodna w punkcie ma wartość równą tangensowi kąta nachylenia wykresu funkcji czyli f'(x₀) = f(x₀) / (x₀ - x₁),
• po przekształceniu otrzymujemy x₁ = x₀ - f(x₀) / f'(x₀) 
•  
  

  Metoda iteracji prostej Banacha

  Wyznaczanie miejsc zerowych funkcji jednej zmiennej (y = x³ - 3x - 20)

  Ogólny opis metody iteracji prostej Banacha

  Metoda iteracji prostej Banacha składa się z dwóch etapów.
• Najpierw równanie f(x) = 0 przekształcamy do równoważnej postaci x = g(x)
  (takie przekształcenie jest zawsze wykonalne (np. dodając x do obu stron równania f(x)=0) i zazwyczaj można je wykonać kilkoma sposobami).
• Wybieramy przybliżenie x₀. Kolejne iteracje obliczamy ze wzoru x_n+1 = g(x_n)

Szczegółowy opis metody iteracji prostej dla rozpatrywanej funkcji

1. Najpierw należy znaleźć funkcję g(x), przekształcając równanie 0 = x³ - 3x - 20 do postaci x = g(x)
2. Po przekształceniach uzyskujemy np. jedną z poniższych funkcji:
   
   
3. Należy wybrać taką funkcję g(x), aby zapewniała zbieżność kolejnych iteracji, tzn. |x_n+1 - x_n| < |x_n - x_n-1| czyli kolejna różnica powinna być mniejsza niż poprzednia.
   Tego warunku nie spełniają funkcje g₁(x), g₂(x) i g₃(x) - pomijamy tu szczegółowe obliczenia. Tylko funkcja g₄(x) gwarantuje spełnienie tego warunku, natomiast sama funkcja g₄(x) ma inne ograniczenia, gdyż nie dla każdego x zachodzi x=g₄(x) <-> f(x)=0. Należy z obliczeń wyłączyć -20/3 < x oraz x < =0. Po wstępnym zbadaniu funkcji f(x) wiemy, że najlepiej poszukiwać jej miejsca zerowego w zakresie x należącym do przedziału [1, 5], zatem można przyjąć funkcję g₄(x) do obliczeń mimo pewnych zastrzeżeń.
4. Należy po wybraniu przybliżenia x₀ policzyć x₁ = g(x₀) ... itd. dla kolejnych x_n, za każdym razem sprawdzając czy |x_n - x_n-1| < epsilon.
   Jeśli ten warunek osiągniemy to przybliżonym miejscem zerowym funkcji f(x) jest ostatnio policzone x_n.

Z czego to wynika?

Metoda iteracji prostej należy do tzw. metod iteracyjnych stałopunktowych. To podejście do wyznaczania miejsca zerowego jest oparte na metodzie Banacha.
Poniżej wykresu wykazuję, że istnieje stała 0 <= L <= 1 np. L=0,6 taka, że
|g(x) - g(y)| <= L * |x -y| dla każdego x, y należącego do D₀ czyli g(x) jest odwzorowaniem zwężającym.
Zgodnie z twierdzeniem Banacha o kontrakcji równanie x=g(x) ma dokładnie jedno rozwiązanie x^*, oraz

dla dowolnego przybliżenia początkowego x₀ należącego do D₀.
Ponieważ równanie x=g(x) jest równaniem równoważnym (tzn. ma te same rozwiązania) co równanie f(x), zaś x dla którego zachodzi równość x=g(x) jest punktem stałym odwzorowania g, więc znajdując rozwiązanie x^* równania x=g(x) znajdujemy jednocześnie rozwiązanie równania f(x)=0.
Na stronie Numerki autor przedstawił ilustację graficzną prezentującą dwa warunki gdy metoda iteracj prostej jest zbieżna i rozbieżna.

Specyfikacja

Dane:

• funkcja g(x) zbieżna dla zadanego punktu x₀, tzn. taka, że |g(g(x)) - g(x)| < |g(x) - x| dla x postaci g^k(x₀), k=0,1,...
• x₀ - początkowe przybliżenie pierwiastka równania x=g(x).
• epsilon - przyjęta dokładność obliczeń (epsilon>0)

Uwaga. g(x) jest równaniem równoważnym otrzymanym po zastąpieniu równania wyjściowego f(x) = 0 przez x=g(x) Wyniki:

• x_bież przybliżona wartość miejsca zerowego funkcji f(x), spełniającego warunek: |x_bież - x_poprz| <= epsilon

Algorytm:

1. Pobierz wartość pierwszego przybliżenia x₀. Podstaw x_poprz=x₀.
2. Oblicz x_bież = g(x_poprz), np. x_bież = sqrt(3+20/x_poprz)
3. Policz różnicę między x_bież i x_poprz.
4. Jeśli |x_bież - x_poprz| <= epsilon zakończ obliczenia - x_bież jest szukanym miejscem zerowym funkcji, wpp podstaw jako nowe x_poprz = x_bież. Wróć do polecenia 2.

Ilustracja graficzna metody iteracji prostej Banacha

Twoja przeglądarka nie obsługuje apletów java

Wybieram D₀ = [2, 10]
Przekształcam nierówność: |pierw(3+20/x)-pierw(3+20/y)| - L * |x-y| <= 0 dla każdego x, y należącego do [2, 10]
Wykonuję obliczenia w arkuszu kalkulacyjnym, zakładając L = 0,6. Wartość maksymalna z całego obszaru obliczeń = 0.
Jeszcze bardziej obrazowo widać na wykresie, wykonanym na podstawie tych obliczeń, spełnienie tej nierówności świadczącej, że dla x należącego do [2, 10] funkcja g(x) jest odwzorowaniem zwężającym.

Można zadać pytanie, ile musimy wykonać kroków (kolejnych przybliżeń), od założonego punktu startowego x₀, aby mieć pewność, że obliczone x_bież rzeczywiście dało rozwiązanie z dokładnością epsilon?
Powołując się na wniosek sformułowany w pracy [5] mówiący: "Przy założeniach twierdzenia Banacha, metoda iteracji prostych jest zbieżna co najmniej liniowo z ilorazem L, tzn.
||x_k - x_bież||<= L^k ||x₀ - x_bież||"
Zakładając, że x_k-x_bież = epsilon oraz przyjmując, że maksymalne oddalenie punktu startowego od rozwiązania (x₀-x_bież) wynosi 8, można dla założonego L=0,6 obliczyć liczbę przybliżeń do rozwiązania z dokładnością epsilon.
epsilon = L^k |x₀ - x_bież|
k = log_L (epsilon/(|x₀ - x_bież|))
Oto fragment obliczeń:
Im bliżej rozwiązania wystartujemy i im L jest mniejsze, tym w mniejszej liczbie przybliżeń osiągniemy wartość x_bież z założoną dokładnością epsilon.