Ringe Kutta 4 Global Truncation error

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

[darova]

Try to check with ode45 also

[James Tursa]

As I stated in your other post, your h values are way too big to get meaningful results.