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

Question

Zane Dobler 2018-12-16

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/436063-how-to-find-error-of-fourth-order-runge-kutta-method

回答： Muhammad Sinan 2020-1-17

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.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Jan 2018-12-17

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/436063-how-to-find-error-of-fourth-order-runge-kutta-method#answer_352769

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