Info

This question is closed. Reopen it to edit or answer.

Numerical solution is off for system of ODEs

1 view (last 30 days)
Wesley Neill
Wesley Neill on 15 Feb 2019
Closed: MATLAB Answer Bot on 20 Aug 2021
Thanks in advance for your help. I'm not looking for an explicit solution to my issue, but rather to be pointed in the right direction.
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 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 and trying to first solve 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 chosed method converges.
ludwig.jpg
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, finally, my question is simply to ask what might be causing the disparity I am seeing in my plots?
Code thus far:
% 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

Answers (1)

SandeepKumar R
SandeepKumar R on 6 Mar 2019
To be consistent use years also in your ode function (better to stick to year). Secondly play with tolerances of ode( Reltol,AbsTol). If that doesn't work change the ode solver (23,15s..etc..)

Tags

Products

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by