Mariusz PYRZ

SIMR (PW), Instytut Pojazdów

Metody numeryczne w mechanice

Rozwiązywanie równań nieliniowych

2.

Rozwiązywanie równania nieliniowego

(z jedną niewiadomą)

Problem:

Dana funkcja f : R

R, regularna w przedziale [a,b].

Wyznaczyć pierwiastek (pierwiastki) równania

f(x)=0

.

Iteracyjna metoda poszukiwania rozwiązania:

2

Iteracyjna metoda poszukiwania rozwiązania:

wyznaczanie i poprawianie kolejnych przybliżeń x

1

, x

2

, x

3

…, zbiegających się do rozwiązania x*.

Numeryczne poszukiwanie rozwiązania:

Jak generować kolejne przybliżenia rozwiązania?
Od jakiego punktu rozpocząć poszukiwanie?
Kiedy zatrzymać algorytm?
Jak szybko zmierzamy do rozwiązania (czy można ulepszyć procedurę)?

M.Pyrz Metody numeryczne w mechanice – Równania nieliniowe 10.2011

Rys. 2.1

3

M.Pyrz Metody numeryczne w mechanice – Równania nieliniowe 10.2011

Lokalizowanie pierwiastków

Do rozdzielenia obszarów występowania pierwiastków można wykorzystać
twierdzenie Rolle’a:

Omawiane metody numeryczne przeszukują przedział [a,b] zakładając, że istnieje tam
jeden pierwiastek równania. Przed przystąpieniem do rozwiązywania problemu w
którym występuje wiele pierwiastków można przeprowadzić próbę lokalizacji i
odseparowania pierwiastków, tzn. określenie przedziałów [a,b] do oddzielnego
przeszukiwania.

4

twierdzenie Rolle’a:

Jeżeli f(x) jest określona i ciągła w [a,b], i jeżeli N oznacza
liczbę pierwiastków f(x)=0 w rozpatrywanym przedziale, wówczas:

Jeżeli f(a)f(b)<0 to f(x) posiada co najmniej jeden pierwiastek w [a,b]
(N jest wówczas liczbą nieparzystą). Jeżeli ponadto f’(x)≠0 dla każdego x z
przedziału [a,b] , to istnieje w nim tylko jeden pierwiastek (N=1).

Jeżeli f(a)f(b)>0 to f(x) nie posiada pierwiastków w [a,b] (N=0), lub w
rozpatrywanym przedziale występuje ich parzysta liczba.

M.Pyrz Metody numeryczne w mechanice – Równania nieliniowe 10.2011

Rys. 2.2

5

M.Pyrz Metody numeryczne w mechanice – Równania nieliniowe 10.2011

Polega na budowaniu ciągu kolejnych
przybliżeń rozwiązania, zaczynając od
wartości początkowej x

0

odpowiednio wybranej w [a,b].

1

0

2

1

1

(

, ...)

(

, ...)

(

, ...)

n

n

x

F x

x

F x

x

F x

−

=

=

=

⋯

Metoda kolejnych przybliżeń

6

M.Pyrz Metody numeryczne w mechanice – Równania nieliniowe 10.2011

Sposób rozwiązywania:

-

określ granice przedziału w którym szukamy pierwiastka (f(a)f(b)<0)

-

wybierz punkt początkowy x

0

-

oblicz kolejne przybliżenia rozwiązania x

n+1

=F(x

n, …

) dla n=0,1,2,…

-

zakończ iteracje po spełnieniu warunku zatrzymania procedury

Metody rekurencyjne

– obliczenia wykonywane w danym kroku procesu

iteracyjnego uwzględniają jedynie rezultaty uzyskane na poprzednim kroku

przykład: metoda stycznych (Newtona)

Klasy metod iteracyjnych

x

F x

n

n

+

=

1

(

)

7

Metody nierekurencyjne

– obliczenia na danym etapie procedury iteracyjnej

biorą pod uwagę wyniki uzyskane na wszystkich poprzednich etapach

przykład: metoda połowienia (bisekcji)

x

F x x

x

x

n

n

n

+

−

=

1

0

1

1

(

,

,

,

,

)

…

M.Pyrz Metody numeryczne w mechanice – Równania nieliniowe 10.2011

Kryteria zatrzymania obliczeń

1

1

2

(

)

n

n

n

x

x

f x

ε

