Wprowadzenie do ODEs w MATLAB-ie

1. Wprowadzenie do ODEs w MATLAB-ie Komputerowe wspomaganie w Inżynierii Materiałowej

2. Równaniem różniczkowym zwyczajnym rzędu pierwszego nazywamy równanie postaci gdzie jest daną funkcją. Rozwiązaniem takiego równania nazywamy każdą funkcję która jest różniczkowalna i spełniania równość Często rozwiązanie będziemy oznaczać także symbolem y(t), więc powyższy warunek zapiszemy jako

3. Przykład Rozważmy równanie Przykładowe rozwiązanie Sprawdzamy przez podstawienie Podane rozwiązanie nie jest jedyne, gdyż na przykład funkcja też spełnia to równanie.

4. Przykład Rozważmy równanie Sprawdzamy, że funkcja jest rozwiązaniem: W ogólnym przypadku każda funkcja postaci jest rozwiązaniem tego równania.

5. Zagadnieniem początkowym (zagadnieniem Cauchy’ego) nazywamy następujące równanie są danymi liczbami (warunek początkowy), a gdzie jest daną funkcją.

6. Przykład Jakie jest rozwiązanie zagadnienia Cauchy’ego Rozwiązanie ogólne równania y’ = t y ma postać Podstawiamy warunek początkowy y(0)=2, co daje C=2. Zatem rozwiązaniem zagadnienia Cauchy’ego jest funkcja

7. Przykład Jakie jest rozwiązanie zagadnienia Cauchy’ego Rozwiązaniem problemu jest funkcja stale równa zero Ale rozwiązaniem jest także funkcja Mamy zatem przykład niejednoznaczności rozwiązania!

8. Metoda rozdzielania zmiennych Rozważmy równanie o rozdzielonych zmiennych Rozwiązywanie możemy symbolicznie opisać następująco

9. Przykład Znaleźć rozwiązanie równania co daje a więc rozwiązaniem jest

10. Ogólne równanie liniowe Rozwiązujemy najpierw równanie jednorodne: Teraz stałą „uzmienniamy”, czyli traktujemy jak funkcję Podstawiamy do wyjściowego równania i uzyskujemy elementarne równanie na C(t).

11. Przykład Rozwiązujemy: Uzmienniane stałej Ostatecznie

12. Synteza bromowodoru z pierwiastków Synteza bromowodoru z pierwiastków jest reakcją złożoną o sumarycznym równaniu W roku 1906 wyznaczono eksperymetalnie następujące równanie kinetyczne tej reakcji Czasami równanie to jest zapisywane równoważnie tak Stałe kinetyczne k1 oraz k2 zależą od warunków przebiegu reakcji (temperatura, ciśnienie itp.). Eksperymetalnie wyznaczono, że w zwykłych warunkach k2≈0,1.

13. Synteza bromowodoru z pierwiastków (c.d.) Wprowadzamy oznaczenie y(t) = [HBr] oraz uwzględniamy bilans masy w równaniu co daje dodatkowe zależności Po podstawieniu do równania kinetycznego na d[HBr]/dt otrzymamy

14. Synteza bromowodoru z pierwiastków (c.d.) Przeprowadzając symulację podanego układu dynamicznego możemy precyzyjnie przewidzieć ewolucję stężenia składników – a w szczególności przewidzieć czas trwania reakcji. Prawdziwa kinetyka syntezy bromowodoru. Gdyby kinetyka syntezy bromowodoru była analogiczna do syntezy chlorowodoru.

15. Układy równań różniczkowych zwyczajnych (ODEs) W ogólnym przypadku możemy mieć n niewiadomych funkcji y1(t),…,yn(t) oraz n równań: z warunkami początkowymi: gdzie liczby są dane.

16. Równania Lotki — jedna reakcja autokatalityczna Rozważmy następującą sekwencję reakcji elementarnych: Powyższy mechanizm opisuje ostatecznie sumaryczną reakcję A  B. Z postaci tego mechanizmu możemy postulować następujący układ równań różniczkowych zwyczajnych:

