How can I solve a system of non-linear equations for positive values only?

Hello,

I wrote a code to solve a system of non-linear equations, but I am not getting the right answer when I run it. I thought the initial guess was the problem, so I used experimental values from a research paper, but the solution is not even close it. Additionally, I get a comment that the solver was stopped prematurely. I am not very experienced in Matlab. Could someone help me find what I am doing wrong?

function f = root(x)
%data
delta = [0.1]*1e-4; QA = 3.57; QB = 20; QC = 60; QD = 1000; ph = [9.86923]*76; pl = [0.986923]*76;               
nF = 277778; yAF = 0.7841; yBF = 0.2084; yCF = 0.0003; yDF = 1-yAF-yBF-yCF; A = 2260000;
% Material Balance
f(1) = x(13)*x(1) + x(14)*x(5) - nF*yAF;
f(2) = x(13)*x(2) + x(14)*x(6) - nF*yBF;
f(3) = x(13)*x(3) + x(14)*x(7) - nF*yCF;
f(4) = x(13)*x(4) + x(14)*x(8) - nF*yDF;
% Material Balance in terms of flux
f(5) = x(9)*A - x(13)*x(1);
f(6) = x(10)*A - x(13)*x(2);
f(7) = x(11)*A - x(13)*x(3);
f(8) = x(12)*A - x(13)*x(4);
% Flux Equations
f(9) = x(9) - (QA/delta)*(ph*x(5) - pl*x(1));
f(10) = x(10) - (QB/delta)*(ph*x(6) - pl*x(2));
f(11) = x(11) - (QC/delta)*(ph*x(7) - pl*x(3));
f(12) = x(12) - (QD/delta)*(ph*x(8) - pl*x(4));
% Component Equations
f(13) = x(5) + x(6) + x(7) + x(8) - 1;
f(14) = x(1) + x(2) + x(3) + x(4) - 1;
end

followed by:

x0 = [0.42, 0.48, 0.05, 0.05, 0.825, 0.075, 0.05, 0.05, 0.00258, 0.00295, 0.00031, 00031, 13888.9, 263889]';

[x, fval] = fsolve(@root, x0) 

The output I get is:

Solver stopped prematurely.

fsolve stopped because it exceeded the function evaluation limit,
options.MaxFunctionEvaluations = 1.400000e+03.

x =

   1.0e+11 *

    0.0000
    0.0000
    0.0000
    0.0000
    0.0000
    0.0000
    0.0000
    0.0000
    0.0000
    0.0000
    0.0000
    0.0000
    1.6223
   -1.9271

Answer 1

[Ayush Gupta]

fsolve doesn’t allow bounds but lsqnonlin does. Change the last piece of code to get answers for positive values only.

x = lsqnonlin(@root,x0,zeros(size(x0))) 

[Walter Roberson]

What is @a2? I would have expected @(x)root(x).^2 instead of @a2

[Julia S]

Thank you so much!