8 nonlinear equations, 8 unknowns - little luck with fsolve/vpasolve

Hi All. I have a system of 8 nonlinear equations that I try to solve with fsolve. If I run the code, I get the solution because I start at the actual solutions. But If I change the guess values (g1,g2,etc) slightly, I get errors, like no solutions found, or fslove stopped prematurely. etc. I tried vpasolve, -- it seems to perform even worse than fslove, returning empty sets. I tried random guess values in fsolve, but it takes too long and I cannot see any solution still. I understand no solver or method is perfect. But perhaps I thought there is an inefficiency in my code? Any help is appreciated. I tried to convert the paramaters to ratios as suggested in earlier threads, -- no luck. Is there a way to solve this problem more efficiently instead of suggesting to give up on this, or research advanced numerical techniques?

[Matt J]

Your TolFun setting of 0.001 looks suspiciously lax. It will encourage the solver to stop iterating very early. I doubt it's the only cause of your problems, but I would go to TolFun=1e-10 or so.

[Alex Tekin]

Dear Matt J, thank you. I have chnaged that setting but to no avail. I received an error: "Solver stopped prematurely. fsolve stopped because it exceeded the function evaluation limit, options.MaxFunctionEvaluations = 100000 (the selected value)."

Answer 1

[Alex Sha]

Hi, Alex Tekin, you may try some solvers or software package with the ability of global optimization, i.e. Baron, Lingo, 1stOpt. the belows are the results obtained from 1stOpt:

1:

x1	4.5153251734112
x2	1.15089340375541
x3	1.02455790981764
x4	0.00214722268127572
x5	0.100000000000004
x6	4.07792365160064
x7	0.0252189443492125
x8	0.330004160606675

2:

x1	2.75653453284178
x2	1.2705558348226
x3	1.54646959816705
x4	0.000816371359806575
x5	-0.894039944060641
x6	1.88605932702964
x7	0.0545283992510777
x8	0.330003462902014

3:

x1	0.817677432070711
x2	-0.734946906055547
x3	-5.97395677713386
x4	-7.07117349295853E-5
x5	-0.894039944060641
x6	0.282420925065193
x7	0.364258018370474
x8	0.329996089923976

4:

x1	5.57459460098411
x2	0.321268890484856
x3	0.205761174232818
x4	0.000833138706775475
x5	0.100000000000003
x6	5.66820754494582
x7	0.018143343701659
x8	0.330004329039833

5:

x1	0.360923161734683
x2	-0.0228328849580109
x3	-0.666063660048423
x4	-6.12135370250361E-7
x5	-0.894039944060643
x6	0.0786951613776357
x7	1.3084212647888
x8	0.329973615235634

6:

x1	2.145899354053
x2	4.9676878398387
x3	14.1405959014085
x4	0.00289857324146092
x5	0.100000000000002
x6	1.27534324528113
x7	0.0806420939672824
x8	0.330002841272532

[Alex Tekin]

Dear Alex Sha, many thanks for your help. By any chance, can you please send me a link to 1stOpt? I suppose it is a software different from MATLAB, it is not a code that can be executed in MATLAB? Also, does it take a long time for 1stOpt to solve a problem like mine, and I suppose the results are reliable? Also, any of these have a trial version? I suppose they are expensive to buy outright. Thank you.

[Alex Sha]

Hi, Tekin, refer the linker below:

http://www.7d-soft.com/

or

http://www.7d-soft.com/en/

[Alex Tekin]

Thank you Alex Sha. I am wondering, those 6 sets of solutions you have sent me, are these the only solutions that 1stOpt would find, or there are potentially many more? Do we have to change like guess values in 1stOpt in order to trace "all" the solutions? Thank you.

[Alex Sha]

Hi, 1stOpt not require end-uesr to guess or supply initial start-values. There are more solutions, for example, for all positive solution:

1:

x1: 1.13746211137498
x2: 3356.08384248914
x3: 25756.3516424922
x4: 0.72631462758202
x5: 0.0999999754181184
x6: 0.473030282935766
x7: 0.217448507176776
x8: 0.329999584648149

2:

x1: 3.50509495101581
x2: 2.69715606884254
x3: 3.56669260501346
x4: 0.00338758240055768
x5: 0.0999999999991526
x6: 2.74524365235143
x7: 0.0374619305147155
x8: 0.330003869164897

...

every run take about 10s, depending in your PC.

[Alex Tekin]

