# Multiple solutions from NSolve

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" 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
• Why not to use FindRoot? Commented Nov 8, 2023 at 14:34
• @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

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 At X>13.5 we can try smaller step like 0.01. Commented Nov 11, 2023 at 4:32
• I just tried smaller steps. No help. Pity.
– Hugh
Commented Nov 11, 2023 at 19:15
• @Hugh We can select the most likely solutions using NMinimize. See Update 1 to my answer. Commented Nov 15, 2023 at 2:46
• Thanks. Another solution.
– Hugh
Commented Nov 16, 2023 at 18:51

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}

• 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?
– Hugh
Commented Nov 8, 2023 at 21:51
• 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. Commented Nov 9, 2023 at 7:59

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}}
*)

• Could be an approach. Any idea if it is mathematically valid? Do sum of errors typically have zero mean? Thanks
– Hugh
Commented Nov 16, 2023 at 14:43
• @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 Commented Nov 16, 2023 at 14:59