Multiple solutions from NSolve

Question

I have a list of polynomial equations which I can solve using NSolve. However, it gives me a large number of similar solutions. The solutions are largely identical except for numerical error. The solution time is also far longer than I would like. Can we find a way of not calculating nearly identical solutions?

The unknowns and the equations are (sorry for the length and Greek)

uk = {A[0], A[1], B[1], A[2], B[2], A[3], B[3]};
eqns = {\[Omega]0^2 A[0] + \[Beta] A[0]^3 + 3/2 \[Beta] A[0] A[1]^2 + 
     3/4 \[Beta] A[1]^2 A[2] + 3/2 \[Beta] A[0] A[2]^2 + 
     3/2 \[Beta] A[1] A[2] A[3] + 3/2 \[Beta] A[0] A[3]^2 + 
     3/2 \[Beta] A[0] B[1]^2 - 3/4 \[Beta] A[2] B[1]^2 + 
     3/2 \[Beta] A[1] B[1] B[2] - 3/2 \[Beta] A[3] B[1] B[2] + 
     3/2 \[Beta] A[0] B[2]^2 + 3/2 \[Beta] A[2] B[1] B[3] + 
     3/2 \[Beta] A[1] B[2] B[3] + 3/2 \[Beta] A[0] B[3]^2 == 
    0, -F - \[Omega]^2 A[1] + \[Omega]0^2 A[1] + 
     3 \[Beta] A[0]^2 A[1] + 3/4 \[Beta] A[1]^3 + 
     3 \[Beta] A[0] A[1] A[2] + 3/2 \[Beta] A[1] A[2]^2 + 
     3/4 \[Beta] A[1]^2 A[3] + 3 \[Beta] A[0] A[2] A[3] + 
     3/4 \[Beta] A[2]^2 A[3] + 3/2 \[Beta] A[1] A[3]^2 + 
     2 \[Zeta] \[Omega] \[Omega]0 B[1] + 3/4 \[Beta] A[1] B[1]^2 - 
     3/4 \[Beta] A[3] B[1]^2 + 3 \[Beta] A[0] B[1] B[2] + 
     3/2 \[Beta] A[1] B[2]^2 - 3/4 \[Beta] A[3] B[2]^2 + 
     3/2 \[Beta] A[1] B[1] B[3] + 3 \[Beta] A[0] B[2] B[3] + 
     3/2 \[Beta] A[2] B[2] B[3] + 3/2 \[Beta] A[1] B[3]^2 == 
    0, -2 \[Zeta] \[Omega] \[Omega]0 A[1] - \[Omega]^2 B[
       1] + \[Omega]0^2 B[1] + 3 \[Beta] A[0]^2 B[1] + 
     3/4 \[Beta] A[1]^2 B[1] - 3 \[Beta] A[0] A[2] B[1] + 
     3/2 \[Beta] A[2]^2 B[1] - 3/2 \[Beta] A[1] A[3] B[1] + 
     3/2 \[Beta] A[3]^2 B[1] + 3/4 \[Beta] B[1]^3 + 
     3 \[Beta] A[0] A[1] B[2] - 3 \[Beta] A[0] A[3] B[2] + 
     3/2 \[Beta] A[2] A[3] B[2] + 3/2 \[Beta] B[1] B[2]^2 + 
     3/4 \[Beta] A[1]^2 B[3] + 3 \[Beta] A[0] A[2] B[3] - 
     3/4 \[Beta] A[2]^2 B[3] - 3/4 \[Beta] B[1]^2 B[3] + 
     3/4 \[Beta] B[2]^2 B[3] + 3/2 \[Beta] B[1] B[3]^2 == 0, 
   3/2 \[Beta] A[0] A[1]^2 - 4 \[Omega]^2 A[2] + \[Omega]0^2 A[2] + 
     3 \[Beta] A[0]^2 A[2] + 3/2 \[Beta] A[1]^2 A[2] + 
     3/4 \[Beta] A[2]^3 + 3 \[Beta] A[0] A[1] A[3] + 
     3/2 \[Beta] A[1] A[2] A[3] + 3/2 \[Beta] A[2] A[3]^2 - 
     3/2 \[Beta] A[0] B[1]^2 + 3/2 \[Beta] A[2] B[1]^2 + 
     4 \[Zeta] \[Omega] \[Omega]0 B[2] + 3/2 \[Beta] A[3] B[1] B[2] + 
     3/4 \[Beta] A[2] B[2]^2 + 3 \[Beta] A[0] B[1] B[3] - 
     3/2 \[Beta] A[2] B[1] B[3] + 3/2 \[Beta] A[1] B[2] B[3] + 
     3/2 \[Beta] A[2] B[3]^2 == 
    0, -4 \[Zeta] \[Omega] \[Omega]0 A[2] + 
     3 \[Beta] A[0] A[1] B[1] - 3 \[Beta] A[0] A[3] B[1] + 
     3/2 \[Beta] A[2] A[3] B[1] - 
     4 \[Omega]^2 B[2] + \[Omega]0^2 B[2] + 3 \[Beta] A[0]^2 B[2] + 
     3/2 \[Beta] A[1]^2 B[2] + 3/4 \[Beta] A[2]^2 B[2] - 
     3/2 \[Beta] A[1] A[3] B[2] + 3/2 \[Beta] A[3]^2 B[2] + 
     3/2 \[Beta] B[1]^2 B[2] + 3/4 \[Beta] B[2]^3 + 
     3 \[Beta] A[0] A[1] B[3] + 3/2 \[Beta] A[1] A[2] B[3] + 
     3/2 \[Beta] B[1] B[2] B[3] + 3/2 \[Beta] B[2] B[3]^2 == 0, 
   1/4 \[Beta] A[1]^3 + 3 \[Beta] A[0] A[1] A[2] + 
     3/4 \[Beta] A[1] A[2]^2 - 9 \[Omega]^2 A[3] + \[Omega]0^2 A[3] + 
     3 \[Beta] A[0]^2 A[3] + 3/2 \[Beta] A[1]^2 A[3] + 
     3/2 \[Beta] A[2]^2 A[3] + 3/4 \[Beta] A[3]^3 - 
     3/4 \[Beta] A[1] B[1]^2 + 3/2 \[Beta] A[3] B[1]^2 - 
     3 \[Beta] A[0] B[1] B[2] + 3/2 \[Beta] A[2] B[1] B[2] - 
     3/4 \[Beta] A[1] B[2]^2 + 3/2 \[Beta] A[3] B[2]^2 + 
     6 \[Zeta] \[Omega] \[Omega]0 B[3] + 3/4 \[Beta] A[3] B[3]^2 == 
    0, -6 \[Zeta] \[Omega] \[Omega]0 A[3] + 3/4 \[Beta] A[1]^2 B[1] + 
     3 \[Beta] A[0] A[2] B[1] - 3/4 \[Beta] A[2]^2 B[1] - 
     1/4 \[Beta] B[1]^3 + 3 \[Beta] A[0] A[1] B[2] + 
     3/2 \[Beta] A[1] A[2] B[2] + 3/4 \[Beta] B[1] B[2]^2 - 
     9 \[Omega]^2 B[3] + \[Omega]0^2 B[3] + 3 \[Beta] A[0]^2 B[3] + 
     3/2 \[Beta] A[1]^2 B[3] + 3/2 \[Beta] A[2]^2 B[3] + 
     3/4 \[Beta] A[3]^2 B[3] + 3/2 \[Beta] B[1]^2 B[3] + 
     3/2 \[Beta] B[2]^2 B[3] + 3/4 \[Beta] B[3]^3 == 0};

