CODEWITHSUNDEEP
matlab
Stralo
Stralo

Reputation: 474

Simulations hang when running ode45

I'm suppose to make a model for an algae population. Here's the code I have so far (all written from examples online). When i run Solve_algaepop, it just hangs for a long time.

Any ideas why? Is there any obvious thing I'm doing wrong? The equations are from a research paper.

This is Solve_algaepop.m. In the equations for r1 and r2, P10 and P20 are supposed to be the values P1 = x(1) and P2 = x(2) defined in algaepop_model.m. I don't know how to access the values when I'm in Solve_algaepop.m 

% Initial conditions
P10 = 560000; %from Chattopadhyay; estimated from graph
P20 = 250000; %same as above
Z0 = 280000; %
N0 = 0.6; %from Edwards

%some variables that the expressions of the parameters use
lambda = .6;
mu = .035;
k = 0.05;

%define parameters (start with estimates from Edwards paper):
r1 = (N0/(.03+N0))*((.2*P10)/(.2 + .4*P10)); 
r2 = (N0/(.03+N0))*((.2*P20)/(.2 + .4*P20));
a = Z0*((lambda*P10^2)/(mu^2 + P10^2));%G1: zooplankton growth function from Edwards paper
% m1 = .15; %r in Edwards paper
m1 = .075; % q in Edwards
m2 = .15;% r in Edwards paper
m3 = .15; % r in Edwards paper
d = 0.5;
cN = k;%*(N-N0);

par = [r1 r2 a m1 m2 m3 d cN]; % Creates vector of parameter values to pass to the ode solver

tspan = 0:1:300; %(Note: can also use the function linspace)

x0 = [P10 P20 Z0 N0]; % Creates vector of initial conditions

[t,x] = ode45(@algaepop_model,tspan,x0,[],par);
plot(t,x)

And here is algaepop_model.m

function dxdt = algaepop_model(t,x,par)

P1 = x(1);
P2 = x(2);
Z = x(3);
N = x(4);

r1 = par(1);
r2 = par(2);
a = par(3);
m1 = par(4);
m2 = par(5);
m3 = par(6);
d = par(7);
cN = par(8);

dxdt = zeros(4,1);

dxdt(1) = r1*N*P1 - m3*P1 - a*P1*Z;
dxdt(2) = r2*N*P2 - a*P2*Z - m2*P2;
dxdt(3) = a*P2*Z + a*P1*Z - m1*Z;
dxdt(4) = d*m2*P2 + d*m1*Z + d*m3*P1 + cN - r2*N*P2 - r1*N*P1;

end

Thanks for the help.

Upvotes: 3

Views: 1005

Answers (1)

horchler
horchler

Reputation: 18484

Let's debug. One of the simplest things that you can do is print out t and x inside of your integration function, algaepop_model. As soon as you do this, you'll probably notice what's happening: ode45 is taking extremely small steps. They're on the order of 1.9e-9. With steps that small, it will take forever to simulate to t = 300 (and even longer if you print stuff out on each step).

This might be caused by a poor choice of initial conditions, poor scaling or dimensionalization, a typo resulting in the wrong equations, or simply that you're using an inappropriate solver (and/or tolerances) for the particular problem. I can't really address the first two situations and must assume that you don't have any errors. Thus, in this case you have what is effectively a stiff system and ode45 is not a particularly good choice in such cases. Simply changing the solver to ode15s results in the following plot almost immediately:

plot

As you can see, there are very large changes over a short period of time in the initial portion of the plot. If you zoom in you'' see that the huge spike happens in the first unit of time (you might output more time points or just let tspan = [0 300]). Some state variables are changing rapidly while others are varying more gradually. Such high frequencies and differences in time scales are the hallmarks of stiff systems. I'd suggest that, in addition to confirming that your code is correct, you also try adjusting the integration tolerances as well via odeset. Make sure that tighter tolerances produce qualitatively similar results. You can also try the other stiff solvers in the ODE suite if you like.

Lastly, it's more efficient and up-to-date to pass your parameters via the function handle itself rather than how you're doing it. Here's how:

[t,x] = ode15s(@(t,x)algaepop_model(t,x,par),tspan,x0);

Upvotes: 1

Related Questions