# How to find error of Fourth Order Runge-Kutta method?

28 views (last 30 days)
Zane Dobler on 16 Dec 2018
I have code which uses fourth order Runge-Kutta to plot a phase diagram of how different initial states reach steady states over time. It involves a system of 2 nonlinear ordinary differential equations:
where x_1 isconcentration in a reactor and x_2 is temperature. The steady states are known to be:
Now I have to calculate the error of the method as a function of
Edit: h is dt in the code
My code for the Runge-Kutta is as follows:
clear
clc
close all
ss = [0.8634, 0.7753;0.8252,0.8121;0.3203,1.2977]
plot(ss(:,1),ss(:,2),"o")
%pause(10)
hold on
%initial values
xi = [ss(2,1)*1.01 ss(2,2)*1.01;
ss(2,1)*0.99 ss(2,2)*0.99;
0.1 0.1;
0.1 0.25;
0.1 0.5;
0.1 0.75;
0.1 1.0;
0.1 1.125;
0.1 1.25;
0.1 1.5;
0.25 0.1;
0.25 1.5;
0.5 0.1;
0.5 1.5;
0.75 0.1;
0.75 1.5;
1.0 0.1;
1.0 1.5;
1.25 0.1;
1.25 1.5;
1.5 0.1;
1.5 0.25;
1.5 0.5;
1.5 0.75;
1.5 0.875;
1.5 1.0;
1.5 1.125;
1.5 1.25;
1.5 1.5];
tf=10;
dt=0.001;
for i=1:length(xi(:,1))
x1i=xi(i,1);
x2i=xi(i,2);
soln=rk_cstr(x1i,x2i,tf,dt);
plot(soln(:,1),soln(:,2),"--")
axis([0 1.5 0 2])
end
hold off
function traj=rk_cstr(x1i,x2i,tf,dt)
a=100;
b=5;
g=149.39;
d=0.553;
f1=@(x1,x2) 1-x1-a*exp(-b/x2)*x1;
f2=@(x1,x2) 1-x2+g*exp(-b/x2)*x1-d*x2;
nt=tf/dt+1;
traj=zeros(nt,2);
x1t=x1i;
x2t=x2i;
traj(1,1)=x1t;
traj(1,2)=x2t;
for i=2:nt
k11=dt*f1(x1t,x2t);
k12=dt*f2(x1t,x2t);
k21=dt*f1(x1t+(1/2)*k11,x2t+(1/2)*k12);
k22=dt*f2(x1t+(1/2)*k11,x2t+(1/2)*k12);
k31=dt*f1(x1t+(1/2)*k21,x2t+(1/2)*k22);
k32=dt*f2(x1t+(1/2)*k21,x2t+(1/2)*k22);
k41=dt*f1(x1t+k31,x2t+k32);
k42=dt*f2(x1t+k31,x2t+k32);
% k21=dt*f1(x1t+(1/2)*dt,x2t+(1/2)*k11);
% k22=dt*f2(x1t+(1/2)*dt,x2t+(1/2)*k12);
% k31=dt*f1(x1t+(1/2)*dt,x2t+(1/2)*k21);
% k32=dt*f2(x1t+(1/2)*dt,x2t+(1/2)*k22);
% k41=dt*f1(x1t+dt,x2t+k31);
% k42=dt*f2(x1t+dt,x2t+k32);
x1t=x1t+((1/6)*(k11+2*k21+2*k31+k41));
x2t=x2t+((1/6)*(k12+2*k22+2*k32+k42));
traj(i,1)=x1t;
traj(i,2)=x2t;
end
end
How on earth do I calculate the error of this method? I've tried to find out how on Google and Wikipedia, but it isn't clear to me how I'm actually supposed to implement those methods into this code.

Jan on 17 Dec 2018
To get the error in terms of h^4 calculate the trajectory with different h and determine the differences in the results. You can compare the results with the known steady values also.
By the way, length(xi(:,1)) wastes time by creating the temporary vector xi(:,1), so better request the size directly: size(xi, 1).

Muhammad Sinan on 17 Jan 2020

### Categories

Find more on Runge Kutta Methods in Help Center and File Exchange

R2018a

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!