Incorrect equilibrium values using fsolve

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

Answer 1

[Teja Muppirala]

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.

[Teja Muppirala]

But to answer your question though ("Are there options I can change or something to make fsolve more accurate?")

Yes. You can use OPTIMSET to set the function tolerance to be smaller.

For example

options = optimset('TolFun',1e-14);

fsolve(@odefun_dynamic_droop, [1, 1, 1],options)

This now yields a much more stringent answer.

[LS]

Yes! Both of these answers were what I was looking for - I need to find both roots (the near zero one isn't helpful for me biologically) but the second one is super helpful! Thank you!!!

And a further question about optimset - I was playing around with the optimset to set the function tolerance smaller, but when I set the function tolerance it looked like it cleared all of the other options. Do those stay at default when you change just TolFun? Or is there a way to set just one of options and keep the rest at default?

Thank you very much for your help!!

[Teja Muppirala]

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')

[LS]

Thank you!!

Answer 2

[Andrew Newell]

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