ε

−

−

<

<

8

M.Pyrz Metody numeryczne w mechanice – Równania nieliniowe 10.2011

Rys. 2.4

Rys. 2.3

Kryteria zatrzymania obliczeń

1

1

2

(

)

n

n

n

x

x

f x

ε

ε

−

−

<

<

9

M.Pyrz Metody numeryczne w mechanice – Równania nieliniowe 10.2011

Rys. 2.6

Rys. 2.5

Inne kryteria zatrzymania obliczeń

1

3

n

n

n

x

x

x

ε

−

−

<

n<N

max

Względna rόżnica

N – numer iteraji

10

M.Pyrz Metody numeryczne w mechanice – Równania nieliniowe 10.2011

Względna rόżnica

kolejnych przybliżeń

N – numer iteraji
N

max

– maksymalna

dopuszczalna liczba iteracji

Rząd metody

Rząd metody informuje o tym, jak szybko ciąg {x

n

} definiowany przez

kolejne przybliżenia x

n

= F(x

n-1

) zbiega się do rozwiązania x* (informuje o

tym, jak szybko maleje błąd przybliżeń w kolejnych iteracjach).

Metoda iteracyjna rozwiązywania równania f(x)=0 jest rzędu p (p

≥1) jeżeli

1

lim

,

0

n

C

C

ε

+

=

≠

*

k

k

x

x

ε

= −

11

C - stała nieujemna (stała błędu asymptotycznego)
p - wykładnik zbieżności metody

Im wyższy rząd metody p tym szybsza zbieżność do x*.

p=1 metoda liniowa

1<p<2 metoda superliniowa

p=2 metoda kwadratowa

lim

,

0

p

n

n

C

C

ε

→∞

=

≠

M.Pyrz Metody numeryczne w mechanice – Równania nieliniowe 10.2011

k

k

x

x

ε

= −

2

3

4

5

1

2

3

10 ,

10 ,

10 ,

10 ,...

k

k

k

k

ε

ε

ε

ε

−

−

−

−

+

+

+

≈

≈

≈

≈

2

4

8

16

1

2

3

10 ,

10 ,

10 ,

10

,...

k

k

k

k

ε

ε

ε

ε

−

−

−

−

+

+

+

≈

≈

≈

≈

Metoda bisekcji (połowienia)

Punkt początkowy x

0

położony jest na środku przedziału początkowego [a,b].

W każdej iteracji przedział przeszukiwania jest dzielony na połowy.
W kolejnym kroku rozpatrywana jest „połowa” zawierająca rozwiązanie (tj.
mająca różne znaki funkcji na granicach przedziału). W n-tej iteracji przedział
początkowy jest podzielony przez 2

n

a rozwiązanie jest ograniczone przez

∆

a b

=

−

12

Można z góry określić maksymalną liczbę iteracji k, niezbędną do uzyskania
dokładności

∆

n

n

=

2

1

2

k

b a

ε

− ≤

M.Pyrz Metody numeryczne w mechanice – Równania nieliniowe 10.2011

1

1

k

k

x

x

ε

−

−

<

2

1

log

b a

k

ε





−

=









k zaokrąglone do liczby całkowitej (od góry)

Metoda bisekcji

- algorytm

1) Oblicz f(a) i f(b) na granicach przedziału [a, b]

powinno być f(a)f(b) < 0

2) Oblicz przybliżenie rozwiązania x=(a+b)/2 ;

Sprawdź kryterium stopu

Jeżeli kryterium stopu spełnione zatrzymaj obliczenia

Jeżeli nie, przejdź do punktu 3.

3) Oblicz f(x)

13

Jeżeli f(a)f(x) > 0 , to rozwiązanie znajduje się w przedziale [x, b];

zawężamy granice przyjmując dla nowego przedziału
a=x, f(a)=f(x) i kontynuujemy algorytm przechodząc do kroku 2)

Jeżeli f(a)f(x) < 0 , to rozwiązanie znajduje się w przedziale [a, x];

zawężamy granice przyjmując dla nowego przedziału
b=x, f(b)=f(x) i kontynuujemy algorytm przechodząc do kroku 2)

Jeżeli f(x) = 0, to rozwiązanie zostało znalezione - zatrzymaj obliczenia.

M.Pyrz Metody numeryczne w mechanice – Równania nieliniowe 10.2011

