How can I solve 3 non linear functions simultaneously?

Hello,

I can´t find how to solve 3 functions which results are interrelated.

In my code, the output "dNO3" of the function "fun_06NO3" is input for the function "fun_07DBO". In the same way, the output "dDBO" of the function "fun_07DBO" is input for the function "fun_08_OD". Finally, the output "dOD" of the function "fun_08_OD" is input for both functions, "fun_06NO3" and "fun_07DBO".

[dNO3,a_NO3,Foxdn,kdn]= fun_06NO3(CNO3, CMNO3, Q, tao, TR, Temp, knNA, dNAm, NO3_kdn, dOD);
[dDBO, a_DBO, kdDBO]  = fun_07DBO(CDBO, CMDBO, Q, tao, TR, Temp, Foxdn, kdn, dNO3, DBO_kd, dOD);
[dOD , a_D  , ka, Os] = fun_08_OD(C_OD, CM_OD, Q, H, tao, TR, Temp, elev, S, vel, knNA, dNAm, kdDBO, dDBO);

All input are single values, for example:

CNO3=0; CMNO3=0; Q=0.0107; tao=0.0217; TR=0.0048; Temp=13.4232; knNA=0.0641; dNAm=0; NO3_kdn=0.1000;
%Supose a dOD value of 7
dOD   = 7;

and the first function is:

function [dNO3, a_NO3, Foxdn, kdn] = fun_06NO3(CNO3, CMNO3, Q, tao, TR, Temp, knNA, dNAm, kdn, dOD)
kn_Na = knNA;
kdn20 = kdn;
tetha = 1.07;
kdn   = kdn20*tetha^(Temp-20);
Foxdn = exp(-kdn*dOD);
%%Main calculation
if dOD >= 0.6
    KNO3 = tao*kn_Na*dNAm;    
    NO3i = (CNO3 + CMNO3)/(Q*86400); 
    dNO3 = NO3i + KNO3 + TR*kn_Na*dNAm;    
    a_NO3 = Q*NO3i/dNO3;   
else
    KNO3 = kn_Na*dNAm*(1-exp(-kdn*Foxdn*tao))/(kdn*Foxdn);   
    NO3i = (CNO3 + CMNO3)/(Q*86400); 
    dNO3 = (NO3i*exp(-kdn*Foxdn*tao) + KNO3 + TR*kn_Na*dNAm)/(1 + TR*kdn*Foxdn);     
    a_NO3 = Q*NO3i*(1 + TR*kdn*Foxdn)/(NO3i*exp(-kdn*Foxdn*tao) + KNO3 + TR*kn_Na*dNAm);
end
end

How can I solve these functions simultaneosly?

Answer 1

[Torsten]

dOD0 = 7.0;
dOD = fzero(@fun,dOD0)

function res = fun(dOD)
  CNO3=0; CMNO3=0; Q=0.0107; tao=0.0217; TR=0.0048; Temp=13.4232; knNA=0.0641; dNAm=0; NO3_kdn=0.1000;
  [dNO3,a_NO3,Foxdn,kdn]= fun_06NO3(CNO3, CMNO3, Q, tao, TR, Temp, knNA, dNAm, NO3_kdn, dOD);
  [dDBO, a_DBO, kdDBO]  = fun_07DBO(CDBO, CMDBO, Q, tao, TR, Temp, Foxdn, kdn, dNO3, DBO_kd, dOD);
  [dOD_iter , a_D  , ka, Os] = fun_08_OD(C_OD, CM_OD, Q, H, tao, TR, Temp, elev, S, vel, knNA, dNAm, kdDBO, dDBO);
  res = dOD - dOD_iter;
end

[Albert Johan Mamani Larico]

Thank you, that help a lot Torsten.

Since al the script are in a loop, the inputs (CNO3, CMNO3, Q...) are changing in each iteration. So I tried to save them outside and then load inside the function, but that process decrease the velocity. Is there another way to do it without lossing performance?

save param_for_OD CNO3 CMNO3 CDBO CMDBO C_OD CM_OD Q tao TR Temp...
            knNA dNAm NO3_kdn DBO_kd H elev S vel
dOD0 = 7.0;
dOD = fzero(@fun,dOD0)
function res = fun(dOD)
load param_for_OD
  [dNO3,a_NO3,Foxdn,kdn]= fun_06NO3(CNO3, CMNO3, Q, tao, TR, Temp, knNA, dNAm, NO3_kdn, dOD);
  [dDBO, a_DBO, kdDBO]  = fun_07DBO(CDBO, CMDBO, Q, tao, TR, Temp, Foxdn, kdn, dNO3, DBO_kd, dOD);
  [dOD_iter , a_D  , ka, Os] = fun_08_OD(C_OD, CM_OD, Q, H, tao, TR, Temp, elev, S, vel, knNA, dNAm, kdDBO, dDBO);
  res = dOD - dOD_iter;
