Dormand/Prince Runge-Kutta with order4 and 5

The Dopri54 is a method of Dormand / Prince Runge-Kutta pairs are formed by Runge-Kutta simply method with order 5 and 4. Therefore a fiveth-order Runge-Kutta method is using each step forward and a fourth-order Runge-Kutta to estimate the error.

We will show at this section, the Matlab code for these method for the initial value problem

y'=y
y(0)=1
We want to know the y value at t = 1. The obvious solution is

y(t)=et, and therefore y(1) = e

Matlab code for the Dormand/Prince method Dopri54: File dopri54.m

% DOPRI54 solves system of ODEs with the Dormand and Prince 54 code.
% INPUT
% funcion - String containing name of user-supplied
% problem description.
% Call: yprime = fun(t,y,rpar) where funcion = 'fun'.
% t - Time (scalar).
% y - Solution column-vector.
% rpar - parameters possibly needed in buinding yprime
% yprime - Returned derivative column-vector; yprime(i)
% = dy(i)/dt.
% t0 - Initial value of t.
% t1 - Final value of t (must be t1 > t0).
% y0 - Initial value column-vector.
% tol - tolerance for the local error
% OUTPUT
% y - computed solution at y=t1.
% stad - vector containing number of steps, rejections and
% function calls
t=t0;
y=y0;
h=tol^(1/5)/4;
step=0;
nrej=0;
fcall=1;
a4=[44/45 -56/15 32/9]';
a5=[19372/6561 -25360/2187 64448/6561 -212/729]';
a6=[9017/3168 -355/33 46732/5247 49/176 -5103/18656]';
a7=[35/384 0 500/1113 125/192 -2187/6784 11/84]';
e=[71/57600 -1/40 -71/16695 71/1920 -17253/339200 22/525]';
k1=feval(funcion,t,y,rpar);

while t < t1
if t+h > t1; h=t1-t; end
k2=feval(funcion,t+h/5,y+h*k1/5,rpar);
k3=feval(funcion,t+3*h/10,y+h*(3*k1+9*k2)/40,rpar);
k4=feval(funcion,t+4*h/5,y+h*(a4(1)*k1+a4(2)*k2+a4(3)*k3),rpar);
k5=feval(funcion,t+8*h/9,y+h*(a5(1)*k1+a5(2)*k2+a5(3)*k3+...
a5(4)*k4),rpar);
k6=feval(funcion,t+h,y+h*(a6(1)*k1+a6(2)*k2+a6(3)*k3+a6(4)*k4+...
a6(5)*k5),rpar);
yt=y+h*(a7(1)*k1+a7(3)*k3+a7(4)*k4+a7(5)*k5+a7(6)*k6);
k2=feval(funcion,t+h,yt,rpar);
est=norm(h*(e(1)*k1+e(2)*k2+e(3)*k3+e(4)*k4+e(5)*k5+...
e(6)*k6),inf);
fcall=fcall+6;
% [t h est]

if est < tol
t=t+h;
k1=k2;
step=step+1;
y=yt;
else
nrej=nrej+1;
end
h=.9*min((tol/(est+eps))^(1/5),10)*h;

end
if nargout > 1
end

file fun.m

This is a file containing the function to evaluate, at this case y' = y

function[y] = fun(x, y, rpar)
u=x;

Running

Keep the two files in the same directory, and from the Matlab console move to that directory
Run the following command from the console
>>y = dopri54c('fun', 0, 1, y0, 0.0001,0);

Now we have the result stored in the variable y, and simply run to see
>& gt; y
and that gives us the result

y = 2.718...

Execute here the Dormand/Prince 5 4 method in our Runge kutta calculator