If I put in some numerical values for the coefficients and solve I get 483 solutions.

eqns1 = eqns /. {\[Omega]0 -> 2 \[Pi] 10, \[Zeta] -> 1/100, \[Beta] ->
      5000, F -> 100., \[Omega] -> 2 \[Pi] 12.};
Timing[sol = NSolve[eqns1, uk, Reals];]

(* {0.171875, Null} *)

The first 49 solutions are identical except for numerical error.

Table[sol[[n]][[All, 2]] - sol[[1]][[All, 2]], {n, 49}] // Chop

(* {{0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0}} *)

Other groups of solutions are also similarly identical. If I put in alternative values for the coefficients e.g. \[Omega] -> 2 \[Pi] 14. I get a different number of solutions. I need to calculate solutions for a range of values of Omega.

DeleteDuplicates does not work but anyway I would like to make NSolve work without producing duplicates. Any ideas?

You have a precision problem. Therefore, increase "WorkingPrecision" — Daniel Huber, Commented Nov 8, 2023 at 13:17
@Nasser well spotted. Bad habit of using NDSolve more than NSolve. Thanks corrected. — Hugh, Commented Nov 8, 2023 at 14:18
@DanielHuber Increasing or decreasing WorkingPrecision results in "Subsystem could not be solved". Also, takes a long time to reach that conclusion. — Hugh, Commented Nov 8, 2023 at 14:32
@AlexTrounev I tried using FindRoot. It needs starting values so I used random ones and then I get a failure due to lack of MachinePrecision. — Hugh, Commented Nov 8, 2023 at 14:47

