Ringe Kutta 4 Global Truncation error

2 次查看（过去 30 天）

Wytse Petrie 2020-1-8

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/499512-ringe-kutta-4-global-truncation-error

评论： James Tursa 2020-1-8

Hi there,

In order to calculate the global truncation error of the Ringe Kutta 4, i calculated y1 and y2 with a dx and y1 and y2 with 2*dx. Then i can subtract them from eachother and devide by (2^p)-1, where p = 4 (because of the Ringe Kutta). The model is a predator prey model.

However the error that i got is to big, can someone help me?

% parameter
a = 3;
c = 3;
alpha = 2*10^(-3);
gamma = 7*10^(-3);
%% k = 0
k = 0
delta_t0 = 0.5/(2^k);
%initial conditions
y10(1)=600;
y20(1)=1000;
t0(1)=0;
t0 = ti:delta_t0:tf;
h0 = delta_t0;
N0 = ceil(tf/h0);
% define the seasonal fishing load factor
W0 = zeros(11,1);
for n =1:N0
W0(n,1)=fishing_load_factor(t0(n));
end
% Define functions
fy1 = @(t,y1,y2,W0) (a-alpha*y2)*y1-W0*y1;
fy2 = @(t,y1,y2) (-c+gamma*y1)*y2;
% Update loop
for i = 1:N0
    %Update time
    t0(1+i)=t0(i)+h0;
    
    % Update y1 and y2
    K11 = fy1(t0(i)      ,y1(i)        , y2(i), fishing_load_factor(t0(i)));
    K12 = fy2(t0(i)      ,y1(i)        , y2(i));
    
    K21 = fy1(t0(i)+h0/2  ,y1(i)+h0/2*K11, y2(i)+h0/2*K12, fishing_load_factor(t0(i)+h0/2));
    K22 = fy2(t0(i)+h0/2  ,y1(i)+h0/2*K11, y2(i)+h0/2*K12);
    
    K31 = fy1(t0(i)+h0/2  ,y1(i)+h0/2*K21, y2(i)+h0/2*K22, fishing_load_factor(t0(i)+h0/2));
    K32 = fy2(t0(i)+h0/2  ,y1(i)+h0/2*K21, y2(i)+h0/2*K22);
 
    K41 = fy1(t0(i)+h0    ,y1(i)+h0*K31   , y2(i)+h0*K32, fishing_load_factor(t0(i)+h0));
    K42 = fy2(t0(i)+h0    ,y1(i)+h0*K31   , y2(i)+h0*K32);
    
    y10(1+i) = y1(i)+h0/6*(K11 + 2*K21 + 2*K31 + K41);
    y20(1+i) = y2(i)+h0/6*(K12 + 2*K22 + 2*K32 + K42);
end
delta_t0f = delta_t0*2;
%initial conditions
y10f(1)=600;
y20f(1)=1000;
t0f(1)=0;
t0f = ti:delta_t0f:tf;
h0f = delta_t0f;
N0f = ceil(tf/h0f);
% define the seasonal fishing load factor
W0f = zeros(N0f,1);
for n =1:N0f
W0f(n,1)=fishing_load_factor(t0f(n));
end
% Define functions
fy1 = @(t,y1,y2,W0f) (a-alpha*y2)*y1-W0f*y1;
fy2 = @(t,y1,y2) (-c+gamma*y1)*y2;
% Update loop
for i = 1:N0f
    %Update time
    t0f(1+i)=t0f(i)+h0f;
    
    % Update y1 and y2
    K11 = fy1(t0f(i)      ,y1(i)        , y2(i), fishing_load_factor(t0f(i)));
    K12 = fy2(t0f(i)      ,y1(i)        , y2(i));
    
    K21 = fy1(t0f(i)+h0f/2  ,y1(i)+h0f/2*K11, y2(i)+h0f/2*K12, fishing_load_factor(t0f(i)+h0f/2));
    K22 = fy2(t0f(i)+h0f/2  ,y1(i)+h0f/2*K11, y2(i)+h0f/2*K12);
    
    K31 = fy1(t0f(i)+h0f/2  ,y1(i)+h0f/2*K21, y2(i)+h0f/2*K22, fishing_load_factor(t0f(i)+h0f/2));
    K32 = fy2(t0f(i)+h0f/2  ,y1(i)+h0f/2*K21, y2(i)+h0f/2*K22);
 
    K41 = fy1(t0f(i)+h0f    ,y1(i)+h0f*K31   , y2(i)+h0f*K32, fishing_load_factor(t0f(i)+h0f));
    K42 = fy2(t0f(i)+h0f    ,y1(i)+h0f*K31   , y2(i)+h0f*K32);
    
    y10f(1+i) = y1(i)+h0f/6*(K11 + 2*K21 + 2*K31 + K41);
    y20f(1+i) = y2(i)+h0f/6*(K12 + 2*K22 + 2*K32 + K42);
end
% calculating errors
p = 4;
error01 = abs((y10(3)-y10f(2))/((2^p)-1))
error02 = abs((y20(3)-y20f(2))/((2^p)-1))