end

[Torsten]

You can pass the parameters needed in "fun" by modifying the call to "fzero" according to

dOD = fzero(@(x)fun(x,CNO3,CMNO3,CDBO,....),dOD0)

and changing "fun" to

function res = fun(dOD,CNO3,CMNO3,CDBO,...)

if that helps.

[Walter Roberson]

http://www.mathworks.com/help/matlab/math/parameterizing-functions.html

[Albert Johan Mamani Larico]

Great, both solutions work quite well. Thank you!

Answer 2

[Walter Roberson]

A different approach is to use the Symbolic Toolbox to arrive at formulas for the three different variables. solve() would probably not end up working for the formulas, but once you have explicit formulas in terms of the inputs, then you could matlabFunction() and pass that into fsolve()

In order to get at the formulas, you would have to rewrite your current functions in order to avoid using if . Here is an example:

function [dNO3, a_NO3, Foxdn, kdn] = fun_06NO3(CNO3, CMNO3, Q, tao, TR, Temp, knNA, dNAm, kdn, dOD)
kn_Na = knNA;
kdn20 = kdn;
tetha = 1.07;
kdn   = kdn20*tetha^(Temp-20);
Foxdn = exp(-kdn*dOD);

NO3i = (CNO3 + CMNO3)/(Q*86400); 

if dOD >= 0.6
    KNO3_high = tao*kn_Na*dNAm;    
    dNO3_high = NO3i + KNO3_high + TR*kn_Na*dNAm;    
    a_NO3_high = Q*NO3i/dNO3_high;   
else
    KNO3_low = kn_Na*dNAm*(1-exp(-kdn*Foxdn*tao))/(kdn*Foxdn);   
    dNO3_low = (NO3i*exp(-kdn*Foxdn*tao) + KNO3_low + TR*kn_Na*dNAm)/(1 + TR*kdn*Foxdn);     
    a_NO3_low = Q*NO3i*(1 + TR*kdn*Foxdn)/(NO3i*exp(-kdn*Foxdn*tao) + KNO3_low + TR*kn_Na*dNAm);
end

KNO3 = piecewise(dOd >= 0.6, KNO3_high, KNO3_low);
a_NO3 = piecewise(dOd >= 0.6, a_NO3_low);

end

You would pass in symbolic dOD.

You would make similar arrangements with the other functions. It will be okay that you have a function that calculates outputs in terms of the other functions: you might end up with complicated formulas, but for this purpose you are more concerned with getting explicit formulas not in terms of function calls -- effectively substituting in variables.

Now you encounter some tricky parts.

If you have formulas involving piecewise() then when you use matlabFunction(), you must use the 'File' option.

With current MATLAB releases (not sure about R2022a but the couple of years immediately before that for sure), if you use matlabFunction() with the 'File' option, I recommend using 'optimize', false which is not the default: optimization has been broken recently.

When you are working with piecewise() formulas, then unless your formulas were built carefully, their derivatives (or second derivatives) might not be continuous. And that is a problem for fsolve(), which assumes continuous derivatives. So you might not be able to use fsolve() on your piecewise() function.

What I would suggest, is that you get your piecewise()-based function, and then you create two versions of it:

assume(dOD >= 0.6)
form_high = simplify( TheFormula_in_dOD );
assume(dOD < 0.6)
form_low = simplify( TheFormula_in_dOD );

This would give you two formulas, neither of which should have any remaining piecewise() in them. You should then be able to matlabFunction() and fsolve() . Do both and choose the one that works.

You do run into the issue that you should be using bounds on dOD for each of those two hypothetical fsolve() calls, but that fsolve() does not explicitly include any way to impose bounds. It is possible to implement bounds with fsolve(), in the sense that you can ask fsolve() to stop when it hits the bounds. It is not obvious that you can do this, but it is possible by using the OutputFcn option, as that has a mechanism to request stopping.

A different way, instead of using fsolve, is to take the three individual formulas, F1, F2, F3, and calculate (F1^2 + F2^2 + F3^2) and matlabFunction() that. This transforms the zero finding task into a minimization task: a perfect set of solutions would generate a squared residue of 0. And that has the advantage that you could use fmincon() -- which supports bounds directly. You would still need to do the two cases, form_low and form_high, as fmincon() also expects continuous derivatives.

Or... you could calculate the residue that I show, but using the formula that still includes the piecewise(), and matlabFunction() 'file' and 'optimize', false that. You would then use ga() on it: ga() does not require continuous derivatives.