Alex Trounev · Accepted Answer · 2023-11-15 02:53:12Z

For parametric research it could be better to use FindRoot as follows

uk = {A[0], A[1], B[1], A[2], B[2], A[3], B[3]};
eqns = {\[Omega]0^2 A[0] + \[Beta] A[0]^3 + 3/2 \[Beta] A[0] A[1]^2 + 
     3/4 \[Beta] A[1]^2 A[2] + 3/2 \[Beta] A[0] A[2]^2 + 
     3/2 \[Beta] A[1] A[2] A[3] + 3/2 \[Beta] A[0] A[3]^2 + 
     3/2 \[Beta] A[0] B[1]^2 - 3/4 \[Beta] A[2] B[1]^2 + 
     3/2 \[Beta] A[1] B[1] B[2] - 3/2 \[Beta] A[3] B[1] B[2] + 
     3/2 \[Beta] A[0] B[2]^2 + 3/2 \[Beta] A[2] B[1] B[3] + 
     3/2 \[Beta] A[1] B[2] B[3] + 3/2 \[Beta] A[0] B[3]^2 == 
    0, -F - \[Omega]^2 A[1] + \[Omega]0^2 A[1] + 
     3 \[Beta] A[0]^2 A[1] + 3/4 \[Beta] A[1]^3 + 
     3 \[Beta] A[0] A[1] A[2] + 3/2 \[Beta] A[1] A[2]^2 + 
     3/4 \[Beta] A[1]^2 A[3] + 3 \[Beta] A[0] A[2] A[3] + 
     3/4 \[Beta] A[2]^2 A[3] + 3/2 \[Beta] A[1] A[3]^2 + 
     2 \[Zeta] \[Omega] \[Omega]0 B[1] + 3/4 \[Beta] A[1] B[1]^2 - 
     3/4 \[Beta] A[3] B[1]^2 + 3 \[Beta] A[0] B[1] B[2] + 
     3/2 \[Beta] A[1] B[2]^2 - 3/4 \[Beta] A[3] B[2]^2 + 
     3/2 \[Beta] A[1] B[1] B[3] + 3 \[Beta] A[0] B[2] B[3] + 
     3/2 \[Beta] A[2] B[2] B[3] + 3/2 \[Beta] A[1] B[3]^2 == 
    0, -2 \[Zeta] \[Omega] \[Omega]0 A[1] - \[Omega]^2 B[
       1] + \[Omega]0^2 B[1] + 3 \[Beta] A[0]^2 B[1] + 
     3/4 \[Beta] A[1]^2 B[1] - 3 \[Beta] A[0] A[2] B[1] + 
     3/2 \[Beta] A[2]^2 B[1] - 3/2 \[Beta] A[1] A[3] B[1] + 
     3/2 \[Beta] A[3]^2 B[1] + 3/4 \[Beta] B[1]^3 + 
     3 \[Beta] A[0] A[1] B[2] - 3 \[Beta] A[0] A[3] B[2] + 
     3/2 \[Beta] A[2] A[3] B[2] + 3/2 \[Beta] B[1] B[2]^2 + 
     3/4 \[Beta] A[1]^2 B[3] + 3 \[Beta] A[0] A[2] B[3] - 
     3/4 \[Beta] A[2]^2 B[3] - 3/4 \[Beta] B[1]^2 B[3] + 
     3/4 \[Beta] B[2]^2 B[3] + 3/2 \[Beta] B[1] B[3]^2 == 0, 
   3/2 \[Beta] A[0] A[1]^2 - 4 \[Omega]^2 A[2] + \[Omega]0^2 A[2] + 
     3 \[Beta] A[0]^2 A[2] + 3/2 \[Beta] A[1]^2 A[2] + 
     3/4 \[Beta] A[2]^3 + 3 \[Beta] A[0] A[1] A[3] + 
     3/2 \[Beta] A[1] A[2] A[3] + 3/2 \[Beta] A[2] A[3]^2 - 
     3/2 \[Beta] A[0] B[1]^2 + 3/2 \[Beta] A[2] B[1]^2 + 
     4 \[Zeta] \[Omega] \[Omega]0 B[2] + 3/2 \[Beta] A[3] B[1] B[2] + 
     3/4 \[Beta] A[2] B[2]^2 + 3 \[Beta] A[0] B[1] B[3] - 
     3/2 \[Beta] A[2] B[1] B[3] + 3/2 \[Beta] A[1] B[2] B[3] + 
     3/2 \[Beta] A[2] B[3]^2 == 
    0, -4 \[Zeta] \[Omega] \[Omega]0 A[2] + 
     3 \[Beta] A[0] A[1] B[1] - 3 \[Beta] A[0] A[3] B[1] + 
     3/2 \[Beta] A[2] A[3] B[1] - 
     4 \[Omega]^2 B[2] + \[Omega]0^2 B[2] + 3 \[Beta] A[0]^2 B[2] + 
     3/2 \[Beta] A[1]^2 B[2] + 3/4 \[Beta] A[2]^2 B[2] - 
     3/2 \[Beta] A[1] A[3] B[2] + 3/2 \[Beta] A[3]^2 B[2] + 
     3/2 \[Beta] B[1]^2 B[2] + 3/4 \[Beta] B[2]^3 + 
     3 \[Beta] A[0] A[1] B[3] + 3/2 \[Beta] A[1] A[2] B[3] + 
     3/2 \[Beta] B[1] B[2] B[3] + 3/2 \[Beta] B[2] B[3]^2 == 0, 
   1/4 \[Beta] A[1]^3 + 3 \[Beta] A[0] A[1] A[2] + 
     3/4 \[Beta] A[1] A[2]^2 - 9 \[Omega]^2 A[3] + \[Omega]0^2 A[3] + 
     3 \[Beta] A[0]^2 A[3] + 3/2 \[Beta] A[1]^2 A[3] + 
     3/2 \[Beta] A[2]^2 A[3] + 3/4 \[Beta] A[3]^3 - 
     3/4 \[Beta] A[1] B[1]^2 + 3/2 \[Beta] A[3] B[1]^2 - 
     3 \[Beta] A[0] B[1] B[2] + 3/2 \[Beta] A[2] B[1] B[2] - 
     3/4 \[Beta] A[1] B[2]^2 + 3/2 \[Beta] A[3] B[2]^2 + 
     6 \[Zeta] \[Omega] \[Omega]0 B[3] + 3/4 \[Beta] A[3] B[3]^2 == 
    0, -6 \[Zeta] \[Omega] \[Omega]0 A[3] + 3/4 \[Beta] A[1]^2 B[1] + 
     3 \[Beta] A[0] A[2] B[1] - 3/4 \[Beta] A[2]^2 B[1] - 
     1/4 \[Beta] B[1]^3 + 3 \[Beta] A[0] A[1] B[2] + 
     3/2 \[Beta] A[1] A[2] B[2] + 3/4 \[Beta] B[1] B[2]^2 - 
     9 \[Omega]^2 B[3] + \[Omega]0^2 B[3] + 3 \[Beta] A[0]^2 B[3] + 
     3/2 \[Beta] A[1]^2 B[3] + 3/2 \[Beta] A[2]^2 B[3] + 
     3/4 \[Beta] A[3]^2 B[3] + 3/2 \[Beta] B[1]^2 B[3] + 
     3/2 \[Beta] B[2]^2 B[3] + 3/4 \[Beta] B[3]^3 == 0};

