Numerical solution of initial value problem using ode45 and Euler method

Using ode45

clearvars; clf

We want to find the solution of the IVP

for

with

f = @(t,y) y-t
f = function_handle with value:
    @(t,y)y-t
t0 = 0; y0 = .5; T = 1;
dirfield(f,0:.1:1,.5:.1:1); hold on

We find the "exact" solution using ode45:

[ts,ys] = ode45(f,[t0,T],y0);

yex = ys(end)

yex = 0.6409

plot(ts,ys,'b')

text(T,yex,' exact')

xlabel('t'); title('Numerical solution of IVP using ode45')

Euler method

The function EulerMethod is defined at the end of this file.

for n = 2.^(0:6)

[ts,ys] = EulerMethod(f,[t0,T],y0,n);

plot(ts,ys,'.-')

Y = ys(end); % approximation for y(T)

fprintf('n = %2g, y(T)≈%8g, error = %8g \n',n,Y,Y-yex)

if n<=4

disp([ts,ys])

end

text(T,Y,sprintf(' n=%g',n))

end

n =  1, y(T)≈       1, error = 0.359141 

         0    0.5000
    1.0000    1.0000

n =  2, y(T)≈   0.875, error = 0.234141 

         0    0.5000
    0.5000    0.7500
    1.0000    0.8750

n =  4, y(T)≈0.779297, error = 0.138438 

         0    0.5000
2500    0.6250
5000    0.7188
7500    0.7734
0000    0.7793

n =  8, y(T)≈0.717108, error = 0.0762487 
n = 16, y(T)≈0.681036, error = 0.0401767 
n = 32, y(T)≈0.661505, error = 0.0206459 
n = 64, y(T)≈0.651328, error = 0.0104684 

hold off

title('Euler Method using 1,2,4,...,64 steps')

Code for Euler method

function [ts,ys] = EulerMethod(f,interval,y,n)
% [ts,ys] = EulerMethod(f,[t0,T],y0,n)
% initial value problem y'=f(t,y), y(t0) = y0
% use n steps of Euler method to approximate solution in interval [t0,T]
t = interval(1); T = interval(2);
h = (T-t)/n;
ts = zeros(n+1,1); ys = zeros(n+1,1);   % vectors for t-values, y-values
ts(1) = t; ys(1) = y;
for k=1:n
  s = f(t,y);
  t = t + h;
  y = y + s*h;
  ts(k+1) = t; ys(k+1) = y;             % store t,y in vectors ts, ys
end
end