As I was reading about the fourth order Runge-Kutta method and how MATLAB’s own ODE45 is a modified version of this, I thought it would be interesting to check it out myself. For details on ODE solvers, take a look at Mark Embree’s notes, which are the best notes on Numerical Analysis I’ve read. Also, have a look at Mathwork’s demo-site on differential equations which shows how to solve systems of differential equations using ODE45.

Included in our comparison is the trusty Euler’s method, and the code is shown below. **DOWNLOAD Euler’s method**.

function [t,y] = euler(yprime, tspan, y0, h)
% function [t,x] = euler(xprime, [t0 tfinal], x0, h)
% Approximate the solution to the ODE: y'(t) = yprime(x,t)
% from t=t0 to t=tfinal with initial condition y(t0)=y0.
% yprime should be a function that can be called like: yprime(t,y).
% h is the step size: reduce h to improve the accuracy.
% The ODE can be a scalar or vector equation.
t0 = tspan(1); tfinal = tspan(end);
% set up the t values at which we will approximate the solution
t=[t0:h:tfinal];
% include tfinal even if h does not evenly divide tfinal-t0
if t(end)~=tfinal, t = [t tfinal]; end
% execute Euler's method
y = [y0 zeros(length(y0),length(t)-1)];
for j=1:length(t)-1
y(:,j+1) = y(:,j) + h*feval(yprime,t(j),y(:,j));
end

The implementation of Runge-Kutta is equaly simple. **Download fourth order Runge-Kutta**.

% Classical fourth-order Runge-Kutta method
function [t,y] = runge_kutta_4order(yprime, tspan, y0, h)
% function function [t,y] = runge_kutta_4order(yprime, [t0 tfinal], y0, h)
% Approximate the solution to the ODE: y'(t) = yprime(x,y)
% from t=t0 to t=tfinal with initial condition y(t0)=y0.
% yprime should be a function that cna be called like xprime(t,y).
% h is the step size: reduce h to improve the accuracy.
% The ODE can be a scalar or vector equation.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Sample code to call runge_kutta_4order.m for the equation x'(t)=sin(t*x)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% yprime = inline('sin(y*t)','y','y');
% Instead of inline we can use yprime=@(t,y) sin(y*t).
% tspan = [0 10]; Integrate from t=0 to t=10.
% y0 = 5;
% h = 0.1; Time step.
% [t,y] = runge_kutta_4order(yprime, tspan, y0, h); Call Runge-Kutta.
% figure(1), clf
% plot(t,y,'b.','markersize',15); Plot output
t0 = tspan(1); tfinal = tspan(end);
% set up the t values at which we will approximate the solution
t = [t0:h:tfinal]';
% include tfinal even if h does not evenly divide tfinal-t0
if t(end)~=tfinal, t=[t tfinal]; end
m = length(t);
y = [y0 zeros(length(y0), length(t)-1)];
for i=1:(m-1)
k1 = feval(yprime,t(i),y(:,i));
k2 = feval(yprime,t(i)+0.5*h, y(:,i)+(h.*k1)/2);
k3 = feval(yprime,t(i)+0.5*h, y(:,i)+(h.*k2)/2);
k4 = feval(yprime,t(i)+h, y(:,i)+h.*k3);
y(:,i+1) = y(:,i)+(h*(k1+2*k2+2*k3+k4))/6;
end

** **We finally write a small script that uses all three solvers and plots the results. **Download ODE comparison script**.

% Comparison of ODE solvers on a system of equations.
% 1. Euler
% 2. Runge-Kutta
% 3. ODE45
clear
clc
yprime = @(t,y) vanderpoldemo(t,y,1); % Use y0=[2 0],tspan=[0 20]
tspan = [0 20];
y0 = [2 0]';
h = 0.1; %Time step.
for i=1:length(y0)
figure(i), clf, hold on, grid on, xlabel('t'), ylabel('y(t)'),...
title('Comparison of ODE solvers on a system of equations')
[teuler,yeuler] = euler(yprime, tspan, y0, h);
plot(teuler,yeuler(i,:),'r.-','markersize',10);
[trunge,yrunge] = runge_kutta_4order(yprime, tspan, y0, h);
plot(trunge,yrunge(i,:),'g.-','markersize',10);
[tode45,yode45] = ode45(yprime, tspan, y0);
plot(tode45,yode45(:,i),'b.-','markersize',10);
legend('Euler','Runge-Kutta','ODE45')
end

Now to the fun part, plots.

As expected, Runge-Kutta is almost identical to ODE45. It should be noted though, that ODE45 uses an adaptive scheme which adjusts the step size when needed. Euler’s is quite a bit off, which was also to be expected. Quite a trivial exercise, but still nice to see theory in practice.