Thank you Alex Sha. Did you have to do any reformulation to the problem? Or you just plugged in the original expressions, and got the bunch of solutions? Because I am trying BARON now; it is perfomring better than fsolbe but struggles a lot. Maybe 1stOpt is simply superior software?

[Alex Sha]

Hi, Tekin, not any reformulation needed for your above equation problem.

Baron is a really good solver espicially on global optimization, however, proviately, 1stOpt is better: easy for using, and often produce much better results compared with others. You may refer a case:

https://www.mathworks.com/matlabcentral/answers/423118-i-want-to-solve-this-equation-please-help-me

[Alex Tekin]

Hello, Alex Sha, many thanks! I referred to 1stOpt people, but they said I need to do some reformulation, which makes things comlicated, like making division as multiplication etc. They sent me some solutions, and I plugged in their numbers -- errors are big. But yours, -- perfect solutions! I've read the thread you sent me, there were 4 equations, and it was like MATLAB would struggle due to double-precision. But you supplied the numbers (a,b,c,d), and they fit perfectly! So, just like that, you plug in as is in 1stOpt and get the solution? Amazing. My problem is driving me crasy. May I ask you please to e-mail me your 1stOpt code for my problem? As I am not an expert, it will help to learn quicker. If this is not too much to ask for. Here is my email: [e-mail] Thank you, that would help me a lot.

[Alex Sha]

Hi, Tekin, have you send the same equations with exacutely same constant values to 1stOpt people? Your equations presented here can be solved easily by 1stOpt, however, if same constant value changed, for instance, "Tl1 = 1.358686751924744" becomes "Tl1 = 0.5", then it will become hard for 1stOpt finding high precision results， one outcome likes below:

Objective Function (Min.): 0.0144534235661664
x1: 0.00889101436256974
x2: 3.65225963270091E-5
x3: 5.04734772446125
x4: 3.35950379013571E-11
x5: 2.66239114539507
x6: 0.000215788866428105
x7: 2133.07291865425
x8: 0.281930090949226

fevl:
0.0299000026051724
0.00479327540684782
-0.0377035006522798
-0.0253917473235021
-0.101604577161084
0.000129045563598817
-4.52553678975676E-8
-0.0338620196508128

So it is important that the constant values your provided are reasonable。

[Alex Tekin]

Hi, I see. Hmmm....so if constant values do change, and the 1stOpt cannot find the precise answer ... What is the solution to this? Does it mean, the solution does not exist? Or it is ill-posed, I mean the problem? Is it dead-end? Does it require reformulation?

[Alex Sha]

Yes, if the original system equations cann't be found with precise answer, reformulation is a good and efficient approaches, better to try.

[Alex Sha]

with the case: "Tl1 = 1.358686751924744" becomes "Tl1 = 0.5", if make reformulation by removing denominator, high precise outcome could be achieved easily:

Objective Function (Min.): 2.18196747888083E-30
x1: 0.0147749379362527
x2: 0.0419344037802394
x3: 2.62626735850057
x4: 4.12662981300126E-6
x5: 0.819015423317213
x6: 0.000533791102422822
x7: 235.879794583329
x8: 0.324422798162672

Fevl:
-1.03231178176173E-36
-1.47708143713109E-15
-1.38777878078145E-17
-1.55040910665427E-16
1.07308910021753E-16
-6.68952740423556E-17
-1.11022302462516E-16
-1.11022302462516E-16

So how to write the formation of each equation are sometime especially important.

[Alex Tekin]

So, for every change of a paramater there shoud be a reformulation of the problem? Or, once the denominator is removed, it should be OK? Thanks.

[Alex Sha]

Hi, it depends on each case。Usually, reformulation by removing denominator is suggested for complicated equations. For example, change your constant values form "Tl1 = 1.358686751924744, TL2 = 0.438403078319232" to "Tl1 = 0.5, TL2 = 0.2", it could be get high precise results without difficult for reformulation equations, but hard for your original equations.

For reformulation one:

Objective Function (Min.): 1.33981412566616E-28
x1: 0.0261587268893943
x2: 0.030658781704042
x3: 85.4185713752472
x4: 0
x5: 0.0984437936656827
x6: 0.001303207052684
x7: 85.3144796649945
x8: 0.327978197428879

Fevl:
7.40821135477289E-31
-1.15750338611737E-14
-3.81639164714898E-17
4.33680868994202E-18
-2.33103467084383E-18
-4.77048955893622E-18
3.33066907387547E-16
4.44089209850063E-16

For original equations

