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

Answer

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.