eqns1 = eqns /. {\[Omega]0 -> 2 \[Pi] 10, \[Zeta] -> 1/100, \[Beta] ->
      5000, F -> 100., \[Omega] -> 2 \[Pi] 12.};

sols[1] = FindRoot[eqns1, Table[{uk[[i]], .1}, {i, Length[uk]}]]

X = Range[2 Pi 12, 2 Pi 13.5, .1]; nn = Length[X] - 1;

Do[sols[n + 1] = 
   FindRoot[
    eqns /. {\[Omega]0 -> 2 \[Pi] 10, \[Zeta] -> 1/100, \[Beta] -> 
       5000, F -> 100., \[Omega] -> X[[n + 1]]}, 
    Table[{uk[[i]], uk[[i]] /. sols[n]}, {i, Length[uk]}], 
    MaxIterations -> 1000, 
    Method -> {"Newton", "StepControl" -> "TrustRegion"}];, {n, nn}]

Visualization

Table[ListPlot[Transpose[{X, Table[uk[[j]] /. sols[i], {i, nn + 1}]}],
   PlotLabel -> uk[[j]], AxesLabel -> {"\[Omega]", ""}, 
  PlotStyle -> Red], {j, Length[uk]}]

Update 1. We can select the most likely solutions using NMinimize as follows

f = eqns[[All, 1]] /. {\[Omega]0 -> 2 \[Pi] 10, \[Zeta] -> 
     1/100, \[Beta] -> 5000, F -> 100, \[Omega] -> y};