Objective Function (Min.): 0.104414953695963
x1: 0.0119428767503493
x2: 0.000253645336492324
x3: 146.99095794656
x4: 0
x5: 8.37423563083838
x6: 0.000293736226647899
x7: 1664.14920244842
x8: 0.29202059266871

Fevl:
0.0444011098139584
0.2375
-0.0848885449713843
-0.0507359335324231
-0.183200927730478
0.000271985114754568
1.37633186375608E-5
-0.0519075239040255

[Alex Tekin]

Thank you very much Alex Sha! So, just to be sure, does it mean that if one simply gets rid of "all" denominators, then it should be fine? Or not? if not, then how to make sure which denominator is causing a problem? Also, why is problem arising to begin with? Is it always with the division? Maybe a negative number is raised to fractional power? Maybe that is the problem?

[Alex Sha]

Hi, Tekin, no common answer for all equation solving problems, it should be treated case by case. The percise of numerical computation has huge infulence some time. So if the system equations in original formulation can be solved correctly, it's OK, otherwise, you may think to make reformulation of your original equation by using multiplication instead of division, if possible. Show you one more case below:

60/(x-0.85)=2664/x;

what's the result of above? the result from 1stOpt or Baron will looks like: x: -1.66332860430034E165 with Fevl = 1.56553551310767E-162

Although the presice is high enough, the result of x is not what we want absolutely.

Make reformulation of above equation to:

60*x=2664*(x-0.85)

we get reasonable outcome: x = 0.869585253456221, this is what we want.

While if the original equation is:

60/(x-0.85)=x/2664

the reformulation, like "60*2664=(x-0.85)*x", will be become unnecessary since the same reaults will be achieved with x = -399.375175869146 or x = 400.225175861622.

Hope helpful.

[Alex Sha]

Here is one case that reformulation will lead to "strange" outcome:

Original equations:

0.946050282299=(86.5312-a)/(b-a)+(75.9074-c)/(d-c);
1.25=(93.4023-a)/(b-a)+(69.2394-c)/(d-c);
1.34=(92.8411-a)/(b-a)+(65.1376-c)/(d-c);
0.95=(77.2-a)/(b-a)+(68.10-c)/(d-c);

There are multi-solutions for above equations:

1:
a: 15.0063003070515
b: 64.3798312420418
c: 55.5697165060803
d: 15.1046370719439

Fevl:
-2.04281036531029E-14
4.70734562441066E-14
-3.79696274421804E-14
1.06581410364015E-14

2:
a: 54.4518156991909
b: 103.825346634182
c: 87.898089288757
d: 47.4330098546204

Fevl:
-1.63202784619898E-14
5.39568389967826E-14
-3.08642000845794E-14
1.25455201782643E-14

Reformulation of equations:

0.946050282299*(b-a)*(d-c)=(86.5312-a)*(d-c)+(75.9074-c)*(b-a);
1.25*(b-a)*(d-c)=(93.4023-a)*(d-c)+(69.2394-c)*(b-a);
1.34*(b-a)*(d-c)=(92.8411-a)*(d-c)+(65.1376-c)*(b-a);
0.95*(b-a)*(d-c)=(77.2-a)*(d-c)+(68.10-c)*(b-a);

results:

1:
a: -0.383222844916884
b: -0.383222844916884
c: 0.415432086157557
d: 0.415432086157557

Fevl:
0
0
0
0

2:
a: 7.54318132141484
b: 7.54318132141484
c: 2.03628568296322
d: 2.03628568296322

Fevl:
0
0
0
0

The results of above are always "a=b" and "c=d", the precise is prefect, however, those results may not what we want, and those results can not be used as the roots of original ones, since will leads to the error of "divide by zero".

So the system equations are usually complicated, they should be treated case by case.

[Alex Tekin]

Hi Alex Sha,

Thanks a lot!

I just used your 1st example in Baron w/o reformulation, and even though I give large range for bounds and start far away from the original solution, it gives me a right answer in no time. I am surprised you Baron requires reformulation for this example. This is my code:

clc

format long

tic

options = baronset('DeltaTerm',0,'MaxTime',500,'filekp',1);

%Objective

fun = @(x) 0;

%Nonlinear Constraints

nlcon = @(x) [60/(x(1)-0.85)-2664/x(1)];

cl = [0];

cu = [0];

%Bounds

lb = [0];

ub = [1000];

% Starting Guess

x0 = [1000];

%Solve

[x,fval,exitflag,info] = baron(fun,[],[],[],lb,ub,nlcon,cl,cu,[],x0,options)

toc

[Alex Sha]

