previous: define/plot functions
next: FFT

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[{diffeq, init}, y[t], t]
gr1 = Plot[ y[t] /. solnB[[1,1]],
    {t,-1.5,5}, PlotRange->{-1,3}];

(3)
Needs["Graphics`PlotField`"]
gr2 = PlotVectorField[ {1, t+Exp[-2t]-3y},
    {t,-1.5,5}, {y,-1,3},
    ScaleFunction->(1&), ScaleFactor->0.2,
    PlotPoints->{25,15}, Axes->True];

(4)
Show[gr1,gr2, PlotRange->{-1,3}];

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( {D(y)(t) + 3.0*y(t) = t + exp(-2.0*t), y(-1.0)=(-0.4)}, y(t));

with(DEtools):
with(plots):
sol := dsolve( {D(y)(t) + 3.0*y(t) = t + exp(-2.0*t), y(-1.0)=(-0.4)}, 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( {gr1,gr2} );

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