No solution found. fsolve stopped because the problem appears to be locally singular.

%equilibrium composition

function F=fsolve1_function(X)

%unpack variables

E2=X(1);

E4=X(2);

E5=X(3);

E6=X(4);

noi=X(5);%initial moles of oxygen

K23=167.86;

K4=4.669E20;

K5=1.29;

K6=3.44E19;

yo2eq=3.75976E-21;%atm

%input equations

eq1=((noi+E4-E2-E6)/(E5-E2))-K23;

eq2=(((2*E4-E5)^2)/((noi+E4-E2-E6)*(noi-E4-E6)))-K4;

eq3=((E5-E2)*(E5-2*E6)/(2*E4-E5)*(2*E6-E5))-K5;

eq4=(((2*E6-E5)^2)*(noi+E4-E2-E6)/((E5-2*E6)^2*(noi-E4-E6)))-K6;

eq5=(noi-E4-E6)/(noi+E4-E2-E6)-yo2eq;

F=[eq1; eq2; eq3; eq4; eq5 ; ];

E=[E2; E4; E5; E6; noi]

end

%%script

X0 = [0.009, 0.008, 0.007, 0.006, 0.01 ];

%Call fsolve function

fhandle = @fsolve1_function;

opts=optimoptions('fsolve','Display','iter','TolFun',1e-20,'TolX',1e-20);

%opts=optimoptions('fsolve','TolFun',1e-6);

[X,fval,exitflag] = fsolve (fhandle, X0,opts);

How do I solve these nonlinear equations with 5 unknons? Is there any other function that I can try?

Answer 1

[John D'Errico]

You have a serious problem here, as I think you are having a problem using double precision arithmetic with this dynamic range of numbers.

syms E2 E4 E5 E6 noi
K23=167.86;
K4=4.669E20;
K5=1.29;
K6=3.44E19;
yo2eq=3.75976E-21;%atm
%input equations
eq1=((noi+E4-E2-E6)/(E5-E2))-K23

eq1 = 

eq2=(((2*E4-E5)^2)/((noi+E4-E2-E6)*(noi-E4-E6)))-K4

eq2 = 

eq3=((E5-E2)*(E5-2*E6)/(2*E4-E5)*(2*E6-E5))-K5

eq3 = 

eq4=(((2*E6-E5)^2)*(noi+E4-E2-E6)/((E5-2*E6)^2*(noi-E4-E6)))-K6

eq4 = 

eq5=(noi-E4-E6)/(noi+E4-E2-E6)-yo2eq

eq5 = 

Do you see that two of those equations are equal to numbers on the order of 10^20? At the same time, another equation wants to be on the order of 10^-21?

I'm sorry, but that is just a recipe for immediate failure. Even if a solution exists, fsolve will NEVER, NEVER, NEVER find it. NEVER. NOT EVER. PERIOD.

Was I firm enough? Double precision arithmetic has a dynamic range of around 16 decimal digits. But your equations span more than 40 powers of 10.

The solution will be one of two options. You must learn to rescale the problem. Learn to use units properly. This is a big reason why we have units like light years, and femotometers. It makes things with huge variations in dynamic range palatable. It turns big or small numbers into things we can deal with.

Or, you can go the symbolic route. Symbolic computations are often very slow, and they have their own issues. But sometimes, maybe it works.

X0 = [0.009, 0.008, 0.007, 0.006, 0.01];
vpasolve([eq1,eq2,eq3,eq4,eq5],[E2 E4 E5 E6 noi],X0)

ans = struct with fields:
     E2: 0
     E4: 0
     E5: 0
     E6: 0
    noi: 0

So vpasolve finds a solution at zero. Which, in fact, if we look at your equations, makes sense. Does a non-trivial solution exist? Maybe, maybe not. You will not find it with fsolve easily, not without some serious effort.

[Walter Roberson]

if you subs the 0s back in you get division by 0.

Answer 2

[Walter Roberson]

If you have the symbolic toolbox,

syms x [1 5]

and calculate through to F. Then

s = solve(F(1:3),X(1:3))

This will give you 4 closed form solution for X1 X2 X3 in terms of X4 X5

Now if you subs those solutions into F(4:5) and simplify() it turns out that you will get expressions that are independent of X4 X5... numeric expressions with no variables. Expressions that are not equal to 0.

Your expressions are rank deficient but also impossible to satisfy.

[Mahe Rukh]

Thanks for your response. Would it be possible to write the commands in more detail? Your expressions are rank deficient but also impossible to satisfy- does it mean equations are not sovable?