How to find error of Fourth Order Runge-Kutta method?
22 次查看(过去 30 天)
显示 更早的评论
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:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/199217/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/199218/image.png)
where x_1 isconcentration in a reactor and x_2 is temperature. The steady states are known to be:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/199219/image.png)
Now I have to calculate the error of the method as a function of ![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/199220/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/199220/image.png)
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 个)
Jan
2018-12-17
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).
0 个评论
Muhammad Sinan
2020-1-17
Hello Dear;
Please study the file using the following link.
Thank you so much!
0 个评论
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!