Euler’s method for y'=-2t y+t

Define the ODE right hand side function:

f[t_,y_]:=-2t y+t;

For this example, the ODE can be solved analytically, with initial condition y(0)=-1:

FullSimplify[DSolve[{y'[t]==f[t,y[t]],y[0]==-1},y[t],t]]

Define the exact solution just obtained as a new function:

yex[t_]:=(1/2)-(3/2)(E^(-(t^2)));

To implement Euler’s method, first define a function that performs one iteration using time-step h:

EulerOneStep[ty_,h_]:=Block[{tn=ty[[1]],yn=ty[[2]]},
{tn+h,yn+h f[tn,yn]}
];

Then use a for-loop to find the approximate solution up to final time (Subscript[t, f]), with initial condition y(0)=-1 as before:

EulerIVP[tf_,h_]:=Block[{},
ty={0,-1};
sol={};
For[i=0,i<tf/h,i++,
ty=EulerOneStep[ty,h];
sol=Append[sol,ty];
];
sol
];

For example, using (Subscript[t, f])=2 and h=0.1, we obtain the approximate solution:

EulerIVP[2,0.1]//MatrixForm

Plot the direction field together with the exact solution (in black), and the approximate solution using Euler’s method (in red), with the step-size h controllable:

Manipulate[Show[VectorPlot[{1,f[t,y]},{t,0,1},{y,-1,0}],ListPlot[EulerIVP[1,h],PlotStyle->Red],Plot[yex[t],{t,0,1},PlotStyle->Black]],{h,0.05,0.1,0.05}]

A note on higher-order methods

To use k-th order method to achieve global error ϵ, we need to take

ϵ∼(h^k)=>h∼(ϵ^(1/k)).

Thus (h^(-1))∼(ϵ^(-1/k)) time steps are needed; each time step costs k operations. Hence the total number of operations is

M[k_]:=k (ϵ^(-1/k));

Take ϵ=0.01, then the number of operations depends on k as

ϵ=0.01;Plot[M[k],{k,1,10},PlotRange->{0,20}]

It can be seen that the number of operations is minimized at k≈4, namely fourth-order methods are the most efficient (at least for this choice of ϵ).