Wprowadzenie do ODEs w MATLAB-ie

Komputerowe wspomaganie w Inżynierii Materiałowej

Równaniem różniczkowym zwyczajnym rzędu pierwszego nazywamy równanie postaci

y

  gdzie

f

:

U



R

jest daną funkcją. Rozwiązaniem takiego równania nazywamy każdą funkcję która jest różniczkowalna i spełniania równość  :

R

 

R

,   ( )   dla

t

 Często rozwiązanie będziemy oznaczać także symbolem warunek zapiszemy jako

y

(

t

), więc powyższy  ( )  dla

t



Przykład

Rozważmy równanie

y t y

, Przykładowe rozwiązanie

e t

1.

Sprawdzamy przez podstawienie   (

e t

1)  1,

t t

1, Podane rozwiązanie nie jest jedyne, gdyż na przykład funkcja 1 też spełnia to równanie.

Przykład

Rozważmy równanie

y

 

y

2 .

Sprawdzamy , że funkcja     2   1 1 

t

jest rozwiązaniem: 1 1 

t

    1 (1 

t

) 2 (1  1

t

) 2 , 1 1 

t

2  1 (1 

t

) 2 , W ogólnym przypadku każda funkcja postaci 

C

1 

t

, jest rozwiązaniem tego równania.

Zagadnieniem początkowym (zagadnieniem Cauchy’ego) nazywamy następujące równanie  

y

 

y t

0 

y

0

,

gdzie

t

0

,

y

0 

R

są danymi liczbami (warunek początkowy), a

f

: ) 

R

jest daną funkcją.

Przykład

Jakie jest rozwiązanie zagadnienia Cauchy’ego  

y

 

y

(0)

t y

,  2.

Rozwiązanie ogólne równania

y’ = t y

ma postać 

Ce

1 2

t

2 .

Podstawiamy warunek początkowy y(0)=2, co daje C=2. Zatem rozwiązaniem zagadnienia Cauchy’ego jest funkcja  2

e

1 2

t

2 .

Przykład

Jakie jest rozwiązanie zagadnienia Cauchy’ego  

y



y

(0)

y

2 / 3 ,  0.

Rozwiązaniem problemu jest funkcja stale równa zero

y t

1 ( )  0.

Ale rozwiązaniem jest także funkcja

y t

2 ( )  1 27

t

3 .

Mamy zatem przykład niejednoznaczności rozwiązania!

Metoda rozdzielania zmiennych

Rozważmy równanie o rozdzielonych zmiennych

y

  Rozwiązywanie możemy symbolicznie opisać następująco

dy



y

 

dt dy

 

dy

  

C

lub 

y

0

dy



t

0

t



Przykład

Znaleźć rozwiązanie równania

y

  2 .

co daje

y dy



dt dy y

2  sin

t dt

, 2 , 

dy y

2   sin

t dt

, cos

t



C

, a więc rozwiązaniem jest  cos

t

1 

C

.

Ogólne równanie liniowe

y

   Rozwiązujemy najpierw równanie jednorodne:

y

  

dy

 

dt



dy y

   

Ce

 

dy

 

y y

   .



const

, Teraz stałą „uzmienniamy”, czyli traktujemy jak funkcję  ( )   .

Podstawiamy do wyjściowego równania i uzyskujemy elementarne równanie na C(t).

Przykład

y ty e



t

2 2 Rozwiązujemy:

y

 0,

y

  

ty

, 

Ce

 

tdt



Ce



t

2 2 .

Uzmienniane stałej

y

  

t

2 2 .



t

2 / 2 

t

2 / 2 

e



t

2 / 2

C

  sin

t



C

  

t

2 / 2 

tCe



t

2 / 2 ( )   

e



t

2 / 2 (cos ) 

t

2 2 .

Ostatecznie 

Ce



t

2 / 2 

e



t

2 / 2

Synteza bromowodoru z pierwiastków

Synteza bromowodoru z pierwiastków jest reakcją złożoną o sumarycznym równaniu H 2  Br 2  2HBr W roku 1906 wyznaczono eksperymetalnie następujące równanie kinetyczne tej reakcji

d

[HBr]

dt



k

1 [H ][Br ] 2 2 3/ 2 [Br ] 2 

k

2 [HBr] .

Czasami równanie to jest zapisywane równoważnie tak

d

[HBr]

dt



k

1 [H ][Br ] 1  2

k

2 2 1/ 2 [HBr] [Br ] 2 .

Stałe kinetyczne

k

1 oraz

k

2 zależą od warunków przebiegu reakcji (temperatura, ciśnienie itp.).

