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
Wyszukiwarka
Podobne podstrony:
minimalizacja funkcji logicznych
02 Minimalizacja funkcji logicznyc (2)
Algorytmy genetyczne w minimalizacji funkcji logicznych 5 część 2
Pomiary wielkosci elektrycznych Minimalizacja funkcji tablica
Minimalizacja funkcji logicznych 3
Algorytmy genetyczne w minimalizacji funkcji logicznych 5 część 1
Algorytmy genetyczne w minimalizacji funkcji logicznych 5 część 3
Minimalizacja funkcji logicznych[1]
Minimalizacja funkcji 3 lub 4 z Nieznany
5w Minimalizacja funkcji 2
5w Minimalizacja funkcji 1
minimalizacja funkcji logicznych
Algorytmy genetyczne w minimalizacji funkcji logicznych Karina Murawko
Pomiary wielkości elektrycznych Instrukcja do ćw 10 Minimalizacja funkcji – tablicami Karnaugh
więcej podobnych podstron