how to use fsolve with more equation

function F = fullfun(S_c,S_d,Ru,T,mu_s0,mu_w0,F,lambda_0,alpha,C_c,C_d,V_cp,r_am,r_cm,b_d,sigma_d,sigma_c,N_cp,h_r,sigma)
F = [T_s+(4*10^(-6)*S_d^(2))-(4*10^(-5)*S_d)-0.96;
    T_w+(4*10^(-5)*S_c^(2))+(1.9*10^(-2)*S_c)-11.2;
    L_s-(min(2*10^(-12)*S_d^(2)-3*10^(-10)*S_d+...
        6*10^(-8),2*10^(-12)*S_c^(2)-3*10^(-10)*S_c+6*10^(-8)));
    L_w-(5*10^(-4)*S_c^(-0.416));
    
    mu_sc-(mu_s0+Ru*T*log(S_c/1000));
    mu_sd-(mu_s0+Ru*T*log(S_d/1000));
    mu_wc-(mu_w0+Ru*T*log((1000-S_c)/1000));
    mu_wd-(mu_w0+Ru*T*log((1000-S_d)/1000));
    
    E_m-(T_s/F*(mu_sc-mu_sd)+T_w/F*(mu_wc-mu_wd));
    E_am-E_m; 
    E_cm-E_m;
    
    lambda_c-lambda_0+alpha*sqrt(C_c);
    lambda_d-lambda_0+alpha*sqrt(C_d);
    lambda_r-((lambda_c + lambda_d)/2); 
    C_r-((C_c+C_d)/2);
    k-(lambda_r * C_r);
    
     I-((V_cp-E_am-E_cm)/(r_am+r_cm+b_d/(sigma_d*lambda_d*C_d)+(b_c/(sigma_c*lambda_c*C_c))+...
         r_cm/N_cp+2*h_r/(sigma*k*N_cp)));
        
    J_s-(T_s*I/F - L_s*(C_c-C_d));
    J_w-(T_w*I/F + L_w*Ru*T*(C_c-C_d)*10^(-5)*2);

Answer 1

[Walter Roberson]

Bundle all of the variables into a single vector for the purposes of fsolve. You can unbundle them inside the function.

vars20init = vector of 20 initial values
sol = fsolve(@fullfun, vars20init)

function result = fullfun(vars20)
   tcell = num2cell(vars20);
   [S_c,S_d,Ru,T,mu_s0,mu_w0,F,lambda_0,alpha,C_c,C_d,V_cp,r_am,r_cm,b_d,sigma_d,sigma_c,N_cp,h_r,sigma] = tcell{:};
   result = [T_s+(4*10^(-6)*S_d^(2))-(4*10^(-5)*S_d)-0.96;
    T_w+(4*10^(-5)*S_c^(2))+(1.9*10^(-2)*S_c)-11.2;
    L_s-(min(2*10^(-12)*S_d^(2)-3*10^(-10)*S_d+...
        6*10^(-8),2*10^(-12)*S_c^(2)-3*10^(-10)*S_c+6*10^(-8)));
    L_w-(5*10^(-4)*S_c^(-0.416));
    
    mu_sc-(mu_s0+Ru*T*log(S_c/1000));
    mu_sd-(mu_s0+Ru*T*log(S_d/1000));
    mu_wc-(mu_w0+Ru*T*log((1000-S_c)/1000));
    mu_wd-(mu_w0+Ru*T*log((1000-S_d)/1000));
    
    E_m-(T_s/F*(mu_sc-mu_sd)+T_w/F*(mu_wc-mu_wd));
    E_am-E_m; 
    E_cm-E_m;
    
    lambda_c-lambda_0+alpha*sqrt(C_c);
    lambda_d-lambda_0+alpha*sqrt(C_d);
    lambda_r-((lambda_c + lambda_d)/2); 
    C_r-((C_c+C_d)/2);
    k-(lambda_r * C_r);
    
     I-((V_cp-E_am-E_cm)/(r_am+r_cm+b_d/(sigma_d*lambda_d*C_d)+(b_c/(sigma_c*lambda_c*C_c))+...
         r_cm/N_cp+2*h_r/(sigma*k*N_cp)));
        
    J_s-(T_s*I/F - L_s*(C_c-C_d));
    J_w-(T_w*I/F + L_w*Ru*T*(C_c-C_d)*10^(-5)*2);]
end

I notice you use F as both input and output. That is unlikely to be a good idea when you are calling fsolve with the function.

[chirag patel]

in m main i have used this :

vars20init = [S_c,S_d,Ru,T,mu_s0,mu_w0,F,lambda_0,alpha,C_c,C_d,V_cp,r_am,r_cm,b_d,sigma_d,sigma_c,N_cp,h_r,sigma];

sol = fsolve(@fullfun, vars20init)

when i use I got following Error :

Unrecognized function or variable 'T_s'.

Error in fullfun (line 4)

result = [T_s+(4*10^(-6)*S_d^(2))-(4*10^(-5)*S_d)-0.96;

Error in fsolve (line 248)

fuser = feval(funfcn{3},x,varargin{:});

Error in main (line 14)

sol = fsolve(@fullfun, vars20init)

Caused by:

Failure in initial objective function evaluation. FSOLVE cannot continue.

[Walter Roberson]

You use T_s in your function, but it is not a parameter to the function. Where do you want your code to find T_s from?

Answer 2

[chirag patel]

I understand your question and modify code as follows :

vars20init = [70,70,1.2875e+03,1.2875e+03,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];
 sol = fsolve(@fullfun, vars20init)

function result = fullfun(vars20)
   tcell = num2cell(vars20);
   [S_c,S_d,C_c,C_d,T_w,T_s,L_s,L_w,mu_sc,mu_sd,mu_wc,mu_wd,E_m,E_am,E_cm,lambda_c,lambda_d,lambda_r,k,C_r,I,J_s,J_w] = tcell{:};
   result = [T_s+(4*10^(-6)*S_d^(2))-(4*10^(-5)*S_d)-0.96;
    T_w+(4*10^(-5)*S_c^(2))+(1.9*10^(-2)*S_c)-11.2;
    L_s-(min(2*10^(-12)*S_d^(2)-3*10^(-10)*S_d+...
        6*10^(-8),2*10^(-12)*S_c^(2)-3*10^(-10)*S_c+6*10^(-8)));
    L_w-(5*10^(-4)*S_c^(-0.416));
    
    mu_sc-(-732*10^3+8.314*298.15*log(S_c/1000));
    mu_sd-(-732*10^3+8.314*298.15*log(S_d/1000));
    mu_wc-(-237*10^3+8.314*298.15*log((1000-S_c)/1000));
    mu_wd-(-237*10^3+8.314*298.15*log((1000-S_d)/1000));
    
    E_m-(T_s/96485*(mu_sc-mu_sd)+T_w/96485*(mu_wc-mu_wd));
    E_am-E_m; 
    E_cm-E_m;
    
    lambda_c-0.012639+2.444*10^-4*sqrt(C_c);
    lambda_d-0.012639+2.444*10^-4*sqrt(C_d);
    lambda_r-((lambda_c + lambda_d)/2); 
    C_r-((C_c+C_d)/2);
    k-(lambda_r * C_r);
    
     I-((0.56-E_am-E_cm)/(3.5*10^(-4)+3.5*10^(-4)+5*10^(-04)/(0.64*lambda_d*C_d)+(5*10^(-4)/(0.64*lambda_c*C_c))+...
         3.5*10^(-4)/1+2*5*10^(-04)/(0.64*k*1)));
        
    J_s-(T_s*I/96485 - L_s*(C_c-C_d));
    J_w-(T_w*I/96485 + L_w*8.314*298.15*(C_c-C_d)*10^(-5)*2);];
end

run but got answer as follows :

Warning: Trust-region-dogleg algorithm of FSOLVE cannot handle non-square systems; using Levenberg-Marquardt algorithm

instead.

> In fsolve (line 316)

In main (line 14)

Equation solved, solver stalled.

fsolve stopped because the relative size of the current step is less than the

value of the step size tolerance and the vector of function values

is near zero as measured by the value of the function tolerance.

<stopping criteria details>

sol =

1.0e+05 *

Columns 1 through 6

0.3139 - 0.0019i 0.8325 - 0.0078i 0.0086 + 0.0000i 0.0086 + 0.0000i -0.4001 + 0.0047i -0.2772 + 0.0052i

Columns 7 through 12

0.0000 - 0.0000i 0.0000 + 0.0000i -7.2346 - 0.0001i -7.2104 - 0.0002i -2.2854 + 0.0777i -2.2607 + 0.0776i

Columns 13 through 18

0.0172 - 0.0003i 0.0172 - 0.0003i 0.0172 - 0.0003i 0.0000 - 0.0000i 0.0000 - 0.0000i 0.0000 - 0.0000i

Columns 19 through 23

-0.0000 - 0.0000i 0.0086 + 0.0000i 0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 + 0.0000i

>>

[Walter Roberson]

Equation solved, solver stalled.

That means fsolve worked and found an answer to you that was within the constraints you specified. Not all of the function outputs were exactly zero, so it kept looking for a better solution, trying smaller and smaller steps near the locations it had arrived at, but was unable to find any improvement by the time it reached the minimum step size you had configured.

Is the solution perfect ? No, it is not. But it is within the tolerances you set.

You should never expect perfect solutions, because of round-off and truncation effects. For example you cannot even solve x^2-2 perfectly because of those effects.

>> 2-sqrt(2)^2

ans =

     -4.44089209850063e-16

[chirag patel]

Atualyy i want to solve following using fsolve :

function [S_c,S_d,mdot_c,mdot_d,I,J_s,J_w,energy_kWhpm3,power]=algo(S_c,S_d,Ru,T,mu_s0,mu_w0,F,lambda_0,...
    alpha,C_c,C_d,V_cp,r_am,r_cm,b_c,b_d,sigma_d,sigma_c,N_cp,h_r,sigma,dA,...
    mdot_c,mdot_d,ndot_sc,ndot_sd,N,MW_s,MW_w,rho_w,A_cp,mdot_cin,mdot_din)

for i=1:N
    T_s(i) = -4*10^(-6)*S_d(i)^(2)+4*10^(-5)*S_d(i)+0.96;
    T_w(i) = -4*10^(-5)*S_c(i)^(2)-1.9*10^(-2)*S_c(i)+11.2;
    L_s(i) = min(2*10^(-12)*S_d(i)^(2)-3*10^(-10)*S_d(i)+...
        6*10^(-8),2*10^(-12)*S_c(i)^(2)-3*10^(-10)*S_c(i)+6*10^(-8));
     L_s(i) = 2*10^(-12)*S_d(i)^(2)-3*10^(-10)*S_d(i)+6*10^(-8);
    L_w(i) = 5*S_c(i)^(-0.416)*10^(-4);
    
    mu_sc(i) = mu_s0+Ru*T*log(S_c(i)/1000); mu_sd(i) = mu_s0+Ru*T*log(S_d(i)/1000);
    mu_wc(i) = mu_w0+Ru*T*log((1000-S_c(i))/1000); mu_wd(i) = mu_w0+Ru*T*log((1000-S_d(i))/1000);
    
    E_m(i) = T_s(i)/F*(mu_sc(i)-mu_sd(i))+T_w(i)/F*(mu_wc(i)-mu_wd(i));
    E_am(i) = E_m(i)/2; E_cm(i) = E_m(i)/2;
    
    lambda_c = lambda_0-alpha*sqrt(C_c(i));

    lambda_d = lambda_0-alpha*sqrt(C_d(i));
    lambda_r = (lambda_c + lambda_d)/2; C_r = (C_c(i)+C_d(i))/2;
    k = lambda_r * C_r;
    
     I(i) = (V_cp-E_am(i)-E_cm(i))/(r_am(i)+r_cm(i)+b_d/(sigma_d*...
         lambda_d*C_d(i))+(b_c/(sigma_c*lambda_c*C_c(i)))+...
         r_cm(i)/N_cp+2*h_r/(sigma*k*N_cp));
    

    
    J_s(i) = T_s(i)*I(i)/F - L_s(i)*(C_c(i)-C_d(i));
    J_w(i) = T_w(i)*I(i)/F + L_w(i)*Ru*T*(C_c(i)-C_d(i))*10^(-5)*2;
    
    ndot_sc(i+1) = ndot_sc(i) + dA*J_s(i);
    ndot_sd(i+1) = ndot_sd(i) - dA*J_s(i);

    
    mdot_c(i+1) = mdot_c(i)+(J_s(i)*dA*MW_s+J_w(i)*dA*MW_w);
    mdot_d(i+1) = mdot_d(i)-(J_s(i)*dA*MW_s+J_w(i)*dA*MW_w);
    S_c(i+1) = ndot_sc(i+1)*1000*MW_s/mdot_c(i+1);
    S_d(i+1) = ndot_sd(i+1)*1000*MW_s/mdot_d(i+1);
    rho_c(i+1) = (1000+S_c(i+1))/1000*rho_w;            % S_c(i+1)/1000*rho_nacl+
    rho_d(i+1) = (1000+S_d(i+1))/1000*rho_w;            % S_d(i+1)/1000*rho_nacl+
    C_c(i+1) = S_c(i+1)*rho_c(i+1)/(1000*MW_s);
    C_d(i+1) = S_d(i+1)*rho_d(i+1)/(1000*MW_s);
    
end

[Walter Roberson]

Okay, go ahead. Follow the same strategy as before: all inputs should go into a single vector, and when you get the vector output from the solver, break it up into meaningful variables. Inside the function to be solved, expect everything as a single input, and break the values out into meaningful variables (this is what the num2cell and {:} is for)