Integrarea ecuatiilor diferentiale ordinare

Integrare ODE grad 1 cu ode()

Exemplul 1

// ---------- O simpla ecuatie de integrat folosind functia ode
//ecuatia de integrat
// dy/dt=y^2-y sin(t)+cos(t), y(0)=0
//care are solutia y=sin(t)
//declararea functiei f
function ydot=f(t,y),ydot=y^2-y*sin(t)+cos(t),endfunction
pas=%pi/50;y0=0;t0=0;  //alegerea pasului si conditiilor initiale
t=0:pas:2*%pi;         //generarea vectorului de timp
y=ode(y0,t0,t,f);      //integrarea ecuatiei
yy=sin(t);             //calculul solutiei analitice
dy=y-yy;        //calculul diferentei dintre solutia numerica si cea analitica
clf();                 //sterg fereastra grafica
plot2d(t,dy,style=color("black"));   //desenez diferenta dintre solutia integrata si cea analitica

//aleg un pas mai fin
pas=%pi/100;y0=0;t0=0;  //alegerea pasului si conditiilor initiale
t=0:pas:2*%pi;         //generarea vectorului de timp
y=ode(y0,t0,t,f);      //integrarea ecuatiei
yy=sin(t);             //calculul solutiei analitice
dy=y-yy;        //calculul diferentei dintre solutia numerica si cea analitica
plot2d(t,dy,style=color("red"));   //desenez diferenta dintre solutia integrata si cea analitica

//aleg un pas foarte fin
pas=%pi/1000;y0=0;t0=0;  //alegerea pasului si conditiilor initiale
t=0:pas:2*%pi;         //generarea vectorului de timp
y=ode(y0,t0,t,f);      //integrarea ecuatiei
yy=sin(t);             //calculul solutiei analitice
dy=y-yy;        //calculul diferentei dintre solutia numerica si cea analitica
plot2d(t,dy,style=color("blue"));   //desenez diferenta dintre solutia integrata si cea analitica
In urma executiei acestui cod se obtine o fereastra grafica precum cea de mai jos.

Daca in programul de mai sus schimbam liniile

t=0:pas:2*%pi;

cu

t=0:pas:nper*2*%pi; //unde nper l-am declarat mai sus ca fiind nper=20

obtinem figura

Aceiasi pasi dar integrare pe 20 de perioade

ceea ce arata ca in general eroarea depinde de pas dar si de timpul pe care se face integrarea

Exemplul 2

Codul

// ---------- O simpla ecuatie de integrat folosind functia ode
//ecuatia de integrat
// dy/dt=y^2-y sin(t)+cos(t), y(0)=0
//care are solutia y=sin(t)
//declararea functiei f
function ydot=f(t,y),ydot=y^2-y*sin(t)+cos(t),endfunction
nper=100;
tip="rkf";//"adams","rk","rkf",
pas=%pi/500;y0=0;t0=0;  //alegerea pasului si conditiilor initiale
t=0:pas:nper*2*%pi;         //generarea vectorului de timp
y=ode(tip,y0,t0,t,f);      //integrarea ecuatiei
yy=sin(t);             //calculul solutiei analitice
dy=y-yy;        //calculul diferentei dintre solutia numerica si cea analitica
clf();                 //sterg fereastra grafica
plot2d(t,dy,style=color("black"));   //desenez diferenta dintre solutia integrata si cea analitica


pas=%pi/400;y0=0;t0=0;  //alegerea pasului si conditiilor initiale
t=0:pas:nper*2*%pi;         //generarea vectorului de timp
y=ode(tip,y0,t0,t,f);      //integrarea ecuatiei
yy=sin(t);             //calculul solutiei analitice
dy=y-yy;        //calculul diferentei dintre solutia numerica si cea analitica
plot2d(t,dy,style=color("red"));   //desenez diferenta dintre solutia integrata si cea analitica

genereaza figura

Eroarea absoluta pe 100 de perioade, integrare tip "rkf", cu 500 si respectiv 400 pasi pe perioada

Integrare ODE grad 2 cu ode()

Aici vom integra ecuatia Mathieu : d2y/dt2=-(a-2*q*cos(2*t))*y.

Integrarea unei ecuatii de gradul 2 de forma

d2y/dt2=f(t,y,dy/dt)

este de fapt echivalenta cu integrarea unui sistem de doua ecuatii de gradul 1:

  1. dy/dt=v
  2. dv/dt=f(t,y,v)

Ecuatia Mathieu rezolvata cu functia ode trebuie sa furnizeze un vector bidimensional de pozitie [y0,v0] la t0 care reprezinta pozitia la momentul initial.

Exemplul 1

// **************************************
// A 2nd order linear diff eq:
// y'' + (a-2*q*cos(2*t))*y = 0.
// **************************************

// Define the equation
function [yprim]=diffeq(t,y)
  yprim(1)=y(2);
  yprim(2)=-(a-2*q*cos(2*t))*y(1);
endfunction
//Parameters
a=0;q=0.92;

// Initial conditions
per=%pi;
t0=0;step=per/50;
nper=15;
tfin=nper*per;

//pozition and velocity
Yini=1;Vini=0;
y(1)=Yini; y(2)=Vini;  
y0=[y(1); y(2)];

// Define time axis
 t=t0:step:tfin; 

// Here the ode gives the solution y,
// row 1: y(1) 501 values
// row 2: y(2) = yprim(1) = dy(1)/dt 501 values
y=ode(y0,t0,t,diffeq);
//Here we compute cos
z=cos(2*t);
// ------------------------------------------------------
//Deseneaza solutia si derivata in subplot

 //aici se traseaza in aceeasi fereastra mai multe axe folosind subplot
subplot(3,1,1) // deseneaza
plot2d(t,y(1,:));xgrid()
xtitle("Pozitia","time","y");

subplot(3,1,2)
plot2d(t,y(2,:));xgrid()
xtitle("Derivata intaia: viteza","time","Vy");

subplot(3,1,3)
plot2d(t,z);xgrid()
xtitle("cos(2*t)","time","cos");

care produce graficul

 

Cornel Mironel Niculae @fizica.unibuc.ro
2008-2009

Last updated: 22 Mar 2009