17. Równania Lotki — jedna reakcja autokatalityczna Symulacje numeryczne Modelu Lotki w MATLAB-ie dla następujących parametrów: Czas symulacji przyjmiemy tend =5·105. Zastosowanie standardowej procedury ode45 (implementujacej metodę Rungego-Kutty 4-tego rzędu) z domyślnymi ustwieniami tolerancji błędów dla przypadku a) daje wyniki:

18. Równania Lotki — przykładowy portret fazowy Poniżej jest przedstawiony portret fazowy układu Lotki dla danych z punktu a). Portret fazowy oznacza, że rysujemy wyniki obliczeń w układzie y1-y2. Tzn. na osi OX odkładane są wartości y1(t), a na osi OY wartości y2(t).

19. Równania Lotki-Volterry (dwie reakcje autokatalityczne) Jest to model podobny do modelu Lotki, ale tym razem występują dwie reakcje autokatalityczne: Sekwencja opisuje sumaryczną reakcję AB. Z postaci powyższego mechanizmu możemy postulować następujący układ równań różniczkowych zwyczajnych

20. Równania Lotki-Volterry (dwie reakcje autokatalityczne – c.d.) W układzie reakcji Lotki-Volterry zakładamy, że stężenie reagenta A jest stałe: [A]=const. Wprowadzając wygodne oznaczenia: [X]=y1(t), [Y]=y2(t), [A]=a możemy układ równań zapisać następująco: Układ ten ma ciekawą własność – występują w nim rozwiązania okresowe. Dokładnej, dla każdej pary warunków początkowych y1(0)=y10 > 0, y2(0)=y20 > 0 rozwiązania y1(t), y2(t) istnieją dla wszystkich t 0 i są funkcjami okresowymi.

21. Układ Lotki-Volterry jako prosty model drapieżnik-ofiara Ten sam układ równań może być wykorzystany do opisu prostego modelu interakcji pomiędzy dwoma populacjami: ofiar i drapieżników. gdzie

22. Układ Lotki-Volterry – przykładowe symulacje Do obliczeń weźmiemy następujące dane:

23. Bruselator Jest to teoretyczny model dla autokatalitycznej reakcji z wszystkimi etapami nieodwracalnymi i takimi samymi stałymi szybkości k1=k2=k3=k4=1. Procesem sumarycznym jest: A+BD+E. Powyższy mechanizm prowadzi do następującego układu, gdy szybkość reakcji jest określona przez postacie reakcji

24. Wprowadzamy wygodniejsze oznaczenia: Parametry a i b są dodatnimi stałymi. Niewiadomymi są funkcje y1=y1(t), y2=y2(t). Równania opisujące Bruselator mają teraz postać: Powyższy układ równań generuje rozwiązania, których jakościowy charakter może istotnie się różnić w zależności od wzajemnej relacji parametrów a i b. W szczególności punkt stacjonarny tego układu staje się niestabilny gdy

25. Przykładowe dane do modelu Bruselator W obu przypadkach czas procesu tend = 80. Jeśli zmieniając stężenie [B] osiągniemy wartość krtytyczną [B]kr=[A]2+1, to następuje tzw. bifurkacja Hopfa. Dotychczasowy pojedynczy stan stacjonarny traci stabilność – w zakresie stężeń [B] > [A]2 + 1 obserwujemy zupełnie inne zachowanie – stabilne oscylacje stężeń [X] i [Y].

26. Przykładowe dane do modelu Bruselator W obu przypadkach czas procesu tend = 120. a) [B]>[B]kr=[A]2+1=2 b) [B]<[B]kr=[A]2+1=2