X = Range[2 Pi 10, 2 Pi 20, .1]; sol = 
 Table[NMinimize[f . f , uk] // Quiet, {y, X}];

To select a solution, we use criteria

dsol = Transpose[{sol[[All, 1]], X, sol[[All, 2, All, 2]]}]

sdsol = Select[dsol, #[[1]] <= 10^-12 &];

Visualization

Table[ListLinePlot[Transpose[{sdsol[[All, 2]], sdsol[[All, 3, i]]}], 
  PlotLabel -> uk[[i]], PlotRange -> All, PlotStyle -> Red, 
  Frame -> True], {i, Length[uk]}]

Note that this is another branch of the solution, not similar to the first.

This is very helpful. Works well for some values but FindRoot runs out of precision for larger values of your X. However, better than NSolve. — Hugh, Commented Nov 10, 2023 at 18:59
@Hugh We can select the most likely solutions using NMinimize. See Update 1 to my answer. — Alex Trounev, Commented Nov 15, 2023 at 2:46

Daniel Huber · Accepted Answer · 2023-11-08 21:10:39Z

1

Setting WorkingPrecision->20, it then takes 85 sec. and gives only 3 solutions, although it gives a warning that the input is less precise than working precision:

Timing[sol = NSolve[eqns1, uk, Reals, WorkingPrecision -> 20];]
sol

{85.8125, Null}

answered Nov 8, 2023 at 21:10

Daniel Huber

56.4k2 gold badges24 silver badges66 bronze badges

$\begingroup$ Thanks. That's a big working precision. I did not try that large. However, the time of solution is now a problem. Also, finding the required precision. Do you know if we can do something to condition the equations beforehand so that multiple solutions are not present? $\endgroup$
– Hugh
Commented Nov 8, 2023 at 21:51
$\begingroup$ Things to try: FullSimplify of the equations before solving. Solve the first equation for B[3] and replace it in the rest of equations, then scond equation for B[2].. Note also, If you work with machine precision, WorkingPrecision is 16, therefore a value of 20 is no very large. $\endgroup$
– Daniel Huber
Commented Nov 9, 2023 at 7:59

Add a comment |

Goofy · Accepted Answer · 2023-11-15 04:13:36Z

1

With the OP's sol, we can gather the similar solutions together, find their mean, and if desired, polish them up with FindRoot:

FindRoot[eqns1, Transpose@{uk, #}] & /@ 
 Mean /@ FindClusters[uk /. sol]
(*
{{A[0] -> 0., 
  A[1] -> 0.522588, B[1] -> 0.46034, A[2] -> 0., 
  B[2] -> 0., A[3] -> -0.00538341, B[3] -> 0.00805524},
 {A[0] -> 0., 
  A[1] -> -0.512465, B[1] -> 0.404015, A[2] -> 0., B[2] -> 0., 
  A[3] -> 0.00335023, B[3] -> 0.00714414},
 {A[0] -> 0., 
  A[1] -> -0.0578141, B[1] -> 0.00317649, A[2] -> 0., B[2] -> 0., 
  A[3] -> -5.06698*10^-6, B[3] -> 8.73366*10^-7}}
*)

answered Nov 15, 2023 at 4:13

Goofy

3,68911 silver badges17 bronze badges

$\begingroup$ Could be an approach. Any idea if it is mathematically valid? Do sum of errors typically have zero mean? Thanks $\endgroup$
– Hugh
Commented Nov 16, 2023 at 14:43
$\begingroup$ @Hugh Depends on the distribution of errors and if they are random. If the distribution mean is zero, the expected value of the mean would be zero. And the standard deviation would be smaller. If the errors are systematic or biased, then no. But honestly, I wouldn't worry about that. Polishing with FindRoot is the cure. NSolve usually produces pretty good candidates for all the solutions. If the problem is ill-conditioned, use FindRoot with a high WorkingPrecision. If using FindRoot, then initializing with the first solution in each cluster is probably just as good as the mean $\endgroup$
– Goofy
Commented Nov 16, 2023 at 14:59

Add a comment |

Stack Exchange Network

Multiple solutions from NSolve

3 Answers 3

Your Answer

Not the answer you're looking for? Browse other questions tagged
equation-solving
numerics
or ask your own question.

Hot Network Questions

Multiple solutions from NSolve

3 Answers 3

Your Answer

Sign up or log in

Post as a guest

Not the answer you're looking for? Browse other questions tagged equation-solvingnumerics or ask your own question.

Related

Hot Network Questions

Not the answer you're looking for? Browse other questions tagged
equation-solving
numerics
or ask your own question.