How many solutions does this system have?

Question

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[%]

one thousand three hundred and fifty-three

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

Answer 1

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)

Answer 2

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.

Answer 3

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 one 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 one thousand six hundred and seventeen 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}} one thousand six hundred and seventeen