Course 8. Ordinary Differential Equations (ODE)

Clear["Global`*"]

Analytic Solutions

Single ODEs of n-order and General solution

Example 1. 2nd order linear ODE with constant coefficients

deq1=y''[t]+ 2 y[t]==Cos[3 t]

DSolve[deq1,y,t]

Example 2. 2nd order linear ODE with non constant coefficients

deq2=y''[t]+  t y[t]==0
DSolve[deq2,y,t]

Example 3. 2nd order linear ODE with non constant coefficients (2nd case)

deq3=y''[t]+  t y[t]==Sin[t]
DSolve[deq3,y,t]

Example 4. 2nd order autonomous non-linear ODE

deq4=y''[t]+ 2 y[t]-y[t]^3==0
DSolve[deq4,y,t]

Example 5. 2nd order non-autonomous, non-linear ODE

deq5=y''[t]+ 2 y[t]-y[t]^3==Sin[t]
DSolve[deq5,y,t]

Initial Value Problem

General solution

deq1=y''[t]+ 2 y[t]==Cos[3 t];
IVP={deq1, y[0]==y0, y'[0]==v0}
DSolve[IVP,y,t]

Partial solutions

IVP={deq1, y[0]==1, y'[0]==0}
sol=DSolve[IVP,y,t]

y1=y[t]/.sol[[1]]   //Simplify

Plot[y1,{t,0,10}]

IVP={deq1, y[0]==0, y'[0]==0.5}
sol=DSolve[IVP,y,t]
y2=y[t]/.sol[[1]]   //Simplify
Plot[y2,{t,0,10}]

Systems of Equations

General Solution

deq1=y'[t]+3 x[t]==0
deq2= x'[t]-y[t]==t

sol=DSolve[{deq1,deq2},{x,y},t]

xx=x[t]/.sol[[1]] //Simplify
yy=y[t]/.sol[[1]] //Simplify

Partial Solution

sol=DSolve[{deq1,deq2, x[0]==1, y[0]==0},{x,y},t]
xx=x[t]/.sol[[1]] //Simplify
yy=y[t]/.sol[[1]] //Simplify

Plot[{xx,yy},{t,0,10},PlotStyle->{Red, Blue}]

ParametricPlot[{xx,yy},{t,0,10}, Frame->True, AspectRatio->1]

Numerical Solutions

Single ODE

deq=y''[t]- 2 y[t]+y[t]^3==Sin[t]
DSolve[deq,y,t]

Partial Numerical Solution

y0=1; dy0=0;   (* initial conditions e.g. at t=0 *)
tmin=-10; tmax=10;
Print["Initial conditions y(0)=",y0,"  dy(0)=",dy0,",    Integration interval from ",tmin, " to ", tmax];

sol=NDSolve[{deq,y[0]==y0, y'[0]==dy0},y,{t,tmin,tmax}]

yy=y[t]/.sol[[1]]

Plot[yy,{t,-10,10}]

numsol=Table[{t,yy},{t,-10,10,0.1}]
ListPlot[numsol]

Systems of Equations

deq1=x'[t]==y[t]+0.5z[t]
deq2=y'[t]==-x[t]
deq3=z'[t]==-z[t]+x[t]^2
x0=1; y0=0; z0=-1;
tmin=0; tmax=100;

sol=NDSolve[{deq1,deq2,deq3,x[0]==x0, y[0]==y0, z[0]==z0},{x,y,z},{t,tmin,tmax}];
xx=x[t]/.sol[[1]]
yy=y[t]/.sol[[1]]
zz= z[t]/.sol[[1]]

Plot[xx,{t,tmin,tmax}, PlotStyle->Red]
Plot[yy,{t,tmin,tmax}, PlotStyle->Green]
Plot[zz,{t,tmin,tmax}, PlotStyle->Blue]

ParametricPlot3D[{xx,yy,zz},{t,tmin,tmax}]

Working Options

deq=y''[t]- 2 y[t]+y[t]^3==(1/4)Sin[t]

tmax=100; (* tmax=1000 *)
sol=NDSolve[{deq,y[0]==1, y'[0]==0},y,{t,0,tmax},MaxSteps->10000];
yy=y[t]/.sol2[[1]]
Plot[yy,{t,0,tmax}]

Solution with default accuracy parameters

tmax=1000; 
sol=NDSolve[{deq,y[0]==1, y'[0]==0},y,{t,0,tmax}];
yy=y[t]/.sol[[1]]
Plot[yy,{t,0,tmax}]
graph1=Plot[yy,{t,300,500}, PlotStyle->Green]

Low Accuracy Solution

tmax=1000;
sol=NDSolve[{deq,y[0]==1, y'[0]==0},y,{t,0,tmax},  AccuracyGoal->5, PrecisionGoal->5];
yy=y[t]/.sol[[1]]
Plot[yy,{t,0,tmax}]
Plot[yy,{t,300,500},PlotStyle->Red]
graph2=Plot[yy,{t,tmax-200,tmax},PlotStyle->Red]

High Precision Solution

tmax=1000; 
sol=NDSolve[{deq,y[0]==1, y'[0]==0},y,{t,0,tmax}, StartingStepSize-> 0.001,AccuracyGoal->13, PrecisionGoal->13];
yy=y[t]/.sol[[1]]
graph3=Plot[yy,{t,tmax-200,tmax},PlotStyle->Blue]

Show[graph1,graph2, graph3]

tmax=1000; 
sol=NDSolve[{deq,y[0]==1, y'[0]==0},y,{t,0,tmax}, Method->"ExplicitRungeKutta"];
yy=y[t]/.sol[[1]]
graph4=Plot[yy,{t,tmax-200,tmax}, PlotStyle->Brown]