Układy równań różniczkowych zwyczajnych Wstęp

(1)

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 t₁( ), y t₂( ) 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

   



  

 ⁽¹⁾

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 LotkiVolterry, 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

  

   

 ⁽²⁾

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 ^dy_{dt y} _y^y_t) 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 y₂, to otrzymalibyśmy

(2)

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 ay₂ 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 y₁ 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 t₁ y t₂( )) 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

  

   

 ⁽⁵⁾

Co ciekawe dokładnie taki sam układ uzyskuje się gdy rozważymy pewien teoretyczny mechanizm reakcji zaproponowany przez Alfreda Lotkę¹ 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

 



(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 molm^-

3 lub moldm^-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 polskoamerykaoski chemik, biofizyk i statystyk, twórca podstaw nauki od dynamice populacji.

(3)

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 LotkiVolterry. Jak widać rozwiązania są okresowe. W opisie chemicznym reakcji AB 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 a_s , a prędkośd styczna v , to z II zasady dynamiki Newtona otrzymujemy zależnośd

sin ,

m   k mg  (8)

czyli m k mgsin0, gdzie: mmasa kulki, długośd linki (pręta), kwspółczynnik oporu ośrodka, gprzyspieszenie grawitacyjne (g9,81 m/s ).² Oznaczmy k m/ ,²  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

(4)

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

  



  



(11)

z warunkami początkowymi y t₁( )₀  y₁⁰, ,y t_n( )₀  y_n⁰, gdzie liczby t₀, y₁⁰, ,y_n⁰ 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

    



     



(12)

W układzie tym funkcje: a_ij a t_ij( ), f_i f t_i( ) :  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 yf t (13)

gdzie występują macierz ( )A t i wektory kolumnowe ( ),y t f t( ) :

(5)

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 a_{i 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 t₁( ), y t₂( ), a współczynniki , , ,a b c d zamiast a₁₁,a₁₂,a₂₁,a₂₂. 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), aybycy0. 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)cb² (a d)adbc0.

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 ²3 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

(6)

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  12 .

Rozwiązanie ogólne ma postad

1 2

1 1 2 2

( ) ,

( ) ( ) ( ) ,

t t

x t C be C be

y t C a e C a e

 

 

    (18)

gdzie C C₁, ₂ są dowolnymi stałymi.

Przypadek b) Δ < 0  12 :    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)

Przypadek c) Δ = 0   ₁ ₂  .

Rozwiązanie ogólne zależy dodatkowo od relacji pomiędzy a b c d, , , . Pojawiają się funkcje liniowe. Na przykład, gdy zachodzi ad, to



²



²

2

1 2

1 2 2

( ) 2 ,

( ) (( ) ( ) ) ,

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)

Zadania

Zad. 1) Zamieo na układ równao podane równania skalarne drugiego rzędu.

a) y y²0. b) y 1 cosy3 .y b) ^y^^²^y^^^{t y}² ^^0. c) ya 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) ²

3

x x y

y x y

  

    

 c) ²

4 6

x x y

y x y

   

   



Zad. 3) Niech ( ),x t y t( ) będzie dowolnym rozwiązanie układu LotkiVolterry. Pokazad, że wtedy funkcja ( )f t cx t( )  dln ( )x t  a y t( )  bln ( )y t jest stała względem czasu.

(7)

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. Przedyskutowad charakter rozwiązao w zależności od relacji pomiędzy parametrami , ,m k g.