Department of Applied Mathematics at the University of Colorado at Boulder
University of Colorado at Boulder Search A to Z Campus Map University of Colorado at BoulderCU Search Links
You are here: Resources » Computing Resources » Lab Policies » Solve ODE
Print this page

Solve ODE

15: Solve & Graph an ODE

Given an ODE (ordinary differential equation), solve it analytically if possible, solve it with a specific initial condition (analytically or numerically); make a graph of its ``direction field'' with little arrows, and plot a specific solution to the ODE + initial condition. For these examples use this ODE for y(t):

y' + 3y = t + e-2t

1. If possible, solve the ODE generally
2. Solve the ODE given initial condition y(-1)=-0.4, and graph the solution
3. graph the direction field of the ODE, for -1.5<t<5 and -1<y<3
4. superimpose the solution (2) on the vector field (3)

Mathematica:

(1)
diffeq = y'[t]+3y[t]==t+Exp[-2t]
solnA=DSolve[ diffeq, y[t], t ] (2)
init = y[-1]==-0.4
solnB=DSolve[, y[t], t]
gr1 = Plot[ y[t] /. solnB[[1,1]],
    , PlotRange->]; (3)
Needs["Graphics`PlotField`"]
gr2 = PlotVectorField[ ,
    , ,
    ScaleFunction->(1&), ScaleFactor->0.2,
    PlotPoints->, Axes->True]; (4)
Show[gr1,gr2, PlotRange->];

 

 

 

 

 Matlab:

(1) (use Mathematica or Maple)

(2)
    make this file "RHS.m":
function f=RHS(t,y)
f=-3*y+t+exp(-2*t);

    now run Matlab:
dt=0.1
timegrid=-1:dt:5;
[T,Y]=ode45(@RHS,timegrid,-0.4);
figure(1)
plot(T,Y);

(3)
dt=0.25;
dy=0.25;
scale=3;
timegrid=-1.5:dt:5;
ygrid=-1:dy:3;
[Tmesh,Ymesh]=meshgrid(timegrid,ygrid);
Dy=RHS(Tmesh,Ymesh)*dt;
Dt=dt*ones(length(ygrid),length(timegrid));
figure(2)
quiver(Tmesh,Ymesh,Dt,Dy,scale);
axis([-1.5 5 -1 3])

(4)
figure(3)
plot(T,Y);
axis([-1.5 5 -1 3])
hold on;
quiver(Tmesh,Ymesh,Dt,Dy,scale);
hold off;

Maple:

 

(1)
dsolve(D(y)(t) + 3*y(t) = t + exp(-2*t), y(t));

(2)
dsolve( , y(t));

with(DEtools):
with(plots):
sol := dsolve( , y(t), type=numeric);
gr1 := odeplot(sol,[t,y(t)],-1.5..5);

(3)
gr2 := dfieldplot( D(y)(t)+3*y(t)=t+exp(-2*t), y(t), t=-1.5..5, y=-1..3 );

(4)
display( );

IDL:

(1) (use Mathematica or Maple)

(2) (numerically, for 0<t<5)
Function dydt, t, y
    return, t + exp(-2*t) - 3*y
End
t0=-1.0 & y0=-0.4 & dydt0=dydt(t0,y0)
h=0.1 & y=[y0] & t=[t0]
For i=1,55 do begin
    y0 = RK4(y0,dydt0,t0,h,'dydt')
    t0 = t0+h
    t=[t,t0] & y=[y,y0]
    dydt0 = dydt(t0,y0)
EndFor
plot, t, y
(3)
tt = (findgen(66)/10-1.5) # (fltarr(41)+1)
yy = (fltarr(66)+1) # (findgen(41)/10-1.0)
dy = dydt(tt,yy)
dt = 0*dy + 1
dl = sqrt(dt*dt + dy*dy)
dy = 0.1 * dy / dl
dt = 0.1 * dt / dl
plot, tt, yy, xstyle=1, ystyle=1, /nodata
arrow, tt, yy, tt+dt, yy+dy, hsize=-0.3, /data (4)
oplot, t, y, thick=3

 

 


^ Top

Previous page: xyN.dat
Next page: Graphing 2D data