# How many solutions does this system have?

Let us consider the system of two transcendental equations x1 - x1*Sin[x1 + 5*x2] - x2*Cos[5 x1 - x2] == 0 && x2 - x2*Sin[5 x1 - 3 x2] + x1*Cos[3 x1 + 5 x2] == 0 and its solutions belonging to the square x1 >= -20 && x1 <= 20 && x2 >= -20 && x2 <= 20. The plot

ContourPlot[{x1 - x1*Sin[x1 + 5*x2] - x2*Cos[5 x1 - x2] == 0,
x2 - x2*Sin[5 x1 - 3 x2] + x1*Cos[3 x1 + 5 x2] == 0}, {x1, -20, 20}, {x2, -20, 20},
PlotPoints -> 50]


suggests the number is big (It's unclear even whether the set of the solutions is finite.). These are my attempts.

NSolve[x1 - x1*Sin[x1 + 5*x2] - x2*Cos[5 x1 - x2] == 0 &&
x2 - x2*Sin[5 x1 - 3 x2] + x1*Cos[3 x1 + 5 x2] == 0 && x1 >= -20 &&
x1 <= 20 && x2 >= -20 && x2 <= 20, {x1, x2}, Reals] // Dimensions


{195, 2}

and the warning "NSolve::incs: Warning: NSolve was unable to prove that the solution set found is complete".

Splitting the square into smaller ones, I obtain

d = 5; Table[Dimensions[NSolve[x1 - x1*Sin[x1 + 5*x2] - x2*Cos[5 x1 - x2] == 0 &&
x2 - x2*Sin[5 x1 - 3 x2] + x1*Cos[3 x1 + 5 x2] == 0 && x1 >= k &&
x1 <= k + d && x2 >= n && x2 <= n + d, {x1, x2}, Reals]][[ 1]], {k, -20, 20 - d, d}, {n, -20, 20 - d, d}] // Flatten


{27, 28, 21, 2, 8, 23, 29, 33, 35, 31, 23, 0, 5, 25, 29, 29, 34, 28, 30, 26, 20, 30, 27, 32, 2, 0, 27, 9, 7, 0, 0, 20, 14, 15, 30, 10, 19, 19, 15, 16, 30, 34, 30, 4, 9, 31, 32, 30, 30, 32, 25, 0, 10, 26, 31, 27, 30, 25, 24, 14, 14, 25, 29, 33}

and the same warnings. The execution lasted a few hours since I can't efficiently use ParallelTable on my comp.

Total[%]


1353

The case d = 2; : the same warnings and 1580 solutions. So the questions arise: how many solutions does this system have? how to deduce it?

• Dealing with a partition function, we need to know its zeros. Commented May 17 at 13:05
• Why not to add option WorkingPrecision -> 30 to NSolve? :) Commented May 17 at 18:22
• @AlexTrounev:NSolve[x1 - x1*Sin[x1 + 5*x2] - x2*Cos[5 x1 - x2] == 0 && x2 - x2*Sin[5 x1 - 3 x2] + x1*Cos[3 x1 + 5 x2] == 0 && x1 >= -20 && x1 <= 20 && x2 >= -20 && x2 <= 20, {x1, x2}, Reals, WorkingPrecision -> 25] // Dimensions results in the same {195, 2}. Commented May 17 at 20:30
• Just a clarification question, do you want to count the root $x1=0, x2=\pi/2$ as a single or double root? I guess for partition function you need to take multiplicities into account, isn't it? Commented May 18 at 21:17
• @yarchik: Thank you for your valuable question. Googling "zeros of partition function", I don't find an answer. Commented May 19 at 3:42

Try ContourPlot and MeshFunctions in order to estimate numerically the number of solutions:

