Getting an InterpolatingFunction from a ContourPlot
This can be done by extracting the point lists that ContourPlot has found.
c = ContourPlot[
D[min20[β][ϵ], ϵ] == 0, {β, 0,
10}, {ϵ, -1, 1}, Evaluated -> True];
lines = Cases[Normal[c], Line[i__] -> i, Infinity];
ParametricPlot[
Evaluate@Map[
Through[Map[ListInterpolation[#, {{0, 1}}] &, Transpose[#]][t]] &,
lines], {t, 0, 1}]
This has found three lines which I'm plotting as parametric plots together in one plot.
You can explicitly assign the interpolating functions to a list as follows:
fns = Map[Map[ListInterpolation[#, {{0, 1}}] &, Transpose[#]] &,
lines]
Then the second function would for example be evaluated like this:
Through[fns[[2]][.4]]
(* ==> {3.1785, 0.912361} *)
The results of ContourPlot are parametric curves - that's why the interpolating functions are also two-component lists (vectors) parametrized by a single variable that I arbitrarily chose to be t between 0 and 1.
Alternatives
Note that because of the parametric nature of the curves, we can't be sure which of the various lists of lines in ContourPlot belong together as one contiguous parametric curve - especially if your plot range happens to cut a single contour line into unnatural segments.
You can tell how the various curves are pieced together from the different colors they have in the ParametricPlot.
If you would prefer something that is not an interpolating function, you can still get somewhere with numerical solutions. For that, I'd simplify the expressions a little by defining the derivative deriv separately:
Clear[minimizeme, Ω, ϵ, β]
minimizeme[Ω_][β_][ϵ_] = ϵ^2 \
Ω -
Log[2 (Cosh[2 β] + Cosh[2 β ϵ])]/(2 β);
deriv[Ω_][β_][ϵ_] =
D[minimizeme[Ω][β][ϵ], ϵ]
(*
==> 2 ϵ Ω - Sinh[2 β ϵ]/(
Cosh[2 β] + Cosh[2 β ϵ])
*)
Then you could define a function that hopefully will yield one of the curves that ContourPlot found - I chose the right-most curve corresponding to the largest β:
g[Ω_][ϵ_] :=
Max[β /.
NSolve[deriv[Ω][β][ϵ] == 0, β, Reals]]
Now try this out:
g[.2][0.1]
Solve::ratnz: Solve was unable to solve the system with inexact coefficients. The answer was obtained by solving a corresponding exact system and numericizing the result.
1.08718
So it works except for the warning message. Now you have an actual function that will give you the values of $\beta$ for a given $\omega$ and $\epsilon\ne 0$ by numerical solution instead of interpolation.
Edit 2: what went wrong in your attempt
You also tried to use NSolve and it didn't work. So I should say why my version works and yours doesn't. The warnings that you mention are the same that I get too, and one can either turn them off or try to rewrite the equations to get rid of them. But then your version
Table[{β,
NSolve[D[min20[β][ϵ], ϵ] == 0, ϵ,
Reals]}, {β, 1.*^-10, 10}]
still doesn't evaluate to numerical solutions. NSolve is picky because its mission is to find all solutions, and when it knows that it can't do it, it doesn't do it. The only thing you really have to change is the added domain specification Reals at the end of NSolve. There apparently are complex solutions that "we know we don't know," so we have to rule them out. That's of course done by ContourPlot as well. What I do then in my solution is to take the list of rules produced by NSolve and select the one corresponding to the desired curve.
In general, if you can't get NSolve to do your bidding for reasons like this, it is usually possible to retreat to FindRoot which is less picky because it only needs to find a root, not all of them. The price that you have to pay is: FindRoot needs a starting value for the search, but it usually works well even if your starting guess is far off.
In your case, FindRoot should be the last resort because you have several different solution curves and FindRoot may jump back and forth between them erratically.
Introduction
I began to answer another question when it was closed as a duplicate of this one. Both questions tend to focus on what is seen in the contour plot, but the underlying question is how to represent a curve defined by an equation either as a function, that is, one variable as a function of the other, or as parameterized curve. I have used a method to answer several similar questions (links at end) that is more precise than simply using the data stored in plots. In general, the aim of plot functions is to produce an accurate picture. I would expect the error to be only somewhat less than a pixel. ContourPlot adaptively subdivides a mesh (see Specific initial sample points for 3D plots for some idea of what happens in other functions, or use Mesh -> All in ContourPlot) and then uses (linear?) interpolation between the mesh points to determine points on a contour.
Here I wish to give two general versions of the method: implicitFunction solves for an implicit function defined by the equation and implicitCurve that finds an arc-length parameterization of a contour. The first method is more appropriate for this question. This question poses a small difficulty due to a singularity in the OP's equation.
The basic idea is to set up a differential-algebraic system whose integral parametrizes the contour of an equation eqn in two variables x, y passing through an initial point {x0, y0}. It will have the form
NDSolve[
{eqn /. {x -> x[t], y -> y[t]},
{x[0], y[0]} == {x0, y0},
<ode for trajectory>}, (* different for the two methods *)
<functions>, (* y or {x, y} *)
{t, tmin, tmax}]
To solve y as a function x, implicitFunction uses an ODE such that x[t] = t. To get an arc-length (or unit-speed) parametrization, we add the constraint x'[t]^2 + y'[t]^2 == 1 and integrate a dummy arc-length variable s[t] satisfying the ODE s'[t] == 1.
It is important that the initial point {x0, y0} be on the contour. This could be programmed into the functions, but there are pathological equations for which one might want greater control over finding the initial point. So in the implicitFunction and implicitCurve, the user is responsible for passing a good starting point.
Usage and examples
Given a curve defined by an equation eqn in terms of variables x and y, a point {x0, y0} on the curve, and a range of parameter values, xmin <= x <= xmax or smin <= arclength <= smax respectively, the following return interpolating functions that describe the curve:
implicitFunction[eqn, {y, y0}, {x, x0, xmin, xmax}, options]
implicitCurve[eqn, {x, x0}, {y, y0}, {smin, smax}, options]
implicitFunction returns a single y -> InterpolatingFunction[<..>] that represents the function y = y[x] implicitly defined by eqn in a neighborhood of {x0, y0}. implicitCurve returns a pair {x -> InterpolatingFunction[<..>], y -> InterpolatingFunction[<..>]} that parameterizes the curve in a neighborhood of {x0, y0}. The options include NDSolve options and an option "Events" that can be used to specify WhenEvent code to be included in the call to NDSolve; if the variables x or y appear in an event, they are translated to the appropriate form x[t] or y[t].
Examples
The contour plot, to the human eye, looks as though it is composed of three curves:
There are singularities along the (horizontal) line ε == 0 where it crosses the other two curves. I will illustrate how to parametrized the rightmost, sideways U-shaped curve. One trick is to combine fractions and set the numerator equal to zero. To get rid of the ε == 0 branch, we need one more, to replace Sinh with Sinc.
minimizeme[Ω_][β_][ϵ_] = ϵ^2 Ω - Log[2 (Cosh[2 β] + Cosh[2 β ϵ])]/(2 β);
min20 = minimizeme[1/5];
num = D[min20[β][ϵ], ϵ] // Together // Numerator;
dmineq = num/ϵ == 0 /. Sinh[x_] :> x Sinc[I x] // Simplify
(* Cosh[2 β] + Cosh[2 β ϵ] == 5 β Sinc[2 I β ϵ] *)
Now we can find β = β[ε]. For an initial point, using ε == 0 causes errors, but another value works. FindRoot will find a point on the equation dmineq. The we'll integrate over -1 <= ε <= 1, but we'll use an event to stop when β == 10, just as in the contour plot. (We cannot integrate to ±1 in any case because of the asymptote.)
ϵ0 = 0.0001;
β0 = β /. FindRoot[dmineq /. ϵ -> ϵ0, {β, 1.08}];
{sol} = implicitFunction[
dmineq, {β, β0}, {ϵ, ϵ0, -1., 1.},
"Events" -> {WhenEvent[β == 10, "StopIntegration"]}]
(* {{β -> InterpolatingFunction[{{-0.977876, 0.977876}}, <>]}} *)
We can plot it over contour plot:
βifn = β /. sol
Show[
ContourPlot[D[min20[β][ϵ], ϵ] == 0, {β, 0, 10}, {ϵ, -1, 1}, Evaluated -> True],
ParametricPlot[{βifn[ϵ], ϵ}, Evaluate@Flatten[{ϵ, βifn["Domain"]}],
PlotStyle -> Red, PlotRange -> All]
]
We can also get an arc-length parametrization, which is given to show an example of the usage of implicitCurve. Here we will only integrate 5 units on either side of the starting point {β0, ε0}.
{sol} = implicitCurve[dmineq, {β, β0}, {ϵ, ϵ0}, {-5., 5.}]
Show[
ContourPlot[D[min20[β][ϵ], ϵ] == 0, {β, 0, 10}, {ϵ, -1, 1}, Evaluated -> True],
ParametricPlot[{β[t], ϵ[t]} /. sol,
Evaluate@Flatten[{t, β["Domain"] /. sol}], PlotStyle -> Red,
PlotRange -> All]
]
Code
The code for variableQ comes from this answer by Mr.Wizard.
variableQ = Quiet@ListQ@Solve[{}, #] &;
ClearAll[implicitFunction];
Options[implicitFunction] = Join[{"Events" -> {}}, Options[NDSolve]];
implicitFunction[
eqn_,
{y_?variableQ, y0_?NumericQ},
{x_?variableQ, x0_?NumericQ, xmin_?NumericQ, xmax_?NumericQ},
opts : OptionsPattern[]] :=
Module[{t, ode},
ode = {x'[t] == 1, (* x[t] == t *)
x[x0] == x0, y[x0] == y0};
NDSolve[
{eqn /. {x -> x[t], y -> y[t]},
ode,
OptionValue["Events"] /. {x -> x[t], y -> y[t]}},
y, {t, xmin, xmax},
FilterRules[{opts}, Options[NDSolve]]]
]
ClearAll[implicitCurve];
Options[implicitCurve] = Join[{"Events" -> {}}, Options[NDSolve]];
implicitCurve[eqn_,
{x_?variableQ, x0_?NumericQ},
{y_?variableQ, y0_?NumericQ},
{tmin_?NumericQ, tmax_?NumericQ},
opts : OptionsPattern[]] :=
Module[{t, ode},
ode = {x'[t]^2 + y'[t]^2 == 1, (* unit speed *)
{x[0], y[0]} == {x0, y0}, (* starting point *)
{x'[0], y'[0]} == Normalize @ (* starting direction *)
Cross[D[eqn /. Equal -> Subtract, {{x, y}}] /. Thread[{x, y} -> {x0, y0}]]};
NDSolve[
{eqn /. {x -> x[t], y -> y[t]},
ode,
OptionValue["Events"] /. {x -> x[t], y -> y[t]}},
{x, y}, {t, tmin, tmax},
FilterRules[{opts}, Options[NDSolve]]]
];
A utility function I sometimes use. It finds a point on eqn starting at {x1, y1} following the gradient (a variation on Newton's method). This can be extended to higher dimensions (see findPoint).
ClearAll[findRootGrad];
findRootGrad[eqn_, {x_?variableQ, x1_?NumericQ}, {y_?variableQ, y1_?NumericQ}] :=
With[{fn = eqn /. Equal -> Subtract,
grad = D[eqn /. Equal -> Subtract, {{x, y}}],
maxstepsize = 0.2},
Module[{fval = 1, norm = 1},
NestWhile[
{x, y} /. vars_ :> (Block[vars,
vars = #;
fval = fn;
norm = Norm[grad];
# -
Sign[fval] Min[Abs[fval/norm], maxstepsize] grad/norm] &),
{x1, y1},
Abs[fval/norm] > 1*^-10 &,
1,
100
]]
];
References
Other questions where I used similar ideas:
- Plotting implicitly-defined space curves
- Solve trigonometric equation
- How to find an approximate solution in a region when Solve or NSolve does not work?