Metoda bisekcji

– interpretacja geometryczna

14

M.Pyrz Metody numeryczne w mechanice – Równania nieliniowe 10.2011

Rys. 2.7

Metoda bisekcji

– przypadek “patologiczny”

15

M.Pyrz Metody numeryczne w mechanice – Równania nieliniowe 10.2011

Rys. 2.8

Metoda bisekcji

– zbieżność

Zatrzymanie obliczeń następuje gdy

- znaleziono rozwiązanie “dokładne”;

- spełniono kryterium zatrzymania (dotyczące np. wielkości przedziału

ograniczającego rozwiązanie lub wartości funkcji f(x))

Zbieżność metody

Metoda bisekcji jest zbieżna liniowo z wykładnikiem p=1 i stałą C=1/2.

16

Metoda bisekcji jest zbieżna liniowo z wykładnikiem p=1 i stałą C=1/2.

Zalety: zbieżność zawsze zapewniona, możliwość oszacowania liczby

iteracji potrzebnych do wygenerowania rozwiązania o określonej
dokładności.

Wady: wolna zbieżność

M.Pyrz Metody numeryczne w mechanice – Równania nieliniowe 10.2011

Metoda cięciw (Regula Falsi)

Nowe przybliżenie pierwiastka jest wyznaczone przez punkt x

p

, w którym prosta

łącząca punkty o współrzędnych (a, f(a)) i (b, f(b)) przecina oś x

Jeżeli f(x ) = 0 to pierwiastek został znaleziony.

( )

( )

,

0

( )

( ) 0

( )

( )

p

p

p

x

a

b

x

af b

bf a

x

f a

f b

f b

f a

−

−

−

=

=

−

−

−

17

Jeżeli f(x

p

) = 0 to pierwiastek został znaleziony.

Jeżeli nie, to zawężamy przedział ograniczający rozwiązanie w zależności
od zmiany znaku funkcji na granicach przedziału, wybierając

przedział [x

p

,b] jeżeli f(x

p

)f(b) < 0 lub przedział [a,x

p

] jeśli f(x

p

)f(b) > 0 .

Rozwiązanie poszukiwane jest iteracyjnie, w każdej kolejnej iteracji przedział zawierający
pierwiastek jest zawężany w zależności od znaku f(x

p

)f(b).

M.Pyrz Metody numeryczne w mechanice – Równania nieliniowe 10.2011

Metoda cięciw

- algorytm

1) Wybierz przedział [a,b] (

w którym funkcja zmienia znak

); oblicz y

a

=f(a) i y

b

=f(b).

2) Oblicz

sprawdź kryterium zatrzymania algorytmu

3) Jeżeli f(x)f(b) > 0 , pierwiastek znajduje się w przedziale [a,x];

( )

( )

( )

( )

a

b

a

a

b

b

a

by

ay

af b

bf a

b a

x

a

y

f b

f a

y

y

y

y

−

−

−

=

=

= −

−

−

−

w celu zmniejszenia

błędów zaokrąglenia

18

3) Jeżeli f(x)f(b) > 0 , pierwiastek znajduje się w przedziale [a,x];

podstaw nową granicę b=x, y

b

=f(x) i wróć do kroku 2)

Jeżeli f(x)f(b) < 0 pierwiastek znajduje się w przedziale [x,b];

podstaw nową granicę a=x, y

a

=f(x) i wróć do kroku 2)

Jeżeli f(x) = 0, znaleziono rozwiązanie, zatrzymaj obliczenia

Obliczenia mogą być zatrzymane po spełnieniu

podobnych kryteriów stopu jak w metodzie bisekcji.

M.Pyrz Metody numeryczne w mechanice – Równania nieliniowe 10.2011

Metoda cięciw

– interpretacja geometryczna

19

M.Pyrz Metody numeryczne w mechanice – Równania nieliniowe 10.2011

Rys. 2.9

Metoda cięciw

– zbieżność

Zbieżność metody

Metoda cięciw jest rzędu

Zalety: pewność, szybsza zbieżność niż metoda bisekcji.

p

= +

≈

1

5

2

1 62

.

20

Zalety: pewność, szybsza zbieżność niż metoda bisekcji.

Wady: zbieżność może być wolna gdy f jest silnie wypukła lub wklęsła.

M.Pyrz Metody numeryczne w mechanice – Równania nieliniowe 10.2011

