fsolve vs. running a model to equilibrium
4 views (last 30 days)
Show older comments
Hey all,
So I've been asking a lot of questions about a set of models I'm trying to find the equilibrium of and I've gotten a lot of great help on here (thank you!) but I'm still having trouble. I'm using fsolve to find the equilibrium values of a series of four differential equations (with the code I've copied below) and using this command in the command window:
>> x = fsolve(@odefun_dynamic, [5, 0.000001, 0.1420, 1e-7], options)
And I'm getting very different equilibrium values from when I just run the model to equilibrium and see what values it ends at. I think I could be hitting other local minima instead of the global, but I've been varying the x0 values for the fsolve command and I still don't get an equilibrium consistent to the one I see when running the model (which is the same set of differential equations). Are there other options I should set besides TolFun? Should I try using a different function than fsolve? Any help would be much appreciated! Thank you!!
function x_prime = odefun_dynamic(x)
options = optimset('TolFun',1e-14);
% Flow rate parameter
F = .01; % flow rate of medium
%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 constsant 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
0 Comments
Answers (0)
See Also
Categories
Find more on Perform Sensitivity Analysis in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!