MATLAB Lotka-Volterra Approximate Solution Error: Trapezoid Method in Matlab

[gpax42]

I'm trying to calculate the error of the Trapezoid Method in my population approximations for wolves and moose over time given certain parameters. Ultimately what i need to do is perform a series of trapezoid approximations using different 'h' values. Each successive 'h' value I use is half the previous value which can be seen in my code below. As my 'h' value gets smaller, the error between successive final population counts of moose and wolves should also be decreasing. Due to the trapezoid method's order of accuracy, i should expect the log-log graph of the 'h-values' vs. 'error' to have a slope of 2. However, mine still has a slope of 1. I know deciphering my code is a pain, but I'd be extremely grateful to anyone that can point me in the right direction... Thanks a million!

%% Trapezoid Error

T = 50;

N = 50;
h = T/N;
w = [20]; %taken from wolf population data
m = [563]; %taken from moose population data

%parameters calculated by least squares for the lotka-volterra equation

c = [-0.0496 -0.0001 -0.1796 -0.0071];
c1 = c(1);
c2 = c(2);
c3 = c(3);
c4 = c(4);

%anonymous function assignments for the two differential equations defined by lotka-volterra

f1 = @(w,m) (c1*w - c2*w*m);
f2 = @(w,m) (-c3*m + c4*w*m);

for i = 1:13
for n = 1:N
w_pred = w(end) + h*f1(w(end),m(end));
m_pred = m(end) + h*f2(w(end),m(end));
w(end+1) = w(end) + h/2*(f1(w(end),m(end)) + f1(w_pred,m_pred));
m(end+1) = m(end) + h/2*(f2(w(end),m(end)) + f2(w_pred,m_pred));
end
wfinal(i) = w(length(w));
mfinal(i) = m(length(m));
hval(i) = h;
h = h/2;

%once the 13th iteration is complete, plot the data, otherwise, perform iterations again with halved 'h' value

if i == 13
TrapError = abs(wfinal(2:end) - wfinal(1:end-1)) + abs(mfinal(2:end) - mfinal(1:end-1));
loglog(hval(2:end),TrapError,'s-');
xlabel('h-value')
ylabel('Trapezoid Error')
title('Error from Trapezoid Approx')
grid on
return
end

w = [20];
m = [563];
end

 

[AlephZero]

Each successive 'h' value I use is half the previous value which can be seen in my code below. As my 'h' value gets smaller, the error between successive final population counts of moose and wolves should also be decreasing. Due to the trapezoid method's order of accuracy, i should expect the log-log graph of the 'h-values' vs. 'error' to have a slope of 2. However, mine still has a slope of 1.

Are you confusing local and global error estimates?

The local error estimate for one step is Ah^2 (where A is a constant).

But the number of steps for the integration is B/h (were B is another constant)

So the global error estimate (Ah^2)(B/h) = ABh.

I haven't read your code in detail, so apologies if this is irrelevant!

 

[gpax42]

Thank you for the helpful response

since I'm plotting the loglog of the difference between my final population values(global error) that resulted from halving my h-value, i believe youre correct in noting that my graph should have a slope of 1.

This, however, begs the question of how to get localized error. My approach once again produced a slope of 1... any further suggestions