Non-linear Algebraic 36 equations unsloved

Hi forum

I have 36 Non-linear Algebraic equations with 36 unknowns to find out.

'solve' comand failed.

'fsolve' comand requires initial choices which I dont have.

Any idea, how to solve? Please help

[Torsten]

An idea ? Without the equations ?

Usually, solving 36 equations without no idea for a good initial guess is hopeless. But why not trying with x0 = ones(36,1) and "fsolve" ?

Answer 1

[John D'Errico]

Actually, a starting value set like ones(36,1) is often a dangerous thing to do. Making all of the starting values the same can sometimes put the solver in a state where it cannot proceed, due to symmetries in the problem. A better choice of simple random numbers is a far better idea, perhaps rand(36,1), or randn(36,1).

Is this likely to be sufficient? Possibly, but that is impossible to know wothout even a clue at your actual equations. It MAY yield a solution, but you need to recognize that without a decent start, a solver can often get lost, and diverge to some bad place.

Worse, even if you use a random start point, this may be insufficient if some of the parameters have very different scalings than the others.

Yes, it is absolutely no surprise an analytical solver like solve failed. In fact, that is far more the expected result than to find a solution from solve on such a problem. (Though I once recall a nonlinear problem with 15 unknowns I posed, where solve found the correct solution, but that itself was a great surprise to me.) As such, you almost certainly will need to use a numerical solver like fsolve. The problem then lies in how to choose a starting set. And if you understand nothing about the parameters in the model to be able to provide an intelligent solution, well, expect random crapola much of the time.

Think of such a problem as if you were to put a blind person down on the face of the earth, and ask them to find the point of lowest altitude, given nothing more than a cane, and an altimeter. (Yes, an altimeter that works underwater, so they may also want scuba gear.) If you start them out in some completely random spot, what are the odds they will find the actual minimum point? Far more likely they may get stuck in a locally sub-optimal point, where there is no direction they can move that will allow them to go downhill anymore.

And I am sorry to say, there is no magical way to automatically choose an intelligent set of starting values. What can you do? You may be left with two alternatives.

1. Use a multi-start method. That has you start a solver like fsolve off with a randomly generated point. Save the result you get. Then repeat. Some of the time, you will get crap, non-convergence, whatever. But save the set of solutions where fsolve thinks it did converge, and choose from among that set of "results", taking the one that makes you happy.
2. Use a solver like GA, or one of the many other tools from the particle swarm family, or a simulated annealing tool, etc. These are tools that are not so heavily tied down to a start point. I'd start with GA, a genetic algorithm, as my personal suggestion.

[MINATI PATRA]

@ Dear John D'Errico

have a try

Answer 2

[Alex Sha]

@MINATI PATRA please giving out the code you used, or the detail description of your 36 equations

[MINATI PATRA]

%% Sorry for the late reply, Actually my Laptop was out of order, Please have a try with this code

tic
syms x N1 N2 N3 Re Peh Pem
N1 = 0.1;N2 = 0.1; N3 = 0.1; Re = 0.1; Peh = 0.1; Pem = 0.1;
A = sym('a', [1 10]);B = sym('b', [1 10]);C = sym('c', [1 8]);D = sym('d', [1 8]); j = 1:10; k = 1:8;
f = sum(A.*x.^(j-1));g = sum(B.*x.^(j-1)); t = sum(C.*x.^(k-1));p = sum(D.*x.^(k-1));
Df = diff(f); D2f = diff(Df); D3f = diff(D2f); D4f = diff(D3f); D1g = diff(g); D2g = diff(D1g); D1p = diff(p); D2p = diff(D1p); D1t = diff(t); D2t = diff(D1t); 
fP = subs(f,x,1); fN = subs(f,x,-1); fdP = subs(Df,x,1); fdN = subs(Df,x,-1);gP = subs(g,x,1); gN = subs(g,x,-1);tP = subs(t,x,1); tN = subs(t,x,-1);pP = subs(p,x,1); pN = subs(p,x,-1); 
F = (1+N1)*D4f - N1*D2g - Re*(f*D3f - Df*D2f); G = N2*D2g - N1*(D2f - 2*g) - N3*Re*(f*D2g - Df*g);T = D2t + Peh*(Df*t - f*D1t); P = D2p + Pem*(Df*p - f*D1p); 
DF = diff(F); D2F = diff(DF); D3F = diff(D2F); D4F = diff(D3F); DG = diff(G); D2G = diff(DG);D3G = diff(D2G); DP = diff(P); D2P = diff(DP); DT = diff(T); D2T = diff(DT); 

FP = subs(F,x,1); FN = subs(F,x,-1); FdP = subs(DF,x,1); FdN = subs(DF,x,-1); Fd2P = subs(D2F,x,1); Fd2N = subs(D2F,x,-1);
GP = subs(G,x,1); GN = subs(G,x,-1); GdP = subs(DG,x,1); GdN = subs(DG,x,-1); Gd2P = subs(D2G,x,1); Gd2N = subs(D2G,x,-1);Gd3P = subs(D3G,x,1); Gd3N = subs(D3G,x,-1);
TP = subs(T,x,1); TN = subs(T,x,-1); TdP = subs(DT,x,1); TdN = subs(DT,x,-1); Td2P = subs(D2T,x,1); Td2N = subs(D2T,x,-1);
PP = subs(P,x,1); PN = subs(P,x,-1); PdP = subs(DP,x,1); PdN = subs(DP,x,-1); Pd2P = subs(D2P,x,1); Pd2N = subs(D2P,x,-1);
e1 = fN == 0; e2 = fP == 0; e3 = fdN == 0; e4 = fdP == -1; e5 = gN == 0; e6 = gP == 1; e7 = pN == 1; e8 = pP == 0; e9 = tN == 1; e10 = tP == 0;
e11 = FN == 0; e12 = FP == 0; e13 = FdN == 0; e14 = FdP == 0;e15 = Fd2N == 0; e16 = Fd2P == 0;
e17 = GN == 0; e18 = GP == 0; e19 = GdN == 0; e20 = GdP == 0;e21 = Gd2N == 0; e22 = Gd2P == 0;e23 = Gd3N == 0; e24 = Gd3P == 0;
e25 = TN == 0; e26 = TP == 0; e27 = TdN == 0; e28 = TdP == 0; e29 = Td2N == 0; e30 = Td2P == 0;
e31 = PN == 0; e32 = PP == 0; e33 = PdN == 0; e34 = PdP == 0;e35 = Pd2N == 0; e36 = Pd2P == 0;

EQN = [e1 e2 e3 e4 e5 e6 e7 e8 e9 e10 e11 e12 e13 e14 e15 e16 e17 e18 e19 e20 e21 e22 e23 e24 e25 e26 e27 e28 e29 e30 e31 e32 e33 e34 e35 e36];
 VARS = [A B C D];  
S = solve(EQN);
toc
 S.a1

Answer 3

[Torsten]

Seems your system has several solutions.

syms x N1 N2 N3 Re Peh Pem
N1 = 0.1;N2 = 0.1; N3 = 0.1; Re = 0.1; Peh = 0.1; Pem = 0.1;
A = sym('a', [1 10]);B = sym('b', [1 10]);C = sym('c', [1 8]);D = sym('d', [1 8]); j = 1:10; k = 1:8;
f = sum(A.*x.^(j-1));g = sum(B.*x.^(j-1)); t = sum(C.*x.^(k-1));p = sum(D.*x.^(k-1));
Df = diff(f); D2f = diff(Df); D3f = diff(D2f); D4f = diff(D3f); D1g = diff(g); D2g = diff(D1g); D1p = diff(p); D2p = diff(D1p); D1t = diff(t); D2t = diff(D1t); 
fP = subs(f,x,1); fN = subs(f,x,-1); fdP = subs(Df,x,1); fdN = subs(Df,x,-1);gP = subs(g,x,1); gN = subs(g,x,-1);tP = subs(t,x,1); tN = subs(t,x,-1);pP = subs(p,x,1); pN = subs(p,x,-1); 
F = (1+N1)*D4f - N1*D2g - Re*(f*D3f - Df*D2f); G = N2*D2g - N1*(D2f - 2*g) - N3*Re*(f*D2g - Df*g);T = D2t + Peh*(Df*t - f*D1t); P = D2p + Pem*(Df*p - f*D1p); 
DF = diff(F); D2F = diff(DF); D3F = diff(D2F); D4F = diff(D3F); DG = diff(G); D2G = diff(DG);D3G = diff(D2G); DP = diff(P); D2P = diff(DP); DT = diff(T); D2T = diff(DT); 
FP = subs(F,x,1); FN = subs(F,x,-1); FdP = subs(DF,x,1); FdN = subs(DF,x,-1); Fd2P = subs(D2F,x,1); Fd2N = subs(D2F,x,-1);
GP = subs(G,x,1); GN = subs(G,x,-1); GdP = subs(DG,x,1); GdN = subs(DG,x,-1); Gd2P = subs(D2G,x,1); Gd2N = subs(D2G,x,-1);Gd3P = subs(D3G,x,1); Gd3N = subs(D3G,x,-1);
TP = subs(T,x,1); TN = subs(T,x,-1); TdP = subs(DT,x,1); TdN = subs(DT,x,-1); Td2P = subs(D2T,x,1); Td2N = subs(D2T,x,-1);
PP = subs(P,x,1); PN = subs(P,x,-1); PdP = subs(DP,x,1); PdN = subs(DP,x,-1); Pd2P = subs(D2P,x,1); Pd2N = subs(D2P,x,-1);
e1 = fN == 0; e2 = fP == 0; e3 = fdN == 0; e4 = fdP == -1; e5 = gN == 0; e6 = gP == 1; e7 = pN == 1; e8 = pP == 0; e9 = tN == 1; e10 = tP == 0;
e11 = FN == 0; e12 = FP == 0; e13 = FdN == 0; e14 = FdP == 0;e15 = Fd2N == 0; e16 = Fd2P == 0;
e17 = GN == 0; e18 = GP == 0; e19 = GdN == 0; e20 = GdP == 0;e21 = Gd2N == 0; e22 = Gd2P == 0;e23 = Gd3N == 0; e24 = Gd3P == 0;
e25 = TN == 0; e26 = TP == 0; e27 = TdN == 0; e28 = TdP == 0; e29 = Td2N == 0; e30 = Td2P == 0;
e31 = PN == 0; e32 = PP == 0; e33 = PdN == 0; e34 = PdP == 0;e35 = Pd2N == 0; e36 = Pd2P == 0;

EQNS = [e1 e2 e3 e4 e5 e6 e7 e8 e9 e10 e11 e12 e13 e14 e15 e16 e17 e18 e19 e20 e21 e22 e23 e24 e25 e26 e27 e28 e29 e30 e31 e32 e33 e34 e35 e36];
EQNS = lhs(EQNS)-rhs(EQNS);
VARS = [A B C D];
EQNS = matlabFunction(EQNS,"Vars",{VARS});
VARS0 = rand(1,36);
VARS = fsolve(EQNS,VARS0)

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

VARS = 1x36

3.6725 0.0906 -7.4061 0.1528 4.0630 -0.3481 -0.3478 0.1256 0.0184 -0.0208 -520.1053 79.3286 666.5327 -148.6066 -177.2206 97.7998 35.9712 -31.6665 -4.6780 3.6447 0.6472 -0.5924 -0.1550 0.1378 -0.0042 -0.0492 0.0120 0.0037 0.6472 -0.5924

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

A = VARS(1:10)

A = 1x10

3.6725 0.0906 -7.4061 0.1528 4.0630 -0.3481 -0.3478 0.1256 0.0184 -0.0208

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

B = VARS(11:20)

B = 1x10

-520.1053 79.3286 666.5327 -148.6066 -177.2206 97.7998 35.9712 -31.6665 -4.6780 3.6447

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

C = VARS(21:28)

C = 1x8

0.6472 -0.5924 -0.1550 0.1378 -0.0042 -0.0492 0.0120 0.0037

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

D = VARS(29:36)

D = 1x8

0.6472 -0.5924 -0.1550 0.1378 -0.0042 -0.0492 0.0120 0.0037

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

norm(EQNS(VARS))

ans = 4.1370e-11

[MINATI PATRA]

@Thanks Torsten, I got the result

[Torsten]

Note that each run with initial conditions changed seems to give a different solution. I don't know if this is to be expected for your underlying problem.

[MINATI PATRA]

Ok, I will take care of that. Thanks for your concern.