Układy równań różniczkowych zwyczajnych
Wstęp
W zastosowaniach bardzo często zamiast pojedynczego równania występują układy równao równao różniczkowych zwyczajnych (ang. System of Ordinary Differential Equations, ODEs). Oznacza to, że szukamy kilku funkcji, których pochodne spełniają pewne związki. W zagadnieniach kinetyki reakcji chemicznych funkcje y t1( ), y t2( ) mogą opisywad stężenia reagujących substancji. Na przykład pewien układ reakcji zwany Brusselatorem może prowadzid do następującego układu równao
2
1 1 2 1
2
2 1 1 2
1 4 ,
3 .
y y y y
y y y y
(1)
Bardzo rzadko się zdarza, aby układy takie jak (1) miały rozwiązania, które można wyrazid „wzorami”
analitycznymi. W tej sytuacji równania takie możemy badad numerycznie, poprzez uzyskiwanie przybliżeo dla konkretnych warunków początkowych. Można też analizowad własności rozwiązao bez wyrażania ich wzorami, ale w oparciu o pewne twierdzenia i matematyczne teorie. Zajmuje się tym dział matematyki zwany teorią układów dynamicznych. Obie metody – numeryczna i jakościową – można łączyd.
Innym obszarem, który dostarcza przykładów układów równao różniczkowych zwyczajnych są różne modele biologiczne. Dobrze znanym jest tzw. układ LotkiVolterry, albo inaczej model drapieżnik- ofiara. Zakładamy, że na danym obszarze występują dwa gatunki: drapieżniki (np. lisy) i ofiary (np.
zające). Drapieżniki żywią się ofiarami. Jeżeli wprowadzimy oznaczenia:
1( )
y t liczebnośd populacji ofiar (można też mówid o gęstości populacji ofiar) w chwili t ,
2( )
y t liczebnośd populacji drapieżników (lub gęstości populacji drapieżników) w chwili t ,
to dośd prosty model opisujący jak zmienia się w czasie populacja ofiar i drapieżników, które na siebie wzajemnie oddziałują, zawarty jest w następującym układzie równao
1 2 1
2 1 2
( ) ,
( ) .
y b ay y y cy d y
(2)
W układzie tym występują dodatnie stałe a b c d, , , charakteryzujące jakośd środowiska oraz możliwości obu gatunków. Jeżeli przepiszemy pierwsze równanie układu (2) następująco
1
2 1
y , y b ay
(3)
to widzimy, że względna zmiana w czasie liczebności ofiar (wszak 1 1
1 1
1 1/ 1 dydt y yyt) y y jest proporcjonalna do b, a zatem parametr ten określa możliwości rozmnażania się ofiar: im więcej jest osobników tym więcej będzie ich przybywad. Gdyby w równości tej pominąd składnik a y2, to otrzymalibyśmy
1
1 1
1
( ) (0) bt,
y b y t y e
y
czyli wzrost wykładniczy zależny od b. Można powiedzied, że parametr ten określa zdolności reprodukcyjne ofiar oraz jakośd środowiska, które są stałe nie zależą od liczebnośd (lub gęstości) żadnej z populacji. Z drugiej strony mamy jednak w równaniu (3) czynnik hamujący wzrost liczby ofiar – są nim drapieżnicy. Mianowicie im więcej drapieżników, tym więcej potrzebują pożywienia czyli ofiar. Czynnik ay2 zawiera parametr a, który charakteryzuje skutecznośd drapieżników w polowaniu na ofiary.
Podobnie można przeanalizowad drugie równanie układu (2) zapisując je w formie
2 1 2
y .
cy d y
(4)
Interpretacja składnika c y1 jest taka, że im więcej ofiar tym szybciej rozmnażają się drapieżniki (jest dużo pożywienia). Dlatego c może oznaczad jakośd tego pożywienia (ofiary) oraz zdolności reprodukcyjne drapieżników.
W przypadku gdy mamy układ dwóch równao to zamiast numerowanie funkcji ( ( ),y t1 y t2( )) używa się często np. oznaczeo ( ),x t y t( ). Wtedy układ (2) zapiszemy następująco
( ) ,
( ) .
x b ay x y cx d y
(5)
Co ciekawe dokładnie taki sam układ uzyskuje się gdy rozważymy pewien teoretyczny mechanizm reakcji zaproponowany przez Alfreda Lotkę1 w 1920 roku, w którym występują dwa autokatalityczne etapy
1
2
3
(i) A X 2X
(ii) X Y 2Y
(iii) Y B
k
k
k
(6)
gdzie A i B to substrat i produkt, X i Y to formy pośrednie, a k1, k2, k3 to stałe szybkości odpowiednich reakcji. Jeżeli symbolem *X+ oznaczymy stężenie molowe (najczęściej wyrażane w jednostkach molm-
3 lub moldm-3), to równania chemiczne (i) oraz (ii) prowadzą do następującego układu równao różniczkowych zwyczajnych
1 2
2 2
[ ] [ ][ ] [ ][ ], [ ] [ ][ ] [ ],
d X k A X k X Y
dt
d Y k X Y k Y dt
(7)
1 Alfred Lotka (1880-1949) urodzony we Lwowie polskoamerykaoski chemik, biofizyk i statystyk, twórca podstaw nauki od dynamice populacji.
gdzie stężenia *X+, *Y+ i *A+ są funkcjami czasu. Gdy stężenie substratu A jest stałe (np. składnik jest w nadmiarze i przez pewien czas ubytek można zaniedbad lub jest on uzupełniany w reaktorze chemicznym), to wtedy k1*A+ = const i układ równao (7) jest formalnie identyczny z układem (5). Inne jest tylko znaczenie symboli: ( ), ( )x t y t to ilości osobników dwóch gatunków w danym terenie, a *X+,
*Y+ to stężenia odpowiednich substancji w reaktorze chemicznym.
Rysunek 1. Przykładowe rozwiązanie układu LotkiVolterry. Jak widać rozwiązania są okresowe. W opisie chemicznym reakcji AB z mechanizmem (6) możemy powiedzieć, że stężenia form pośrednich X i Y wykazują oscylacje o stałej amplitudzie.
Układy równao różniczkowych zwyczajnych pierwszego rzędu (czyli takie, w których występują tylko pochodne pierwszego rzędu) pojawiają się także gdy przekształcamy równania różniczkowe wyższych rzędów. Rozważamy dla przykładu znane nam już równanie wahadła matematycznego w ośrodku stawiającym opór (Błąd! Nie można odnaleźd źródła odwołania.). Jeżeli ( ) t oznacza kąt wychylenia (liczony względem pionu), to z rysunku widad, że w kierunku stycznym do toru (który jest okręgiem o promieniu ) działają dwie siły: mg sin (składowa siły ciężkości m g w kierunku stycznym) oraz
k v (siła oporu ośrodka – zakładamy, że jest proporcjonalna do prędkości i oczywiście przeciwnie skierowana). Ponieważ przyspieszenie styczne as , a prędkośd styczna v , to z II zasady dynamiki Newtona otrzymujemy zależnośd
sin ,
m k mg (8)
czyli m k mgsin0, gdzie: mmasa kulki, długośd linki (pręta), kwspółczynnik oporu ośrodka, gprzyspieszenie grawitacyjne (g9,81 m/s ).2 Oznaczmy k m/ ,2 g/ i przepisujemy równanie (8)
2sin .
(9)
Jeżeli wprowadzimy teraz oznaczenia: y1, y2, to powyższe równanie można sprowadzid do następującego układu równao pierwszego rzędu
1 2
2
2 2 1
,
sin .
y y
y y y
(10)
Rysunek 2. Wahadło matematyczne o długości .
W ogólnym przypadku może występowad n niewiadomych funkcji y1(t),…, yn(t), których ewolucje czasowe są powiązane. Oznacza to, że danych jest n równao różniczkowych zwyczajnych (czyli układ równań):
1 1 1
1
( , , , ),
( , , , ),
n
n n n
y f t y y
y f t y y
(11)
z warunkami początkowymi y t1( )0 y10, ,y tn( )0 yn0, gdzie liczby t0, y10, ,yn0 są dane.
Układy równań liniowych
W ogólnym przypadku układ liniowy pierwszego rzędu ma postad
1 11 1 1 1
1 1
( ) ( ) ( ),
( ) ( ) ( ).
n n
n n nn n n
y a t y a t y f t
y a t y a t y f t
(12)
W układzie tym funkcje: aij a tij( ), fi f ti( ) : J są dane i określone na odcinku .J W praktyce najczęściej zakładamy, że są one przynajmniej ciągłe. Wtedy problem początkowy dla układu (12) ma jednoznaczne rozwiązanie dla każdych warunków początkowych.
Układ (12) można zapisad w sposób zwarty używając notacji macierzowej (algebra liniowa):
( ) ( ),
y A t yf t (13)
gdzie występują macierz ( )A t i wektory kolumnowe ( ),y t f t( ) :
11 1 1 1
1
( ) ( ) ( ) ( )
( ) , ( ) , ( ) .
( ) ( ) ( ) ( )
n
n nn n n
a t a t y t f t
A t y t f t
a t a t y t f t
(14)
Dokładniej zajmiemy się szczególnym przypadkiem: dwa równania (n=2), stałe współczynniki (tzn.
funkcje ai j( )t nie zależą od czasu). Zaczniemy od równania jednorodnego ( ( ) 0).f t Dla wygody funkcje niewiadome będziemy oznaczad ( ), ( )x t y t zamiast y t1( ), y t2( ), a współczynniki , , ,a b c d zamiast a11,a12,a21,a22. Mamy zatem układ
, , x ax by y cx dy
(15)
którego macierzą jest a b .
A c d
W przypadku układu (15) można zawsze wyrazid rozwiązanie prostymi wzorami. Sytuacja jest tutaj bardzo podobna do równania skalarnego liniowego drugiego rzędu ze stałymi współczynnikami (które było rozważane wcześniej), aybycy0. Musimy tym razem znaleźd wartości własne macierzy układu, czyli rozwiązad równanie
det a b 0,
c d
(16)
czyli (a)(d)cb2 (a d)adbc0.
Oczywiście (16) jest równaniem kwadratowym na , które ma zawsze rozwiązanie (w dziedzinie zespolonej, ). Równanie (16) nazywamy równaniem charakterystycznym (macierzy, układu równao), a sam wielomian zmiennej – wielomianem charakterystycznym. Pierwiastki wielomianu charakterystycznego nazywamy wartościami własnymi (macierzy, układu). Na przykład dla układ
2 3 , ,
x x y
y x y
(17)
ma macierz I równanie charakterystyczne:
2 3 2 3
, det (2 )(1 ) 3 0.
1 1 1 1
A
Skąd (2)(1) 3 23 5 0, czyli pierwiastki charakterystyczne układu (17) wynoszą
1 2 1 2
3 11 3 11 3 11 3 11
, , .
2 2 2 i 2 2 i 2
Aby podad rozwiązanie układu (15) można obliczyd wektory własne macierzy układu dla każdej wartości własnej, lub bezpośrednio rozpatrzyd trzy zasadnicze przypadki: (a) 0 (dwa rzeczywiste
pierwiastki), (b) 0 (dwa zespolone pierwiastki, (c) 0 (jeden rzeczywisty pierwiastek podwójny (musi byd rzeczywisty, bo współczynniki równania charakterystycznego są liczbami rzeczywistymi)).
Przypadek a) Δ > 0 12 .
Rozwiązanie ogólne ma postad
1 2
1 2
1 2
1 1 2 2
( ) ,
( ) ( ) ( ) ,
t t
t t
x t C be C be
y t C a e C a e
(18)
gdzie C C1, 2 są dowolnymi stałymi.
Przypadek b) Δ < 0 12 : 1 i , 2 i . Rozwiązanie ogólne ma postad
1 2
1 2 1 2
( ) ( sin cos ),
( ) ((( ) )sin ( ( ) ) cos ),
t t
x t be C t C t
y t e a C C t C a C t
(19)
gdzie C C1, 2 są dowolnymi stałymi.
Przypadek c) Δ = 0 1 2 .
Rozwiązanie ogólne zależy dodatkowo od relacji pomiędzy a b c d, , , . Pojawiają się funkcje liniowe. Na przykład, gdy zachodzi ad, to
2
22
1 2
1 2 2
( ) 2 ,
( ) (( ) ( ) ) ,
a d
a d
t C
a d
t
x t b C C t e
y t d a C C d a C t e
(20)
gdzie C C1, 2 są dowolnymi stałymi.
Zadania
Zad. 1) Zamieo na układ równao podane równania skalarne drugiego rzędu.
a) y y20. b) y 1 cosy3 .y b) y2yt y2 0. c) ya t y( ) b t y( ) c t( ).
Zad. 2) Podaj rozwiązanie ogólne układów ze stałymi współczynnikami.
a) x x y y x y
b) 2
3
x x y
y x y
c) 2
4 6
x x y
y x y
Zad. 3) Niech ( ),x t y t( ) będzie dowolnym rozwiązanie układu LotkiVolterry. Pokazad, że wtedy funkcja ( )f t cx t( ) dln ( )x t a y t( ) bln ( )y t jest stała względem czasu.
Zad. 4) Podaj rozwiązanie problemu wahadła matematycznego przyjmując założenie małych wychyleo. Oznacza to, że przybliżamy sinus dla małych kątów znanym wzorem (szereg Taylora!):
sin . Tak więc zamiast (8) rozważamy równanie
0.
m k mg
Dla uproszczenia rozważyd następujący warunek początkowy: (0) 0, (0)0. Przedyskutowad charakter rozwiązao w zależności od relacji pomiędzy parametrami , ,m k g.