ODE solver for stiff system with 1 variable close to 0

Hi,

I'm solving a system of equations to model Daphnia and algae in a chemostat. I've already asked one question about this system (see: http://mathworks.com/matlabcentral/answers/6826-giant-spikes-in-model-from-equilibrium-when-run-for-relatively-long-time-periods), but I have another one. The code is the same from that question (but I'll copy it again here below). I'm having trouble running the system of equations, probably because it's very stiff and one variable (R - concentration of available resource) is very close to 0. The biggest problem is that the values for Daphnia and/or algae (N_A and N_d) run negative, which is not biologically realistic. Any advice y'all can give, probably on a different solver and/or solver options (absolute tolerance?), would be much appreciated! Thank you!

function algae_daphnia_NOrecycling()%(figure_handle) for multiple graphs

% Time Variable T = [0 500];

% Solver Options options = odeset('RelTol', 1E-12, 'MaxStep', 1E-3);

%state variables N_Ainit = 5; %N_Ainit = 1; %Q_Ainit = 1; Q_Ainit = 0.000001; N_dinit = 0.1420; R_init = 1.030400E-09; x = [N_Ainit; Q_Ainit; N_dinit; R_init];

% Solver [T x] = ode23s(@odefun_dynamic, T, x, options);

%Plot results - four graphs figure() subplot(2,2,1), plot(T, x(:,1)); xlabel('time (days)'); ylabel('Algae biomass (mg-C/L)'); subplot(2,2,2), plot(T, x(:,2)); xlabel('time (days)'); ylabel('Algae quota (Mol-P/mg-C)'); subplot(2,2,3), plot(T, x(:,3)); xlabel('time (days)'); ylabel('Daphnia biomass (mg-C/L)'); subplot(2,2,4), plot(T, x(:,4)); xlabel('time (days)'); ylabel('Available resource (Mol-P/L)');

end

function x_prime = odefun_dynamic(T, x)

% flow rate parameter F = 0.5;

% set parameters S = 1e-7; %initial substrate concentration p_max = 3.38e-6; %max intake rate of nutrient (P) by algae K_p = 1.29e-8; %half saturation constant of nutrient intake by algae u_Amax = 1; %apparent max growth rate Q_Amin = (2.5E-7); %subsistence quote for algae q_max = 1; %max uptake rate of algae by Daphnia K_q = 0.164; %half saturation constant of uptake of algae by Daphnia gamma = 0.5; %fraction of biomass of algae assimilated by Daphnia m_d = 0.02; %mortality of Daphnia

% Get state variables from x N_A = x(1); Q_A = x(2); N_d = x(3); R = x(4);

%functions

M_A = F * N_A;

p = (p_max * R) / (K_p + R);

I_R = N_A * p;

r_A = (u_Amax * N_A) * ((Q_A - Q_Amin) / (Q_A));

M_d = m_d * N_d;

I_A = (N_d * q_max * N_A) / (K_q + N_A);

r_d = gamma * I_A;

% differential equations

N_A_prime = r_A - M_A - I_A;

Q_A_prime = p - ((Q_A - Q_Amin) * u_Amax);

N_d_prime = r_d - M_d;

R_prime = (F * S) - (F * R) - I_R;

% Return derivatives in a single column vector x_prime = [N_A_prime; Q_A_prime; N_d_prime; R_prime]

end

Answer 1

[Laura Proctor]

Try using ODESET with the NonNegative property set.