7
$\begingroup$

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?

Improve this question
$\endgroup$
3
  • $\begingroup$ 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]. $\endgroup$ Commented Feb 7, 2014 at 13:12
  • 1
    $\begingroup$ 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. $\endgroup$ Commented Mar 9, 2014 at 15:16
  • $\begingroup$ You've accepted none of the answers to your questions as yet. Please remember to accept the answer, if any, that best solves your problem, by clicking the checkmark sign! And please leave feedback if none of them do. $\endgroup$ Commented Oct 26, 2014 at 17:04

4 Answers 4

Reset to default
4
$\begingroup$

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]

enter image description here

Improve this answer
$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$ Commented Mar 9, 2014 at 23:16
3
$\begingroup$

One can use NDSolve to construct interpolating functions that parametrize implicit curves. Below I'll show how to get x as a function of y and y as function of x (over appropriate domains, of course).

Unless the denominators of an equation are important, it can help to simplify the equation to eliminate them. In particular, the identity Sinc[u] == u / Sin[u] is extremely helpful in this equation as it removes the discontinuties. (Note in explanation: The curious looking replacement t_Cot :> t, which appears to do nothing, does this: It matches the Tan expressions from eqn, which are automatically converted to Cot expressions. This prevents the second pattern u : x | y from matching the x and y inside the Cot expressions, since they matched the first rule. A cute trick I felt forced to resort to.)

eqn = x/Tan[π x] == y/Tan[π y]; denoms = Sinc[π x] Sinc[π y]; eqn2 = denoms * (eqn /. {t_Cot :> t, u : x | y :> Sin[π u]/Sinc[π u], Equal -> Subtract}) == 0 // Simplify (* Cos[π y] Sinc[π x] == Cos[π x] Sinc[π y] *) ContourPlot[Evaluate@eqn2, {x, -4, 4}, {y, -4, 4}]

Mathematica graphics

We can use the plot to find initial points.

ivReg = RegionNearest[DiscretizeGraphics @ ContourPlot[Evaluate@eqn2, {x, -4., 4.}, {y, -4., 4.}]]

We can find a point on the curve using ivReg to find a starting point for FindRoot. We can then use the point to set up a simple DAE for NDSolve in which y[t] == t and x[t] essentially is equivalent to x as a function y. (In another region, one might want y as a function of x, x as a function of y is better in this part, judging from the contour plot.)

p0 = ivReg[{1, 0}]; y0 = Last[p0]; x0 = x /. FindRoot[eqn2 /. y -> y0, {x, First[p0]}]; xfn = First @ NDSolve[ {eqn2 /. {x -> x[t], y -> y[t]}, y'[t] == 1, x[y0] == x0, y[y0] == y0}, x, {t, -3, 3}]

The interpolating function in xfn has a domain from -3 to 3 but I'll plot just half of it.

ParametricPlot[{x[t], t} /. xfn // Evaluate, {t, 0, 3}]

Mathematica graphics

For those who demand x as a function of y, then here it is:

p0 = ivReg[{1.4, 0.1}]; (* avoid the singularity at y == 0 *) x0 = First[p0]; y0 = y /. FindRoot[eqn2 /. x -> First[p0], {y, Last[p0]}]; yfn = NDSolve[ {eqn2 /. {x -> x[t], y -> y[t]}, x'[t] == 1, x[x0] == x0, y[x0] == y0}, y, {t, 1.4, 3}, "ExtrapolationHandler" -> {Indeterminate &} (* optional *) ]

One gets a NDSolve::ndsz: warning, of course, when the integration hits the x axis. All it means is that the integration has reached the lower limit of the domain. One could Quiet it, I suppose, I prefer to see it, since it confirms my understanding of the problem. Since I let it find the domain automatically, I don't know what it is. I could either extract it from the InterpolatingFunction(see Interpolating Function Anatomy, or, as I've done here, I can make extrapolated values all return Indeterminate. This works nicely with Plot. One can also turn off the extrapolation messages with "ExtrapolationHandler" -> {Indeterminate &, "WarningMessage" -> False}.

ParametricPlot[{t, y[t]} /. yfn // Evaluate, {t, 1.43, 3}]

InterpolatingFunction::dmval: Input value {1.43003} lies outside the range of data in the interpolating function. Extrapolation will be used. >>

Mathematica graphics

Improve this answer
$\endgroup$
2
$\begingroup$

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[ {grad, grad0, t, contourPoint}, grad = D[f, {pt}]; grad0 = grad /. Thread[pt -> pt0]; contourPoint = grad0 t + pt0 /. First@FindRoot[f /. Thread[pt -> grad0 t + pt0], {t, 0}]; Sow[contourPoint]; grad = Normalize[grad /. Thread[pt -> contourPoint]]; contourPoint + rot90[grad] resolution ]

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.

Improve this answer
$\endgroup$
0
$\begingroup$

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}]

enter image description here

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

Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.