I. MINIMALIZACJA FUNKCJI.
4.1 Ekstrema funkcji. Punkty siodłowe.
Funkcja f (x) jednej zmiennej ma w punkcie x
0
minimum (maksimum) lokalne, jeśli istnieje
takie otoczenie punktu x
0
, Ŝe dla dowolnego punktu x z tego otoczenia zachodzi f(x) > f(x
0
)
(f(x) < f(x
0
) ). Minima i maksima lokalne to ogólnie ekstrema lokalne.
Minimum (maksimum) globalnym danej funkcji nazywamy najmniejszą (największą) wartość
przyjmowaną przez funkcję w całej dziedzinie.
Twierdzenie Fermata mówi, Ŝe w punkcie x
0
, w którym funkcja ma ekstremum lokalne znika
jej pierwsza pochodna: f’(x
0
) = 0. ZauwaŜmy, Ŝe zerowanie się pierwszej pochodnej to
warunek konieczny, ale niewystarczający istnienia ekstremum lokalnego.
Sprawdzenie istnienia bądź nie ekstremum oraz jego typu (minimum czy maksimum) funkcji
jednej zmiennej w danym punkcie moŜna oprzeć na badaniu zmian znaku pierwszej
pochodnej w otoczeniu tego punktu.
Skupmy się jednak na teście drugiej pochodnej:
Jeśli w punkcie x
0
pierwsza pochodna funkcji przyjmuje wartość zero f’(x
0
) = 0 a druga
pochodna jest róŜna od zera f’’(x
0
)
≠
0, to
jeśli f’’(x
0
) > 0
funkcja ma w x
0
minimum lokalne
jeśli f’’(x
0
) < 0
funkcja ma w x
0
maksimum lokalne
Jeśli natomiast f’’(x
0
) = 0, to musimy zbadać wyŜsze pochodne.
Uzasadnienie powyŜszego twierdzenia moŜna oprzeć na wzorze Taylora.
K
+
+
+
=
+
)
(
''
2
1
)
(
'
)
(
)
(
0
2
0
0
0
x
f
h
x
hf
x
f
h
x
f
Jeśli h jest dostatecznie małe, to człony z pochodnymi wyŜszymi niŜ druga moŜemy
zaniedbać i mamy (uwzględniając, Ŝe f’(x
0
) = 0):
)
(
''
2
1
)
(
)
(
0
2
0
0
x
f
h
x
f
h
x
f
+
≈
+
Gdy f’’(x
0
) > 0, to dla dowolnego dostatecznie małego h mamy f(x
0
+ h) > f(x
0
) a to oznacza
istnienie w x
0
minimum lokalnego. Analogicznie z faktu f’’(x
0
) < 0 wynika f(x
0
+ h) < f(x
0
),
czyli maksimum lokalne w x
0
. Gdy f’’(x
0
) = 0, musimy sięgnąć do dalszych członów we
wzorze Taylora.
Analogicznie do ekstremów funkcji jednej zmiennej określa się ekstrema funkcji n zmiennych
(zaleŜnej od wektora x = (x
1
, x
2
, ..., x
n
)).
Gdy f(x) > f(x
0
) dla kaŜdego x w pewnym otoczeniu x
0
, to funkcja ma w x
0
minimum lokalne.
Gdy f(x) < f(x
0
) dla kaŜdego x w pewnym otoczeniu x
0
, to funkcja ma w x
0
maksimum
lokalne.
Podobnie jak w przypadku jednej zmiennej, warunkiem koniecznym istnienia w punkcie x
0
ekstremum lokalnego jest zerowanie się w nim wszystkich pierwszych pochodnych
cząstkowych:
0
2
1
=
∂
∂
=
=
∂
∂
=
∂
∂
n
x
f
x
f
x
f
L
, wszystkie pochodne liczone w punkcie x
0
Wektor o składowych równych pochodnym cząstkowym funkcji nazywamy gradientem
funkcji a operację tworzenia gradientu zaznaczamy przez operator
∇
∂
∂
∂
∂
=
∇
n
x
f
x
f
f
M
1
Twierdzenie Fermata dla funkcji wielu zmiennych mówi więc, Ŝe w punkcie ekstremum
lokalnego gradient funkcji jest równy zero.
Gradient funkcji w danym punkcie określa kierunek najszybszego wzrostu funkcji.
Punkty, w których gradient się zeruje, to punkty krytyczne funkcji.
Macierz drugich pochodnych cząstkowych funkcji wielu zmiennych nazywamy hesjanem
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
=
n
n
n
n
n
n
x
x
f
x
x
f
x
x
f
x
x
f
x
x
f
x
x
f
x
x
f
x
x
f
x
x
f
2
2
2
1
2
2
2
2
2
2
1
2
2
1
2
2
1
2
1
1
2
)
(
L
M
M
M
M
L
L
x
H
,
pochodne są tu liczone w punkcie x.
Z wykorzystaniem tych wielkości moŜemy zapisać wzór Taylora do członu z drugą pochodną
f(x
0
+ h) = f(x
0
) + h
T
·
f
∇
(x
0
)+
2
1
h
T
H(x
0
)h + ...
PoniewaŜ w punkcie krytycznym
f
∇
(x
0
)= 0, więc drugi człon po prawej stronie wzoru znika
i charakter punktu krytycznego zaleŜy od znaku formy kwadratowej h
T
H(x
0
)h. Jeśli macierz
H(x
0
) jest dodatnio określona, to dla dowolnego dostatecznie małego h mamy f(x
0
+ h) >
f(x
0
), czyli w x
0
jest minimum lokalne funkcji
f. Analogicznie, jeśli H(x
0
) jest określona
ujemnie, to w x
0
mamy maksimum lokalne
f.
Określoność macierzy H moŜemy zbadać sprowadzając ją do postaci diagonalnej. Istotnie,
przy diagonalnej macierzy mamy
h
T
Hh =
∑
n
i
i
i
h
2
'
)
(
λ
,
gdzie
λ
i
to wartości własne H. Znakiem ’ zaznaczyliśmy, Ŝe po diagonalizacji macierzy
mamy nowe składowe h określone jako kombinacje starych wartości przy pomocy macierzy
wektorów własnych. PoniewaŜ przy badaniu ekstremum h ma być dowolne, ta transformacja
nie ma znaczenia dla dalszych rozwaŜań.
Widzimy zatem, Ŝe jeśli wszystkie wartości własne macierzy hesjanu są dodatnie, to macierz
ta jest dodatnio określona i f ma w x
0
minimum lokalne. Jeśli wszystkie wartości własne
hesjanu są ujemne, to jest on ujemnie określony i w x
0
jest maksimum lokalne. Jeśli część
wartości własnych hesjanu jest dodatnia a część ujemna, to w x
0
jest punkt siodłowy.
W szczególności, jeśli jedna wartość własna jest ujemna a pozostałe są dodatnie, to w x
0
jest
punkt siodłowy pierwszego rzędu.
Dla funkcji dwu zmiennych jedna wartość dodatnia i jedna ujemna hesjanu oznacza istnienie
siodła. MoŜna zauwaŜyć, Ŝe w zaleŜności od wyboru wektora h przy zbliŜaniu się do punktu
siodłowego w jednym kierunku funkcja rośnie (podchodzimy do przełęczy po zboczu) a w
innym maleje (schodzimy ku przełęczy po grani).
Znajdowanie minimów lub maksimów funkcji jest przedmiotem optymalizacji. W fizyce
zwykle interesują nas minima ze względu na dąŜenie przez układ do minimum energii
potencjalnej i róŜne zasady oparte na minimalizacji pewnej wielkości.
PoniewaŜ znajdowanie maksimów funkcji f jest równowaŜne znajdowaniu minimów funkcji
–f, zajmujemy się jedynie sposobami wyznaczania minimum funkcji.
4.2 Znajdowanie minimum funkcji jednej zmiennej.
ZałóŜmy, Ŝe mamy dane trzy punkty a, b i c uporządkowane tak, Ŝe
a < c < b
(czyli c znajduje się wewnątrz przedziału (a, b)).
Jeśli
f(a) > f(c),
f(c) < f(b),
to w przedziale (a, b) znajduje się punkt, w którym funkcja f ma minimum lokalne. MoŜemy
teraz obliczyć wartość funkcji f(d) w kolejnym punkcie d wewnątrz przedziału i po
porównaniu wartości wybrać nowy, krótszy przedział (w zaleŜności od wartości funkcji
odrzucając punkt a lub b) zawierający minimum. Powtarzając tę procedurę moŜemy coraz
bardziej zawęŜać przedział znajdując coraz lepsze ograniczenie dla minimum.
Metoda ta jest podobna to znajdowaniu miejsca zerowego funkcji metodą połowienia
przedziału, z tą róŜnicą, Ŝe do określenia przedziału zawierającego miejsce zerowe wystarcza
znajomość funkcji w dwu punktach, do wyznaczenia przedziału z minimum musimy zaś
wyznaczyć f w trzech punktach.
Jeśli wstępnie wyznaczymy przedział, w którym znajduje się tylko jedno minimum, to dzieląc
przedział na mniejsze będziemy w stanie znaleźć połoŜenie tego ekstremum.
Stosowana zwykle metoda wyboru kolejnego punktu w przedziale opiera się na złotym
podziale odcinka. Przy złotym podziale stosunek części, na które podzielono odcinek jest taki
sam jak stosunek dłuŜszej części do całości. Rozwiązując odpowiednią proporcję łatwo
stwierdzić, Ŝe stosunek ten jest równy
2
1
5
−
=
ϕ
Stosując metodę złotego podziału odcinka w pierwszym kroku wyznaczamy w przedziale
(a, b) dwa punkty c, d takie, Ŝe
c
= a +
ϕ
(b – a)
d
= b –
ϕ
(b – a)
Po obliczeniu wartości funkcji odrzucamy ten z końców odcinka (a, b), w którym funkcja
przyjmuje większą wartość i znajdujemy punkt dzielący w złotym stosunku dłuŜszy z
pozostałych dwu przedziałów. Kroki te następnie powtarzamy.
Metodą tą moŜemy takŜe znaleźć punkt, w którym na zadanej prostej funkcja wielu
zmiennych przyjmuje wartość minimalną.
Jeśli znajdujemy się w punkcie x
0
dostatecznie blisko minimum funkcji, moŜemy spróbować
przybliŜyć funkcję formą kwadratową
f
(x) = f(x
0
) – b(x – x
0
) +
2
1
A (x – x
0
)
2
,
gdzie b = –f’(x
0
) a A = f’’(x
0
); zaleŜność ta wynika oczywiście z wzoru Taylora.
RóŜniczkując tę zaleŜność stronami dostajemy
)
(
A
)
(
0
x
x
b
dx
x
df
−
+
−
=
W punkcie x
m
, w którym funkcja ma minimum jej pierwsza pochodna się zeruje, czyli
–b + A(x
m
– x
0
) = 0
skąd znajdujemy przybliŜenie połoŜenia minimum
A
0
b
x
x
m
+
=
W znalezionym punkcie moŜemy operacje powtórzyć dostając kolejne, lepsze przybliŜenie.
4.3 Minima funkcji wielu zmiennych – metody gradientowe.
Wiele metod wyznaczania połoŜenia minimum funkcji wielu zmiennych wykorzystuje
informacje o gradiencie badanej funkcji.
Po pierwsze, zauwaŜmy, Ŝe gradient
f
∇
w danym punkcie wyznacza kierunek najszybszego
wzrostu funkcji. Oznacza to, Ŝe w kierunku przeciwnym
f
∇
−
mamy do czynienia z
najszybszym spadkiem. MoŜemy zatem przesunąć się w tym kierunku o odległość zaleŜną od
wartości gradientu, albo poszukać punktu na prostej wyznaczonej przez gradient, w którym
funkcja ma najmniejszą wartość. W ten sposób znajdziemy się w nowym punkcie, w którym
znowu wyznaczamy gradient funkcji i podąŜamy w kierunku przeciwnym do niego. PoniewaŜ
z kaŜdego punktu przesuwamy się w kierunku
f
∇
−
, czyli największego spadku funkcji
metoda ta nazywa się metodą najszybszego spadku.
Jak się okazuje, metoda najszybszego spadku powoduje zygzakowanie przy zbliŜaniu się do
minimum. Związane jest to z faktem, iŜ przy wyborze nowego kierunku zatracamy informację
zdobytą wcześniej.
Metodą o lepszej zbieŜności do połoŜenia minimum jest metoda kierunków sprzęŜonych o
gradientu funkcji.
Kierunek p sprzęŜony do danego kierunku q względem macierzy A zadany jest warunkiem
p
T
Aq = 0
Przy minimalizacji metodą kierunków sprzęŜonych do gradientu wyznacza się kierunki
sprzęŜone względem macierzy hesjanu, posuwając się kolejno wzdłuŜ kierunku sprzęŜonego
do gradientu w danym punkcie a potem kierunków sprzęŜonych do niego.
Kolejna grupa metod opiera się na przybliŜeniu funkcji formą kwadratową i poszukiwaniu
rozwiązania sposobem podobnym do opisanego w 4.2.
f
(x) = f(x
0
) + (x – x
0
)
f
∇
(x
0
)+
2
1
(x – x
0
)H(x
0
) (x –x
0
)
RóŜniczkując stronami i przyrównując pochodną do zera znajdziemy
0 = f
∇
(x
0
)+ H (x
m
– x
0
)
a stąd
x
m
= x
0
– H
-1
· f
∇
(x
0
)
Istnieją metody, w których nie oblicza się hesjanu, natomiast przybliŜony hesjan tworzy się
iteracyjnie. Są to np. metody Davidona-Fletchera-Powella lub Broydena-Fletchera-Shanno.
ZauwaŜmy jeszcze, Ŝe gradient funkcji zeruje się nie tylko w punkcie, w którym funkcja
osiąga minimum, ale teŜ w punktach siodłowych i w maksimach. Zatem jeśli metoda bazująca
na gradientach przypadkiem trafi w któryś z tych punktów, to w nim pozostanie.
4.4 Metoda sympleksów
Jako przykład metody, która nie wykorzystuje informacji o gradiencie a jedynie obliczanie
wartości funkcji w róŜnych punktach, naszkicujemy metodę sympleksów (Nelder-Meada).
Sympleksem w przestrzeni n wymiarowej nazywamy układ n + 1 punktów (wierzchołków
sympleksu). Na płaszczyźnie sympleksami są trójkąty.
W metodzie sympleksów oblicza się wartości minimalizowanej funkcji w wierzchołkach
sympleksu a następnie wykonuje na sympleksie jedną z operacji (pokazanych na przykładzie
sympleksów w przestrzeni dwuwymiarowej, stary sympleks rysowany jest linią pogrubioną):
odbicie
rozciągnięcie
spłaszczenie
kontrakcja
W przypadku pierwszych trzech operacji (odbicie, rozciągnięcie, spłaszczenie) zmienia się
jedynie połoŜenie tego wierzchołka sympleksu, w którym funkcja przyjmowała największą
wartość (przemieszcza się on wzdłuŜ prostej wyznaczonej przezeń oraz środek cięŜkości
pozostałych wierzchołków). Przy kontrakcji przemieszczają się wszystkie punkty z wyjątkiem
tego, w którym funkcja przyjmuje wartość najmniejszą.
Algorytm Nelder-Meada określa (w zaleŜności od wartości funkcji w wierzchołkach
sympleksu), którą operację i z jakim współczynnikiem naleŜy wykonać.
Procedurę przerywa się kiedy róŜnica najmniejszej i największej wartości funkcji w
wierzchołkach sympleksu stanie się dostatecznie mała i połoŜenie wierzchołka z najmniejszą
wartością określa lokalizację minimum funkcji.
4.5 Minimalizacja funkcji a chemia kwantowa.
Minimalizacja odgrywa istotną rolę w badaniu struktury cząsteczek metodami chemii
kwantowej. Znajdując optymalną geometrię cząsteczki minimalizujemy energię potencjalną
(jako funkcję połoŜenia atomów) obliczaną którąś z metod kwantowochemicznych.
Wykorzystuje się w tym celu metody gradientowe, poszukując punktu stacjonarnego, tzn.
takiego, którym gradient energii potencjalnej jest równy zero. Pomocna jest w tym dostępność
w uŜywanym programie analitycznych wyraŜeń na gradienty energii dla danej metody
kwantowochemicznej. Jeśli ich nie ma, gradienty muszą być obliczane numerycznie, co
wydłuŜa czas obliczeń przy jednoczesnym zmniejszeniu dokładności.
Zwróćmy jednak uwagę, Ŝe metoda gradientowa moŜe przypadkiem doprowadzić nas nie do
minimum, ale do punktu siodłowego. Dlatego chcąc się upewnić, Ŝe znaleziona geometria
odpowiada minimum energii powinniśmy zbadać wartości własne macierzy hesjanu w
badanym punkcie (przydatna jest tu dostępność analitycznych wyraŜeń na hesjan, by uniknąć
numerycznego róŜniczkowania gradientu). W praktyce wyznacza się wartości własne hesjanu
waŜonego przez masy atomów, co daje nam częstości drgań normalnych cząsteczki
(wyznaczone wartości własne są kwadratami częstości drgań normalnych). Chcąc zatem
potwierdzić charakter znalezionego punktu stacjonarnego, wyznaczamy drgania normalne
cząsteczki. Jeśli wszystkie wartości własne hesjanu (waŜonego masami) są dodatnie, to
znaleźliśmy geometrię odpowiadającą minimum energii potencjalnej (a pierwiastki z tych
wartości to częstości drgań). Jeśli oprócz dodatnich wartości własnych występują takŜe
ujemne (co formalnie odpowiada urojonym częstościom drgań), to znaleziono geometrię
punktu siodłowego.
Standardowe metody minimalizacji prowadzą do punktu siodłowego jedynie przypadkiem.
Istnieją jednak przypadki, kiedy zaleŜy nam nie na wyznaczeniu geometrii odpowiadającej
minimum energii potencjalnej a właśnie na punkcie siodłowym. OtóŜ punkt siodłowy
pierwszego rzędu odpowiada stanowi przejściowemu reakcji – na mapie energii potencjalnej
dolina produktów jest oddzielona przełęczą (punktem siodłowym) od doliny substratów.
Badając zatem przebieg reakcji wzdłuŜ pewnej współrzędnej poszukujemy punktu
siodłowego. Do takich potrzeb stworzono specjalne algorytmy optymalizacji przystosowane
do poszukiwania punktów siodłowych.
ZauwaŜmy jeszcze, Ŝe większość algorytmów minimalizacji prowadzi nas do wyznaczenia
minimum lokalnego. Tymczasem chemika interesują nie tylko róŜne konformacje cząsteczki
odpowiadające minimom lokalnym, ale teŜ znalezienie konformacji najniŜej energetycznej –
minimum globalnego. Liczba minimów lokalnych energii potencjalnej dla interesujących
cząsteczek moŜe łatwo przekraczać dziesiątki i setki tysięcy, znajdowanie minimum
globalnego jest zatem bardzo trudnym i wciąŜ intensywnie badanym zadaniem. Pomocne są w
tym algorytmy takie, jak symulowane wyŜarzanie czy algorytmy genetyczne. Ich wspólną
cechą jest to, iŜ przy poszukiwaniu minimum z pewnym prawdopodobieństwem akceptują
przejściowo gorsze rozwiązania (posuwają się w stronę większych wartości funkcji) – czasem
bowiem jedyna droga z jednej kotliny do innej, głębszej, wiedzie przez przełęcz, czyli pod
górę.
II. RÓWNANIA RÓśNICZKOWE ZWYCZAJNE.
4.6 Terminologia. Istnienie rozwiązań.
Równanie róŜniczkowe zawiera pochodne nieznanej funkcji. Jeśli szukana funkcja zaleŜy
tylko od jednej zmiennej (w równaniu nie występują pochodne cząstkowe), to równanie
nazywamy równaniem róŜniczkowym zwyczajnym.
Rzędem równania róŜniczkowego nazywamy rząd najwyŜszej występującej w nim pochodnej.
Stopniem równania róŜniczkowego to potęga, w jakiej występuje najwyŜsza pochodna
szukanej funkcji, jeśli równanie da się zapisać w postaci wielomianu od występujących w nim
pochodnych. Np. równanie róŜniczkowe 1-szego rzędu i 1-szego stopnia da się zapisać jako
)
,
(
y
x
f
dx
dy
=
Funkcja y(x) jest rozwiązaniem równania róŜniczkowego, jeśli w pewnym przedziale
zmiennej x spełnia to równanie toŜsamościowo po podstawieniu doń funkcji y i jej
pochodnych.
Podręczniki matematyki traktujące o równaniach róŜniczkowych zajmują się znajdowaniem
funkcji będących rozwiązaniami dla róŜnych typów równań. W praktyce numerycznej
musimy zadowolić się wyznaczeniem tabeli wartości poszukiwanej funkcji y dla róŜnych
wartości x (później moŜemy ewentualnie uŜyć tej tabeli w interpolacji).
Teoria równań róŜniczkowych pozwala nam stwierdzić, Ŝe ze względu na stałe całkowania
dla znalezienia jednoznacznego rozwiązania musimy podać jeszcze dodatkowe informacje o
szukanej funkcji y, np. wartość funkcji i jej pochodnych (w zaleŜności od rzędu równania) w
zadanym punkcie: y(x
0
), y’(x
0
) .... Są to tzw. warunki początkowe. Znalezienie rozwiązania
równania róŜniczkowego zwyczajnego z zadanymi warunkami początkowymi w metodach
numerycznych zwane jest rozwiązaniem zagadnienia początkowego. Będziemy się nim
zajmować i dalszej części tego wykładu. (Alternatywnie zamiast wartości funkcji i jej
pochodnych w jednym punkcie moŜna zadać wartości funkcji w róŜnych punktach – są to
warunki brzegowe.)
O jednoznaczności rozwiązania zagadnienia początkowego mówi nam twierdzenie.
Jeśli f(x,y) i
y
f
∂
∂
/
są funkcjami rzeczywistymi, ograniczonymi, jednowartościowymi i
ciągłymi w pewnym otoczeniu punktu (x
0
, y
0
), to w przedziale [x
0
– h, x
0
+ h] wewnątrz tego
otoczenia istnieje dokładnie jedno rozwiązanie równania pierwszego rzędu
)
,
(
y
x
f
dx
dy
=
z warunkiem początkowym y(x
0
) = y
0
.
Szukając numerycznie rozwiązania zagadnienia początkowego, czyli wartości funkcji y w
interesującym nas punkcie x wyznaczamy krok po kroku wartości funkcji y w punktach
pośrednich między x
0
(warunki początkowe) a x. W kaŜdym takim kroku popełniamy pewien
błąd wynikający z przyjętych przybliŜeń oraz zaokrągleń. Jest to błąd lokalny. Błędy lokalne
kumulują się w błąd globalny – róŜnicę między obliczoną przez nas wartością szukanej
funkcji w interesującym nas punkcie a dokładnym rozwiązaniem równania w tym punkcie.
Oczywiście zaleŜny nam na metodach dających jak najmniejsze błędy.
4.7 Układy równań róŜniczkowych rzędu pierwszego. Równania wyŜszych rzędów.
W zastosowaniach praktycznych w naukach ścisłych często pojawiają się układy równań
róŜniczkowych pierwszego rzędu. Układ taki ma postać
)
,
,
,
,
(
2
1
'
n
i
i
y
y
y
x
f
y
K
=
i
= 1, ..., n
Niewiadomymi są tu funkcje y
i
zaleŜne od argumentu x. Dodatkowo potrzebne są warunki
początkowe y
i
(x
0
) określające wartości poszczególnych funkcji w punkcie x
0
.
Układ taki moŜna wyrazić w notacji wektorowej. Niech y będzie wektorem kolumnowym o
składowych y
i
a f wektorem kolumnowym o składowych f
i
. Wtedy układ da się zapisać jako
y’ = f(x, y)
y(x
0
) = y
0
Przy takim zapisie znaczna część metod rozwiązywania jednego równania pierwszego rzędu
przenosi się na rozwiązywanie układu, zatem dalej będziemy zajmować się rozwiązywaniem
pojedynczego równania.
Równanie róŜniczkowe rzędu wyŜszego niŜ pierwszy
y
(n)
= f(x, y, y’, y’’, ... y
(n-1)
)
moŜna przekształcić na układ równań rzędu pierwszego wprowadzając nowe zmienne równe
kolejnym pochodnym funkcji y.
Oznaczając bowiem:
)
1
(
'
1
'
1
2
1
'
−
−
=
=
=
=
=
n
n
n
y
y
y
η
η
η
η
η
L
dostajemy z wyjściowego równania układ równań rzędu pierwszego:
)
,
,
,
,
(
2
1
'
'
1
3
'
2
2
'
1
n
n
n
n
x
f
η
η
η
η
η
η
η
η
η
η
K
L
=
=
=
=
−
Podobnie moŜna przekształcić układ równań rzędu wyŜszego niŜ 1.
PoniewaŜ równania róŜniczkowe rzędu wyŜszego niŜ 1 da się sprowadzić do układu równań
rzędu pierwszego, ograniczymy się do rozwiązywania równań rzędu 1.
Czasem powyŜszy układ równań zapisuje się w postaci, w której zmienna x nie występuje
jawnie – jest to układ autonomiczny. Osiągamy to zastępując zmienną x dodatkową zmienną
0
η
i dokładając do układu równanie
1
'
0
=
η
.
4.8 Zastosowanie wzoru Taylora. Metoda Eulera.
Rozwiązanie równania
y’ = f(x,y)
moŜemy oprzeć na wzorze Taylora
K
+
+
+
+
=
+
)
(
!
3
1
)
(
''
2
1
)
(
'
)
(
)
(
)
3
(
3
2
x
y
h
x
y
h
x
hy
x
y
h
x
y
Zatrzymujemy składniki do pewnego rzędu. Potrzebne pochodne funkcji y obliczamy
róŜniczkując względem x obie strony wyjściowego równania:
y’’ = f’(x,y)
y
(3)
= f’’(x,y)
...
Mając warunek początkowy y(x
0
) obliczamy korzystając z powyŜszego wzoru wartość funkcji
y(x + h) w punkcie x + h, następnie wartość funkcji w x + 2h, itd.
Rząd metody określa rząd zachowanej w rozwinięciu pochodnej.
Najprostszą metodą jest metoda rzędu pierwszego, w której zachowaliśmy tylko pierwszą
pochodną. Jest to metoda Eulera, opisana wzorem
y(x + h)
≈
y(x) + hf(x, y)
Unikamy w niej obliczania pochodnej funkcji f(x,y), ale musimy stosować mały krok h.
4.9 Metody Rungego-Kutty.
Metody Rungego-Kutty pozwalają na podstawie znajomości szukanej funkcji w pewnym
punkcie y
n
= y(x
n
) znaleźć jej przybliŜenie w punkcie odległym o h: y
n+1
= y(x
n
+ h). W tym
celu obliczamy wartości funkcji f w specjalnie dobranych punktach i tworzymy kombinacje
tych wartości, aby uzyskać jak najlepsze przybliŜenie funkcji w x
n+1
.
Metoda Rungego-Kutty rzędu drugiego (zwana metodą Heuna) opisana jest wzorami
)
(
)
,
(
)
,
(
2
1
2
1
1
1
2
1
F
F
y
y
F
y
h
x
hf
F
y
x
hf
F
n
n
n
n
n
n
+
+
=
+
+
=
=
+
Najbardziej znana jest metoda Rungego-Kutty rzędu czwartego, której wzory są bardziej
skomplikowane, ale łatwe do zaprogramowania:
)
2
2
(
)
,
(
)
,
(
)
,
(
)
,
(
4
3
2
1
6
1
1
3
4
2
2
1
2
1
3
1
2
1
2
1
2
1
F
F
F
F
y
y
F
y
h
x
hf
F
F
y
h
x
hf
F
F
y
h
x
hf
F
y
x
hf
F
n
n
n
n
n
n
n
n
n
n
+
+
+
+
=
+
+
=
+
+
=
+
+
=
=
+
Podobnie jak w metodzie Eulera, nie musimy obliczać wartości pochodnych funkcji f.
Metoda Eulera jak i metody Rungego-Kutty pozwalają na wykorzystanie ekstrapolacji
Richardsona do poprawienia dokładności wyniku.
4.10 Metody wielokrokowe. Metody ekstrapolacyjno-interpolacyjne.
Metody Eulera i Rungego-Kutty są jednokrokowe: aby znaleźć rozwiązanie y(x + h) w
punkcie x + h musimy znać rozwiązanie y(x) tylko w jednym punkcie. W metodach
wielokrokowych uŜywamy do tego wartości w kilku punktach. Ogólna metoda
wyprowadzania takich rozwiązań jest następująca.
Szukamy rozwiązania zagadnienia początkowego
y’ = f(x, y),
y(x
0
) = y
0
w punktach x
1
, x
2
, ... . Przez scałkowanie obu stron powyŜszego równania róŜniczkowego
otrzymujemy
∫
+
+
=
+
1
))
(
,
(
)
(
)
(
1
n
n
x
x
n
n
dx
x
y
x
f
x
y
x
y
Do obliczenia całki moŜemy wykorzystać kwadraturę z węzłami x
i
.
W ten sposób moŜemy otrzymać wzory Adamsa-Bashforta
y
n+1
= y
n
+ af
n
+ bf
n-1
+ cf
n-2
+ ...
(*)
Oznaczyliśmy f
i
= f(x
i
,y
i
), gdzie y
i
jest przybliŜoną wartością rozwiązania w punkcie x
i
.
Na przykład, przy równo rozmieszczonych punktach x
i
= x
0
+ ih, uŜywając pięciu węzłów dla
wyznaczenia współczynników kwadratury otrzymać moŜemy wzór
)
251
1274
2616
2774
1901
(
4
3
2
1
720
1
1
−
−
−
−
+
+
−
+
−
+
=
n
n
n
n
n
n
n
f
f
f
f
f
h
y
y
Wyprowadzając wzory kwadratury w metodzie wielokrokowej moŜemy teŜ załoŜyć, Ŝe
wartość całki po prawej stronie zaleŜy teŜ od y
n+1
. Dostajemy wtedy wzory Adamsa-Moultona
y
n+1
= y
n
+ af
n+1
+ bf
n
+ cf
n-1
+ ...
(**)
np. dla punktów rozmieszczonych jak w przykładzie wyŜej, wzór
).
19
106
264
646
251
(
3
2
1
1
720
1
1
−
−
−
+
+
−
+
−
+
+
=
n
n
n
n
n
n
n
f
f
f
f
f
h
y
y
ZauwaŜmy, Ŝe szukane y
n+1
występuje takŜe po prawej stronie powyŜszych wzorów – musimy
je uŜyć do obliczenia f
n+1
. Jest to przykład metody niejawnej.
Równanie (**) moŜna rozwiązywać iteracyjnie podstawiając uzyskane y
n+1
do prawej strony
wzoru. Wstępne przybliŜenie moŜna otrzymać uŜywając wzoru Adamsa-Bashforta (*) do
obliczenia przybliŜonej wartości y
n+1
(wzór wstępny) a następnie podstawiając ją do prawej
strony wzoru Adamsa-Moultona (**) (wzór korekcyjny). Jest to metoda ekstrapolacyjno-
interpolacyjna (ang.: predictor-corrector). Iterację przerywa się, jeśli obliczone wartości
róŜnią się dostatecznie mało, albo po wykonaniu określonej liczby kroków.
ZauwaŜmy jeszcze, Ŝe stosując metody wielokrokowe musimy na początku obliczeń (gdy
znamy tylko y
0
= y(x
0
) ) wyznaczyć y
1
, y
2
, ... . MoŜemy do tego uŜyć np. metody typu
Rungego-Kutty.
III. RÓWNANIA RÓśNICZKOWE CZĄSTKOWE.
4.11 Definicje i klasyfikacja.
W równaniach róŜniczkowych cząstkowych występują pochodne cząstkowe szukanej funkcji
u po zmiennych przestrzennych (x, y, z) i czasie t.
Równania róŜniczkowe cząstkowe opisują wiele waŜnych procesów fizycznych:
przewodnictwo ciepła, dyfuzję, rozkłady potencjałów, ruch drgający; do równań
róŜniczkowych cząstkowych naleŜy teŜ równanie Schrödingera.
Znalezienie jednoznacznego rozwiązania wymaga zadania warunków granicznych. Oprócz
warunków początkowych potrzebujemy warunków brzegowych. Mogą nimi być:
- warunki brzegowe Dirichleta: zadana jest wartość funkcji na brzegu badanego obszaru
- warunki brzegowe von Neumanna: na brzegu badanego obszaru zadana jest wartość
pochodnej normalnej szukanej funkcji
- warunki brzegowe Cauchy’ego: w t = 0 podana jest wartość funkcji i jej pochodnej
normalnej.
Równania róŜniczkowe cząstkowe klasyfikuje się ze względu na współczynniki przy drugich
pochodnych.
Dla równania
0
2
2
2
2
2
=
+
∂
∂
+
∂
∂
+
∂
∂
d
y
u
c
y
dx
u
b
x
u
a
,
gdzie d nie zaleŜy od drugich pochodnych obliczamy wartość b
2
– 4ac.
Równanie nazywamy:
hiperbolicznym,
gdy b
2
– 4ac > 0
parabolicznym,
gdy b
2
– 4ac = 0
eliptycznym,
gdy b
2
– 4ac < 0.
W podręcznikach matematyki omawia się typowe sposoby rozwiązywania równań
róŜniczkowych cząstkowych: metodę rozdzielenia zmiennych i transformaty całkowe (np.
szeregi Fouriera). W następnych punktach przedstawimy jeden ze sposobów numerycznego
podejścia (metodę róŜnic skończonych) do dwu typowych zagadnień.
4.12 Równanie przepływu ciepła w jednym wymiarze.
Równanie przewodnictwa cieplnego ma postać
t
T
z
T
y
T
x
T
∂
∂
=
∂
∂
+
∂
∂
+
∂
∂
2
2
2
2
2
2
2
α
,
gdzie α
2
jest współczynnikiem przewodności cieplnej. Dla przypadku jednowymiarowego
(np. przewodzenie ciepła wzdłuŜ jednorodnego pręta) przyjmuje ono postać
t
T
x
T
∂
∂
=
∂
∂
2
2
2
α
.
Dla uproszczenia dalszych zapisów wprowadźmy oznaczenia u
xx
i u
t
na odpowiednie
pochodne cząstkowe i przyjmijmy, Ŝe uŜywamy jednostek, w których α
2
= 1. Niech pręt ma
końce w punktach x = 0 i x = 1. Mamy wtedy
u
xx
= u
t
Równanie jest drugiego rzędu po współrzędnej przestrzennej i pierwszego rzędu ze względu
na pochodną po czasie.
Potrzebujemy zatem jednego warunku początkowego i dwu warunków
brzegowych. MoŜemy przyjąć, Ŝe znamy początkowy rozkład temperatury wzdłuŜ pręta a
takŜe, Ŝe znamy temperaturę na końcach pręta w kaŜdej chwili czasu. Nasze równanie z
warunkami początkowymi i brzegowymi moŜemy zatem zapisać jako
u
xx
= u
t
u(x, 0) = g(x)
u(0, t) = a(t)
u(1, t) = b(t)
Mamy dwie zmienne x i t, rozwiązania szukamy na siatce punktów (x
i
, t
j
):
t
j
= jk , j = 0, 1, ...
x
i
= ih, i = 0, 1, ..., n+1;
h = 1/(n + 1)
Długości kroku k i h dla czasu i zmiennej przestrzennej mogą być róŜne.
Pochodne występujące w równaniu zastępujemy wyraŜeniami przybliŜonymi:
[
]
)
(
)
(
1
)
(
'
x
f
h
x
f
h
x
f
−
+
≈
[
]
)
(
)
(
2
)
(
1
)
(
''
2
h
x
f
x
f
h
x
f
h
x
f
−
+
−
+
≈
UŜywając ich w naszym równaniu dostajemy
[
] [
]
)
,
(
)
,
(
1
)
,
(
)
,
(
2
)
,
(
1
2
t
x
u
k
t
x
u
k
t
h
x
u
t
x
u
t
h
x
u
h
−
+
=
−
+
−
+
(*)
Co dla punków siatki moŜemy zapisać
)
(
1
)
2
(
1
1
,
,
1
,
1
2
ij
j
i
j
i
ij
j
i
u
u
k
u
u
u
h
−
=
+
−
+
−
+
,
gdzie u
ij
= u(x
i
, t
j
). Przeprowadziliśmy w ten sposób dyskretyzację rozwiązania.
Równość (*) przekształcamy do postaci
ij
j
i
ij
j
i
j
i
u
u
u
u
h
k
u
+
+
−
=
−
+
+
)
2
(
,
1
,
1
2
1
,
,
czyli
j
i
ij
j
i
j
i
su
u
s
su
u
,
1
,
1
1
,
)
2
1
(
+
−
+
+
−
+
=
,
gdzie
2
h
k
s
=
Znając rozwiązanie w trzech punktach u
i-1,j
, u
ij
oraz u
i+1,j
siatki moŜemy je wyznaczyć w
punkcie u
i,j+1
:
t
x
u
i,j+1
u
i+1,j
u
ij
u
i-1,j
PoniewaŜ znamy wartości na osi x (warunek początkowy) oraz na lewym i prawym brzegu
siatki (warunki brzegowe) nie mamy z tym kłopotu i zaczynając od wartości dla j = 0
obliczamy u
i1
a potem dla coraz większych czasów. PoniewaŜ nowa wartość u
i,j+1
wyraŜa się
przez wartości juŜ znane, jest to metoda jawna.
ZauwaŜmy, Ŝe oznaczając U
j
= (u
1j,
..., u
nj
) wektor wartości rozwiązania dla czasu t = jk mamy
U
j+1
= AU
j
gdzie A jest macierzą stopnia n
−
−
−
−
=
s
s
s
s
s
s
s
s
2
1
0
0
0
0
2
1
0
0
2
1
0
0
2
1
L
M
M
M
M
M
L
L
L
A
MoŜna pokazać, Ŝe metoda ta jest stabilna, gdy k/h
2
≤
1/2. Oznacza to konieczność
stosowania bardzo małego kroku czasowego k.
UŜywając innego wzoru na pierwszą pochodną:
[
]
)
(
)
(
1
)
(
'
h
x
f
x
f
h
x
f
−
−
≈
Dostajemy dla naszego zagadnienia zaleŜność
[
] [
]
)
,
(
)
,
(
1
)
,
(
)
,
(
2
)
,
(
1
2
k
t
x
u
t
x
u
k
t
h
x
u
t
x
u
t
h
x
u
h
−
−
=
−
+
−
+
,
czyli na punktach siatki
)
(
1
)
2
(
1
1
,
,
1
,
1
2
−
−
+
−
=
+
−
j
i
ij
j
i
ij
j
i
u
u
k
u
u
u
h
Tym razem równanie wiąŜe jedną znaną wartość dla j-1 z trzema nieznanymi dla j (jest to
więc metoda niejawna), zatem jeśli znamy wszystkie u
i,j-1
to znalezienie wszystkich u
ij
wymaga rozwiązania układu równań. Po przekształceniu do postaci
1
,
,
1
,
1
)
2
1
(
−
+
−
=
−
+
+
−
j
i
j
i
ij
j
i
u
su
u
s
su
dostajemy
AU
j
= U
j-1
gdzie
+
+
−
−
+
−
−
+
=
s
s
s
s
s
s
s
s
2
1
0
0
0
0
2
1
0
0
2
1
0
0
2
1
L
M
M
M
M
M
L
L
L
A
.
Jak widać, macierz A jest trójdiagonalna, więc rozwiązując powyŜszy układ równań
korzystamy z procedury dostosowanej do układu równań liniowych z macierzą
trójprzekątniową.
Jakkolwiek metoda niejawna wymaga większego nakładu pracy, jej zaletą jest lepsza
stabilność w porównaniu z metodą jawną.
4.13 Równanie Laplace’a w dwu wymiarach.
Równanie Laplace’a jest przykładem równania róŜniczkowego cząstkowego, w którym czas
nie występuje. Dla dwu wymiarów ma ono postać
0
2
2
2
2
=
∂
∂
+
∂
∂
y
u
x
u
Typowy warunek brzegowy to warunek Dirichleta – na brzegu obszaru zadane są wartości
funkcji u.
RozwaŜmy to zagadnienie w kwadracie 0 < x < 1, 0 < y < 1
Mamy
u
xx
+ u
yy
= 0
u(x, y) = g(x, y) na brzegu kwadratu
Podobnie jak poprzednio dyskretyzujemy zadanie na siatce (x
i
, y
i
) pokrywającej kwadrat
równomiernie
(x
i
, y
i
) = (ih, jh),
i = 0, ..., n+1; j = 0, ..., n+1;
h = 1/(n+1)
Stosując wzór róŜnicowy
[
]
)
(
)
(
2
)
(
1
)
(
''
2
h
x
f
x
f
h
x
f
h
x
f
−
+
−
+
≈
dostajemy
0
)
2
(
1
)
2
(
1
1
,
1
,
2
,
1
,
1
2
=
+
−
+
+
−
+
−
+
−
j
i
ij
j
i
j
i
ij
j
i
u
u
u
h
u
u
u
h
czyli
0
4
1
,
1
,
,
1
,
1
=
−
−
−
−
+
−
+
−
j
i
j
i
j
i
j
i
ij
u
u
u
u
u
Wartości u
ij
znamy na brzegach kwadratu, natomiast w punktach wewnętrznych
(zaznaczonych kółeczkami) musimy je wyznaczyć rozwiązując układ n
2
równań liniowych.
ZauwaŜmy, Ŝe kaŜde równanie wiąŜe ze sobą pięć współczynników. Oznacza to, Ŝe w
kaŜdym wierszu macierzy układu równań jedynie kilka elementów jest róŜnych od zera.
Macierz takiej postaci to macierz rzadka i do rozwiązywania układu równań stosuje się w tym
przypadku specjalne metody.
y
x