Eksperymetalnie wyznaczono, że w zwykłych warunkach k 2 ≈0,1.

Synteza bromowodoru z pierwiastków (c.d.)

Wprowadzamy oznaczenie

y

(

t

) = [HBr] oraz uwzględniamy bilans masy w równaniu H 2  Br 2 co daje dodatkowe zależności  2 0 2  2 0   1 2 1 2  2HBr   2 0   1 2 1 2

y

,

y

.

Po podstawieniu do równania kinetycznego na d[HBr]/dt otrzymamy

dy dt



k

1 ([H ] 2 0  [Br ] 2 0

y

2 0  (

k

2   0.5)

y

3/ 2 .

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.

d

[HBr]

dt



k

1 [H ][Br ] 2 2 3/ 2 [Br ] 2 

k

2 [HBr] .

4 3 2 1 0 0 10 20 30 40 t, czas 50 60 70 80 90 Gdyby kinetyka syntezy bromowodoru była analogiczna do syntezy chlorowodoru.

d

[HBr] 

k

1 [H ][Br ].

dt

5 2 1 4 3 0 0 10 20 30 40 t, czas 50 60 70 80 90

Układy równań różniczkowych zwyczajnych (ODEs)

W ogólnym przypadku możemy mieć n niewiadomych funkcji

y

1 (t),…,

y

n (t) oraz n równań:   

y

1  

y n

( , 1 ,

n

( , 1 , ,

y n

), ,

y n

), z warunkami początkowymi: 1 ( ) 0 

y

1 0 , ,

n

( ) 0 

y n

0 , gdzie liczby

t

0 ,

y

1 0 , ,

y n

0 

R

są dane.

Równania Lotki — jedna reakcja autokatalityczna

Rozważmy następującą sekwencję reakcji elementarnych:

A X



Y Y B X

2

Y

(produkcja X) (autokatalityczna produkcja Y) (rozkład Y) 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:     

dt dt

  [ ][ ], 

k X Y

2 [ ][ ] 

k Y

3 [ ].

Równania Lotki — jedna reakcja autokatalityczna

Symulacje numeryczne Modelu Lotki w MATLAB-ie dla następujących parametrów: ) 1  0.9,

k

2  9

k

3    4

A



X

0 

Y

0  4 ) 1  1.0,

k

2  1.0,

k

3 

A



X

0 

Y

0  2.0.

Czas symulacji przyjmiemy tend =5·10 5 . 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: 14 x 10 4 5 x 10 4 12 10 8 6 4 2 0 0 1 2 3 4 x 10 5 5 4.5

4 3.5

3 2.5

2 1.5

1 0.5

0 0 1 2 3 4 x 10 5 5

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

y

1 (t), a na osi OY wartości

y

2 (t).

2.5

2 1.5

1 0.5

5 x 10 4 4.5

4 3.5

3 0 0 2 4 6 8 10 12 14 x 10 4

Równania Lotki-Volterry (dwie reakcje autokatalityczne)

Jest to model podobny do modelu Lotki, ale tym razem występują dwie reakcje autokatalityczne:

A



X

2

X

(autokatalityczna produkcja X)

k X Y

2

Y

(autokatalityczna produkcja Y)

Y B

(zanik Y) Sekwencja opisuje sumaryczną reakcję A  B.

Z postaci powyższego mechanizmu możemy postulować następujący układ równań różniczkowych zwyczajnych     

dt dt

  [ ][ ], 

k X Y

2 [ ][ ] 

k Y

3

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]=

y

1 (t), [Y]=

y

2 (t), [A]=a możemy układ równań zapisać następująco:     

dy

1

dt dy

2

dt

  1 

k y y

2 1 2 1 2 

k y

3 2 .

, Układ ten ma ciekawą własność – występują w nim rozwiązania okresowe.

Dokładnej, dla każdej pary warunków początkowych 0 rozwiązania okresowymi.

y

1 (t),

y

2 (t) istnieją dla wszystkich

y t

1 

y

10 20 > 0 i są funkcjami

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.

 

y y

