Metody Numeryczne w Fizyce 1 Minima i maksima funkcji jednej zmiennej 1
Znajdowanie miejsc zerowych oraz minimów i maksimów funkcji
W grze Space Quest 2 była scenka ze skokiem z rozkołysanej liny. Ciekawostka gracz powinien puścić się w chwili
maksymalnego wychylenia, po czym następował lot (sic !) jeszcze dalej. W związku z tym zajmiemy się takim zadaniem:
Nale\y podać wartość kąta wychylenia liny ą
L
(mierzonego od pionu), przy którym zasięg skoku
H
będzie maksymalny. Zakładamy, \e znane są:
długość liny L, wysokość punktu zaczepienia liny
H oraz amplituda wahań ą0.
Analityczne uzyskanie wzoru na zasięg w funkcji
kąta nie jest zadaniem prostym. Jak widać potrzebne
będą metody numeryczne znajdowania minimum funkcji
(je\eli potrzebujemy tak jak w tym przypadku
maksimum mo\emy pomno\yć funkcję przez 1, albo
w prosty sposób zmodyfikować przedstawiane metody).
Metody znajdowania minimów funkcji jednej zmiennej
Zakładamy, \e badana funkcja ma minimum w przedziale (a,b ) (dla x = ą), przy czym w przedziale
jest malejąca,
a w przedziale <ą,b> rosnąca. W takiej sytuacji wystarczy znalezć wartości funkcji dla dwóch punktów le\ących
wewnątrz przedziału (oznaczmy je t1 i t2, przy czym t1 < t2). Jeśli f (t1) < f (t2) to ą 0 , a jeśli f (t1) > f (t2) to ą 0 .
Warto podkreślić, \e wszystkie omawiane metody znajdowania minimów i miejsc zerowych wymagają jakiejś początkowej
znajomości badanej funkcji trzeba znać chocia\ w przybli\eniu przedział, w którym znajduje się minimum bądz miejsce
zerowe.
Nale\y zwrócić uwagę na wyznaczenie przedziału (a,b ) musi on być tak wybrany, aby minimum znajdowało się
wewnątrz przedziału.
Podczas poszukiwań nale\y z góry zało\yć po\ądaną dokładność obliczeń umo\liwi to sformułowanie kryterium
przerwania obliczeń gdy długość badanego przedziału (wewnątrz którego znajduje się minimum) będzie dwa razy
większa od zało\onej dokładności. Jako wynik podamy poło\enie środka przedziału ą zało\ony błąd.
Metoda podziału na trzy części
Punkty wewnętrzne obliczamy jako: t1 = b a + a b
oraz t2 = a a + b b (czyli dzielimy przedział na trzy
równe części. W zale\ności od ró\nicy f (t1) i f (t2) punkt
a przenosimy do t1 lub punkt b przenosimy do t2. Widać,
\e ka\da kolejna iteracja zmniejsza przedział o a jego
długości. W ka\dej kolejnej iteracji mamy ju\ obliczone
wartości funkcji w punktach brzegowych przedziału, ale
drugi punkt wewnętrzny "marnuje się" musimy
obliczać dwie następne wartości wewnętrzne.
Jako przykład posłu\y nam funkcja f (x) = (x 2)2+1,
dla której wybierzemy przedział początkowy (0,3). Zało-
\ymy te\, \e zadowala nas błąd względny 0,001
wiersz wytłuszczony odpowiada spełnieniu tego
a t1 t2 b
warunku.
f(t1) > f(t2)
Andrzej Brozi, Instytut Fizyki Politechniki Aódzkiej
Metody Numeryczne w Fizyce 1 Minima i maksima funkcji jednej zmiennej 2
Metoda podziału na 3 części. a = 0 długość przedziału = 3
Funkcja: y = (x-2)2 + 1 b = 3 dokładna wartość pierwiastka = 2
N aN bN błąd błąd
0 t1 f(t1) t2 f(t2) 0 3 pierw. szacow. dokł.
1 1 2 2 1 1 3 2 0.5 0
2 1.666667 1.111111 2.333333 1.111111 1.666667 3 2.333333 0.333333 0.16667
3 2.111111 1.012346 2.555556 1.308642 1.666667 2.555556 2.111111 0.222222 0.05556
4 1.962963 1.001372 2.259259 1.067215 1.666667 2.259259 1.962963 0.148148 -0.0185
5 1.864198 1.018442 2.061728 1.00381 1.864198 2.259259 2.061728 0.098765 0.03086
6 1.995885 1.000017 2.127572 1.016275 1.864198 2.127572 1.995885 0.065844 -0.0021
7 1.951989 1.002305 2.039781 1.001582 1.951989 2.127572 2.039781 0.043896 0.01989
8 2.010517 1.000111 2.069044 1.004767 1.951989 2.069044 2.010517 0.029264 0.00526
9 1.991007 1.000081 2.030026 1.000902 1.951989 2.030026 1.991007 0.019509 -0.0045
10 1.978001 1.000484 2.004014 1.000016 1.978001 2.030026 2.004014 0.013006 0.00201
11 1.995343 1.000022 2.012684 1.000161 1.978001 2.012684 1.995343 0.008671 -0.0023
12 1.989562 1.000109 2.001123 1.000001 1.989562 2.012684 2.001123 0.005781 0.00056
13 1.99727 1.000007 2.004977 1.000025 1.989562 2.004977 1.99727 0.003854 -0.0014
14 1.994701 1.000028 1.999839 1 1.994701 2.004977 1.999839 0.002569 -8e-05
15 1.998126 1.000004 2.001552 1.000002 1.998126 2.004977 2.001552 0.001713 0.00078
16 2.00041 1 2.002693 1.000007 1.998126 2.002693 2.00041 0.001142 0.0002
17 1.999648 1 2.001171 1.000001 1.998126 2.001171 1.999648 0.000761 -0.0002
18 1.999141 1.000001 2.000156 1 1.999141 2.001171 2.000156 0.000507 0.00008
Wykresy poni\ej ilustrują przebieg zwę\ania przedziału podczas rozwiązywania przykładowego problemu oraz przebieg
zmniejszania się błędu rozwiązania.
Poszukiwanie minimum funkcji metod podziau na trzy cz Ńci.
Ń
Ń
Ń
Przebieg zaw ó u w metodzie podziau na trzy.
óania przedzia
ó
ó
Zadana dok Ńł
adnoŃł 0.100000. Przedzia pocztkowy 0.000000 - 3.000000.
Ńł
Ńł
5 3
rozw. dok.
a
4.5
b
2.5
4 rozw. num.
3.5
2
3
2.5
1.5
2
1
1.5
1
0.5
0.5
0
a a1 a3 a4 a6 b7 b5 b2 b
0
0 0.5 1 1.5 2 2.5 3 1 2 3 4 5 6 7 8
Wynik programu 1.960219 +- 0.087791, liczba kroków 7.
Numer kroku
Wynik dokadny: 2.000000, faktyczny b
d bezwzgl dny -0.039781, wzgl dny -0.019890.
Andrzej Brozi, Instytut Fizyki Politechniki Aódzkiej
Metody Numeryczne w Fizyce 1 Minima i maksima funkcji jednej zmiennej 3
Metoda połowienia (podział na cztery części)
W tym przypadku jako punkty wewnętrzne przyjmujemy: t1 = a + ź b, t2 = a + b oraz t3 = ź a + b. W ten
sposób mamy pięć punktów: a, b, t1, t2 i t3. Nowymi granicami przedziału są punkty sąsiadujące z tym z nich, w którym
funkcja ma najmniejszą wartość. Warto zwrócić uwagę na to, \e po wybraniu nowego przedziału zawsze mamy jego punkty
brzegowe oraz punkt środkowy trzeba więc obliczyć jeszcze dwa tak samo jak przy podziale na trzy części, ale
długość kolejnego przedziału jest połową długości początkowej.
Pisząc program wykorzystujący tę metodę warto
zauwa\yć, \e istnieją szczególne przypadki (dla funkcji
takiej jak na rysunku: f (t1) $ f (a) i f (t3) $ f (b)), w
których na podstawie zało\enia o monotoniczności
funkcji wokół minimum mo\na stwierdzić, \e minimum
znajduje się w przedziale (a,t1) lub (t3,b) w takim
korzystnym przypadku przedział maleje nie dwu- ale
czterokrotnie, trzeba jednak wówczas wewnątrz nowego
przedziału obliczyć trzy, a nie dwa punkty zysk na
czasie jest więc nieco mniejszy.
W celu ilustracji posłu\ymy się tą samą funkcją
a t1 t2 t3 b
przykładową co poprzednio.
a t1 t2 t3 b
Tak samo jak poprzednio wytłuszczony wiersz w
tabeli odpowiada spełnieniu warunku dokładności jak
widać w metodzie podziału na trzy części warunek ten został spełniony w kroku 17, a w przypadku metody połowienia ju\
w kroku 10.
Metoda połowienia (podział na 4 części) a =0 długość przedziału = 3
funkcja:y = (x-2)2 + 1 b =3 dokł. wart. pierwiastka = 2
N a' b' pierwias- błąd błąd
t1 f(t1) t2 f(t2) t3 f(t3)
tek szacow. dokł.
0 0 3
1 0.75 2.5625 1.5 1.25 2.25 1.0625 1.5 3 2.25 0.375 0.125
2 1.875 1.0156 2.25 1.0625 2.625 1.3906 1.5 2.25 1.875 0.1875 -0.06
3 1.6875 1.0977 1.875 1.0156 2.0625 1.0039 1.875 2.25 2.0625 0.0938 0.031
4 1.9688 1.001 2.0625 1.0039 2.1563 1.0244 1.875 2.0625 1.9688 0.0469 -0.02
5 1.9219 1.0061 1.9688 1.001 2.0156 1.0002 1.9688 2.0625 2.0156 0.0234 0.008
6 1.9922 1.0001 2.0156 1.0002 2.0391 1.0015 1.9688 2.0156 1.9922 0.0117 0
7 1.9805 1.0004 1.9922 1.0001 2.0039 1 1.9922 2.0156 2.0039 0.0059 0.002
8 1.998 1 2.0039 1 2.0098 1.0001 1.9922 2.0039 1.998 0.0029 0
9 1.9951 1 1.998 1 2.001 1 1.998 2.0039 2.001 0.0015 5e-04
10 1.9995 1 2.001 1 2.0024 1 1.998 2.001 1.9995 0.0007 0
11 1.9988 1 1.9995 1 2.0002 1 1.9995 2.001 2.0002 0.0004 1e-04
12 1.9999 1 2.0002 1 2.0006 1 1.9995 2.0002 1.9999 0.0002 0
13 1.9997 1 1.9999 1 2.0001 1 1.9999 2.0002 2.0001 9e-05 3e-05
Andrzej Brozi, Instytut Fizyki Politechniki Aódzkiej
Metody Numeryczne w Fizyce 1 Minima i maksima funkcji jednej zmiennej 4
14 2 1 2.0001 1 2.0002 1 1.9999 2.0001 2 5e-05 0
15 1.9999 1 2 1 2 1 2 2.0001 2 2e-05 8e-06
16 2 1 2 1 2 1 2 2 2 1e-05 0
17 2 1 2 1 2 1 2 2 2 6e-06 2e-06
18 2 1 2 1 2 1 2 2 2 3e-06 0
Poszukiwanie minimum funkcji metod podzia Ń
Ń
u na cztery cz Ńci.
Ń
Przebieg zaw ó
óania przedziau w metodzie podziau na cztery.
ó
ó
Zadana dok Ńł
adnoŃł 0.100000. Przedzia pocztkowy 0.000000 - 3.000000.
Ńł
Ńł
5 3
rozw. dok.
a
4.5
b
2.5
4 rozw. num.
3.5
2
3
2.5
1.5
2
1
1.5
1
0.5
0.5
0
a a1 a3 b4 b2 b
0
0 0.5 1 1.5 2 2.5 3 1 2 3 4 5
Wynik programu 1.968750 +- 0.093750, liczba kroków 4.
Numer kroku
Wynik dokadny: 2.000000, faktyczny b
d bezwzgl dny -0.031250, wzgl dny -0.015625.
Metoda Johnsona (optymalnych podziałów)
W metodzie tej do określenia poło\enia dwóch punktów wewnętrznych wykorzystuje się tzw. liczby Fibonacciego. Jest
to ciąg liczb zdefiniowanych wzorem rekurencyjnym:
F0 = F1 = 1
Fi = Fi 1 + Fi 2 i = 2,3,...
czyli: Fi = 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, ...
We wzorach poni\ej indeksy górne w nawiasach będą oznaczać numer iteracji: a(0), b(0) są więc wartościami
początkowymi granic. Wymaganą dokładność (błąd bezwzględny) oznaczymy przez ". Obliczamy wielkość c =
= (b(0) a(0))/", po czym znajdujemy taką liczbę N aby FN 1 < c # FN. Następnie kolejne punkty podziału obliczamy jako:
FN&i &1
(i
t1 ) ' b&a %a
FN&i %1
i ' 1,2,,N&2
FN&i
(i
t2 ) ' b&a %a
FN&i %1
(i (i (i
Je\eli f (t1 )) # f (t2 )) to wartość a pozostaje bez zmiany, a za b podstawiamy t2 ) ;
(i (i (i
Jeśli zaś f (t1 )) > f (t2 )) to nie zmieniamy b, natomiast a = t1 ).
Liczbę iteracji znamy z góry poniewa\ do wyznaczenia N wykorzystaliśmy dopuszczalny błąd bezwzględny wyniku.
Mo\na pokazać, \e w metodzie tej wykorzystywane są oba znalezione punkty podziału jeden z nich staje się nową
granicą przedziału, a drugi pokrywa się z punktem kolejnego podziału dzięki temu w ka\dej kolejnej iteracji konieczne
jest obliczenie tylko jednej nowej wartości funkcji.
Stosując tę metodę do naszego przykładu otrzymamy:
Andrzej Brozi, Instytut Fizyki Politechniki Aódzkiej
Metody Numeryczne w Fizyce 1 Minima i maksima funkcji jednej zmiennej 5
Metoda optymalnych podziałów a = 0 błąd = 0.001 c = 3000
funkcja:y=(x-2)2+1 b = 3 przedział = 3 N = 18
N Fi
t1 f (t1) t2 f (t2) aN bN x0 faktyczny błąd sza-
blad cowany
0 1
1 1 1.999226 2.000387 1.999807 -0.0001
2 2 2.000387 1 2.000387 1 1.999226 2.001548 2.000387 0.000193 0.00058
3 3 2.000387 1 2.001548 1.000002 1.999226 2.002709 2.000967 0.000484 0.000871
4 5 1.999226 1.000001 2.000387 1 1.996904 2.002709 1.999807 -0.0001 0.001451
5 8 2.000387 1 2.002709 1.000007 1.996904 2.006192 2.001548 0.000774 0.002322
6 13 1.996904 1.00001 2.000387 1 1.991099 2.006192 1.998646 -0.00068 0.003773
7 21 1.991099 1.000079 1.996904 1.00001 1.981811 2.006192 1.994002 -0.003 0.006095
8 34 1.996904 1.00001 2.006192 1.000038 1.981811 2.021285 2.001548 0.000774 0.009868
9 55 1.981811 1.000331 1.996904 1.00001 1.95743 2.021285 1.989358 -0.00532 0.015964
10 89 1.996904 1.00001 2.021285 1.000453 1.95743 2.060759 2.009094 0.004547 0.025832
11 144 2.021285 1.000453 2.060759 1.003692 1.95743 2.124613 2.041022 0.020511 0.041796
12 233 1.95743 1.001812 2.021285 1.000453 1.854102 2.124613 1.989358 -0.00532 0.067628
13 377 2.021285 1.000453 2.124613 1.015528 1.854102 2.291796 2.072949 0.036474 0.109423
14 610 1.854102 1.021286 2.021285 1.000453 1.583591 2.291796 1.937693 -0.03115 0.177051
15 987 1.583591 1.173396 1.854102 1.021286 1.145898 2.291796 1.718847 -0.14058 0.286474
16 1597 1.854102 1.021286 2.291796 1.085145 1.145898 3 2.072949 0.036474 0.463526
17 2584 1.145898 1.729491 1.854102 1.021286 0 3 0.75
18 4181
Interpretując powy\szą tabelę trzeba pamiętać, \e do obliczania kolejnych przedziałów wykorzystujemy najpierw
końcowe wartości liczb z ciągu Fibonacciego z tej przyczyny prawe kolumny zapisywane są od dołu. Stanowi to
równocześnie pewien niedostatek tej metody jeśli zmienimy zdanie co do po\ądanej dokładności obliczenia musimy
zacząć od początku.
Przebieg zaw ó Przebieg zaw óania przedziau w metodzie optymalnych podziaów.
óania przedziau w metodzie optymalnych podziaów. ó
ó ó
ó ó
2.5 2.01
wartoŃł dokadna wartoŃł dokadna
Ńł Ńł
Ńł Ńł
Ńł Ńł
wartoŃł obliczona wartoŃł obliczona
Ńł Ńł
Ńł Ńł
Ńł Ńł
2.008
zakresy zakresy
przedziaów przedziaów
2.006
2.004
2
2.002
2
1.998
1.5
1.996
1.994
1.992
1 1.99
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 8 9 10 11 12 13 14 15 16
numer kroku numer kroku
Andrzej Brozi, Instytut Fizyki Politechniki Aódzkiej
Metody Numeryczne w Fizyce 1 Minima i maksima funkcji jednej zmiennej 6
Drugi rysunek jest powiększeniem końcowego fragmentu pierwszego, dopiero tu widać, \e obliczona wartość xmin ró\ni
się od wartości faktycznej.
Metoda złotego podziału
Metoda ta jest dość podobna do poprzedniej, równie\ wykorzystuje oba punkty poprzedniego podziału. Jednak zamiast
5& 1
liczb z ciągu Fibonacciego wykorzystuje się liczbę ( ' . 0,62 ).
2
Algorytm wygląda następująco:
1
t1 ' a % 1& b&a
1
t2 ' b & 1& b&a
(i (i
Je\eli f (t1 )) # f (t2 )) to wartość a nie zmienia się, natomiast:
(i
b ' t2 )
(i (i
t2 %1) ' t1 )
(i
t1 %1) ' a % (1&)(b&a)
(i (i
Je\eli natomiast f (t1 )) > f (t2 )) to nie zmienia się wartość b, zaś:
(i
a ' t1 )
(i (i
t1 %1) ' t2 )
(i
t2 %1) ' b & (1&)(b&a)
Udowodniono, \e podczas poszukiwania minimum funkcji z identyczną dokładnością metodami złotego podziału i
Johnsona, metoda złotego podziału będzie wymagać obliczenia wartości funkcji co najwy\ej o jeden raz więcej.
W przykładzie wykorzystuję funkcję identyczną jak poprzednio.
Metoda złotego podziału a = 0 = 0.61803399
funkcja: y=(x-2)2+1 b = 3 przedział = 3 x0 dokł.= 2
N
t1 f(t1) t2 f(t2) aN bN x0 wzgl. błąd wzgl. błąd
dokł. szac.
0
1 1.145898 1.72949 1.854102 1.021286 0 3 1.5 -0.25 0.75
2 1.854102 1.021286 2.291796 1.085145 1.145898 3 2.072949 0.036475 0.463525
3 1.583592 1.173396 1.854102 1.021286 1.145898 2.291796 1.718847 -0.14058 0.286475
4 1.854102 1.021286 2.021286 1.000453 1.583592 2.291796 1.937694 -0.03115 0.177051
5 2.021286 1.000453 2.124612 1.015528 1.854102 2.291796 2.072949 0.036475 0.109424
6 1.957428 1.001812 2.021286 1.000453 1.854102 2.124612 1.989357 -0.00532 0.067627
7 2.021286 1.000453 2.060753 1.003691 1.957428 2.124612 2.04102 0.02051 0.041796
8 1.996894 1.00001 2.021286 1.000453 1.957428 2.060753 2.00909 0.004545 0.025831
9 1.981819 1.000331 1.996894 1.00001 1.957428 2.021286 1.989357 -0.00532 0.015965
10 1.996894 1.00001 2.006211 1.000039 1.981819 2.021286 2.001553 0.000776 0.009867
11 1.991136 1.000079 1.996894 1.00001 1.981819 2.006211 1.994015 -0.00299 0.006098
Andrzej Brozi, Instytut Fizyki Politechniki Aódzkiej
Metody Numeryczne w Fizyce 1 Minima i maksima funkcji jednej zmiennej 7
12 1.996894 1.00001 2.000453 1 1.991136 2.006211 1.998674 -0.00066 0.003769
13 2.000453 1 2.002653 1.000007 1.996894 2.006211 2.001553 0.000776 0.002329
14 1.999094 1.000001 2.000453 1 1.996894 2.002653 1.999773 -0.00011 0.00144
15 2.000453 1 2.001293 1.000002 1.999094 2.002653 2.000873 0.000437 0.00089
16 1.999934 1 2.000453 1 1.999094 2.001293 2.000193 0.000097 0.00055
17 1.999613 1 1.999934 1 1.999094 2.000453 1.999773 -0.00011 0.00034
18 1.999934 1 2.000132 1 1.999613 2.000453 2.000033 0.000017 0.00021
Przewaga tej metody nad metodą optymalnego podziału polega na tym, \e w razie potrzeby mo\na przedłu\yć
obliczenia w celu uzyskania wy\szej dokładności.
Poszukiwanie minimum funkcji metod z
otego podziau.
Zadana dok Ńł
adnoŃł 0.100000. Przedzia pocztkowy 0.000000 - 3.000000.
Ńł
Ńł
Przebieg zaw ó otego podziau.
óania przedziau w metodzie z
ó
ó
5
3
rozw. dok.
4.5
a
b
4 2.5
rozw. num.
3.5
2
3
2.5
1.5
2
1.5
1
1
0.5 0.5
0
a a1 a3 a4 a6 b5 b2 b
0 0.5 1 1.5 2 2.5 3 0
1 2 3 4 5 6 7
Wynik programu 2.041020 +- 0.083592, liczba kroków 6.
Wynik dokadny: 2.000000, faktyczny bd bezwzgl dny 0.041020, wzgl dny 0.020510. Numer kroku
Aatwo mo\na sobie wyobrazić sytuację, gdy nie będzie nam chodzić o uzyskanie maksymalnej długości skoku, ale o
skok o dokładnie zadanej długości (np. s). Numerycznie sprowadzi się to do zadania polegającego na znalezieniu miejsca
zerowego funkcji zasięg(ą) s, albo jeśli uda się wyznaczyć pochodną zasięgu będziemy szukać jej miejsca zerowego.
Metody znajdowania miejsc zerowych funkcji jednej zmiennej
Zało\ymy, \e funkcja f(x) jest dana w taki sposób, \e dla ka\dego x mo\emy wyliczyć f(x). Przedstawiane dalej metody
wymagają dodatkowo abyśmy znali dwie wartości x1 i x2 takie, \e znaki f(x1) i f(x2) są przeciwne wyznaczają więc one
przedział, w którym musi zawierać się pierwiastek. Warto zauwa\yć, \e nie mamy tu \adnej mo\liwości znalezienia
pierwiastków, w których pochodna funkcji jest równa zeru a więc np. dla funkcji y = x2.
Andrzej Brozi, Instytut Fizyki Politechniki Aódzkiej
Metody Numeryczne w Fizyce 1 Minima i maksima funkcji jednej zmiennej 8
Metoda połowienia (bisekcji)
x2
Jest to chyba najprostsza z metod. Jeśli znamy dwie wartości
f (x)
x1 i x2 to obliczamy wartość funkcji f dla x=(x2+x1)/2. Jeśli znak
funkcji w punkcie środkowym jest zgodny z x1 to podstawiamy
x1=x, a jeśli zgodny z x2 to x2=x i powtarzamy całą procedurę.
(Oczywiście jeśli f(x) wynosi zero to mamy wynik.) Problemem
jest dość powolna zbie\ność tej metody. Ka\de przeprowadzenie
x1
tej procedury zmniejsza o połowę przedział, w którym musi xN
N
N
x N
O
znajdować się szukany pierwiastek tzn. dokładność wyniku xO x
O
O
poprawia się o czynnik 2 (albo inaczej o 1 rząd, ale w układzie pierwiastek
dwójkowym). Aby więc uzyskać poprawę dokładności o rząd
wielkości (ale w układzie dziesiętnym) musimy przeprowadzić
procedurę ok. 3.3-krotnie.
Jako przykład rozwa\ymy znalezienie miejsca zero- Wykres funkcji badanej w przykadach poszukiwania miejsc zerowych
1
lnx
wego funkcji f x ' sin . Zajmiemy się pier- 0.8
x
0.6
wiastkiem x0 = 1 (znamy jego dokładną wartość
0.4
dzięki czemu łatwo będzie wyznaczyć faktycznie po-
0.2
pełniany błąd).
0
We wszystkich następnych przykładach jako prze-
dział początkowy będę wybierać przedział (0.98,1.1) i -0.2
\ądać dokładności 10 6.
-0.4
-0.6
-0.8
-1
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1
Metoda bisekcji a =0.98 b =1.1
błąd wynik
a f(a) x0 f(x0) b f(b)
obliczony
0 0.98 -0.0206135 1.04 0.03770329 1.1 0.08653724
1 0.98 -0.0206135 1.01 0.00985165 1.04 0.03770329
2 0.98 -0.0206135 0.995 -0.0050377 1.01 0.00985165
3 0.995 -0.0050377 1.0025 0.00249065 1.01 0.00985165
4 0.995 -0.0050377 0.99875 -0.0012523 1.0025 0.00249065
5 0.99875 -0.0012523 1.000625 0.00062441 1.0025 0.00249065
6 0.99875 -0.0012523 0.9996875 -0.0003126 1.000625 0.00062441
7 0.9996875 -0.0003126 1.00015625 0.00015621 1.000625 0.00062441
8 0.9996875 -0.0003126 0.99992188 -0.0000781 1.00015625 0.00015621
9 0.99992188 -0.0000781 1.00003906 0.00003906 1.00015625 0.00015621
10 0.99992188 -0.0000781 0.99998047 -0.0000195 1.00003906 0.00003906
11 0.99998047 -0.0000195 1.00000977 0.00000977 1.00003906 0.00003906
12 0.99998047 -0.0000195 0.99999512 -0.0000049 1.00000977 0.00000977
Andrzej Brozi, Instytut Fizyki Politechniki Aódzkiej
Metody Numeryczne w Fizyce 1 Minima i maksima funkcji jednej zmiennej 9
13 0.99999512 -0.0000049 1.00000244 0.00000244 1.00000977 0.00000977
14 0.99999512 -0.0000049 0.99999878 -0.0000012 1.00000244 0.00000244
15 0.99999878 -0.0000012 1.00000061 0.00000061 1.00000244 0.00000244
16 0.99999878 -0.0000012 0.99999969 -0.0000003 1.00000061 0.00000061 + 0.00000092 0.99999969
Jak widać zadowalającą zbie\ność osiągnęliśmy dopiero po wykonaniu 16 iteracji. Warto zauwa\yć, \e obliczony błąd
rozwiązania jest w tym przypadku trzykrotnie większy ni\ faktycznie popełniony.
Poszukiwanie miejsca zerowego funkcji metod bisekcji.
Przebieg zaw ó
óania przedziau w metodzie bisekcji.
ó
ó
Zadana dokadnoŃł 0.050000. Przedzia pocztkowy 0.500000 - 2.000000.
Ńł
Ńł
Ńł
0.5 2
rozw. dok.
a
b
rozw. num.
0 1.5
a a2 a4 b5 b3 b1 b
-0.5 1
-1 0.5
0.5 1 1.5 2 1 2 3 4 5 6
Wynik programu 0.992188 +- 0.023438, liczba kroków 5.
Numer kroku
Wynik dokadny: 1.000000, faktyczny bd bezwzgl dny -0.007813, wzgl dny -0.007813.
Poszukiwanie miejsca zerowego funkcji metod regula falsi.
Zadana dokadnoŃł 1.0e-002. Przedzia pocztkowy 0.50 - 2.00.
Ńł
Ńł
Ńł
0.5
Regula falsi
Metoda zwana regula falsi jest swego rodzaju mo-
0
a b8 b4 b3 b2 b1 b
b7b5
b6
dyfikacją wy\ej opisanej metody połowienia (bisekcji).
Zamiast wybierać za ka\dym razem nowy argument w po-
łowie odległości między poprzednią parą wybiera się go w
miejscu przecięcia osi x z prostą przechodzącą przez
-0.5
punkty (x1,f(x1)) i (x2,f(x2)) czyli przez poprzednią parę
punktów.
x1f(x2)&x2f(x1)
x '
-1
f(x2)&f(x1)
0.5 1 1.5 2
Wynik programu 1.007 +- 7.1e-003, liczba kroków 8.
Wynik dokadny: 1.000, faktyczny bd bezwzgl dny 6.8e-003, wzgl dny 6.8e-003.
Przebieg poszukiwania pierwiastka w metodzie regula falsi.
Zasada wyboru nowej granicy przedziału jest
2
rozw. dok.
identyczna jak poprzednio do punktu x przenosimy tę
a
granicę, w której znak funkcji jest zgodny ze znakiem w b
rozw. num.
nowym punkcie.
Osobną kwestią pozostaje określenie dokładności z jaką
1.5
wyznaczono miejsce zerowe. Jak łatwo się domyśleć
mo\liwa jest sytuacja, w której szukane miejsce zerowe
le\y znacznie bli\ej jednego z punktów brzegowych
przedziału. Jeśli kolejne punkty le\ą dostatecznie blisko
1
siebie (a "dostateczność" tej bliskości zale\y zwłaszcza od
szybkości zmian pochodnej funkcji) to jako wynik osta-
teczny mo\na podać ostatni wyznaczony punkt (będący
0.5
jedną z granic przedziału), a jego błąd mo\na określić jako
1 2 3 4 5 6 7 8 9
Numer kroku
stosunek wartości funkcji w tym punkcie do wartości jej
Andrzej Brozi, Instytut Fizyki Politechniki Aódzkiej
Metody Numeryczne w Fizyce 1 Minima i maksima funkcji jednej zmiennej 10
pochodnej (wyznaczonej na podstawie punktu ostatniego i przedostatniego).
f (xk) xk&xk&1
/ / @ (xk )/
" / . .
/ /
0 0
0 0
/pierwiastek&xk/ /f
) 0 0
0 0 f (xk)&f (xk&1)
f (xk) 0 0
0 0
0 0
0 0
Regula falsi a =0.98 b =1.1 błąd wynik błąd
(poch. (poch.
a f(a) x0 f(x0) b f(b)
num.) dokł.)
0 0.98 -0.0206135 1.00308546 0.00307123 1.1 0.08653724
1 0.98 -0.0206135 1.00009195 0.00009193 1.00308546 0.00307123
2 0.98 -0.0206135 1.00000274 0.00000274 1.00009195 0.00009193
3 0.98 -0.0206135 1.00000008 0.00000008 1.00000274 0.00000274 + 8e 8 1+8e 8 8e 8
Metoda siecznych
Poszukiwanie miejsca zerowego funkcji metod siecznych (od lewej).
W zasadzie jest to ta sama metoda (regula falsi)
Zadana dokadnoŃł 1.0e-002. Przedzia pocztkowy 0.40 - 2.00.
Ńł
Ńł
Ńł
0.5
pozbawiona jednak \ądania, aby na obu krańcach
kolejnego przedziału funkcja miała przeciwne znaki
(\ądanie to pozostaje jednak w mocy w przypadku krańców
przedziału początkowego dając gwarancję, \e wewnątrz
0
a 3 65 4 2 1 b
7
przedziału w ogóle jest czego szukać). Jako punkty
brzegowe przedziału słu\ą zawsze dwa ostatnio
wyznaczone punkty niezale\nie od tego czy funkcja ma
w nich znaki przeciwne czy jednakowe. Na rysunku obok
-0.5
pierwsza sieczna została oczywiście poprowadzona przez
punkty a i b i wyznaczyła punkt x1. Drugą sieczną (wyzna-
czającą punkt x2) poprowadzono przez punkty x1 i b
dlatego metoda została określona jako od lewej . Trzecią -1
0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Wynik programu 1.000 +- 1.9e-003, liczba kroków 7.
sieczną (wyznaczającą punkt x3) poprowadzono
Wynik dokadny: 1.000, faktyczny bd bezwzgl dny 5.1e-005, wzgl dny 5.1e-005.
oczywiście przez punkty x1 i x2.
W momencie kreślenia drugiej siecznej mamy swobodę
Przebieg poszukiwania pierwiastka w metodzie siecznych.
wyboru mo\emy poprowadzić ją przez punkty (x1,b)
1.6
rozw. dok
.
albo przez (a,x1) (w dalszych obliczeniach zawsze jest ju\
rozw. od lewej
1.5
jasne, które punkty są "dwoma ostatnimi"). W zale\ności
1.4
od tego wyboru dalsze przybli\enia mogą być inne jak
widać na dalszych ilustracjach.
1.3
Metoda ta jest szybciej zbie\na ni\ poprzednia (regula
1.2
falsi), ale te\ ma większe wymagania dotyczące "jakości"
1.1
przybli\enia początkowego, tzn. początkowy przedział
musi być dostatecznie wąski. Objawem niewłaściwego do- 1
boru przedziału początkowego jest szybki wzrost ró\nicy
0.9
pomiędzy kolejnymi przybli\eniami (metoda się "roz-
0.8
biega"), w programach realizujących tę metodę nale\y
0.7
zawrzeć procedury kontrolujące zachowanie się kolejnych
1 2 3 4 5 6 7 8
Numer kroku
przybli\eń wartości funkcji.
Sposób oszacowania błędu jest taki sam jak opisany
powy\ej (z wykorzystaniem wartości pochodnej badanej funkcji obliczonej dla ostatniego przedziału i wartości funkcji w
ostatnim obliczonym punkcie).
Andrzej Brozi, Instytut Fizyki Politechniki Aódzkiej
Metody Numeryczne w Fizyce 1 Minima i maksima funkcji jednej zmiennej 11
W przypadku gdy xn . xn 1 oraz fnfn 1 > 0 mogą wystąpić błędy związane z wzajemnym znoszeniem się składników.
Nale\y wówczas korzystać z wzoru:
xn&xn&1
xn%1 ' xn % hn gdzie: hn ' &fn fn & fn&1
fn&fn&1
Wykazano, \e chocia\ iloraz w definicji poprawki hn mo\e mieć małą dokładność, to jednak dominujący wpływ na błąd
poprawki ma błąd wartości f (xn) i wzór ten daje lepsze wyniki ni\ podany wy\ej.
Metoda siecznych a =0.98 b =1.1 błąd błąd
wynik
od lewej x0 f(x0) (poch.num.) (poch.dokł.)
0 1.00308546 0.00307123
1 1.00009195 0.00009193
2 0.99999957 -0.0000004 + -0.0000004 0.99999957 -0.0000004
Metoda siecznych a =0.98 b =1.1 błąd błąd
wynik
od prawej x0 f(x0) (poch.num.) (poch.dokł.)
0 1.00308546 0.00307123
1 0.99951938 -0.000481
2 1.00000223 0.00000223
3 1 1.6050e-09 + 1.6050e-09 1 1.6050e-09
Poszukiwanie miejsca zerowego funkcji metod siecznych (od prawej).
Poszukiwanie miejsca zerowego funkcji metod siecznych (od lewej).
Zadana dokadnoŃł 1.0e-002. Przedzia pocztkowy 0.80 - 2.00.
Ńł
Ńł
Ńł
Zadana dokadnoŃł 1.0e-002. Przedzia pocztkowy 0.80 - 2.00.
Ńł
Ńł
Ńł
1
0.5
0.8
0.6
0.4
0
a 3 5 2 1 b
4
0.2
0
2 a 3 6 8 5 14 b
7
-0.2
-0.5 -0.4
-0.6
-0.8
-1
-1 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
0.8 1 1.2 1.4 1.6 1.8 2
Wynik programu 1.000 +- 3.9e-003, liczba kroków 8.
Wynik programu 1.001 +- 8.4e-003, liczba kroków 5.
Wynik dokadny: 1.000, faktyczny bd bezwzgl dny 2.5e-004, wzgl dny 2.5e-004.
Wynik dokadny: 1.000, faktyczny bd bezwzgl dny 8.1e-004, wzgl dny 8.1e-004.
Poszukiwanie miejsca zerowego funkcji metod stycznych.
Znana postał pochodnej. Zadana dokadnoŃł 1.0e-004. Punkt pocztkowy 0.500.
ł Ńł
ł Ńł
ł Ńł
0.5
Przebieg poszukiwania pierwiastka w metodzie siecznych.
1.4
rozw. dok.
rozw. od lewej
1.2 rozw. od prawej
0
1
a 2 3 6 1
4
5
0.8
0.6
-0.5
0.4
0.2
-1
0.5 1 1.5 2
0
Wynik programu 1.000 +- 3.7e-006, liczba kroków 6.
1 2 3 4 5 6 7 8 9
Wynik dokadny: 1.000, faktyczny bd bezwzgl dny -2.0e-011, wzgl dny -2.0e-011.
Numer kroku
Andrzej Brozi, Instytut Fizyki Politechniki Aódzkiej
Metody Numeryczne w Fizyce 1 Minima i maksima funkcji jednej zmiennej 12
Metoda Newtona (metoda stycznych)
Metoda ta (znana te\ jako metoda Newtona-Raphsona) wymaga abyśmy znali przedział początkowy, w którym zawiera
się pierwiastek i w którym fN i fO mają stały znak. Z tego końca przedziału, w którym znak fO jest taki sam jak znak f
prowadzimy styczną do wykresu funkcji, a jej punkt przecięcia z osią x wyznacza kolejne przybli\enie szukanego pier-
wiastka przy takim wyborze punktu początkowego i jeśli spełnione są powy\sze zało\enia mamy gwarancję, \e znale-
ziony punkt nadal znajduje się wewnątrz początkowego przedziału.
Poszukiwanie miejsca zerowego funkcji metod stycznych.
Działanie tej metody szczególnie silnie zale\y od wy-
Ńł
Pochodna liczona numerycznie (krok 1.0e-002). Zadana dokadnoŃł 1.0e-004. Punkt pocztkowy 0.500.
Ńł
Ńł
0.5
boru punktu początkowego. Jej ideę mo\na zinterpretować
jako wykorzystanie rozwinięcia badanej funcji w szereg
Taylora i odrzucenie wyrazów rozwinięcia zawierających
pochodne wy\sze ni\ pierwsza. Takie postępowanie jest
0
a 2 3 6 1
4
5
dopuszczalne pod warunkiem, \e odległość od punktu po-
czątkowego jest bardzo mała wówczas wy\sze potęgi
małej wielkości "x "gwarantują" słuszność takiego
zało\enia. Jeśli jednak wybrane przez nas punkty począt-
-0.5
kowe są zbyt odległe od szukanego pierwiastka zachowanie
metody będzie zale\eć od konkretnej postaci funkcji.
Przykład takiej właśnie sytuacji znajduje się na rysunku
ni\ej kolejne punkty mogą nawet oddalać się od szu- -1
0.5 1 1.5 2
Wynik programu 1.000 +- 5.9e-006, liczba kroków 6.
kanego pierwiastka (w zilustrowanym przypadku zbiegają
Wynik dokadny: 1.000, faktyczny bd bezwzgl dny -1.0e-009, wzgl dny -1.0e-009.
się do pierwiastka poło\onego dalej). Mo\na te\ łatwo
wyobrazić sobie sytuację, w której omawiana metoda "zapętla" się nie zbli\ając się w kolejnych iteracjach do pierwiastka
ani od niego nie oddalając (np. gdy mamy do czynienia z funkcją o wykresie symetrycznym względem pierwiastka).
Wzór iteracyjny dla kolejnych wartości x ma postać: Przebieg poszukiwania pierwiastka w metodzie stycznych.
1.3
rozw. dok
.
f xi
1.25 rozw. numeryczne dokadne
xi%1 ' xi &
rozw. numeryczne przybl.
)
1.2
f xi
1.15
Błąd określamy w taki sam sposób jak w metodzie
1.1
regula falsi, jest to tym wygodniejsze \e z samej natury
metody Newtona musimy obliczać wartość pochodnej 1.05
funkcji w poprzednim punkcie.
1
Zastosowanie omawianej metody mo\e nie być
0.95
mo\liwe jeśli nie potrafimy z zadowalającą dokładnością
0.9
określić granic przedziału wokół szukanego pierwiastka,
0.85
jeśli jednak znajdziemy się ju\ dostatecznie blisko to do-
0.8
kładność przybli\enia podwaja się z ka\dą iteracją, w wielu
1 2 3 4 5 6 7
Numer kroku
wypadkach zatem korzystne mo\e być rozpoczęcie
poszukiwań pierwiastka którąś z metod mniej wra\liwych
na szerokość przedziału początkowego i udokładnienie uzyskanego wyniku za pomocą kilku (dosłownie) przebiegów
metody Newtona-Raphsona.
Andrzej Brozi, Instytut Fizyki Politechniki Aódzkiej
Metody Numeryczne w Fizyce 1 Minima i maksima funkcji jednej zmiennej 13
Poszukiwanie miejsca zerowego funkcji metod Poszukiwanie miejsca zerowego funkcji metod
stycznych. stycznych.
Znana postał pochodnej. Zadana dokadnoŃł 1.0e-007. Punkt pocztkowy 0.500. Pochodna liczona numerycznie (krok 5.0e-002). Zadana dokadnoŃł 1.0e-007. Punkt pocztkowy 0.500.
ł Ńł Ńł
ł Ńł Ńł
ł Ńł Ńł
0.5 0.5
0 0
a 2 3 7 1 210 a 1
4 7
5 8
6 9
4
5
6
3
-0.5 -0.5
-1 -1
0.5 1 1.5 2 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Wynik programu 1.000 +- NaN, liczba kroków 7. Wynik programu 0.342 +- 2.0e-008, liczba kroków 10.
B Ń B Ń
d wartoŃci zerowej 0.000e+000 d wartoŃci zerowej 4.344e-008
Ń Ń
Ń Ń
Przebieg poszukiwania pierwiastka w metodzie stycznych. Dwa rysunki powy\ej ró\nią się tylko tym, \e do
Zadana dokadnoŃł 1.0e-007. Punkt pocztkowy 0.500.
Ńł
Ńł
Ńł
1.6
obliczeń zilustrowanych po lewej stronie wykorzystano
rozw. numeryczne dokadne
rozw. numeryczne przybl.
znaną postać pochodnej (a więc jej wartości były liczone
1.4
dokładnie), a po prawej stronie pochodną liczono nume-
rycznie. Ró\nice wartości dokładnych i numerycznych
1.2
były wystarczające by, zupełnie poprawnie działająca
1
metoda, znajdowała dwa ró\ne (faktycznie istniejące)
pierwiastki.
0.8
Istnieje jeszcze jedno ograniczenie stosowalności tej
metody: musimy być w stanie określić w dowolnym
0.6
punkcie (nale\ącym do wybranego przedziału) nie tylko
0.4
wartość funkcji ale i jej pochodnej. Mo\na czywiście
liczyć pochodną numeryczne, ale wymaga to dwukrotnego
0.2
1 2 3 4 5 6 7 8 9 10 11
w ka\dym cyklu obliczania wartości funkcji w
Numer kroku
przypadku skomplikowanych funkcji wydłu\a to czas
obliczeń, a poza tym wprowadza ograniczenie
dopuszczalnych wartości kroku dx: w przypadku zbyt małych kroków wartości funkcji na krańcach przedziału ró\niłyby
się zbyt mało i błędy obcięcia uniemo\liwiłyby sensowne prowadzenie obliczeń, zaś w przypadku kroków zbyt du\ych
obliczenia mo\e uniemo\liwić zbytnia niedokładność obliczanej numerycznie wartości pochodnej.
Inną zaletą metody Newtona-Raphsona jest fakt, \e dość łatwo mo\na zaadaptować ją do przypadku funkcji większej
liczby zmiennych (o rozwiązywaniu układów równań będzie pózniej).
Metoda a =0.98 b =1.1
błąd błąd
Newtona wynik
(poch.num.) (poch.dokł.)
f(a) x0 f(x0) f'(x0) f'(a)
z pktu a
0 -0.02061 0.999409 -0.000591 1.001774 1.062043
1 0.999999 -5e-07 1.000002 + 5.23e-07 0.999999 5.23e-07
Metoda a =0.98 b =1.1
błąd błąd
Newtona wynik
(poch.num.) (poch.dokł.)
f(b) x0 f(x0) f'(x0) f'(b)
z pktu b
0 0.086537 0.983823 -0.016577 1.049863 0.744873
1 0.999612 -0.000388 1.001164
2 0.999999 -2e-07 1.0000007 + 2.25e-07 0.999999 2.25e-07
Andrzej Brozi, Instytut Fizyki Politechniki Aódzkiej
Metody Numeryczne w Fizyce 1 Minima i maksima funkcji jednej zmiennej 14
Metoda mieszana
Jak łatwo mo\na było zauwa\yć w metodzie Newtona i regula falsi do szukanego pierwiastka dochodzi się z
"przeciwnych stron" mo\na więc posłu\yć się obiema równocześnie. W metodzie stycznych (Newtona) ostatnia
obliczona wartość xk jest większa od pierwiastka, natomiast w metodzie siecznych jest ona od tego samego pierwiastka
mniejsza. Za pierwiastek mo\na więc przyjąć punkt w połowie między oboma przybli\eniami, a za jego błąd połowę ich
odległości.
Porównanie przedstawionych metod
Przedstawione metody były stosowane do znajdowania tego samego pierwiastka w takich samych warunkach. W
przypadku metody Newtona punkty brzegowe przedziału początkowego a i b były punktami początkowymi.
a = 0.98 błąd zało\ony = 1e-06
b = 1.1
Metoda bisekcji: Regula falsi:
podany wynik 9.999997e-01 podany wynik 1.00000008146
błąd obliczony 9.155273e-07 błąd obliczony 7.90342e-08
numer kroku 16 numer kroku 3
Metoda stycznych "od lewej" Metoda stycznych "od prawej"
podany wynik 9.999996e-01 podany wynik 1.0000000016
błąd obl.num. -4.261246e-07 błąd obl.num. 1.604988e-09
błąd obl.dokł. -4.260666e-07 błąd obl.dokł. 1.604983e-09
numer kroku 2 numer kroku 3
Metoda Newtona z pktu a Metoda Newtona z pktu b
podany wynik 0.999999694824 podany wynik 0.999999774702
błąd obliczony 5.230831e-07 błąd obliczony 2.252981e-07
numer kroku 1 numer kroku 2
Błędy w przedstawionych metodach
Warto zdawać sobie sprawę z faktu, \e zwłaszcza podczas obliczania pochodnych wykonujemy dzielenie dwóch liczb
będących ró\nicami. W chwili gdy zbli\ymy się ju\ bardzo do faktycznej wartości szukanego pierwiastka ró\nice pomiędzy
kolejnymi współrzędnymi x i odpowiadającymi im wartościami funkcji f (x ) są bardzo małe. W przypadku odejmowania
od siebie liczb o bardzo bliskich wartościach gwałtownie rośnie błąd jakim obarczony jest wynik.
Dla przykładu popatrzmy na wynik działania Turbo Pascala 7.0. W programie testowym zdefiniowane zostały trzy pary
stałych typu real, single, double i extended o wartościach:
1.1111111111211111111111 i 1.1111111111311111111111
1.1111111111121111111111 i 1.1111111111131111111111
1.1111111111112111111111 i 1.1111111111113111111111
po czym przeprowadzono odejmowanie mniejszej od większej dla poszczególnych typów. Poniewa\ dokładnie wiadomo
ile powinno wyjść obliczone zostały równie\ błędy bezwzględne i względne. Wyniki znajdują się poni\ej (błędy
Andrzej Brozi, Instytut Fizyki Politechniki Aódzkiej
Metody Numeryczne w Fizyce 1 Minima i maksima funkcji jednej zmiennej 15
bezwzględne z lewej, względne z prawej):
--------real--------
1.11111111113132210E+0000 B 1.11111111112040817E+0000 = 1.09139364212751389EB0011
B9.13936421277035EB0013
B9.13936421277413EB0002
1.11111111111313221E+0000 B 1.11111111111131322E+0000 = 1.81898940354585648EB0012
B8.18989403546377EB0013
B8.18989403546766EB0001
1.11111111111131322E+0000 B 1.11111111111131322E+0000 = 0.00000000000000000E+0000
1.00000000000031EB0013 1.00000000000000E+0000
BBBBBBBBsingleBBBBBBBB
1.11111116409301758E+0000 B 1.11111116409301758E+0000 = 0.00000000000000000E+0000
9.99999996004197EB0012 1.00000000000000E+0000
1.11111116409301758E+0000 B 1.11111116409301758E+0000 = 0.00000000000000000E+0000
9.99999996004197EB0013 1.00000000000000E+0000
1.11111116409301758E+0000 B 1.11111116409301758E+0000 = 0.00000000000000000E+0000
9.99999982451670EB0014 1.00000000000000E+0000
BBBBBBBBdoubleBBBBBBBB
1.11111111113111116E+0000 B 1.11111111112111116E+0000 = 1.00000008274037100EB0011
B8.27403710595934EB0019
B8.27403710595934EB0008
1.11111111111311112E+0000 B 1.11111111111211103E+0000 = 1.00008890058234101EB0012
B8.89005823410317EB0017
B8.89005823410317EB0005
1.11111111111131100E+0000 B 1.11111111111121108E+0000 = 9.99200722162640886EB0014
7.99277837359144EB0017 7.99277837359144EB0004
BBBBBBBBextendedBBBBBBBB
1.11111111113111111E+0000 B 1.11111111112111111E+0000 = 9.99999996004197200EB0012
3.99580279983147EB0020 3.99580279983147EB0009
1.11111111111311111E+0000 B 1.11111111111211111E+0000 = 9.99999996004197200EB0013
3.99580279985119EB0021 3.99580279985119EB0009
1.11111111111131111E+0000 B 1.11111111111121111E+0000 = 9.99999779163762703EB0014
2.20836237296903EB0020 2.20836237296903EB0007
Jak widać błędy względne mogą być bardzo du\e, a wynik podzielenia lub pomno\enia przez siebie dwóch liczb
obarczonych błędami sam jest obarczony błędem względnym będącym sumą błędów względnych czynników mo\e więc
być zupełnie bezu\yteczny.
Metoda kolejnych podstawień (metoda iteracyjna)
Rozwiązywane równanie, które poprzednio zapisywaliśmy
y=x
w postaci f(x) = 0 mo\na przekształcić do postaci: x = F(x).
Oczywiście zawsze gdy f(x) = 0 zachodzić będzie równie\
x = F(x).
Jeśli zało\ymy jakąś początkową wartość przybli\enia
F(x)
szukanego pierwiastka x0 to mo\emy zastosować zale\ność
rekursyjną:
xi%1 ' F (xi)
i mieć nadzieję, \e będzie się ona zbiegać do jego dokładnej
wartości.
Zbie\ność takiej metody jest mo\liwa tylko wówczas, gdy
dla ka\dego x takiego \e *x ą* # *x1 ą*, spełniony jest
warunek:
F(x)&F(ą) # x&ą
gdzie ą jest dokładną wartością pierwiastka, a jest liczbą z
x4 x3 x2 x1 x0 x
przedziału (0,1). Mo\na pokazać, \e:
x2&ą ' F(x1)&ą ' F(x1)&F(ą) # x1&ą
Andrzej Brozi, Instytut Fizyki Politechniki Aódzkiej
Metody Numeryczne w Fizyce 1 Minima i maksima funkcji jednej zmiennej 16
i dalej:
x3&ą ' F(x2)&F(ą) # x2&ą # 2 x1&ą
a zatem:
xi&ą # i&1 x1&ą
i xi dą\y do wartości pierwiastka ą (oczywiście gdy jest mniejsze od 1).
Tymczasem z pierwszego z zapisanych wzorów wynika, \e
F(x)&F(ą)
)
' F (x) # < 1
x&ą
czyli pochodna funkcji F dla ka\dego x le\ącego wewnątrz wybranego przedziału musi być mniejsza co do wartości
bezwzględnej od jedności.
Kolejne rysunki ilustrują przypadki gdy pochodna w otoczeniu szukanego pierwiastka nie spełnia tego warunku.
y=x
y=x
F(x)
F(x)
x0 x0 x
x4 x1 x0 x2 x
Co ciekawe (i bardzo po\yteczne) prawdziwe jest równie\ twierdzenie odwrotne:
Jeśli w pewnym przedziale pochodna funkcji F(x) jest mniejsza od jedności i ka\da wartość xi nale\y
do tego przedziału to w badanym przedziale istnieje pierwiastek.
Jako przykład rozpatrzymy równanie: Wykres funkcji badanej
4
3
f x ' x & 3x % 1 ' 0
Aatwo zauwa\yć, \e trzy pierwiastki tego równania
3
spełniają warunki: 2<ą1< 1, 0<ą2<1 i 1<ą3<2.
Rozwa\ymy kilka wersji algorytmu wykorzystują- 2
cych następujące przekształcenia:
1 3
1
3
F1 x ' x %1 ; F2 x ' x &2x%1
3
&1
3 0
F3 x ' ; F4 x ' 3x&1 -1.879 0.347 1.532
2
x &3
-1
Funkcje te mają następujące pochodne:
) )
2 2
F1 x ' x ; F2 x ' 3x &2
-2
-2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5
3
&
f = x3-3x+1
) 2x )
2
F3 x ' ; F4 x ' 3x &1
2
2
x &3
Andrzej Brozi, Instytut Fizyki Politechniki Aódzkiej
Metody Numeryczne w Fizyce 1 Minima i maksima funkcji jednej zmiennej 17
Pierwsza funkcja pomocnicza i jej pochodna Druga funkcja pomocnicza i jej pochodna
4 4
funkcja funkcja
pochodna
pochodna
3 3
2 2
1 1
0 0
-1.879 0.347 1.532 -1.879 0.347 1.532
-1 -1
-2 -2
-2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5
F1 = 1/3(x3+1), Fp1 = x2 F2 = x3-2x+1, Fp2 = 3x2-2
Trzecia funkcja pomocnicza i jej pochodna Czwarta funkcja pomocnicza i jej pochodna
4 4
funkcja
funkcja
pochodna
pochodna
3 3
2 2
1 1
0 0
-1.879 0.347 1.532 -1.879 0.347 1.532
-1 -1
-2 -2
-2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5
F3 = -1/(x2-3), Fp3 = 2x/(x2-3)2 F4 = (3x-1)(1/3), Fp4 = (3x-1)(-3/2)
Jak widać pochodna funkcji F1 : F1N(x ) = x 2 jest mniejsza od jedności dla wszystkich x z zakresu ( 1,1). Znaczy to,
\e funkcją tą mo\na posłu\yć się do znalezienia pierwiastka ą2, ale nie ą1 ani ą3. Próba wycelowania w pierwiastek ą1
lub ą3 spowoduje albo rozbiegnięcie się rozwiązania albo znalezienie pierwiastka ą2.
Funkcja F2 jest jeszcze gorsza w okolicy 0 wartość bezwzględna jej pochodnej równie\ przekracza 1. Oprócz
rozbiegania się (czyli nieograniczonego wzrostu kolejnych otrzymywanych wartości) mo\liwe jest równie\ zachowanie
chaotyczne, gdy w kolejnych cyklach otrzymujemy bardzo ró\niące się wartości.
Funkcja F3 będzie nadawać się do wyznaczenia pierwiastka ą2. Nawet wybranie wartości początkowej dokładnie równej
ą1 lub ą3 spowoduje otrzymanie wyniku ą2.
Zmodyfikowana pierwsza funkcja pomocnicza (z "k") i jej pochodna
4 Funkcja F4 jest odpowiednia do wyznaczenia ą1 i ą3.
funkcja
pochodna Odwrotnie ni\ poprzednio wybranie wartości początkowej
3
bliskiej ą2 doprowadzi nas do ą1 lub ą3 (zale\nie od tego
czy wybrana wartość jest większa czy mniejsza od ą2).
2
Funkcje F (x) mo\na jednak wybrać na jeszcze więcej
sposobów. Jeśli zauwa\ymy, \e gdy zachodzi ą = F (ą) to
1
musi zachodzić równie\ ą = (1 k )ą + k@ F (ą). Oka\e się
więc, \e równie dobrze mo\na rozwa\ać równanie
0
-1.879 0.347 1.532 rekurencyjne postaci:
xi%1 ' 1&k xi % k F xi
-1
gdzie F (x ) jest dowolną spośród mo\liwych funkcji, a
-2
stałą k mo\emy dobrać dowolnie.
-2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5
F1k = 3/2 x - (x3+1)/6, Fp1k = 3/2 - x2/2
Andrzej Brozi, Instytut Fizyki Politechniki Aódzkiej
Metody Numeryczne w Fizyce 1 Minima i maksima funkcji jednej zmiennej 18
W przypadku funkcji F1(x ) mo\emy zastosować k = i wzór rekurencyjny przyjmie postać:
3 1
xi%1 ' xi & xi3%1
2 6
zaś pochodna funkcji F1k(x ) będzie wynosić:
2
) 3 x
F1k x ' &
2 2
)
Zatem w przedziałach 1 < x < 5 zachodzi F x <1 , co umo\liwia obliczenie pierwiastków ą1 i ą3.
Poni\ej znajdują się wykresy ilustrujące przebieg obliczania pierwiastków ą1, ą2 i ą3 z wykorzystaniem poszczególnych
funkcji (F1, F2, F3 i F4).
Porównanie przedstawionych metod iteracyjnych
Na poni\szych wykresach, przedstawiających rozwiązania z wykorzystaniem poszczególnych funkcji kólkiem na osi
x zaznaczono punkt początkowy iteracji.
Zastosowanie pierwszej funkcji pomocniczej z obu stron dochodzimy do pierwiastka drugiego.
Ilustracja metody iteracyjnej - funkcja pierwsza: F = (x3+1)/3 Ilustracja metody iteracyjnej - funkcja pierwsza: F = (x3+1)/3
1 1
4 4
y=x
y=x
funkcja 1
funkcja 1
3 3
2 2
1 1
0 0
-1.80 1.50
-1 -1
-2 -2
-2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5
Znaleziono pierwiastek 2: 0.347 +- 2.1e-005. Bd faktyczny -2.8e-006. Wykonano 9 p Znaleziono pierwiastek 2: 0.347 +- 8.4e-005. Bd faktyczny 1.2e-005. Wykonano 11 p
tli. tli.
Funkcja druga nie uzyskujemy zbie\ności (prowadzenie obliczeń przez dłu\szy czas tylko zagmatwałoby rysunek).
Ilustracja metody iteracyjnej - funkcja druga: F = x3-2x+1 Ilustracja metody iteracyjnej - funkcja druga: F = x3-2x+1
2 2
4 4
y=x
y=x
funkcja 2
funkcja 2
3 3
2 2
1 1
0 0
0.50 1.20
-1 -1
-2 -2
-2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5
Program ROZBIEGA SI ! Ostatni wynik 0.739. Wykonano 15 p tli. Program ROZBIEG
SI ! Ostatni wynik 1.133. Wykonano 15 p tli.
Andrzej Brozi, Instytut Fizyki Politechniki Aódzkiej
Metody Numeryczne w Fizyce 1 Minima i maksima funkcji jednej zmiennej 19
Wykorzystanie funkcji trzeciej przedstawiłem tylko na jednym przykładzie. Obok porównanie wyników obliczeń
wykorzystujących poszczególne funkcje. Jak widać \adna z nich nie znalazła pierwiastka nr 3.
Ilustracja metody iteracyjnej - funkcja trzecia: F = -1/(x2-3)
3 Przebieg poszukiwania pierwiastka w metodzie iteracyjnej
4 2
y=x
funkcja 3
3
1
2
0
1
-1
-1.80
0
-2
-1
pierwiastek 1
pierwiastek 2
-2 -3
pierwiastek 3
funkcja 1
funkcja 2
-3
funkcja 3
-4
funkcja 4
funkcja 1k
-4
-5
-4 -3 -2 -1 0 1 2 0 5 10 15 20 25
Znaleziono pierwiastek 2: 0.347 +- 8.4e-005. Bd faktyczny -7.7e-006. Wykonano 7 p
tli.
Dwa przykłady obliczeń z wykorzystaniem funkcji czwartej funkcja uparcie ucieka od pierwiastka drugiego.
Ilustracja metody iteracyjnej - funkcja czwarta: F = (3x-1)(1/3) Ilustracja metody iteracyjnej - funkcja czwarta: F = (3x-1)(1/3)
4 4
4 4
y=x y=x
funkcja 4 funkcja 4
3 3
2 2
1 1
0 0
0.35 -0.50
-1 -1
-2 -2
-2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5
Znaleziono pierwiastek 3: 1.531 +- 8.8e-004. Bd faktyczny -6.5e-004. Wykonano 13 p tli. Znaleziono pierwiastek 1: -1.879 +- 7.7e-004. Bd faktyczny 3.0e-004. Wykonano 8 p
tli.
Ilustracja metody iteracyjnej - zmodyfikowana funkcja pierwsza: F
k = 3/2 x - (x3+1)/6
1
4
y=x
Rysunek obok przedstawia wykorzystanie
funkcja 1k
zmodyfikowanej funkcji pierwszej (F1k). W podpisie
3
zaznaczono, \e program wykonał 9 pętli trudno je
zauwa\yć na wykresie.
2
1
-1.80
0
-1
-2
-2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5
Znaleziono pierwiastek 1: -1.879 +- 7.1e-006. B tli.
d faktyczny 1.5e-006. Wykonano 9 p
Andrzej Brozi, Instytut Fizyki Politechniki Aódzkiej
Metody Numeryczne w Fizyce 1 Minima i maksima funkcji jednej zmiennej 20
Rząd metody iteracyjnej
Załó\my, \e badamy rozwiązanie iteracyjne postaci xn+1 = F (xn). Jeśli ą jest dokładną wartością pierwiastka równania
f (x ) = 0 to oczywiście musi zachodzić ą = F (ą). Oznaczmy błąd bezwzględny n-tego przybli\enia przez "n: "n = xn ą.
Jeśli tylko funkcja F ma wystarczającą liczbę pochodnych będziemy mogli zapisać:
2
2
xn&ą
dF d F
xn%1 ' F xn ' F ą % xn&ą % % '
2
dx 2
dx
"2 d 2F
n
' F ą % "n dF % %
2
dx 2
dx
Jeśli od obu stron tak zapisanego równania odejmiemy ą = F(ą) otrzymamy:
1
) ))
"n%1 ' "nF % "2F %
n
2
praktycznie mo\emy zało\yć, \e błąd przybli\enia n+1 jest równy pierwszemu niezerowemu składnikowi ostatniego
równania. Jeśli pierwszy taki składnik zawiera ju\ pierwszą pochodną to mówimy, \e metoda jest pierwszego rzędu, jeśli
drugą to mamy do czynienia z metodą rzędu drugiego itd..
W przypadku metody rzędu drugiego jeśli błąd w danym kroku jest rzędu 10 k to po kolejnej iteracji będzie rzędu 10 2k
(razy jakaś stała wartość drugiej pochodnej funkcji F w punkcie ą), zatem ka\da iteracja (w metodzie rzędu drugiego)
podwaja liczbę cyfr dokładnych rozwiązania.
Znajdowanie miejsc zerowych funkcji wielu zmiennych
Jeśli mamy do czynienia z pojedynczą funkcją więcej ni\ jednej zmiennej to oczywiście mamy zbyt mało danych, aby
znalezć jej miejsce zerowe. Jeśli natomiast dany jest układ równań odpowiednio wysokiego rzędu:
f x ' 0
gdzie: f ' f1,f2,,fn i x ' x1,x2,,xn , to jest on równowa\ny układowi n równań z n zmiennymi (ka\dy z wierszy
takiego układu mo\e być jedną funkcją n zmiennych). Zagadnienie to bedzie więc równowa\ne rozwiązaniu układu równań
nieliniowych zajmiemy się tym w następnym wykładzie.
Znajdowanie minimów funkcji wielu zmiennych
Ta kwestia bedzie omówiona w następnym wykładzie.
Andrzej Brozi, Instytut Fizyki Politechniki Aódzkiej
Metody Numeryczne w Fizyce 1 Minima i maksima funkcji jednej zmiennej 21
Rozwiązanie analityczne rzutu z wahadła
L
Na zakończenie powróćmy do rozpoczynającego ten
H
rozdział zadania ze skokiem z wahadła.
Ogólne rozwiązanie równania dla rzutu ukośnego:
x ' x0%v0t cosą
2
gt
y ' y0%v0t siną&
2
Zale\ności x0, y0 i v0 od kąta ą:
x0 ą ' L siną
y0 ą ' H &L cosą
v0 ą ' 2gL cosą &cosą0
Wyrugowanie czasu i równanie toru:
x&x0
t '
v0cosą
2
x&x0
g
y ' y0% x&x0 tgą& ' 0
2
2
v0 cos2ą
Podstawienie z zamiast (x x0), równanie kwadratowe i delta równania kwadratowego:
z 6 x&x0
2
g z
y0 % z tgą & ' 0
2
2
v0 cos2ą
2gy0
" ' tg2ą %
2
v0 cos2ą
Rozwiązanie równania kwadratowego względem z:
2gy0
&tgą & tg2ą %
2
v0 cos2ą
z '
g
&
2
v0 cos2ą
2
v0 2gy0
z ' cosą siną % sin2ą%
2
g
v0
Rozwiązanie względem x:
2
v0 2gy0
x ' x0 % cosą sin2ą% % siną
2
g
v0
H
& cosą
L
x ' Lsiną % 2L cosą&cosą0 cosą siną% sin2ą%
cosą&cosą0
Andrzej Brozi, Instytut Fizyki Politechniki Aódzkiej
Metody Numeryczne w Fizyce 1 Minima i maksima funkcji jednej zmiennej 22
Poni\ej znajduje się wykres trajektorii lotu po skoku z wahadła dla kilku kątów wyskoku (do obliczeń przyjęto h = 2 m,
L = 1 m, ąo = 45E). Obok zale\ność zasięgu od kąta wyskoku.
Trajektorie skoku z wahada ZaleónoŃł zasi gu od kta wyskoku
1.4
15
20
1.45
25
1.2
30
35
1 40
1.425
0.8
0.6
1.4
0.4
0.2
1.375
0
-0.2
1.35
0.2 0.4 0.6 0.8 1 1.2 1.4 1.6
x 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
kt ą [w stopniach]
Andrzej Brozi, Instytut Fizyki Politechniki Aódzkiej
Wyszukiwarka
Podobne podstrony:
04a
04a
04a
04a
MNF 05
04a
04a
04a
04a?5 Power Supply and Bus Systems
Wykład 04a Bayes
04a
04a
04a
04a
04a
04a
MNF 04b
więcej podobnych podstron