( ) ( ) Całkowanie równań różniczkowych I rzędu

Całkowanie równań różniczkowych I rzędu

Omówimy tutaj tylko całkowanie równań różniczkowych zwyczajnych I rzędu, czyli postaci:

dy

dx = F ( x , y)

i oczywiście szukamy rozwiązania, czyli funkcji:

y= f (x )+C

gdzie C oznacza dowolną stałą, co w efekcie oznacza, że rozwiązaniem równania różniczkowego jest rodzina funkcji różniących się o tę stałą. Ponieważ w trakcie całkowania numerycznego rozwiązanie podajemy za pomocą jego wartości w poszczególnych punktach na osi 0X, więc musimy ustalić, któ- rą z powyższych funkcji obliczamy. Jednoznacznie ustala nam to podanie wartości początkowej y0

w punkcie początkowym x0 – takie podanie punktu początkowego jednoznacznie ustala wartość sta- łej C.

Najbardziej popularną metodą całkowania równań różniczkowych zwyczajnych I rzędu jest metoda Rungego-Kutty 4 rzędu (RK4rz). Formalnie metodami RK są również metoda punktu środko- wego (RK2rz) i metoda Eulera (RK1rz), pierwotnie wcześniej opracowane.

Poniżej podaję wzory na obliczenie wartości funkcji w punkcie xi+1, gdy znamy wartość tej funkcji w punkcie poprzednim xi dla tych trzech metod.

Oznaczenia:

h — krok całkowania równy xi+1 – xi,

yi — wyliczona wartość funkcji w punkcie xi, yi+1 — wyliczona wartość funkcji w punkcie xi+1,

I. Metoda Eulera

k

= F (x

, y

) y

= y

+ h⋅k

II. Metoda punktu pośredniego

k

= F (x

, y

)

k

= F ( ^x

+ h , y

+ h⋅k

)

y

= y

+ h

2 ( k

+k

)

III. Metoda Rungego-Kutty 4 rzędu

k

= F (x

, y

)

k

= F ( ^x

^{+ 2, h y}

^{+ h 2 ⋅ k}

)

k

= F ( ^x

^{+ 2, h y}

^{+ h 2 ⋅ k}

)

k

=F ( ^x

+ h , y

+ h⋅k

)

y

= y

+ h

6 ( k

+ 2⋅k

+2⋅k

+k

)

Realizację wszystkich trzech metod pokażę w arkuszu kalkulacyjnym na przykładzie prostego rów- nania różniczkowego:

2

y    x y, czyli dla ^{F ,}

 

Metoda Eulera

1. Na arkusz wprowadziłem kolumnę wartości x (na arkuszu krok h wynosi 0,01) oraz wpisałem wartość y0 do pierwszej komórki w kolumnie y (na arkuszu w komórce C2). Jak wprowadza się szereg liczb, pokazałem w temacie Funkcje uwikłane. Wartość początkowa w kolumnie y (komórka B2) jest wpisywana (nie liczona).

2. Korzystając z wartości x i y obliczyłem w pierwszym wierszu tabeli (komórka B2) wartość k1

według podanego wzoru – jest to prawa strona równania.

3. W komórce C3 poniżej wartości początkowej y0 obliczyłem wartość następnego y, według poda- nego wzoru, czyli y poprzednie (adres komórki wyżej) plus h·k1.

4. Skopiowałem wzór na k1 do następnego wiersza — w ten sposób dostałem wiersz wypełniony formułami dla wszystkich kolumn tabeli.

5. W ostatnim etapie skopiowałem formuły z drugiego wiersza wzdłuż całej tabeli.

Metoda punktu pośredniego

1. Na arkusz wprowadziłem kolumnę wartości x (na arkuszu krok h wynosi 0,01) oraz wpisałem wartość y0 do pierwszej komórki w kolumnie y (na arkuszu w komórce D2).

2. Korzystając z wartości x i y obliczyłem w pierwszym wierszu tabeli wartość k1 (komórka B2) i k2

(komórka C2)według podanych wzorów – dla k1 jest to (jak w metodzie Eulera) prawa strona równania, a dla k2 modyfikuję prawą stronę równania tak, że zamiast każdego x podstawiam w równaniu x+h, a zamiast każdego y podstawiam w równaniu y+h·k1, gdzie jako wartość k1

