- #1
gpax42
- 25
- 0
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
%% 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