Metoda siecznych

Niech [a,b] będzie przedziałem zawierającym pierwiastek równania f(x)=0 a x

0

wartością

rzędnej [a,b] dobraną tak aby f(x

0

)≠0 . Wybierzmy inny punkt startowy o rzędnej x

1

. Sieczna

łącząca punkty (x

0

,f(x

0

)) i (x

1

,f(x

1

)) przecina oś x w punkcie x

2

. Iteracyjnie kontynuujemy

procedurę przyjmując punkty (x

1

,f(x

1

)) et (x

2

,f(x

2

)) w celu wyznaczenia kolejnego przybliżenia

rozwiązania o odciętej x

2

, itd. Otrzymujemy w ten sposób ciąg {x

i

} zdefiniowany przez

.

1

1

(

) (

)

:

(

)

n

n

n

n

n

n

f x

x

x

x

x

F x

−

+

−

= −

=

−

21

.

Zbieżność jest zagwarantowana jedynie jeżeli początkowa wartość x

0

wybrana jest w bliskim

sąsiedztwie poszukiwanego pierwiastka który spełnia warunek

.

1

1

:

(

)

(

)

(

)

n

n

n

n

n

x

x

F x

f x

f x

+

−

= −

=

−

( )

1

F x

′

<

M.Pyrz Metody numeryczne w mechanice – Równania nieliniowe 10.2011

Metoda siecznych

- algorytm

Wybierz dwa punkty startowe x

0

i x

1

z przedziału [a,b]

Oblicz kolejne przybliżenie rozwiązania (n=1,2,…)

Przerwij obliczenia jeżeli

1

1

1

1

(

) (

)

;

(

);

1

(

)

(

)

n

n

n

n

n

n

n

n

f x

x

x

x

x

f x

n

n

f x

f x

−

+

+

−

−

= −

= +

−

22

- znaleziono rozwiązanie problemu f(x

n+1

)=0

- lub spełniono warunek zatrzymania algorytmu, np.:

Uwaga : Początkowe wzory metody siecznych są podobne do metody cięciw ale sieczne

buduje się bez sprawdzania warunku f(x

1

)f(x

2

) < 0 i wykorzystuje się zawsze dwa

ostatnio wyznaczone punkty x

n

i x

n-1

.

x

x

n

n

+

−

≤

1

1

ε

M.Pyrz Metody numeryczne w mechanice – Równania nieliniowe 10.2011

Metoda siecznych

– interpretacja geometryczna

23

M.Pyrz Metody numeryczne w mechanice – Równania nieliniowe 10.2011

Rys. 2.10

Metoda stycznych (Newtona-Raphsona)

Załóżmy, ze f(x) jest klasy C

2

w otoczeniu pierwiastka x*.

Rozwinięcie f(x) w szereg Taylora wokół oszacowania x

n

wartości x* :

Dla x

n

wystarczająco bliskiego x*, można pominąć składniki rzędu (p

≥0)

Ponieważ x* jest pierwiastkiem f(x*) =0 , otrzymujemy przybliżenie

f x

f x

x

x

f

x

x

x

f

x

n

n

n

n

n

(

)

(

)

(

)

(

)

(

)

!

(

)

*

*

*

=

+

−

′

+

−

′′

+

2

2

…

ε

n

p

24

Ponieważ x* jest pierwiastkiem f(x*) =0 , otrzymujemy przybliżenie

Przybliżenie błędu związanego z zastąpieniem x* przez x

n

.

Można zatem zbudować

następującą procedurę iteracyjną

f x

x

x

f

x

n

n

n

(

)

(

)

(

)

*

+

−

′

≈

0

ε

n

n

n

n

x

x

f x

f

x

=

−

≈ −

′

*

(

)

(

)

x

x

x

x

f x

f

x

F x

n

n

n

n

n

n

n

n

+

+

=

+

=

−

′

=

1

1

ε

(

)

(

)

(

)

ε

n

M.Pyrz Metody numeryczne w mechanice – Równania nieliniowe 10.2011

Metoda stycznych

- algorytm

1) Wybierz punkt startowy

, oblicz f(x

0

), f '(x

0

), n=0

2) Oblicz przybliżenie rozwiązania x

n+1

3) Jeżeli kryterium zatrzymania nie jest spełnione, podstaw n=n+1

