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
0 个评论
采纳的回答
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
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
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
0 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Special Values 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!