1 2 ( (

b cy

1 2 ) , 1



2 .

gdzie

y t

1 ( )  ofiary (np. zające),

y t

2 

y

1 

y

1 względna zmiana populacji zajęcy 

b



ay

2 ,

y

2 

y

2 względna zmiana populacji lisów 

cy

1 

d

.

Układ Lotki-Volterry – przykładowe symulacje

Do obliczeń weźmiemy następujące dane:

t k

1   4 1.5 10 ,

A



X

0

k

2    5

end

 8

Y

0

k

3   10000;  4 Wykres y1(t) i y2(t) Portret fazowy 15000 10000 5000 16000 14000 12000 10000 8000 6000 4000 0 0 0.5

1 1.5

2 2.5

3 3.5

4 4.5

x 10 5 2000 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 11000

Bruselator

Jest to teoretyczny model dla autokatalitycznej reakcji z wszystkimi etapami nieodwracalnymi i takimi samymi stałymi szybkości k 1 =k 2 =k 3 =k 4 =1.

A



X

2

X B



X

3

X



E X D

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     

dt dt



A



X

2

Y

 

B X



X

2

Y B



X

Wprowadzamy wygodniejsze oznaczenia:

y

1 

y

2 

a



b

 Parametry

a

i

b

są dodatnimi stałymi. Niewiadomymi są funkcje Równania opisujące Bruselator mają teraz postać:

y

1 =

y

1 (t),

y

2 =

y

2 (t).



y y



1 2

 

by

1

y y



1 2

y y

1 2 2 2 .

y

1 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

b



a

2



1.

Przykładowe dane do modelu Bruselator

a A



B



X

0 

Y

0 

b A



B



X

0 

Y

0  W obu przypadkach czas procesu t end = 80.

X

0 

X

0   4.0;  4.0; 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].

Przykładowe dane do modelu Bruselator



X

0 

Y



X

0 

A



B



X

0 

Y



X

0  W obu przypadkach czas procesu t end = 120.

 4.0;  4.0; 0 2.5

2 1.5

1 0 1.4

a) [B]>[B]

kr

=[A]

2

+1=2

1.2

1 0.8

20 40 60 t 20 40 60 t 80 80 100 120 100 120 3

b) [B]<[B]

kr

=[A]

2

+1=2

3 2 1 0 2 1 0 0 4 20 40 20 40 60 t 60 t 80 80 100 120 100 120

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.

Model „Oregonator”

Formalna postać Oregonatora przedstawia się następująco:

 

A Y X Y k

1

k

2

X

2

P P

  2

X X k

4

Z k

5



Y k

3

2

X A



Z

Związek z mechanizmem FKN dla reakcji BŻ jest następujący:

A X  



BrO , P

3

HBrO , Y

2

Z  2Ce(IV).

  HOBr,



Br ,

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.

      

dt dt dt

 

k A Y

1 

k A X

3 [ ][ ] 

k X Y

2 [ ][ ] 

k A X

3 [ ][ ]   

k A Y

1 

k X Y

2 [ ][ ] 

k Z

5 [ ],

k Z

5 [ ].

k X

4 2 [ ] ,  0.06;

k

1  2.1,

k

2

t X

0

end

 10  100.

 10   9

Y

0 

k

3  7  4 10 , 0

k

4     8 10 ; 7

k

5  1;

Wprowadzając oznaczenia:

y

1 

y

2 

y

3  otrzymujemy       

dy

1

dt dy dt dy

3

dt

2 

k ay

1 2    3

k ay

1

k ay

1  2

k y y

 2 1 2 

k ay

3 1 

k y

4 1 2 ,

k y y

2 1 2 

k y

5 3 .



k y

5 3 ,

a

  0.06;

k

1  2.1,

k

2

t X

0

end

 10  100.

 10   9

Y

0 

k

3  7  4 10 ,

Z

0

k

4     8 10 ; 7

k

5  1,

f

 1;

Oregonator

W ten sposób otrzymujemy słynny przykładem reakcji chemicznej z cyklem granicznym w

R

3 . Jest układ reakcji pomiędzy HBrO 2 , Br  oraz Ce(IV).

Równania kinetyki mają postać: 

y

1  

y y

2 3 77, 27(

y

2 

y

1 1 77, 27 (

y

3 0,161(

y

1 

y

3 ),

y y

1 ) 2 ),  6

y

1 

y

2 )), z przykładowymi warunkami początkowymi

y

1

y

2 (0)  2,

y

3 (0)  3.

Następujący układ reakcji chemicznych został zaproponowany przez H. Robertsona w 1966 roku:

A B



B

 7

B C



B

(powolna) (bardzo szybka) (szybka) Prowadzi on do układu równań  0,04 

y

1  

y

2

y

3   0,04 3 10 7

y

1

y

1

y

2 2   10 10 4 4

y y

2

y y

2 3 3   7

y

2 2 Typowe testy przeprowadza się dla następujących warunków początkowych:

y

1

y

2 (0)  0,

y

3 (0)  0, oraz dla czasów t end = 10, 10 2 , 10 3 , ..., 10 11 .