podaję adres komórki B2.

3. W komórce D3 poniżej wartości początkowej y0 obliczyłem wartość następnego y, według poda- nego wzoru – czyli y poprzednie plus krok h razy średnia z k1 i k2.

4. Skopiowałem wzory na k1 i k2 do następnego wiersza — w ten sposób dostałem wiersz wypełnio- ny formułami dla wszystkich kolumn tabeli.

5. W ostatnim etapie skopiowałem formuły z drugiego wiersza wzdłuż całej tabeli.

Metoda Rungego-Kutty

1. Na arkusz wprowadziłem kolumnę wartości x (na arkuszu krok h wynosi 0,01) oraz wpisałem wartość y0 do pierwszej komórki w kolumnie y (na arkuszu w komórce F2).

2. Korzystając z wartości x i y obliczyłem w pierwszym wierszu tabeli wartość k1 (komórka B2), k2

(komórka C2), k3 (komórka D2) i k4 (komórka E2), według podanych wzorów – dla k1 jest to (jak w metodzie Eulera) prawa strona równania, dla k2 modyfikuję prawą stronę równania tak, że zamiast każdego x podstawiam w równaniu x+h/2, a zamiast każdego y podstawiam w równa- niu y+h/2·k1, gdzie jako wartość k1 podaję adres komórki B2, dla k3 modyfikuję prawą stronę równania tak, że zamiast każdego x podstawiam w równaniu x+h/2, a zamiast każdego y podstawiam w równaniu y+h/2·k2, gdzie jako wartość k2 podaję adres komórki C2, ostatecznie dla k4 modyfikuję prawą stronę równania tak, że zamiast każdego x podstawiam w równaniu x+h, a zamiast każdego y podstawiam w równaniu y+h·k3, gdzie jako wartość k3 podaję adres komórki D2.

3. W komórce F3 poniżej wartości początkowej y0 obliczyłem wartość następnego y, według poda- nego wzoru – czyli y poprzednie plus krok h razy średnia ważona z k1, k2, k3 i k4, z wagami 1, 2, 2, 1.

4. Skopiowałem wzory na k1, k2, k3 i k4 do następnego wiersza — w ten sposób dostałem wiersz wypełniony formułami dla wszystkich kolumn tabeli.

5. W ostatnim etapie skopiowałem formuły z drugiego wiersza wzdłuż całej tabeli.

Przykłady

Spadek w ośrodku płynnym

W zakładce Skoczek pokazałem przykład rozwiązania równania spadku w ośrodku płynnym, gdy opory ruch są proporcjonalne do kwadratu prędkości. Jest tak w przypadku, gdy ruch nie jest lami- narny.

Na ciało spadające działa siła grawitacji Fg=m·g oraz siła oporu proporcjonalna do kwadratu prędkości Fo=k·v². Przy spadku swobodnym siły te mają przeciwne zwroty, więc siła wypadkowa Fw=Fg−Fo=m·g− k·v². Dzieląc siłę wypadkową przez masę ciał spadającego, dostajemy przyspiesze- nie aw=g− (k/m)·v², a ponieważ przyspieszenie jest pochodną prędkości po czasie, więc dostajemy równanie różniczkowe na prędkość podczas spadania:

dv

dt =g− k m v

Warto zauważyć, że prawa strona równania różniczkowego jest zależna tylko od zmiennej zależnej v (odpowiednik y w poprzednich przykładach). Zmienna niezależna t (odpowiednik x w poprzednich przykładach) po prawej stronie nie występuje. Wspólczynnik k w tym równaniu zależy od pola prze- kroju poprzecznego (względem ruchu) ciała które spada – A, gęstości ośrodka, w którym odbywa się ruch – ρ i czynnika kształtu – f:

k = A ρ f

W arkuszu rozwiązywałem problem skoczka spadochronowego, który najpierw spada swobodnie aż

2

do zrównoważenia siły grawitacji przez siłę, a potem otwiera spadochron. Ponieważ w tym wypadku czynnik kstałtu się zmienia, podzieliłem problem na dwa etapy. Przyjąłem masę skoczka m=80 kg.

