Numeryczne rozwiązywanie równań nieliniowych
W wielu zagadnieniach o znaczeniu technicznym i naukowym napotykamy
problem poszukiwania miejsc zerowych pewnej funkcji f(x). Zagadnienie to jest
równoważne rozwiązaniu równania
f(x) = 0.
Dotyczy to także równań algebraicznych stopnia wyższego niż 4, dla których nie
można sformułować analitycznych wzorów określających pierwiastki. W takich
przypadkach stosowane są metody przybliżonego wyznaczania miejsc zerowych dla
narzuconego przedziału [a,b]. Matematycznie, zagadnienie to sprowadza się do
wyznaczenie takiej wartości zmiennej x*, dla której zachodzi f(x*) = 0. Większość
metod stosować można jedynie gdy w przedziale [a,b] występuje pojedynczy
pierwiastek. Do określenia wielkości przedziału zaleca się wstępną separację
pierwiastków, którą można dokonać przez tabelaryzację wartości funkcji w
wyjściowym przedziale [a,b] z pewnym krokiem Dx. W tabeli zawierającej wartości
pary xk,f(xk) poszukujemy takich podprzedziałów [xm,xm+1] dla których następuje
zmiana znaku wartości funkcji. Jeśli Dx zostałnie dobrane odpowiednio małe,
można oczekiwać, że w każdym z podprzedziałów występuję 1 pierwiastek. Inną
metodą określenia przedziału może być wykorzystanie wykresu funkcji f(x).
Dla wielomianów istnieją analityczne metody separacji pierwiastków.
Większość metod poszukiwania miejsc zerowych to metody iteracyjne polegające
na poprawianiu kolejnych przybliżeń pierwiastka xn+1= xn+dx przy dx �� 0. Miarą
szybkości zbieżności metody jest rząd metody
Mówimy, że metoda jest rzędu p, jeśli istnieje pewna stała K taka, że dla dwóch
kolejnych przybliżeń pierwiastka xk oraz xk+1 zachodzi
| xk+1 -� x | Ł� K | xk -� x | p
Metody o rzędzie zbieżności 1 nazywamy liniowymi, przy rzędzie zbieżności 2
nazywamy kwadratowymi, itd.
Miarą efektywności metody w pobliżu dokładnego pierwiastka jest wskaz
nik
efektywności,
1
E =� pk
gdzie k oznacza koszt pojedynczej iteracji danej metody.
Zwykle koszt pojedynczej iteracji zależy głownie od kosztu wszystkich
potrzebnych (w jednej iteracji-kroku) obliczeń wartości funkcji f(x) i jej
pochodnych f (x), natomiast mało zależy od operacji arytmetycznych
wykowywanych na tych wartościach. Dlatego dla porównywania efektywności
różnych metod, można założyć koszt obliczeń f(x) równy 1.
Metoda bisekcji (połowienia przedziału)
Założenia: funkcja f(x) jest ciągła w przedziale [a,b], oraz zachodzi
sgn(f(a)) ą� sgn(f(b)) ,
czyli funkcja różni się znakiem na obu końcach przedziału.
Z twierdzenia o wartości średniej wynika, że w przedziale [a,b] istnieje taki punkt
x*, dla którego f(x*)=0.
Uwaga! Z tego twierdzenia nie wynika wcale, że istnieje tylko jeden taki punkt, może być ich wiele, ale
tylko nieparzysta ilość (tzn. 1, 3, 5 itp.). Załóżmy jednak że istnieje jeden taki punkt.
Generalnie, metody przybliżone rozwiązania takich zagadnień sprowadzają się do
generowania ciągu kolejnych przybliżeń pierwiastka równania, który jest zbieżny
do dokładnego rozwiązania zagadnienia
{xk} x*
��
k��Ą�
x0 =� A , xk =� F(xk -�1) k =�1,2,3,...
Dwa pierwsze kroki metody bisekcji (każdy krok wedle takiego samego schematu)
f (A)�� f (B) <� 0
f(x)=0 x*=?
f(B)
f(x)
f(C)
D=(A+C)/2
A x*
B
x3 x1
x2
x
C=(A+B)/2
f(A)
Procedura generowania kolejnych elementów ciągu ci musi być przerwana po
skończonej liczbie kroków w oparciu o jedno z 3 możliwych kryteriów:
Kryterium 1: | ci  ci-1 | < e�
Kryterium 2: | ci  ci-1 | / | ci | < e� gdy | ci | ą� 0
Kryterium 3: | f(ci) | < e�
gdzie e� - przyjęta dokładność wyznaczenia pierwiastka (e� > 0)
Wybór kryterium zależy od własności problemu, ale preferowane jest kryterium 2.
W każdym kroku algorytmu bisekcji, przedział (czyli błąd bezwzględny
oszacowania pierwiastka) zmniejszany jest o połowę. Zatem ciąg błędów
bezwzględnych jest ciągiem geometrycznym o ilorazie � co możemy zapisać
następująco
| ci  x* | Ł� (b  a)/2i gdzie i  numer iteracji
Wartości a oraz b są znane. To oznacza, że przy stosowaniu kryterium 1
i narzuconej dokładności obliczeń e�, możemy z góry oszacować liczbę kroków
iteracji
(b  a)/2N Ł� e� czyli (b  a)/e� Ł� 2N czyli log2 [ (b  a)/e� ] Ł� log2 2N = N
N ł� log2 [ (b  a)/e� ]
Zaletą algorytmu są bardzo słabe wymagania dotyczące funkcji f(x), i duża
pewność zbieżności generowanego ciągu bowiem w każdej iteracji szukany
pierwiastek znajduje się w przedziale [ai,bi]. Przyjęta dokładność wyznaczania
pierwiastka e� powinna być liczbą odpowiednio małą ale większą od
maksymalnej granicznej dokładności obliczeń maszyny.
Rząd zbieżności metody wynosi
p = 1
(czyli zbiezność liniowa) zaś wskaz
nik efektywności (przy założeniu że koszt
obliczenia wartości funkcji f(x) wynosi k = 1 )
E = 11/1 = 1
Regula falsi (fałszywa linia)
Założenia: - funkcja f(x) jest ciągła w przedziale [a,b],
- sgn(f(a)) ą� sgn(f(b)) ,
- w przedziale [a,b] występuje jeden pierwiastek pojedynczy,
- funkcja f(x) jest klasy C2 w przedziale [a,b],
- pierwsza i druga pochodna mają stały znak w przedziale [a,b].
Dwa ostatnie założenia nie są warunkami koniecznymi zbieżności, zapewniają
tylko ustalenie stałego punktu iteracji i są pomocne przy oszacowaniu błędu.
Dodatkowo przy takich założeniach funkcja f(x) nie ma ekstremów w przedziale
[a,b] i posiada stałą wypukłość.
f(a) f(a) f(b) f(b)
b b a
a
a a
b b
f(b) f(b) f(a) f(a)
f(b)
x1 x2
a=x0
x
b
f(a)
Ciąg przybliżeń pierwiastka konstruujemy zgodnie ze wzorem
(linie sieczne z jednego punktu stałego)
f (xk )
xk +�1 =� xk -� (�b -� xk )�
x0 =� a
f (b) -� f (xk )
Rząd metody wynosi 1 (czyli liniowa), a efektywność metody E=1 (przy założeniu
że koszt obliczenia wartości funkcji f(x) wynosi k=1 ). Oznacza to że regula falsi i
metoda połowienia przedziału są równoważne pod względem szybkości
zbieżności nakładu obliczeń.
Metoda ta jest zbieżna dla dowolnej funkcji ciągłej w [a,b], (f(a)*f(b) < 0 ) jeśli
pierwsza pochodna f (x) jest ograniczona i różna od zera w otoczeniu pierwiastka.
Błąd metody, przy założeniu że punkt xk+1 znajduje się w pobliżu szukanego
pierwiastka x* można oszacować wg wzoru
xk+�1 -� xk
x* -� xk+�1 � �� f (xk+�1)
f (xk+�1) -� f (xk )
Metoda siecznych
Metodę regula falsi, można znacznie przyspieszyć gdy zrezygnuje się z żądania
by punkty odcinka w każdym kroku charakteryzowały się różnym znakiem
wartości funkcji. Wówczas przy wyznaczania kolejnego przybliżenia pierwiastka
opieramy się na kolejnych przybliżeniach xk-1 oraz xk .
f (xk )
xk +�1 =� xk -� (�xk -� xk -�1)�
f (xk ) -� f (xk -�1)
Zalety:
Rząd zbieżności p @� 1.62
wskaz
nik efektywności E = 1.62
Wady:
- Metoda może nie być zbieżna, gdy początkowe przybliżenie x0 nie leży
odpowiednio blisko pierwiastka x*.
- Często konieczne jest wprowadzenie dodatkowego kryterium zbieżności ze
względu na ograniczoną dokładność obliczeń maszyny , bowiem gdy różnica
(xk - xk-1) jest tego samego rzędu co mianownik f(xk) - f(xk-1) , nowa wartość xk+1
może być obarczona dużym błędem.
Metoda Newtona (metoda stycznych)
Założenia:
- funkcja zmienia znak na krańcach
f(x0) przedziału [a,b] w którym jest jeden
pojedynczy pierwiastek,
- pierwsza i druga pochodna mają w nim
stały znak (brak ekstremów i punktów
przegięcia)
a
Rząd zbieżności wynosi p = 2
x1
x2
x0 x
Wskaz
nik efektywności E= 21/(1+y)
gdzie y jest kosztem obliczenia 1
f(x1)
pochodnej f(x). Jeżeli y < 0.44 dla
kosztu obliczenia wartości funkcji
f (x0) f (x0)
=1, bardziej opłaca się stosować
tana =� =� f '(x0) ��
x1 =� x0 -�
x0 -� x1 f '(x0)
metodę Newton a aniżeli siecznych.
f (xk )
xk +�1 =� xk -�
f '(xk )
Punkt startowy w metodzie stycznych
ma istotne znaczenie. Dla funkcji
pokazanej na rysunku procedura jest
a
zbieżna w przedziale [a,b] natomiast nie
b
zbieżna poza nim.
Metoda Newtona jest szeroko stosowana w obliczeniach komputerowych np. z niej
wywodzą się wzory pozwalające iteracyjnie wyznaczać pierwiastek kwadratowy
x =� a
dowolnej liczby. Przyjmując, że f(x) = x2- a , której rozwiązaniem jest , po
analitycznym wyznaczaniu pochodnej dostaniemy zależność
ć� ��
1 a
��
xk+�1 =�
��xk +� xk ��
��
2
Ł� ł�
Analogicznie pierwiastek stopnia 3 wyznaczany jest iteracyjnie wg wzoru
��
1ć�2x +� a
�� ��
xk+�1 =�
2
3�� k xk ��
Ł� ł�
Oba powyższe procesy iteracyjnie są szybko zbieżne a wynik dostajemy po kilku
krokach.
Iteracyjne wyznaczanie pierwiastka 2 i 3 stopnia liczby 10
z dokładnością e� = 10-12
��
ć� ��
1 a
1ć�2x +� a
Nr.
�� ��
��
xk+�1 =�
xk+�1 =�
��xk +� xk ��
��
2
Iteracji
2
3�� k xk ��
Ł� ł�
Ł� ł�
1 10,00000000000 10,00000000000
2 5,50000000000 6,70000000000
3 3,65909090909 4,54092225440
4 3,19600508187 3,18893705034
5 3,16245562280 2,45374135571
6 3,16227766518 2,18945956482
7 3,16227766017 2,15499199671
8 2,15443483415
9 2,15443469003
Metody dla pierwiastków wielokrotnych
poprawiona metoda Newtona
f (xn )
xn+�1 =� xn -� r r  krotność pierwiastka
f '(xn )
zmodyfikowana metoda Newtona
f (x)
u(x) =�
w metodzie tej szukamy miejsca zerowego funkcji
f '(x)
która ma takie same ale pojedyncze pierwiastki
u(xn )
xn+�1 =� xn -�
u'(xn )
gdzie
f ''(xn )
u'(xn ) =� 1-� u(xn )
f '(xn )
lub po rozpisaniu
f (xn ) f '(xn )
xn+�1 =� xn -�
2
(� f '(xn ))� -� f ''(xn ) f (xn )
 Szybkość zbieżności metod Newtona dla funkcji f(x) = x4 -2x3 +2x2 - 2x +1
poprawiona metoda zmodyfikowana
metoda Newtona
Newtona metoda Newtona
3 3 3
2,375 1,75 0,575758
1,91414 1,18314 0,947395
1,58177 1,01517 0,998622
1,35141 1,00011 0,999999
1,20098 1 1
1,10953
1,05759
1,0296
1,01502
1,00756
1,0038
1,0019
1,00095
1,00048
1,00024
1,00012
1,00006
1,00003
1,00001
1,00001
Metoda iteracji prostej
Najprostsza metodą rozwiązania równania postaci f(x)=0 jest przekształcenie
równania do postaci x = g(x) (co zawsze można uzyskać), oraz założenie, że
iteracyjne poszukiwanie kolejnych przybliżeń pierwiastka prowadzić będziemy wg
zależności
xk+1 = g(xk) , startując z założonego x0 .
W matematyce rozwiązanie równania x = g(x) nazywamy punktem stałym funkcji g.
Warunek wystarczający istnienia punktu stałego
TW. Jeżeli funkcja g(x) jest ciągła w [a,b], oraz wartości funkcji także należą do
przedziału [a,b] dla wszystkich x z przedziału [a,b] to funkcja g(x) ma punkt stały w
przedziale [a,b].
Warunek unikalności punktu stacjonarnego
TW. Jeśli spełnione są założenia poprzedniego twierdzenia dla funkcji g(x), oraz
dodatkowo pochodna g (x) spełnia nierówność dla każdego x z przedziału [a,b]
| g (x) | < k<1
to funkcja g(x) ma jeden punkt stały w [a,b].
Zatem jeśli spełnione są warunki podane w obu twierdzeniach dla funkcji g(x),
możemy być pewni że istnieje pojedynczy pierwiastek równania x=g(x) w [a,b].
Podane warunki nie zapewniają jednak zbieżności procedury iteracyjnego
poszukiwania pierwiastka f(x)=0. Również sposób przekształcenia f(x)=0 do postaci
x=g(x) wpływa na zbieżność procedury.
Przykład. Funkcja f(x) = x3+4x2-10 ma w przedziale [1,2] pojedynczy pierwiastek.
Można ją przekształcić do postaci x = g(x) na wiele sposobów:
a b c d
a)
x =� x -� x3 -� 4x2 -�10
0 1,5 1,5 1,5 1,5
1 -0,875 1 1,2869538 1,34839972
2 6,73242188 1,8171206 1,4025408 1,36737637
3
x =� 10-� 4x2 3 -469,72001 -1,474795 1,3454584 1,36495702
b)
4 102754555 1,0913702 1,3751703 1,36526475
5 1,7364277 1,3600942 1,36522559
1
6 -1,2725456 1,367847 1,36523058
c)
x =� 10 -� x3
7 1,5215426 1,363887 1,36522994
2
8 0,9043545 1,3659167 1,36523002
9 1,8878796 1,3648782 1,36523001
10
10 -1,6206132 1,3654101 1,36523001
d)
x =�
11 -0,7966259 1,3651378
x +� 4
12 1,9540829 1,3652772
13 -1,7406313 1,3652059
16 14 -1,2844679 1,3652424
15 1,5037784 1,3652237
12
16 0,9846323 1,3652333
8 17 1,8293538 1,3652284
18 -1,5016488 1,3652309
4
19 0,9933573 1,3652296
0
20 1,8224519 1,3652302
-1 -0,5 0 0,5 1 1,5 2
21 -1,4865951 1,3652299
-4
22 1,0507598 1,3652301
-8
-12
nie zbieżne zbieżne
Dlaczego tak jest? Sprawdzimy, jakie własności posiada w każdym z przypadków
funkcja g(x)
6
przypadek a
4
przypadek b
przypadek c
przypadek d
2
0
-2
-4
-6
1.0 1.2 1.4 1.6 1.8 2.0
x
g(x)
Powiększając skalę rysunku widzimy oczywistą przyczynę : warunek
wystarczający istnienia punktu stałego (jego spełnienie!)
2
1
przypadek a
0
przypadek b
przypadek c
przypadek d
-1
-2
1.0 1.2 1.4 1.6 1.8 2.0
x
g(x)