27. Oregonator Oregonator jest najprostszym realistycznym modelem opisującym dynamikę reakcji oscylacyjnych typu Biełousowa-Żabotyńskiego (BŻ). Model ten jest jednak na tyle bogaty, że pozwala symulować wiele ciekawych i złożonych zachowań takich, jak rozwiązania okresowe, cykle graniczne czy bifurkacje. Twórcami tego modelu są R. J. Field i R. M. Noyes z Uniwersytetu stanu Oregon, USA. W latach 50-tych XX wieku B. Biełousow, chemik pracujący w tajnych laboratoriach wojskowych Związku Sowieckiego, badając m.in. wpływ gazów trujących i promieniowania na organizm ludzki zauważył, że mieszanina bromianu potasu, kwasu cytrynowego i jonów ceru wykazuje okresowe zmiany koloru. Oznaczało to, że w reakcji nie następowała monotoniczna zmiana stężenia reagentów jak to ma miejsce w większości „typowych” reakcji chemicznych – tylko były to zmiany okresowe. Około dekady później A. Żabotyński dokładniej zbadał tę reakcję (zmienił też kwas cytrynowy na kwas malonowy) i zasugerował, że przyczyną zmian koloru mogą być jonu ceru na różnych stopniach utlenienia i rozpropagował temat na konferencji w Pradze, 1968. Następnie w 1972 R.J. Field, R. M. Noyes i E. Körös zaproponowali pierwszy zbliżony do realnego mechanizm reakcji BŻ. W literaturze określa się go jako mechanizm FKN. W oparciu o ten mechanizm reakcji powstał modelowy ciąg reakcji z trzema zmiennymi stężeniami – Oregonator. Określenie to oznacza zarówno mechanizm reakcji jak i układ równań różniczkowych, który je opisuje.

28. Model „Oregonator” Formalna postać Oregonatora przedstawia się następująco: Związek z mechanizmem FKN dla reakcji BŻ jest następujący: Przy założeniu, że [A]=const zmiennymi zależnymi od czasu stają się stężenia [X], [Y], [Z]. Stężenie [P] dotyczy tylko produktu, który nie występuje jako substrat.

30. otrzymujemy Wprowadzając oznaczenia:

31. Oregonator W ten sposób otrzymujemy słynny przykładem reakcji chemicznej z cyklem granicznym w R3. Jest układ reakcji pomiędzy HBrO2, Br oraz Ce(IV). Równania kinetyki mają postać: z przykładowymi warunkami początkowymi

32. Następujący układ reakcji chemicznych został zaproponowany przez H. Robertsona w 1966 roku: Prowadzi on do układu równań Typowe testy przeprowadza się dla następujących warunków początkowych: oraz dla czasów tend = 10, 102, 103, ..., 1011.

33. Równania drugiego rzędu Rozważmy równanie ruchu wahadła matematycznego ((t)=kąt wychylenia): Wprowadzając w2=g/l otrzymujemy Zamieniamy na układ wprowadzając:

34. Prosta implementacja w MATLAB-ie wraz z wykresem pokazującym zależność kata od czasu. Dla wygody kąt jest na wykresie podawany w stopniach. Rozwiązanie jest okresowe. Dla małych wychyleń jest to prawie sinusoida. W ogólnym przypadku tak nie jest. W szczególności okres drgań zależy do amplitudy.

35. Równanie van der Pola Równanie to zostało pierwotnie zaproponowane do opisu pewnego układu elektronicznego składającego się z obwodu rezonansowego RLC sprzężonego indukcyjnie z cewką. Przez cewkę płynie prąd zależny do napięcia nieliniowo. Dokładniej była to zależność trzeciego stopnia: Stosując prawa Kirchoffa do tego układu i wprowadzając przeskalowane nowe zmienne otrzymujemy standardowe równanie różniczkowe zwyczajne: gdzie y=y(t) jest funkcją czasu, a stała   0 charakteryzuje wielkość tłumienia. Równoważny układ równań pierwszego rzędu to:

36. Ewolucja czasowa i portrety fazowe Przykładowe rozwiązania y1(t), y2(t) oraz portret fazowy dla układu Lotki-Volterry.

37. Numeryczne rozwiązywanie problemu początkowego dla równań różniczkowych w MATLAB-ie Środowisko MATLAB oferuje wiele numerycznych procedur rozwiązywania układów równań różniczkowych zwyczajnych. Do większości problemów – w szczególności takich, które nie są sztywne – wystarczy używać funkcji ode23 lub ode45. Obie funkcje implementują odpowiednie wersje metody Rungego-Kutty. W pierwszym przypadku są wykorzystywane metody RK drugiego i trzeciego rzędu, a w drugim – metody czwartego i piątego rzędu. W obu przypadkach stosowany jest zmienny krok całkowania (adaptacja kroku w zależności od lokalnie oszacowanego błędu). Dla problemów sztywnych mozna stosować np. Procedurę ode23s. Układy równań mogą mieć postać jawną lub postać uwikłaną liniowo gdzie