Gdy w równaniach różniczkowych występują jakieś parametry, korzystnie jest zapisać ich wartości w komórkach obok i w równaniach odwoływać się do adresów ich wartości adresami bezwzględnymi (na arkuszu komórki F1:F7 dla pierwszego etapu i M1:M7 dla drugiego etapu).

W etapie swobodnego spadku przyjąłem przekroj poprzeczny A odpowiadający mniej więcej powie- rzchni dorosłego mężczyzny z drobnym naddatkiem na ubranie i rozpostate ramiona (oczywiście wszystko to jest szacunkowe), czynnik kształtu f (współczynnik oporu) przyjąłęm podawany w litera- turze dla wyprostowanego człowieka, a gęstość ośrodka, to gęstość powietrza w temperaturze 20°C pod normalnym ciśnieniem; oczywiście g to przyspieszenie grawitacyjne na powierzchni Ziemi, któ- re nie zmienia się znacząco na odległości przelotu samolotu.

Krok w pierwszym etapie przyjąłem równy 1 s i to wystarczyło, by rozwiązanie przebiegało w rozsąd- nych granicach. I tutaj:

UWAGA

W przypadku całkowania równań różniczkowych, bardzo ważny jest dobór kroku całkowania. Jeśli przyjmiemy za duży krok, to rozwiązanie może nie mieć nic wspólnego z rzeczywistością. Najlepiej patrzeć na przyrost zmian w wartości zmiennej zależnej. Jeśli są one zbyt duże w stosunku do war- tości z poprzedniego kroku, to jest to sygnał, żeby zmniejszyć krok. W zasadzie nic nie stoi na przeszkodzie, żeby ten krok zmieniać płynnie, wtedy trzeba krok w równaniach wpisywać jako różnicę przyrost wartości zmiennej niezależnej (wcześniej x, tutaj t) pomiędzy krokiem poprzednim a obecnym, zamiast sztywnej wartości, jak to jest w arkuszu. Zawsze lepiej dać mniejszy krok, bo płacimy za to jedynie czasem obliczeń, a potem ewentualnie go wydłużać, sprawdzając, jak to wpływa na wyniki. Oczywiście ne zawsze tak się da.

Sledząc wynik rozwiazania, widzimy, że w pewnym momencie wartość się zeruje i dalej nie ma zmian w prędkości. Jest to tzw. Prędkość graniczna albo wysycenia. Jak widać jest ona dosyć duża, bo aż 41 m/s czyli 147 km/h.

Tę prędkość wstawiłem jako wartość początkową w dugim etapie, kiedy badam prędkośc po otwo- rzeniu spadochronu. Zmieniają się wtedy wspólczynnik oporu (kształtu) i powierzchnia przekroju po- przecznego, bo teraz głownym źródłem oporu jest kształt czaszy spadochronu i jej przekrój poprze- czny. Przyjąłem podawaną w zaleceniach powierzchnię 9 m² na 10 kg masy, czyli , w tym przypadku 72 m².

W tym wypadku musiałem radykalnie zmniejszyć krok, bo na początku, że względu na dużą prędkość zmiany są gwałtowne i zbyt duże zmiany powodowały niestabilność rozwiązania. Chętnym zalecam przyjęcie kroku, jak w pierszym etapie i prześledzenie na wykresie tych niestabilnych rozwiązań i porównanie ich z tymi, kiedy krok jest dostatecznie mały.

W rozwiązaniu widać, że na tym etapie następuje gwałtowne hamowanie, wręcz szarpnięcie, i końcowa prędkość ustala się na 3,6 m/s, co mniej więcej odpowiada zalecanej prędkości końcowej w skokach spadochronowych.

Zachęcam do sporządzenia sobie wykresów zależności prędkości v od czasu t, a także wstawienia do równań swoich parametrów i dobrania odpowiedniej czaszy spadochronu.

Równanie Streetera-Phelpsa

Oryginalnie równanie Streetera-Phelpsa opisuje dynamikę zmian deficytu tlenowego w wodzie pod wplywe procesów deaeracji wskutek zużycia tlenu procesami biochemicznymi i aeracji wskutek dostarczania tlenu do środowiska wodnego. Zmieniając znak prawej strony, otrzymamy równanie dynamiki zmian zawartości tlenu w wodzie.

