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

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.

Answer 1

[Jan]

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).

Answer 2

[Muhammad Sinan]

Hello Dear;

Please study the file using the following link.

https://www.researchgate.net/publication/31135539_Error_estimate_for_Runge-Kutta_type_solutions_to_systems_of_ordinary_differential_equations

Thank you so much!