MS-C1650 - Numeerinen analyysi
Lecture 7: Numerical Methods for ODEs

Harri Hakula, 2019

Setting

To begin, consider an initial value problem for a linear first-order ODE.

This is a linear first-order ODE:

linearequation = y'[t]-3*t*y[t]==1;

Notice that the general solution is a linear function of the arbitrary constant (RefLink[C, paclet:ref/C])[1]:

generalsolution = DSolve[linearequation, y[t],t]

This finds a particular solution for the initial condition y[0]==4:

particularsolution = DSolve[{linearequation, y[0]==4}, y,t]

This verifies that the solution satisfies both the equation and the initial condition:

linearequation/.particularsolution[[1]]

y[0]/.particularsolution[[1]]

Here is the solution to the same problem with the general initial condition y[0]==a:

particularsolution = DSolve[{linearequation, y[0]==a}, y,t]

This plots several integral curves of the equation for different values of a. The plot shows that the solutions have an inflection point if the parameter a lies between -1 and 1, while a global maximum or minimum arises for other values of a:

base=Plot[Evaluate[Table[y[t]/.particularsolution[[1]]/.a-> i, {i,-2,2,0.6}]], {t,-1.8,1.8}, PlotRange -> {-4,4}]

Challenge

Any numerical method has to be able to track the particular solution fixed by the initial value.

Detail

base=Plot[Evaluate[Table[y[t]/.particularsolution[[1]]/.a-> i, {i,-1.5,0,0.2}]], {t,-1.8,1.8}, PlotRange -> {-4,4},ImageSize->Full]

Euler’s Method

h=0.4

f[t_,y_]:=1+3*t*y

t0=0;
y0=-0.8;
trajectory={{t0,y0}};
y=y0;
Do[
y=y+h f[(k-1)h,y];
AppendTo[trajectory,{k h,y}],{k,Ceiling[2/h]}];

path=ListPlot[trajectory,Joined->True,Axes->False,PlotStyle->{Thick,Black},PlotRange->All];

Show[base,path]

Euler’s Method

h=0.02

f[t_,y_]:=1+3*t*y

t0=0;
y0=-0.8;
trajectory={{t0,y0}};
y=y0;
Do[
y=y+h f[(k-1)h,y];
AppendTo[trajectory,{k h,y}],{k,Ceiling[2/h]}];

path=ListPlot[trajectory,Joined->True,Axes->False,PlotStyle->{Thick,Black},PlotRange->All];

Show[base,path]