dC

dt =−k

L

e

−k

(C

−C )

gdzie C to koncentracja tlenu w wodzie, Cs to koncentracja tlenu w stanie nasycenia, L0 to biochemiczne zapotrzebowanie na tlen (BZT), kb to stała szybkości usuwania BZT, a ka to stała reaeracji. Na arkuszu Streeter-Phelps umieściłem tabelkę ze zmierzonymi w pewnym zbiorniku wodnym wartościami C (O_pocz) Cs (O_nasyc) i L0 (BZT). Jeśli się Państwo przyjrzycie, to stężenie nasycenia się zmienia, co nie jest dziwne, bo rozpuszczalność gazów w wodzie zależy silnie od

temperatury, natomiast, co już może dziwić, czasami widać, że zarejestrowane stężenie tlenu jest wyższe niż stężenie nasycenia, ale i to jest normalne i zdarza się, chociażby wskutek zmian tem- peratury. W obliczeniach przyjąłem ostatnie wartości z tabelki. Spotykane wartości Kb i Ka podałem obok rubryk, które zawierają wartości przyjęte w rozwiązaniu: dla kb to są wartości 0,05 – 0,5, a dla ka to 0,4 – 1,5. Kolejności tworzenia tabelki z rozwiązaniami nie będę tu omawiał, bo opisywałem to na wcześniejszych przykładach. Proszę sobie ją przeanalizować.

Zadania do wykonania

1. W przykładzie Skoczek zastosować metodę punktu środkowego i sprawdzić, czy nie da się dzięki temu zwiększyć krok w przypadku otworzenia spadochronu.

2. W przykładzie Streeter-Phelps w tabelce parametrów wstawić inne wartości BZT, Onasyc i Opocz

i zobaczyć na wykresie zmiany zawartósci tlenu rozpuszczonego w wodzie. Wartości dobierać z tablicy zamieszczonej na arkuszu.

3. Zmiana w czasie t wysokości h napełnienia zbiornika, z którego wypływa woda przez otwór w dnie jest dana wzorem:

dh

dt =− μ A

A √ 2 g h

gdzie μ to wspólczynnik wydatku otworu, A0 powierzchnia otworu odpływowego, A powierzchnia lustra wody w zbiorniku i g przyspieszenie grawitacyjne. Jeśli zbiornik ma symetrię obrotową (jest walcem, kulą stożkiem, elipsoidą obrotową itp.) i ma otwór w dnie dokładnie w osi, to przyjmując r0 jako promień otworu, a R(z) jako promień zbiornika na wysokości z, otrzymamy:

dz

dt =− μ π r

π (R(z ))

√ 2 g z=−μ r

(R( z))

√ 2 g z

Rozwiązać to równanie dla trzech przypadków:

1) zbiornik jest walcem o promieniu R=1 m i wysokości 1 m, czyli:

R(z)=1 czyli dz

dt =−μ r

√ ^{2 g z}

Zadanie rozwiązać metodą Eulera w granicach od 0 do 400 sekund z krokiem 0,4 s. Wspólczyn- nik μ=0,5, a r0=0,05 m.

2) zbiornik jest odwróconym stożkiem ściętym (kształt donicy) o promieniu górnym 1 m, promie- niu dolnym 0,5 m i wysokości 1 m, czyli:

R(z)= 1

2 (1+z ) czyli dz

dt =− μ r

4 ⋅ √ 2 g z (1+ z)

Zadanie rozwiązać metodą Eulera w granicach od 0 do 200 sekund z krokiem 0,2 s. Wspólczyn- nik μ=0,5, a r0=0,05 m.

3) zbiornik jest elipsoidą obrotową o wysokości 1 m, której promień jest dany wzorem:

R(z)= z

+r

czyli dz

dt =−μ r

⋅ √ 2 g z ( z

+ r

)

Zadanie rozwiązać metodą Eulera w granicach od 0 do 50 sekund z krokiem 0,05 s. Wspólczyn- nik μ=0,5, a r0=0,05 m.