Równania drugiego rzędu

Rozważmy równanie ruchu wahadła matematycznego (  (t)=kąt wychylenia):    

g

Wprowadzając  2 =

g

/

l

otrzymujemy 2 sin



 0.

Zamieniamy na układ wprowadzając:

y

1 



,

y

2 



 .

 

y y

1 2  

y

2 ,



2 sin

y

1 .

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.

-5 -10 -15 0 15 10 5 0 5 10 15 20 czas 25 30 35 40

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: 

gu

(1 

u

2 / 3

U

2 ).

a

Stosując prawa Kirchoffa do tego układu i wprowadzając przeskalowane nowe zmienne otrzymujemy standardowe równanie różniczkowe zwyczajne: 2

d y

  (1 

dt

2

y

2 )

dy dt

0, 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:     

dy

1

dt dy dt

2  

y

2 ,  (1 

y

1 2 )

y

2 

y

1 .

Ewolucja czasowa i portrety fazowe

Przykładowe rozwiązania

y

1 (

t

),

y

2 (

t

) oraz portret fazowy dla układu Lotki-Volterry.

  

y y

1 2  

y

1 (0) (

a

(

cy

1  

by y

 20, 2

y

) 2 1 2 (0)

a

 1,

b

 0.01,

c

  20 0.02,

d

  1

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ą

y

  lub postać uwikłaną liniowo gdzie

y



y

1     ,

M

    

m

11

m n

1

m

1

n m nn

    .

Przykład

Numeryczne całkowanie pojedynczego równania różniczkowego zwyczajnego Rozważmy zagadnienie początkowe

y

 

y y

(0)  (2

y

0  .

y

), 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.

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].

Drugi sposób definiowanie prawej strony wpisaniu definicji funkcji do odrębnego pliku .m i wykorzystaniu go w procedurze

ode23

(lub

ode45

zagadnienie początkowego polega na ). 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:

Ewolucja populacji z uwzględnieniem dyfuzji – model Fishera

Założenia:

1) Występuje tylko jeden gatunek.

2) Problem jest jednowymiarowy (stworzenia żyją na linii prostej) ograniczony do skończonego odcinka.

3) Gdyby nie było tłoku, to układ rozwijał by się zgdonie z modelem logistycznym.

4) Na skutek tłoku osobniki mają tendencję do migracji z obszarów o dużym zagęszczeniu do obszarów o mniejszym zagęszczeniu.

5) 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 

c

   

t c x



D

 2

c



x

2 



 0, 

c



x

 ),  0.

Matematyczny opis modelu Fishera

Niech

c

(

x

,

t

) oznacza gęstość osobników w punkcie Wtedy możemy rozważać model opisany równaniem

x

w momencie czasu

t

.

z warunkami brzgowymi 

c



t

  2

c D



x

2 





c



x

 0, 

c



x

 ),  0.

Ponadto zakładamy, że znamy początkowy rozkład populacji

c

0 (

x

):  0 ( ), 0 .

Dyskretyzacja równania Fishera



c



t x t k



D

 2

c



x

2

x t k





( , )(

k

Stąd mamy

dc k dt



D c k

 1

t



c t k

( 

x

) 2 

c k

 1 



k

( , )), dla

k k

 0,1,

k

( )), dla

k

 1, ,

n

 1.

Zatem

c t

1 ( )  

x c t

0 ( )  0,

c n

 1  

x n

 0

c t

1 ( ) 

c t

0 ( ),

c n

 1 

n

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:         

dc

1

dt dc k dt dc n dt

 1  

D



D D c

2 (  

x c

) 2 1

c k

 1 

c n

 (  1  2

c k



( 

x

)  

x

) 2 2  1

c n

 2 (

c k

 1  





c

1



n

),  1

k

( 

k

), 

n

 1 ), który rozwiązujemy z warunkami początkowymi

k

 1

k

 2,

k

1 ,

n

 2

c k

(0)  0 (

k

),

k

 1, ,

n

 1.

Funkcja c 0 (x) jest dana i oznacza początkowy rozkład populacji.

W powyższym równaniu nie występują wartości w węzłach x 0 mocy warunków brzegowych mamy c 0 =c 1 oraz c n =c n-1 .

i x n , gdyż na

Ciągłe układy dynamiczne w środowisku Mathematica

W środowisku Mathematica do numerycznego rozwiązywania dynamicznych służy funkcja NDSolve. Podstawowa składnia ma postać: układów NDSolve [{równania, warunki}, {t, t

min , t max

}]

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}]