i powróć do kroku 2)

x

a b

0

∈

[ , ]

1

(

)

(

)

n

n

n

n

f x

x

x

f x

+

= −

′

25

i powróć do kroku 2)

Przerwij obliczenia jeżeli

- znaleziono rozwiązanie problemu f(x

n+1

)=0

- lub spełniono warunek zatrzymania algorytmu, np.

x

x

1

0

1

−

≤

ε

M.Pyrz Metody numeryczne w mechanice – Równania nieliniowe 10.2011

max

n

N

≤

W praktyce dodaje się też warunek

Metoda stycznych

– interpretacja geometryczna

26

M.Pyrz Metody numeryczne w mechanice – Równania nieliniowe 10.2011

Rys. 2.11

Metoda stycznych

- zbieżność

Metoda nie zawsze jest zbieżna

Jeżeli
Oraz f’(x) i f’’(x) nie zmieniają znaku w [a,b], to przyjmując x

0

tak żeby f’(x

0

)f’’(x

0

)>0

zapewnimy zbieżność metody z wykładnikiem p=2 (do pierwiastka pojedynczego)

2

( )

([ , ])

( ) ( )

0

f x

C

a b

f a f b

∈

<

27

Zalety: bardzo szybka zbieżność (kwadratowa p=2)

Wady: trzeba obliczać pierwszą pochodną funkcji f(x) (znajomość analityczna

lub przybliżenie numeryczne) i dobrze wybrać punkt startowy x

0

.

M.Pyrz Metody numeryczne w mechanice – Równania nieliniowe 10.2011

Uwagi

Metoda siecznych

.

Metoda stycznych

1

1

1

1

1

(

) (

)

(

)

(

)

(

)

(

)

(

)

(

)

n

n

n

n

n

n

n

n

n

n

n

n

n

f x

x

x

f x

x

x

x

f x

f x

f x

f x

x

x

−

+

−

−

−

−

= −

= −

−

−

−

1

(

)

(

)

n

n

n

n

f x

x

x

f x

+

= −

′

Przybliżenie pochodnej

28

Strategia rozwiązywania trudnych przypadków:

najpierw wolno zbieżna ale pewna metoda bisekcji, następnie szybka metoda siecznych

.

M.Pyrz Metody numeryczne w mechanice – Równania nieliniowe 10.2011

(

)

n

f x

′

Metody szukania punktu stałego

Równanie f(x)=0 może być zapisane w równoważnej postaci

x=F(x)

,

gdzie F : R R jest pewną funkcją zmiennej rzeczywistej x.

Poszukiwanie x spełniającego f(x)=0 jest równoważne szukaniu takiego x, aby na

przykład spełnić równanie

x = x +a f(x)

(stała a ≠ 0).

29

Można zatem przyjąć F(x)= x +a f(x).

Pierwiastek x* równania x=F(x) nazywany jest

„punktem stałym" funkcji F(x).

F(x) może posiadać wiele punktów stałych.

M.Pyrz Metody numeryczne w mechanice – Równania nieliniowe 10.2011

1

0

2

1

1

(

)

(

)

(

)

n

n

x

F x

x

F x

x

F x

−

=

=

=

⋯

algorytm

Metody szukania punktu stałego -przykład

30

M.Pyrz Metody numeryczne w mechanice – Równania nieliniowe 10.2011

Rys. 2.12

Metody szukania punktu stałego - zbieżność

Warunek wystarczający zbieżności

Jeżeli F(x) : R R jest ciągła w przedziale [a,b], oraz

to metoda punktu stałego jest zbieżna.

( )

1

[ , ]

F x

x

a b

′

<

∀ ∈

31

M.Pyrz Metody numeryczne w mechanice – Równania nieliniowe 10.2011

Metody przyspieszania zbieżności

Przyspieszanie zbieżności metod pierwszego rzędu polega na budowaniu na

podstawie ciągu wartości przybliżonych {x

i

} nowego ciągu wartości {x°

i

}, który

zbiega się do pierwiastka x* szybciej niż ciąg początkowy.

Metoda Aitken’a

32

Metoda Aitken’a

Metoda Steffensen’a

Rozwiazywanie układu równań nieliniowych

- patrz koniec rozdziału " Układy równań liniowych"

M.Pyrz Metody numeryczne w mechanice – Równania nieliniowe 10.2011