Incorrect equilibrium values using fsolve

2 次查看(过去 30 天)
Hi!
I'm using fsolve to find equilibrium values of a series of three differential equations ( the Droop model). I've been using fsolve and getting weird equilibrium values, so I went back and solved the equations by hand and fsolve has been giving me incorrect values. The script I've been using is below and I've been using this command in the Command Window:
fsolve(@odefun_dynamic_droop, [1, 1, 1])
I've tried various x0 values but that doesn't seem to change things. For the parameter values in the script, fsolve finds zeros at
2.6605e-07 3.6300e-06 2.5009e-06
when I calculate
2.5253e-07 2.2529e-09 4.9704e-02
Are there options I can change or something to make fsolve more accurate? I haven't entered the parameters or equations wrong.
Thank you!
Code:
function x_prime = odefun_dynamic_droop(x)
%flow rate parameter
F = 0.01; %flow rate
%fixed parameters
S = 2.5e-6; %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
% Get state variables from x
N_A = x(1);
Q_A = x(2);
R = x(3);
%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);
% differential equations
N_A_prime = r_A - M_A;
Q_A_prime = p - ((Q_A - Q_Amin) * u_Amax);
R_prime = (F * S) - (F * R) - I_R;
% Return derivatives in a single column vector
x_prime = [N_A_prime; Q_A_prime; R_prime];
end

采纳的回答

Teja Muppirala
Teja Muppirala 2011-5-5
If you're going to do it by hand, then at least let MATLAB help you check your answers (using symbolic variables):
%flow rate parameter
F = 0.01; %flow rate
%fixed parameters
S = 2.5e-6; %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
% Set up symbolic variables to be solved for
N_A = sym('N_A');
Q_A = sym('Q_A');
R = sym('R');
%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);
% differential equations
N_A_prime = r_A - M_A;
Q_A_prime = p - ((Q_A - Q_Amin) * u_Amax);
R_prime = (F * S) - (F * R) - I_R;
% Return derivatives in a single column vector
x_prime = [N_A_prime; Q_A_prime; R_prime];
sol = solve(x_prime);
sol = [sol.N_A sol.Q_A sol.R];
x_solution1 = double(sol(1,:))
x_solution2 = double(sol(2,:))
Which yields two solutions.
x_solution1 =
1.0e-005 *
0 0.361264873254009 0.250000000000000
x_solution2 =
9.899961805784013 0.000000252525253 0.000000000009645
Looks like FSOLVE didn't do so bad.
  4 个评论
Teja Muppirala
Teja Muppirala 2011-5-5
Yes. The empty values (the ones you did not explicitly set) retain their default values. You can view the default values for the FSOLVE function using
optimset('fsolve')

请先登录,再进行评论。

更多回答(1 个)

Andrew Newell
Andrew Newell 2011-5-5
The answer provided by fsolve looks better than yours:
>> odefun_dynamic_droop([2.6605e-07 3.6300e-06 2.5009e-06])
ans =
1.0e-06 *
0.2451
-0.0173
-0.0000
>> odefun_dynamic_droop([2.5253e-07 2.2529e-09 4.9704e-02])
ans =
1.0e-03 *
-0.0278
0.0036
-0.4970

标签

Community Treasure Hunt

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

Start Hunting!

Translated by