Solving a system of nonlinear equations 2

Question

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}]]

Answer 1

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. [Edit] 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.

Answer 2

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.