# Solving a system of nonlinear equations 2

Regarding this question: Solving a system of nonlinear equations When I try to solve the same system (with the code proposed in Daniel Lichtblau's answer) with different Right Hand Side values, Mathematica does not return any answer (while "fsolve" in Matlab gives me the exact answer as I am anticipating, that is: α1=3, α2=2, α3=1, β1=1, β2=1, β3=1) I used the following code in Mathematica. Please let me know what is the reason?

M = {{α1, β3, β2}, {β3, α2, β1}, {β2, β1, α3}};

Timing[
sol2 = NSolve[{
β1 - (β2 β3)/α1 == 2/3,
α2 - (β3 β3)/α1 == 5/3,
α3 - (β2 β2)/α1 == 2/3,
α1 - (β2 β2)/α3 == 2,
α2 - (β1 β1)/α3 == 1,
β3 - (β1 β2)/α3 == 0,
(4/3*Pi)^2 == 2.96^2*(CharacteristicPolynomial[M, x] /. x -> 0)},
{α1, α2, α3, β1, β2, β3}]]

• Is there a typo?: 1) The Matlab solution is not a solution to your system. Try Block[{α1 = 3, α2 = 2, α3 = 1, β1 = 0.2, β2 = 0.2, β3 = 0.2}, equations /. Equal -> Subtract] on your equations. The first three equations have sides that differ by -0.0000333333. 2) You have more equations than unknowns. Getting them to be consistent with machine numbers might be tricky. You might have to solve the first six and then select the solution(s) that satisfy the seventh within a suitably chosen tolerance. – Michael E2 Oct 18 '14 at 12:27
• if I show the coefficients in Right Hand Side(RHS) with X and Z as follows: ,β1 - (β2 β3)/α1 =X1, α2 - (β3 β3)/α1 =X2, α3 - (β2 β2)/α1 =X3, α1 - (β2 β2)/α3 =Z1, α2 - (β1 β1)/α3 =Z2, β3 - (β1 β2)/α3 =Z3, then there is a dependency between first 6 equations as following: (Z3^2/Z1)-Z2=(X1^2/X3)-X2 I think I have actually 5 equation instead of 6 and that's why I have added the last equation (it is a conservation of mass) and ...(please go to the next comment) – vorujak Oct 18 '14 at 12:46
• the above system has no problem in solving many cases (like: mathematica.stackexchange.com/questions/60861/…) but, I have problem with these specific above mentioned values for RHS (as mentioned in this recent question above). I do not know, what should I do in this situation? – vorujak Oct 18 '14 at 12:47
• Very minor improvement on @MichaelE2 approach: First create expressions by subtracting right sides from left sides of the equations. Then can do NSolve[Flatten[{Thread[Rest[exprs] == 0], First[exprs]^2 < 10^(-4)}], vars]. Advantage is it filters out the "parasite" solutions automatically. – Daniel Lichtblau Oct 18 '14 at 18:09
• My thought is that for now the method of @MichaelE2 is the best I can think of. Solve for all but the first equation (so it's neither over nor under determined), then use the first equation to rule out "bad" solutions. The problem is the system is not what's called a "complete intersection", and NSolve is not able to discern that it is consistent within a certain tolerance. If I get time I'll look into that but in any case NSolve is not now able to handle the full system. – Daniel Lichtblau Oct 20 '14 at 21:54

A system with approximate (Real) coefficients sometimes has only approximate solutions. Minimizing the norm of the residuals, approach 3 below, may be the best way to approximate the solutions. In this case, we have seven equations in six unknowns.

equations = {
β1 - (β2 β3)/α1 == 0.1867,
α2 - (β3 β3)/α1 == 1.9867,
α3 - (β2 β2)/α1 == 0.9867,
α1 - (β2 β2)/α3 == 2.96,
α2 - (β1 β1)/α3 == 1.96,
β3 - (β1 β2)/α3 == 0.16,
(4/3*Pi)^2 == 1.743^2*(CharacteristicPolynomial[M, x] /. x -> 0)};
forms = equations /.  Equal -> Subtract; (* differences between the sides of the equations *)
variables = Variables[forms]
(*
{α1, α2, α3, β1, β2, β3}
*)


The OP mentioned in a comment that the first six are dependent, but this is true only approximately so. It seems that NSolve thinks that they are independent and, as a result, inconsistent. (The functions in forms were introduced for purposes that become clear below; the system of equations is equivalent to forms == 0.)

A few approaches come to mind:

1. Try to increase the tolerance so that NSolve solves the system as intended.
2. Omit one of the equations, solve the complementary system, and select solutions that are approximate solutions of the omitted equation.
3. Minimize the distance between the two sides of the equations.

Approach 1

I was unsuccessful. There are several avenues (e.g., precision, the system option "NSolveOptions" -> {"Tolerance" -> tol}, Internal\$EqualTolerance), but I could find no combination of them that worked.

Approach 2

One can drop an equation with Drop. It turns out that the first equation is dependent on the rest and may be dropped.  Following Daniel Lichtblau's advice in a comment, we can add a condition, Abs@forms[[1]] < 0.0001, that the first equation be satisfied within a certain tolerance, say, 0.0001. Then we get two solutions:

sols = NSolve[Append[Drop[equations, 1], Abs@forms[[1]] < 0.0001], variables]
(*
{{α1 -> 2.99959, α2 -> 1.99187, α3 -> 0.999897, β1 -> 0.178501, β2 -> -0.198961, β3 -> 0.124482},
{α1 -> 2.99959, α2 -> 2.00001, α3 -> 0.999897, β1 -> 0.20001, β2 -> 0.198961, β3 -> 0.199798}}
*)


One drawback is that the chosen equation to be dropped will be approximately satisfied while the rest are satisfied exactly. Indeed all the error is forced on the chosen equation. Minimizing the norm of the residuals is probably to be preferred, since it shares out the error. This is done below in approach 3 by processing the results of this section.

Approach 3

Here we want to minimize the distance between the two sides. Thus our objective function could be

forms^2 // Total


Or perhaps better, we could scale forms by the magnitude of the gradients at the solutions, given by

df = ComplexExpand /@ Norm /@ D[forms, {variables}];


So that the objective function would be each of the following (for each respective solution):

(forms^2).(1/df /. sols2[[1]])
(forms^2).(1/df /. sols2[[2]])


The best way to proceed is to start with the approximate solutions found in Approach 2. One might use NMinimize, but that would turn out to be less satisfactory (see below). Instead let's use FindArgMin. We can use each solution found above as a starting point:

FindMinimum[(forms^2).(1/df /. #), List @@@ #] & /@ sols2
(*
{{2.99957, 1.99186, 0.999892, 0.178447, -0.198922, 0.124498},
{2.99957, 2., 0.999892, 0.199952, 0.198922, 0.199778}}
*)


Or if you want a Rule:

Thread /@ (variables -> # &) /@ %
(*
{{α1 -> 2.99957, α2 -> 1.99186, α3 -> 0.999892, β1 -> 0.178447, β2 -> -0.198922, β3 -> 0.124498},
{α1 -> 2.99957, α2 -> 2., α3 -> 0.999892, β1 -> 0.199952, β2 -> 0.198922, β3 -> 0.199778}}
*)


Remark: NMinimize might seem like a good approach but it returns only one result. With luck, one might coax it to return different results by tweaking the methods, their parameters, or by using the "RandomSearch" method with different random seeds. But one would not know when to stop except by analyzing the system as in Approach 2.

• Very nice (I think I am the first upvote). At some point I may need to investigate why the Tolerance method is coming up short. (I had of course also tried that, and got nowhere.) – Daniel Lichtblau Oct 18 '14 at 18:12
• @DanielLichtblau Thanks. I was hoping you would know how to use the tolerance option. Ah, well... – Michael E2 Oct 18 '14 at 18:13
• I think I know how to use it. I just can't beat a better result out of it. Alas. – Daniel Lichtblau Oct 18 '14 at 19:26
• The suggestions there seem fine. A third camera or some other independent measurement can be used to exclude one of the two solutions. – Daniel Lichtblau Oct 26 '14 at 22:27
• @Hossein As I said, I think the best way to proceed is to start with the approximate solutions found in approach 2 and apply approach 3. Unless you're sure the approximate solutions of 2 are close enough, in which case you can save time by skipping 3. In other words, I think 3, which is basically to minimize a certain error measure, improves the results of 2. You could also use your own error measure that reflects the sensitivity of your objective to the parameters. E.g., if the result is highly sensitive to small changes in α1, then weight α1 more. I hope that makes sense. :) – Michael E2 Nov 4 '14 at 15:14

I finally got a chance to look into the tolerancing. It seems that there is a fairly narrow window for which NSolve, with a nondefault tolerance setting, will handle this.

mat = {{a1, b3, b2}, {b3, a2, b1}, {b2, b1, a3}};
eqns = {
b1 - (b2 b3)/a1 == 0.1867,
a2 - (b3 b3)/a1 == 1.9867,
a3 - (b2 b2)/a1 == 0.9867,
a1 - (b2 b2)/a3 == 2.96,
a2 - (b1 b1)/a3 == 1.96,
b3 - (b1 b2)/a3 == 0.16,
(4/3*Pi)^2 == 1.743^2*(CharacteristicPolynomial[mat, x] /. x -> 0)};

SetSystemOptions["NSolveOptions"->{"Tolerance"->{10^(-4),10^(-2)}}];


Roughly, this will clip individual polynomial terms that are made a factor of 10^(-4) smaller during the reduction process (in GroebnerBasis, to be technical). It will also entirely remove polynomials when ALL terms have been reduced by a factor of 10^(-2). With this setting:

NSolve[eqns, {a1, a2, a3, b1, b2, b3}]

(* Out[176]//InputForm=
{{a1 -> 2.9995896242394378, a2 -> 1.9999796638931548, a3 -> 0.9998969872422477,
b1 -> 0.19992449742212776, b2 -> 0.19896116707304107,
b3 -> 0.19978673137863223}, {a1 -> 2.9995896242394378,
a2 -> 1.9918501376365532, a3 -> 0.9998969872422477,
b1 -> 0.17846968017618514, b2 -> -0.19896116707304115,
b3 -> 0.1244932279986043}} *)
`

I think it would be correct to say that this is not the best method for the situation at hand, as it relies too heavily on a difficult-to-predict interplay between tolerance settings that only seem to be useful in a very narrow range. In hindsight we have some idea of what range to look in, since, a posteriori, we have an idea of the sizes if residuals from the approximate solutions. But that's not of much help when one does not already know the answer. In that case about the best to be done is to base tolerances on some measure of the estimated experimental error that gave rise to the equation coefficients.

• Dear Daniel, Thanks for your helpful answer. As you know, the Right Hand Side (RHS) values in question are just an example here and they may change in next situations (for next droplets frames). Is there any way that I can calculate that based on my RHS values, which Tolerance is suitable to use in order to converge to the correct answer(s)? If not, what is your suggestion for me to be able to solve a system while RHS values vary? Thanks again – vorujak Nov 4 '14 at 5:07
• I think the method of Approach 2 proposed by @Michael E2 is likely to be the most reliable in general. – Daniel Lichtblau Nov 4 '14 at 14:50