# How to find an approximate solution in a region when Solve or NSolve does not work?

Consider the relation

x/Tan[Pi x] == y/Tan[Pi y]


Solve or NSolve can't give a solution for it. But if you choose a particular region say $1.3 \leq x \leq 1.9$ and $.1 \leq y \leq .7$, there can be an approximate solution of the form $y=f(x)$. One way to get that would be to use ContourPlot, extract the data and then use Interpolate.

c = ContourPlot[x/Tan[Pi x] == y/Tan[Pi y], {x, 1.3, 1.9}, {y, 0.1, 0.7}];
pts = First@Cases[Normal[c], Line[i__] -> i, Infinity];
f = Interpolation[pts];


Is there any other way to get an approximate solution ($y=f(x)$) in a given region for a given relation?

Regarding this, I would be interested to know how ContourPlot works. Does it use any kind of root finding method to find all the points satisfy the relation and plot them? It that case will it be possible to use the same method to find an approximate solution as well?

-
You could use a series expansion in the region you are interested in and then solve a polynomial equation, like in Solve[Normal[ Series[x/Tan[Pi x] - y/Tan[Pi y], {x, 1.6, 1}, {y, 0.4, 1}]] == 0, y]. –  b.gatessucks Feb 7 at 13:12
You could say that ContourPlot uses a root finding method in a sense, but it's geared towards visualization, not precision. That's the biggest problem with it. There's no reliable precision (i.e. something similar to PrecisionGoal) control other than MaxRecursion. You can see the actual points it uses with the option Mesh -> All. –  Szabolcs Mar 9 at 15:16

You can use Solve.

f[x_] :=  With[
{sol = y /. Quiet[
Solve[x/Tan[Pi x] == y/Tan[Pi y] && Floor[x] - 1 < y < Floor[x], y],
Solve::ratnz]},
If[Length[sol] > 0, First[sol], Null]
];
Plot[f[x], {x, 1.4, 3},
PlotStyle -> Thick,
AspectRatio -> Automatic]


-
To anyone reading this: it's also worth looking at blog.wolfram.com/2008/12/18/… Generally, this will work if some bounds are specified for the roots, e.g. Floor[x] - 1 < y < Floor[x] in this case. –  Szabolcs Mar 9 at 23:16

You may use FindRoot with a list of x values. But I am not sure if this solution is better than yours though. Your question is very interesting in that how ContourPlot finds the exact line.

finder = Function[x, {x,
y /.FindRoot[x/Tan[Pi x] == y/Tan[Pi y], {y, x - 1.2}, AccuracyGoal -> 17]}]
f = Interpolation[finder /@ Range[1.3, 1.9, 0.01]]

Plot[Evaluate@f[x], {x, 1.3, 1.9}]


The plot is not regular. But at the range of your contour plot, I think it should have consistent values.

-

There was a very similar question on StackOverflow. I recommend you read my answer there.

I'm going to reproduce part of that answer here:

ContourPlot is designed for plotting. It does do something similar to root finding, but it's designed for visualization so the curves extracted form a ContourPlot won't be of a very high precision.

We can do better in two ways:

1. Extract the curves from ContourPlot and improve their precision using FindRoot. See the link for that.

2. Trace the contour ourselves, and again use FindRoot to obtain high precision solutions.

Let's do the second one now as it'll show you a method of tracing the contours:

(* just rotate a vector by 90 degrees *)
rot90[{x_, y_}] := {y, -x}

(* find a point where f[x,y]==0 by descending along f's gradient from {x0,y0} *)
step[f_, pt : {x_, y_}, pt0 : {x0_, y0_}, resolution_] :=
Module[
contourPoint =
Sow[contourPoint];
]


This step function takes an expression in terms of x,y and a starting point. If will move along the gradient of the expression to find a point where that expression is zero. Then it will take a step in a perpendicular direction to obtain a point which can be used a new starting point when tracing the contour.

You equation can written in the following implicit form: x Tan[Pi y] - y Tan[Pi x] == 0. We can use ContourPlot to visually obtain a starting point for our algorithm. {1.4, 0.1} looks like it's close enough to an interesting solution.

Then run:

result = Reap@NestList[step[x Tan[Pi y] - y Tan[Pi x], {x, y}, #, .08] &, {1.4, 0.1}, 50];


And plot the result:

ListPlot[{result[[1]], result[[-1, 1]]}, PlotStyle -> {Red, Black},
Joined -> True, AspectRatio -> Automatic, PlotMarkers -> Automatic]


The red points are used for tracing the contour but they're not precisely on it. The black points are precisely on the countour and are obtained with the red ones as starting points.

Note: ContourPlot doesn't work like this. I don't know for certain what it does, but I believe ContourPlot[f[x,y]==0, ...] will first compute the values of f[x,y] on a regular x,y grid to obtain a rectangular mesh. Then it will locally refine this mesh along the edges where f[x,y] changed sign. If you evaluate ContourPlot[f[x,y], ... , Mesh -> All] instead, you'll be able to see this mesh.

-