Hi, Tekin, you may download a software package GAMS trial version from https://www.gams.com, Baron solver is included, the code for first case (slight changed： “-2664/x” -> "-2664/(x-0.01)", avoided devision by zero) looks like:

variables obj, x;
equations defobj, con1;
defobj.. obj =e= 1;
con1.. 60/(x-0.85) =e= -2664/(x-0.01);
model m / all /;
option nlp=baron;
solve m minimizing obj using nlp;
option decimals=8;
display obj.l, x.l;

It gives result:

VARIABLE obj.L                 =   1.00000000  
VARIABLE x.L                   =  -1.0000E+51

the x value obtained by Baron is -1.0000E+51

Reformulation code:

variables obj, x;
equations defobj, con1;
defobj.. obj =e= 1;
con1.. 60*(x-0.01) =e= -2664*(x-0.85);
model m / all /;
option nlp=baron;
solve m minimizing obj using nlp;
option decimals=8;
display obj.l, x.l;

result

VARIABLE obj.L                 =   1.00000000  
VARIABLE x.L                   =   0.83149780 

the x value is 0.83149780

Answer 2

[John D'Errico]

This same question gets asked on literally a daily basis. The problem is not fsolve. It lies in understanding nonlinear equations, and things like basins of attraction. Finally, part of the problem lies in understanding numerical computation in floating point arithmetic, because you have lots of exponents in there. Tiny changes in a parameter will result in huge differences, perhaps in some terms suddenly going into the complex domain.

Without looking seriously at your equations, except to note they form an incredible mess of terms, perhaps even a hundred lines or so in length. I see no trig functions or other periodic functions, so I'll predict there are problably only a few dozen or so real solutions. There will almost certainly be multiple solutions, if any solutions exist at all. (The presence of periodic functions will often turn a problem with finitely many solutions into usually one with infinitely many solutions.)

For example, here is ONE of the 8 equations posed to be solved:

F(1) = - (x(6)^(4-4*n)*((x(1)^2*n^2*x(6)^(2*n-2)*(r-1)^2*(c-1)^2*(x(1)*x(6)^n*(x(1)*x(4)*n*x(6)^(n-1)*(r-1)*(c-1) - x(1)*x(6)^n*(x(2)*(x(5) + x(5)^2 + x(1)*n*x(6)^(n-1)*(r-1)*(c-1) + x(1)*n*x(6)^(n-1)*x(5)*(r-1)*(c-1)+1) + x(1)^2*x(3)*n*x(6)^n*x(6)^(n-1)*u*(n-1)*(r-1)*(c-1))*(x(5)^2+1)*(n-1)*(f-1))*(n-1)*(f-1) + x(1)*x(4)*x(6)^n*x(5)*(x(5) + x(1)*n*x(6)^(n-1)*(r-1)*(c-1))*(x(5)^2+1)*(n-1)*(f-1)) + x(1)*x(4)*x(6)^n*(n-1)*(f-1)*(2*x(5)^2 + 2*x(5)^4 + x(5)^6 + x(1)^2*n^2*x(6)^(2*n-2)*x(5)^2*(x(5)^2 + 2)*(r-1)^2*(c-1)^2 + 2*x(1)*n*x(6)^(n-1)*x(5)*(x(5)^2+1)^2*(r-1)*(c-1)+1))*(1/(exp(-x(6)*x(7))+1)^(1/x(8))-1) - x(1)^2*x(4)*n^2*x(6)^(2*n-2)*(x(1)*x(6)^n*(n-1)*(f-1) - x(1)^2*n*x(6)^n*x(6)^(n-1)*(n-1)*(f-1)*(r-1)*(c-1))*(r-1)^2*(c-1)^2 + x(1)^3*n^2*x(6)^n*x(6)^(2*n-2)*(x(5)^2+1)*(x(4)*x(5)*(x(5) + x(1)*n*x(6)^(n-1)*(r-1)*(c-1)) + x(1)^2*n^2*x(6)^(2*n-2)*(x(4) - x(1)*x(2)*x(6)^n*(n-1)*(f-1))*(r-1)^2*(c-1)^2)*(n-1)*(f-1)*(r-1)^2*(c-1)^2*(1/(exp(-x(6)*x(7))+1)^(1/x(8))-1)))/(x(1)^7*x(3)*n^4*x(6)^(3*n)*(x(5)^2+1)*(n-1)^3*(f-1)^3*(r-1)^4*(c-1)^4*(1/(exp(-x(6)*x(7))+1)^(1/x(8))-1)) - (x(4)*x(6)^(4-4*n)*(x(1)*x(6)^n*(n-1) - x(1)*x(6)^n*f*(n-1) + x(1)^2*n*x(6)^n*x(6)^(n-1)*(n-1)*(f-1)*(r-1)*(c-1))*(x(5)^2 + x(1)*n*x(6)^(n-1)*x(5)*(r-1)*(c-1)+1)^2)/(x(1)^7*x(3)*n^4*x(6)^(3*n)*(exp(-x(6)*x(7))+1)^(1/x(8))*(x(5)^2+1)*(n-1)^3*(f-1)^3*(r-1)^4*(c-1)^4*(1/(exp(-x(6)*x(7))+1)^(1/x(8))-1)) - a;

I should probably have split it into multiple lines to make it possible to read without scrolling too far to the right.

ANY numerical solver is sensitive to initial guesses. Change the start point, get a completely different result. This is the concept of a basin of attraction, thus the set of start points that will result in any given solution. A basin of attaction can sometimes be a very strangely shaped set, especially so in 8 dimensions.

As I said, your problem is not in fsolve. It lies in assuming that you can just throw any complete mess of computations at a computer and expect a reasonable result, without also understanding the methods that must be used to solve the problem.

[Alex Tekin]

Thank you, this is very useful to know, much appreciated. Is there a way, for this system to find all of these multiple solutions? Can you suggest a different method to tackle the problem?

[John D'Errico]

I answered the same question a day or so ago, but the question, and my lengthy answer and followup comments got deleted with it.

Is there a way to find all solutions of a nonlinear multi-variate probem?

No, in general this is not possible. In fact, it is provably true that you cannot do so for many actually simple looking provblems. Nothing stops you from using multi-start methods with different random starting points to hope to gain a solution. But at best you can only hope to find perhaps many of the solutions. To know that you have found them all would require bounds on the derivatives of your function within regions, which could allow you to exclude a possible solution in that region. (I recall seeing work on global solvers some years ago, try to achieve exactly this sort of thing. Those ideas won't go anywhere on such a nasty problem. 8 dimensional spaces are big, and the derivatives of your expressions will be really nasty.)

Worse, the very fact that fsolve fails when given a start point away from the solution tends to imply that your function has nasty derivatives. Exponentials, and lots of exponents in general make that just expected.

Can I suggest a different method?

No. Too often I see people create massive complicated systems like this, then expecting their powerful computer to make short work of it. Not gonna happen most of the time. Computers are not yet that powerful. At least some of the time, you cannot even work with the expressions involved using double precision arithmetic. The problem is though, using higher precision is usually a route to failure too, since high precision computation is so slow to work with.

At best, you could consider using vpasolve, using different start points. So a numerical solver in the symbolic toolbox. Such an approach will be highly problematic, and even if you make any progress at all, it will be incredibly slow. Again, high precision symbolic computations are terribly slow.

My suggestion, if anything would work is to step back, to reformulate your problem. I presume this problem comes from some real world system, and is an approximation thereof. You need to reduce the complexity of your approximation, eliminating terms that are of little significance. This is a difficult thing to do, and requires sometimes great skill and true understanding of the problem. But done well, it can make a once impossible problem solvable.

Or, you can keep trying to hack this problem to death, hopefully before it does the same to you.

[Alex Tekin]

Thank you for your reply, which is informative and useful.

Unfortunately, making a system less complicated is not an option for me, so either I have to keep on trying to solve it (or get as many of those solutions as possible), or I just have to totally give up.

By the way, I tried vpasolve for the above problem, and it performs even worse than fsolve. Even when I insert the initial guesses that are extremely close to the solutions, vpasolve returns an empty answer while fsolve does not. I have checked manually, for the first three equations vpasolve encounters no problem. Second I introduce the fourth one, it returns empty set, even though there is cleary a solution.

I do understand that no numerical technique is meant to be “perfect”, and computers are not meant to be all powerful. However, the truth is, fsolve and vpasolve in MATLAB will not solve fairly complicated problems, and occasionally even if one starts quite close to the solutions. This is certainly obvious for more experienced users, but not for those who are fairly new to programming and numerical techniques.

I am not technically equipped to deeply research numerical approximation techniques, so I would like to request further help, or some reference to the codes that can be used in MATLAB (in place of fsolve vpasolve) that can more successfully tackle the above problem. I know there has been some good research on generic algorithms etc but I could not find a straightforward example, or user-supplied code to apply them to this particular problem. Thanks.