ode45 converges to correct curve shape, but with wrong solution

Question

Thanks in advance for your help. I'm not looking for an explicit solution to my problem, but rather to have my probably obvious errors pointed out.

I have been plugging away at solving a system of non-linear, first order ODEs in MATLAB. The system was solved numerically in this study: http://web.math.ku.dk/~moller/e04/bio/ludwig78.pdf

I have been following the documentation for ode45, and have code that runs.

I have done all of the work to understand and recreate the model from scratch. I presented the qualitative part for a class project. What I am doing now is taking that project a step farther by solving the system in MATLAB with runge-kutta (or any method that works). Finally, I want to dive into the theory behind the numerical analysis to find out why the chosen method converges.

Here is a plot of the numerically solved system, which I am trying to re-create:

I have found that I can create a plot with roughly the same shape, but there are several problems:

1. The time-scale over which the change occurs is three times that of the above plot.
2. The range of function values is is vastly wrong.
3. The desired shapes only occur if I tweak the initial conditions to be significantly different than what is shown near t=0 above.

So what I'm looking for is a reason for these discrepancies. I've checked my system of ODEs and parameter values so many times my eyes are blurry. Perhaps I am missing something conceptually?

Code:

% System Parameters:
r_b = 1.52;
k_b = 355;
alph = 1.11;
bet = 43200;
r_e = 0.92;
k_e = 1;
p = 0.00195;
r_s = 0.095;
k_s = 25440;

tspan = [0 200];
init = [1 1 1];

[t, Y] = ode45(@(t,y) odefcn(t, y, r_b, k_b, alph, bet, r_e, k_e, p, r_s, k_s), tspan, init);
subplot(3,1,1);
plot(t,Y(:,1),'b');
title('Budworm Density');

subplot(3,1,2)
plot(t,Y(:,2),'g');
title('Branch Density');

subplot(3,1,3);
plot(t,Y(:,3),'r');
title('Foliage Condition');

function dydt = odefcn(t, y, r_b, k_b, alph, bet, r_e, k_e, p, r_s, k_s)
dydt = [ r_b*y(1)*(1 - y(1)/(k_b*y(2))) - bet*(y(1)^2/((alph*y(2))^2 + y(1)^2));
         r_s*y(2)*(1 - (y(2)*k_e)/(k_s*y(3)));
         r_e*y(3)*(1 - (y(3)/k_e)) - p*y(1)/y(2)
        ];

end

Answer 1

I don't see anything wrong with your code as such. But I think there are some subtleties involved in producing the figure which are not well explained in the paper.

1) The S axis is scaled (it says 'relative' in the label). I believe they've scaled S by k_s. I think you also need to scale the parameter p (set p = p*k_s) else the final term in the equation for E will be tiny and the E population won't decrease over the required timescales.

2) I think they must have enforced some lower limit on E, to avoid dividing by 0. You can see in the figure that E->0 first, but in your equation for S, if this happened then you would be dividing by 0 and the solver wouldn't converge.

Putting these together, the following slight modification of your code produces a result more similar to that in the paper:

% System Parameters:
r_b = 1.52;
k_b = 355;
alph = 1.11;
bet = 43200;
r_e = 0.92;
k_e = 1;
p = 0.00195;
r_s = 0.095;
k_s = 25440;

% Scale p with k_s
p = p*k_s;

tspan = [0 50]; % [0 200];
init = [1e-16 0.075*k_s 1]; % [1 1 1];

[t, Y] = ode45(@(t,y) odefcn(t, y, r_b, k_b, alph, bet, r_e, k_e, p, r_s, k_s), tspan, init);

% To scale before plotting, so everything fits on a 0->1 y axis.
maxB = 500;
S_scale = k_s;

figure('Position', [200 200 1000 600]);
hold on;
plot(t,Y(:,1)/maxB,'b');
plot(t,Y(:,2)/(S_scale),'g');
plot(t,Y(:,3),'r');

ylim([0, 1]);

hold off;

box on;

legend({['Budworm Density, B / ', num2str(maxB)], 'Branch Density, S / 0.75', 'Foliage Condition, E'}, ...
    'Location', 'eastoutside')

function dydt = odefcn(t, y, r_b, k_b, alph, bet, r_e, k_e, p, r_s, k_s)

% Place lower limit on E
E = max(y(3), 1e-5);

dydt = [ r_b*y(1)*(1 - y(1)/(k_b*y(2))) - bet*(y(1)^2/((alph*y(2))^2 + y(1)^2));
         r_s*y(2)*(1 - (y(2)*k_e)/(k_s*E));
         r_e*E*(1 - (E/k_e)) - p*y(1)/y(2)
        ];

end

There is a lot of sensitivity to the initial conditions.

A further tweak gets you closer still to the original figure, but I'm not sure if this is just a hack: in the first equation, replace k_b*y(2) with just k_b. Without this, the Budworm density becomes too big before decreasing. The new plot is below.