38. Przykład Numeryczne całkowanie pojedynczego równania różniczkowego zwyczajnego Rozważmy zagadnienie początkowe Aby móc użyć procedurę ode23 lub ode45 należy zdefiniować prawą stronę równania, czyli f(t,y)=y(2-y). Można to zrobić w postaci tzw. funkcji anonimowej lub zawierając definicję funkcji w osobnym pliku. Przykład z funkcją anonimową i komendą w wierszu interpretera poleceń podany jest poniżej. Przedział całkowania to [0, 10], a warunek początkowy to y(0)=1.2.

39. Po wykonaniu poleceń w takiej formie jak widzimy na obrazie, MATLAB wyświetli rozwiązanie na wykresie, czyli zbiór punktów (t, y(t)) dla t[0, 10].

40. Drugi sposób definiowanie prawej strony zagadnienie początkowego polega na wpisaniu definicji funkcji do odrębnego pliku .m i wykorzystaniu go w procedurze ode23 (lub ode45). W oknie MATLAB-a klikamy ikonkę New Script i wpisujemy tam definicję funkcji, a następnie zapisujemy w bieżącym (lub wybranym przez nas) katalogu. Domyślna nazwa pliku jest taka sama jak nazwa funkcji plus rozszerzenie .m. W naszym przypadku będzie to f.m. Po zapisaniu prawej strony układu w pliku f.m możemy uruchomić procedurę numerycznego całkowania równania:

41. Ewolucja populacji z uwzględnieniem dyfuzji – model Fishera Założenia: Występuje tylko jeden gatunek. Problem jest jednowymiarowy (stworzenia żyją na linii prostej) ograniczony do skończonego odcinka. Gdyby nie było tłoku, to układ rozwijał by się zgdonie z modelem logistycznym. Na skutek tłoku osobniki mają tendencję do migracji z obszarów o dużym zagęszczeniu do obszarów o mniejszym zagęszczeniu. Tempo migracji jest proporcjonalne do gradientu gęstości populacji. Współczynnik tego tempa nazywamy współczynnikiem dyfuzji i oznaczamy przez D. Niech c(x,t) oznacza gęstość osobników w punkcie x w momencie czasu t. Wtedy możemy rozważać model opisany równaniem

42. Matematyczny opis modelu Fishera Niech c(x,t) oznacza gęstość osobników w punkcie x w momencie czasu t. Wtedy możemy rozważać model opisany równaniem z warunkami brzgowymi Ponadto zakładamy, że znamy początkowy rozkład populacji c0(x):

43. Dyskretyzacja równania Fishera Stąd mamy Zatem

44. Uwzględniając warunki brzegowe otrzymujemy ostatecznie następujący układ równań różniczkowych zwyczajnych (metoda linii) aproksymujący wyjściowe równanie Fishera: który rozwiązujemy z warunkami początkowymi Funkcja c0(x) jest dana i oznacza początkowy rozkład populacji. W powyższym równaniu nie występują wartości w węzłach x0 i xn, gdyż na mocy warunków brzegowych mamy c0=c1oraz cn=cn-1.

45. Ciągłe układy dynamiczne w środowisku Mathematica W środowisku Mathematica do numerycznego rozwiązywania układów dynamicznych służy funkcja NDSolve. Podstawowa składnia ma postać: NDSolve[{równania, warunki}, {t, tmin, tmax}] Przykład 1: sol=NDSolve[{y'[t]==y[t]Cos[t+y[t]],y[0]==1},y,{t,0,30}] Plot[Evaluate[y[t]/.sol], {t, 0, 30}, PlotRange->All] Przykład 2: sol=NDSolve[{x'[t]==-y[t]-x[t]^2,y'[t]==2x[t]-y[t]^3,x[0]==y[0]==1},{x,y},{t,20}] ParametricPlot[Evaluate[{x[t], y[t]} /. sol], t, 0, 20}]