I. MINIMALIZACJA FUNKCJI.
4.1 Ekstrema funkcji. Punkty siodłowe.
Funkcja f (x) jednej zmiennej ma w punkcie x0 minimum (maksimum) lokalne, jeśli istnieje
takie otoczenie punktu x0, \
e dla dowolnego punktu x z tego otoczenia zachodzi f(x) > f(x0)
(f(x) < f(x0) ). 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 x0, w którym funkcja ma ekstremum lokalne znika
jej pierwsza pochodna: f (x0) = 0. Zauwa\
my, \
e zerowanie się pierwszej pochodnej to
warunek konieczny, ale niewystarczający istnienia ekstremum lokalnego.
Sprawdzenie istnienia bądz
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 x0 pierwsza pochodna funkcji przyjmuje wartość zero f (x0) = 0 a druga
pochodna jest ró\
na od zera f  (x0) `" 0, to
jeśli f  (x0) > 0 funkcja ma w x0 minimum lokalne
jeśli f  (x0) < 0 funkcja ma w x0 maksimum lokalne
Jeśli natomiast f  (x0) = 0, to musimy zbadać wy\
sze pochodne.
Uzasadnienie powy\
szego twierdzenia mo\
na oprzeć na wzorze Taylora.
1
f (x0 + h) = f (x0 ) + hf '(x0 ) + h2 f ''(x0 ) +K
2
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 (x0) = 0):
1
f (x0 + h) H" f (x0 ) + h2 f ''(x0 )
2
Gdy f  (x0) > 0, to dla dowolnego dostatecznie małego h mamy f(x0 + h) > f(x0) a to oznacza
istnienie w x0 minimum lokalnego. Analogicznie z faktu f  (x0) < 0 wynika f(x0 + h) < f(x0),
czyli maksimum lokalne w x0. Gdy f  (x0) = 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 = (x1, x2, ..., xn)).
Gdy f(x) > f(x0) dla ka\
dego x w pewnym otoczeniu x0, to funkcja ma w x0 minimum lokalne.
Gdy f(x) < f(x0) dla ka\
dego x w pewnym otoczeniu x0, to funkcja ma w x0 maksimum
lokalne.
Podobnie jak w przypadku jednej zmiennej, warunkiem koniecznym istnienia w punkcie x0
ekstremum lokalnego jest zerowanie się w nim wszystkich pierwszych pochodnych
cząstkowych:
"f "f "f
= = L = = 0 , wszystkie pochodne liczone w punkcie x0
"x1 "x2 "xn
Wektor o składowych równych pochodnym cząstkowym funkcji nazywamy gradientem
funkcji a operację tworzenia gradientu zaznaczamy przez operator"
"f
�ł �ł
�ł �ł
"x1
�ł �ł
"f = M
�ł �ł
"f
�ł �ł
�ł �ł
"xn
�ł łł
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
�ł łł
"2 f "2 f "2 f
L
�ł
"x1"x1 "x1"x2 "x1"xn śł
�ł śł
"2 f "2 f "2 f
�ł śł
�ł"x2"x1 "x2"x2 L "x2"xn śł
H(x) = ,
�ł śł
M M M M
�ł śł
"2 f "2 f "2 f
�ł śł
L
�ł"xn"x1 "xn"x2
"xn"xn śł
�ł �ł
pochodne są tu liczone w punkcie x.
Z wykorzystaniem tych wielkości mo\
emy zapisać wzór Taylora do członu z drugą pochodną
1
f(x0 + h) = f(x0) + hT�"f (x0)+ hTH(x0)h + ...
2
Poniewa\
w punkcie krytycznym "f (x0)= 0, więc drugi człon po prawej stronie wzoru znika
i charakter punktu krytycznego zale\
y od znaku formy kwadratowej hTH(x0)h. Jeśli macierz
H(x0) jest dodatnio określona, to dla dowolnego dostatecznie małego h mamy f(x0 + h) >
f(x0), czyli w x0 jest minimum lokalne funkcji f. Analogicznie, jeśli H(x0) jest określona
ujemnie, to w x0 mamy maksimum lokalne f.
Określoność macierzy H mo\
emy zbadać sprowadzając ją do postaci diagonalnej. Istotnie,
przy diagonalnej macierzy mamy
n
hTHh = (hi' )2 ,
"i
i
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 x0 minimum lokalne. Jeśli wszystkie wartości własne
hesjanu są ujemne, to jest on ujemnie określony i w x0 jest maksimum lokalne. Jeśli część
wartości własnych hesjanu jest dodatnia a część ujemna, to w x0 jest punkt siodłowy.
W szczególności, jeśli jedna wartość własna jest ujemna a pozostałe są dodatnie, to w x0 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 znalez
ć 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
5 -1
� =
2
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 znalez
ć punkt, w którym na zadanej prostej funkcja wielu
zmiennych przyjmuje wartość minimalną.
Jeśli znajdujemy się w punkcie x0 dostatecznie blisko minimum funkcji, mo\
emy spróbować
przybli\
yć funkcję formą kwadratową
1
f(x) = f(x0)  b(x  x0) + A (x  x0)2 ,
2
gdzie b =  f (x0) a A = f  (x0); zale\
ność ta wynika oczywiście z wzoru Taylora.
Ró\
niczkując tę zale\
ność stronami dostajemy
df (x)
= -b + A(x - x0 )
dx
W punkcie xm, w którym funkcja ma minimum jej pierwsza pochodna się zeruje, czyli
 b + A(xm  x0) = 0
skąd znajdujemy przybli\
enie poło\
enia minimum
b
xm = x0 +
A
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
pTAq = 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.
1
f(x) = f(x0) + (x  x0)"f (x0)+ (x  x0)H(x0) (x  x0)
2
Ró\
niczkując stronami i przyrównując pochodną do zera znajdziemy
0 = "f (x0)+ H (xm  x0)
a stąd
xm = x0  H-1�"f (x0)
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łaszczyz
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
znalez
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
dy
= f (x, y)
dx
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óz
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(x0), y (x0) .... 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 "f / "y są funkcjami rzeczywistymi, ograniczonymi, jednowartościowymi i
ciągłymi w pewnym otoczeniu punktu (x0, y0), to w przedziale [x0  h, x0 + h] wewnątrz tego
otoczenia istnieje dokładnie jedno rozwiązanie równania pierwszego rzędu
dy
= f (x, y)
dx
z warunkiem początkowym y(x0) = y0.
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 x0 (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ć
yi' = fi (x, y1, y2 ,K, yn ) i = 1, ..., n
Niewiadomymi są tu funkcje yi zale\
ne od argumentu x. Dodatkowo potrzebne są warunki
początkowe yi(x0) określające wartości poszczególnych funkcji w punkcie x0.
Układ taki mo\
na wyrazić w notacji wektorowej. Niech y będzie wektorem kolumnowym o
składowych yi a f wektorem kolumnowym o składowych fi. Wtedy układ da się zapisać jako
y = f(x, y) y(x0) = y0
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 = y
'
�2 = �1 = y'
L
'
�n = �n-1 = y(n-1)
dostajemy z wyjściowego równania układ równań rzędu pierwszego:
'
�1 = �2
'
�2 = �3
L
'
�n-1 = �n
'
�n = f (x,�1,�2 ,K,�n )
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 �0 = 1.
4.8 Zastosowanie wzoru Taylora. Metoda Eulera.
Rozwiązanie równania
y = f(x,y)
mo\
emy oprzeć na wzorze Taylora
1 1
y(x + h) = y(x) + hy'(x) + h2 y''(x) + h3 y(3) (x) +K
2 3!
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(x0) 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) 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 yn = y(xn) znalez
ć jej przybli\
enie w punkcie odległym o h: yn+1 = y(xn + 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 xn+1.
Metoda Rungego-Kutty rzędu drugiego (zwana metodą Heuna) opisana jest wzorami
F1 = hf (xn , yn )
F2 = hf (xn + h, yn + F1)
1
yn+1 = yn + (F1 + F2 )
2
Najbardziej znana jest metoda Rungego-Kutty rzędu czwartego, której wzory są bardziej
skomplikowane, ale łatwe do zaprogramowania:
F1 = hf (xn , yn )
1 1
F2 = hf (xn + h, yn + F1)
2 2
1 1
F3 = hf (xn + h, yn + F2 )
2 2
F4 = hf (xn + h, yn + F3 )
1
yn+1 = yn + (F1 + 2F2 + 2F3 + F4 )
6
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 znalez
ć 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(x0) = y0
w punktach x1, x2, ... . Przez scałkowanie obu stron powy\
szego równania ró\
niczkowego
otrzymujemy
xn+1
y(xn+1) = y(xn ) + f (x, y(x))dx
+"
xn
Do obliczenia całki mo\
emy wykorzystać kwadraturę z węzłami xi.
W ten sposób mo\
emy otrzymać wzory Adamsa-Bashforta
yn+1 = yn + afn + bfn-1 + cfn-2 + ... (*)
Oznaczyliśmy fi = f(xi,yi), gdzie yi jest przybli\
oną wartością rozwiązania w punkcie xi.
Na przykład, przy równo rozmieszczonych punktach xi = x0 + ih, u\
ywając pięciu węzłów dla
wyznaczenia współczynników kwadratury otrzymać mo\
emy wzór
1
yn+1 = yn + h(1901 fn - 2774 fn-1 + 2616 fn-2 -1274 fn-3 + 251 fn-4 )
720
Wyprowadzając wzory kwadratury w metodzie wielokrokowej mo\
emy te\
zało\
yć, \
e
wartość całki po prawej stronie zale\
y te\
od yn+1. Dostajemy wtedy wzory Adamsa-Moultona
yn+1 = yn + afn+1 + bfn + cfn-1 + ... (**)
np. dla punktów rozmieszczonych jak w przykładzie wy\
ej, wzór
1
yn+1 = yn + h(251 fn+1 + 646 fn - 264 fn-1 +106 fn-2 -19 fn-3).
720
Zauwa\
my, \
e szukane yn+1 występuje tak\
e po prawej stronie powy\
szych wzorów  musimy
je u\
yć do obliczenia fn+1. Jest to przykład metody niejawnej.
Równanie (**) mo\
na rozwiązywać iteracyjnie podstawiając uzyskane yn+1 do prawej strony
wzoru. Wstępne przybli\
enie mo\
na otrzymać u\
ywając wzoru Adamsa-Bashforta (*) do
obliczenia przybli\
onej wartości yn+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 y0 = y(x0) ) wyznaczyć y1, y2, ... . 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
"2u "2u "2u
a + b + c + d = 0 ,
"x2 dx"y "y2
gdzie d nie zale\
y od drugich pochodnych obliczamy wartość b2  4ac.
Równanie nazywamy:
hiperbolicznym, gdy b2  4ac > 0
parabolicznym, gdy b2  4ac = 0
eliptycznym, gdy b2  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ć
�ł �ł
"2T "2T "2T "T
2
ą �ł �ł
+ + = ,
�ł 2 �ł
"x2 "y2 "z "t
�ł łł
gdzie ą2 jest współczynnikiem przewodności cieplnej. Dla przypadku jednowymiarowego
(np. przewodzenie ciepła wzdłu\
jednorodnego pręta) przyjmuje ono postać
"2T "T
2
ą = .
"x2 "t
Dla uproszczenia dalszych zapisów wprowadz
my oznaczenia uxx i ut 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
uxx = ut
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
uxx = ut
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 (xi, tj):
tj = jk , j = 0, 1, ...
xi = 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
f ' (x) H" [f (x + h) - f (x)]
h
1
f ''(x) H" [f (x + h) - 2 f (x) + f (x - h)]
h2
U\
ywając ich w naszym równaniu dostajemy
1 1
[u(x + h,t) - 2u(x,t) + u(x - h,t)]= [u(x,t + k) - u(x,t)] (*)
h2 k
Co dla punków siatki mo\
emy zapisać
1 1
(ui+1, - 2uij + ui-1, j ) = (ui, j+1 - uij ) ,
h2 j k
gdzie uij = u(xi, tj). Przeprowadziliśmy w ten sposób dyskretyzację rozwiązania.
Równość (*) przekształcamy do postaci
k
ui, j+1 = (ui+1, j - 2uij + ui-1, j ) + uij ,
h2
czyli
k
ui, j+1 = sui-1, j + (1- 2s)uij + sui+1, j , gdzie s =
h2
Znając rozwiązanie w trzech punktach ui-1,j, uij oraz ui+1,j siatki mo\
emy je wyznaczyć w
punkcie ui,j+1:
t
ui,j+1
ui-1,j uij ui+1,j
x
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 ui1 a potem dla coraz większych czasów. Poniewa\
nowa wartość ui,j+1 wyra\
a się
przez wartości ju\
znane, jest to metoda jawna.
Zauwa\
my, \
e oznaczając Uj = (u1j, ..., unj) wektor wartości rozwiązania dla czasu t = jk mamy
Uj+1 = AUj
gdzie A jest macierzą stopnia n
1- 2s s 0 L 0
�ł łł
�ł śł
s 1- 2s s L 0
�ł śł
�ł śł
A = 0 s 1- 2s L 0
�ł śł
M M M M M
�ł śł
�ł
0 0 0 L 1- 2sśł
�ł �ł
Mo\
na pokazać, \
e metoda ta jest stabilna, gdy k/h2 d" 1/2. Oznacza to konieczność
stosowania bardzo małego kroku czasowego k.
U\
ywając innego wzoru na pierwszą pochodną:
1
f '(x) H" [f (x) - f (x - h)]
h
Dostajemy dla naszego zagadnienia zale\
ność
1 1
[u(x + h,t) - 2u(x,t) + u(x - h,t)]= [u(x,t) - u(x,t - k)],
h2 k
czyli na punktach siatki
1 1
(ui+1, - 2uij + ui-1, j ) = (uij - ui, j-1)
h2 j k
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 ui,j-1 to znalezienie wszystkich uij
wymaga rozwiązania układu równań. Po przekształceniu do postaci
- sui-1, j + (1+ 2s)uij - sui+1, j = ui, j-1
dostajemy
AUj = Uj-1
gdzie
1+ 2s - s 0 L 0
�ł łł
�ł śł
- s 1+ 2s - s L 0
�ł śł
�ł śł
A = 0 - s 1+ 2s L 0 .
�ł śł
M M M M M
�ł śł
�ł
0 0 0 L 1+ 2sśł
�ł �ł
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ć
"2u "2u
+ = 0
"x2 "y2
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
uxx + uyy = 0
u(x, y) = g(x, y) na brzegu kwadratu
Podobnie jak poprzednio dyskretyzujemy zadanie na siatce (xi, yi) pokrywającej kwadrat
równomiernie
(xi, yi) = (ih, jh), i = 0, ..., n+1; j = 0, ..., n+1; h = 1/(n+1)
Stosując wzór ró\
nicowy
1
f ''(x) H" [f (x + h) - 2 f (x) + f (x - h)]
h2
dostajemy
1 1
(ui-1, - 2uij + ui+1, j ) + (ui, j-1 - 2uij + ui, j+1) = 0
h2 j h2
czyli
4uij - ui-1, j - ui+1, j - ui, j-1 - ui, j+1 = 0
Wartości uij znamy na brzegach kwadratu, natomiast w punktach wewnętrznych
(zaznaczonych kółeczkami) musimy je wyznaczyć rozwiązując układ n2 równań liniowych.
y
x
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.