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
function [y,stad]=dopri54c(funcion,t0,t1,y0,tol,rpar)
% 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
stad=[step nrej fcall];
end
function[y] = fun(x, y, rpar)
u=x;
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
function [y,stad]=dopri54c(funcion,t0,t1,y0,tol,rpar)
% 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
stad=[step nrej fcall];
end
file fun.m
This is a file containing the function to evaluate, at this case y' = yfunction[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
Was useful? want add anything?
Post here
Post from other users
yeshambel:
2020-02-13 23:46:47Is it possible to apply Dorman-Prince method to solve systems of non-linear differential equation in a model? How to compare Dorman-Prince method and Nonstandard finite difference method by plotting in the same axis? Thank you!!
Carlos:
2020-02-20 23:46:47Hi yeshambel. First of all i want to apology for delay in answer to you. The first you question you asked is: of course. Dormand-Prince (as the Runge-Kutta methods) method can be applied to non linear system. In fact, almost every ODE problem problems solved thu numeric methods costumates to be non linear (but can be linear too). The second question is quite bite complexer to answer, i will try to help you but it admits to more more arged Dormand-Prince (DP) Method is applied to initial value problems in the form: y'=f(x,y,y') y(a)=b I only used Finite diferences (DF) to solve partial derivates equation problem, wich are different. Depending on what (DP) method you use, you can have have much more accuracy than (DF) if using the same number of nodes, because order of (DP) is (usually) greater than (DF) case. I hope i've been helpful. Carlos.
Post here