f[x1_, x2_] := x1 - x1*Sin[x1 + 5*x2] - x2*Cos[5 x1 - x2]
g[x1_, x2_] := x2 - x2*Sin[5 x1 - 3 x2] + x1*Cos[3 x1 + 5 x2]
picS = ContourPlot[{f[x1, x2] == 0(*, g[x1, x2] == 0*)}, {x1, - 20,20}, {x2, - 20, 20},
MeshFunctions -> {f[#1, #2] &, f[#1, #2] - g[#1, #2] &},Mesh -> {{0}}, MeshStyle -> Black,PlotPoints->200 ] (* needs ~20s to evaluate*)


The number of Point' s found in picS gives the number of intersections(solutions)

punkte = picS[[1, 1]][[1]];
pS = Cases[picS, Point[p_] :> punkte[[p]], -1][[1]];
Length[pS] (*1613 solutions*)

Show[ContourPlot[{f[x1, x2] == 0 , g[x1, x2] == 0 }, {x1, - 20,20}, {x2, - 20, 20}, PlotPoints -> 200],
ListPlot[pS , PlotStyle -> {  PointSize[Medium], Black}]]


The result depends on the accuracy of ContourPlot (options PlotPoints, MaxRecursion)

• +1. Thank you. This is helpful, though does not answer the question. Commented May 17 at 12:41
• It gives an numerical estimation of the number of solutions, not more! Commented May 17 at 12:54
• It does not seem like there should be 3226 points. I guess it is too much. Commented May 17 at 13:05
• 1617 is a very good estimation I think. Commented May 17 at 17:22
• Actually ContourPlot it is not a method to compute roots :) Commented May 17 at 19:21

We can try FindRoot as follows

h=1/4;grid = Range[-20, 20, h]; s =
DeleteDuplicates@
Flatten[ParallelTable[{x1, x2} /.
Quiet@FindRoot[
x1 - x1*Sin[x1 + 5*x2] - x2*Cos[5  x1 - x2] == 0 &&
x2 - x2*Sin[5  x1 - 3  x2] + x1*Cos[3  x1 + 5  x2] == 0, {{x1,
grid[[i]]}, {x2, grid[[j]]}}, WorkingPrecision -> 30], {i,
Length[grid]}, {j, Length[grid]}], 1]


In this example we have

s // Dimensions

Out[]= {10313, 2}


But some of the roots lie outside the region Rectangle[{-20, -20}, {20, 20}]. Then we can decries step size from h=1/4 to h=1/8 and increase from 1/4 to 1/2 and look if there is convergence. For h=1/2 we have {3855, 2}, and for h=1/8 we have {34175, 2} (it takes about 12m on my Silver Pentium). As we can see there is no convergence at all. Therefore, we need some deep analytical research, since number of roots increases with h decreases - see Figure below.

Using Select we can take roots from Rectangle[{-20, -20}, {20, 20}] only,

s0 = Select[s, Abs[#[[1]]] <= 20 && Abs[#[[2]]] <= 20 &]


For this number we have

s0 // Dimensions

Out[]= {9085, 2}


For h=1/2,1/8 numbers are {3423, 2}, and {30484, 2} consequently. Roots shown below

Update 1. To find all branches and exclude some of them we use Reduce as follows

 Reduce[{x1 - x1*Sin[a] - x2*Cos[b] == 0,
x2 - x2*Sin[c] + x1*Cos[d] == 0}, {x1, x2}]

Out[]= (Sin[c] == 1 && x1 == 0 && x2 == 0 &&
Cos[b] Cos[d] != 0) || (Sin[c] == 1 && Cos[b] == 0 && x1 == 0 &&
Cos[d] != 0) || (Sin[c] == 1 && Cos[d] == 0 && Cos[b] != 0 &&
x2 == Sec[b] (x1 - x1 Sin[a])) || (Sin[c] == 1 && Sin[a] == 1 &&
Cos[d] == 0 && Cos[b] == 0) || (Cos[d] != 0 &&
Cos[b] ==
Sec[d] (-1 + Sin[a] + Sin[c] - Sin[a] Sin[c]) && -1 + Sin[c] !=
0 && x2 == (x1 Cos[d])/(-1 + Sin[c])) || (-1 + Sin[c] != 0 &&
Sin[a] == 1 && Cos[d] == 0 &&
x2 == 0) || (1 + Cos[b] Cos[d] - Sin[a] - Sin[c] + Sin[a] Sin[c] !=
0 && x1 == 0 && -1 + Sin[c] != 0 && x2 == 0) || (Sin[c] == 1 &&
Cos[d] == 0 && Cos[b] == 0 && -1 + Sin[a] != 0 && x1 == 0)


Where {a = x1 + 5*x2, b = 5 x1 - x2, c = 5 x1 - 3 x2, d = 3 x1 + 5 x2}. First 4 branches describe trivial solution $$x1=0, x2=0$$. We interested in nontrivial solution that described by branch 5

(Cos[d] != 0 &&
Cos[b] ==
Sec[d] (-1 + Sin[a] + Sin[c] - Sin[a] Sin[c]) && -1 + Sin[c] != 0 &&
x2 == (x1 Cos[d])/(-1 + Sin[c]))


As we can see it is a nontrivial solution of homogenous linear system with zero discriminant that we can represent with

{vec, mat} =
CoefficientArrays[{x1 - x1*Sin[a] - x2*Cos[b] == 0,
x2 - x2*Sin[c] + x1*Cos[d] == 0}, {x1, x2}];
mat // Normal

(*Out[]= {{1 - Sin[a], -Cos[b]}, {Cos[d], 1 - Sin[c]}}*)
With[{a = x1 + 5*x2, b = 5 x1 - x2, c = 5 x1 - 3 x2,
d = 3 x1 + 5 x2}, Det[{{1 - Sin[a], -Cos[b]}, {Cos[d], 1 - Sin[c]}}]]

Out[]= 1 + Cos[5 x1 - x2] Cos[3 x1 + 5 x2] - Sin[5 x1 - 3 x2] -
Sin[x1 + 5 x2] + Sin[5 x1 - 3 x2] Sin[x1 + 5 x2]


Finally we can compute nontrivial roots with FindRoot using norm of equations to select right solutions

grid1 = Range[-20, 20, 1]; snorm1 =
Flatten[ParallelTable[{x1, x2,
Norm[{1 + Cos[5 x1 - x2] Cos[3 x1 + 5 x2] - Sin[5 x1 - 3 x2] -
Sin[x1 + 5 x2] + Sin[5 x1 - 3 x2] Sin[x1 + 5 x2],
x2 - (x1 Cos[3 x1 + 5 x2])/(-1 + Sin[5 x1 - 3 x2])}]} /.
Quiet@FindRoot[{1 + Cos[5 x1 - x2] Cos[3 x1 + 5 x2] -
Sin[5 x1 - 3 x2] - Sin[x1 + 5 x2] +
Sin[5 x1 - 3 x2] Sin[x1 + 5 x2] == 0,
x2 - (x1 Cos[3 x1 + 5 x2])/(-1 + Sin[5 x1 - 3 x2]) == 0}, {{x1,
grid1[[i]]}, {x2, grid1[[j]]}}, WorkingPrecision -> 30], {i,
Length[grid1]}, {j, Length[grid1]}], 1]


Selection

sol1 = DeleteDuplicates@Select[snorm1, #[[3]] < 10^-16 &];

Length[sol1]

(*Out[]= 425*)


Using 4 grids with h=1,1/2,1/4,1/8 we plot results as

As we can see number of roots increasing with h decreasing even we select the reliable solutions only with norm of equations <10^-16.

Update 2 Ignoring small differences we can select roots as follows

grid = Range[-20, 20, #] & /@ {1/5, 1/10};
snorm[k_] :=
Flatten[ParallelTable[{x1, x2,
Norm[{1 + Cos[5  x1 - x2]  Cos[3  x1 + 5  x2] -
Sin[5  x1 - 3  x2] - Sin[x1 + 5  x2] +
Sin[5  x1 - 3  x2]  Sin[x1 + 5  x2],
x2 - (x1  Cos[3  x1 + 5  x2])/(-1 + Sin[5  x1 - 3  x2])}]} /.
Quiet@FindRoot[{1 + Cos[5  x1 - x2]  Cos[3  x1 + 5  x2] -
Sin[5  x1 - 3  x2] - Sin[x1 + 5  x2] +
Sin[5  x1 - 3  x2]  Sin[x1 + 5  x2] == 0,
x2 - (x1  Cos[3  x1 + 5  x2])/(-1 + Sin[5  x1 - 3  x2]) ==
0}, {{x1, grid[[k, i]]}, {x2, grid[[k, j]]}},
WorkingPrecision -> 30], {i, Length[grid[[k]]]}, {j,
Length[grid[[k]]]}], 1]

sol = snorm[#] & /@ {1, 2}; // AbsoluteTiming

sol1 =
DeleteDuplicates@
Select[sol[[1]], #[[3]] < 10^-16 && Abs[#[[1]]] <= 20 &&
Abs[#[[2]]] <= 20 &]; Length[sol1]

(*Out[]= 2462*)

sol2 =
DeleteDuplicates@
Select[sol[[2]], #[[3]] < 10^-16 && Abs[#[[1]]] <= 20 &&
Abs[#[[2]]] <= 20 &]; Length[sol2]

(*Out[]= 5522*)


Now we apply filter to select unique roots and so ignoring multiple roots

sol11 = DeleteDuplicates[sol1, Norm[#1 - #2] < 10^-15 &]

sol21 = DeleteDuplicates[sol2, Norm[#1 - #2] < 10^-15 &]
Length[#] & /@ {sol11, sol21}

Out[]= {1385, 1421}


These numbers look as reasonable and maybe we can reach number 1617 computed with other methods.

• I tried FindRoot, and I have found some issues that you have overlooked. First, some solutions provided by FindRoot are not the roots of the original equation, you have to check them. Second, DeleteDuplicates  cannot delete the same root because they could be near but not the same. Instead, I use DeleteDuplicates[#, Norm[#1 - #2] < 10^-10 &] &. Commented May 17 at 16:48
• @AlexTrounev So what is your conclusion? How many roots are there? Your method is not correct as you use Quiet that just ignores warning messages. With each warning message you introduce another incorrect root, so there is no surprise that number of your found "roots" increases . But that was already pointed out in the comment by Jie Zhu. Commented May 19 at 11:59
• @azerbajdzan You don't pay attention that I used norm of equations to test solutions, therefore, with my method we compute all roots with error less than 10^-16, it is much better then that with NSolve as in your answer or with ContourPlot . My conclusion is that we can't numerically define number of roots, so we need some analytical research. :) Commented May 19 at 21:00
• @AlexTrounev Again, have you read the comment by Jie Zhu? Your DeleteDuplicates in last section "Selection" is there probably for decoration because it does nothing. Do you understand how DeleteDuplicates works? Do it properly by the code provided in Jie Zhu's comment. Your final result is full of duplicates that were not properly deleted. For example with h=1 the root sol1[[6]] is same as root sol1[[25]], yet you take it as two separate roots. Commented May 19 at 21:37
• @AlexTrounev Oh boy, you really do not understand what precision means. Evaluate FindRoot[{x1 - x1*Sin[x1 + 5*x2] - x2*Cos[5 x1 - x2] == 0, x2 - x2*Sin[5 x1 - 3 x2] + x1*Cos[3 x1 + 5 x2] == 0}, Evaluate[Sequence @@ Transpose[{{x1, x2}, sol1[[6, {1, 2}]]}]], WorkingPrecision -> 1000] FindRoot[{x1 - x1*Sin[x1 + 5*x2] - x2*Cos[5 x1 - x2] == 0, x2 - x2*Sin[5 x1 - 3 x2] + x1*Cos[3 x1 + 5 x2] == 0}, Evaluate[Sequence @@ Transpose[{{x1, x2}, sol1[[25, {1, 2}]]}]], WorkingPrecision -> 1000] to see that roots sol1[[6]] and sol1[[25]] are the same. Commented May 20 at 12:58

This is a exact and "error/warning-messages-free" method.

Firstly, there are roots that can be expressed exactly as {x1, x2} -> {(4 m - n) Pi/2, (1 + n) Pi/2}. I am not going into detail how I found them as it is not necessary for the answer to the question. But we can verify it by:

FullSimplify[{x1 - x1*Sin[x1 + 5*x2] - x2*Cos[5 x1 - x2] == 0,
x2 - x2*Sin[5 x1 - 3 x2] + x1*Cos[3 x1 + 5 x2] == 0} /.
Thread[{x1, x2} -> {(4 m - n) Pi/2, (1 + n) Pi/2}],
Assumptions -> Element[{m, n}, Integers]]


{True, True}


Some of those exact roots are singular, namely those with x1==0 or x2==0. These singular roots are the cause of problems when trying to find a root around these singular points.

The next code searches for roots in squares m <= x1 < m + 1, n <= x2 < n + 1 if it takes longer than 1 second the region is added to fail if not the found roots are added to roots.

Those failed region are regions that contain singular roots so then we search again those regions but avoiding x1==0 or x2==0. Those roots are added to rootsavoid.

Finally we search for singular roots in failed regions by explicitly setting x1==0 or x2==0. These roots are added to rootssingular.

Then all found roots are joined in allroots = Join[roots, rootsavoid, rootssingular].

The length of allroots is 1617 surprisingly exactly same as @UlrichNeumann's answer (when PlotPoints -> 220 is used).

The code monitors {m,n}. At output we see failed regions and total number of roots.

Clear[m, n]
roots = {};
fail = {};
Monitor[Do[
TimeConstrained[
AppendTo[roots,
NSolve[{x1 - x1*Sin[x1 + 5*x2] - x2*Cos[5 x1 - x2] == 0,
x2 - x2*Sin[5 x1 - 3 x2] + x1*Cos[3 x1 + 5 x2] == 0,
m <= x1 < m + 1, n <= x2 < n + 1}]], 1,
AppendTo[fail, {m, n}]], {m, -20, 19}, {n, -20, 19}], {m, n}]
roots = Values[Flatten[roots, 1]];
fail

avoidzero = fail /. {-1 -> -1 - 1/1000, 0 -> 0 + 1/1000};
rootsavoid =
Flatten[NSolve[{x1 - x1*Sin[x1 + 5*x2] - x2*Cos[5 x1 - x2] == 0,
x2 - x2*Sin[5 x1 - 3 x2] + x1*Cos[3 x1 + 5 x2] == 0,
m <= x1 < m + 1, n <= x2 < n + 1} /. Thread[{m, n} -> #]] & /@
avoidzero, 1];

singular = fail /. {-1 -> 0} // Union;
rootssingular =
Flatten[Join[
NSolve[{x1 - x1*Sin[x1 + 5*x2] - x2*Cos[5 x1 - x2] == 0,
x2 - x2*Sin[5 x1 - 3 x2] + x1*Cos[3 x1 + 5 x2] == 0,
x1 == 0, n <= x2 < n + 1} /. Thread[{m, n} -> #]] & /@
Cases[singular, {0, _}],
NSolve[{x1 - x1*Sin[x1 + 5*x2] - x2*Cos[5 x1 - x2] == 0,
x2 - x2*Sin[5 x1 - 3 x2] + x1*Cos[3 x1 + 5 x2] == 0,
m <= x1 < m + 1, x2 == 0} /. Thread[{m, n} -> #]] & /@
Cases[singular, {_, 0}]], 1] // Values;

allroots = Join[roots, rootsavoid, rootssingular];
Length[allroots]


{{-18, -1}, {-18, 0}, {-11, -1}, {-11, 0}, {-5, -1}, {-5,
0}, {-1, -18}, {-1, -11}, {-1, -5}, {-1, 1}, {-1, 7}, {-1,
14}, {0, -18}, {0, -11}, {0, -5}, {0, 1}, {0, 7}, {0,
14}, {1, -1}, {1, 0}, {7, -1}, {7, 0}, {14, -1}, {14, 0}}

1617

• Thank you for your work. As I understand it, some solutions are of the form {(4 m - n) Pi/2, (1 + n) Pi/2} with Element[{m, n}, Integers] and you find other solutions numerically. Can you kindly elaborate your "Some of those exact roots are singular, namely those with x1==0 or x2==0. These singular roots are the cause of problems when trying to find a root around these singular points", giving us details? Concerning the comparison with the UlrichNeumann's answer, if PlotPoints -> 240 are taken, I obtain 1613 solutions by his code and this is noticed in my comment to his answer. Commented May 18 at 19:56
• @azerbajdzan Nice attempt (+1), but I think that we can't compute all roots with NSolve. See Update 1 to my answer with FindRoot application on 4 grids. Commented May 19 at 11:47
• @AlexTrounev Surely we can when used wisely. Give me a root that is missing in my list of 1617 roots stored in the variable allroots. Commented May 19 at 12:02
• @azerbajdzan: You wrote "Yes, the set is complete - you do not seem to understand the code when you are asking". Sorry, I prefer arguments over empty words. Commented May 19 at 16:01
• @azerbajdzan The grid form your "guessed" analytical solution is much smaller than the grid you use in your following calculation ...,{m,-20,19},{n,-20,19}] , I don't understand why that should be correct Commented May 19 at 16:39