Getting an InterpolatingFunction from a ContourPlot

Question

I have a function, say

minimizeme[Ω_][β_][ϵ_] = ϵ^2 Ω - Log[2 (Cosh[2 β] + Cosh[2 β ϵ])]/(2 β);

I want to find its critical points in $\epsilon$ for a given value of $\omega$, varying $\beta$, within certain ranges, i.e. $\omega\geq 0$, $\beta\geq 0$ and $-1\leq \epsilon \leq 1$. When I pick, say, $\omega$ equal to $0.2$, e.g. I define

min20 = minimizeme[0.20];

I can make a ContourPlot of the locus of couples $(\beta, \epsilon)$ for which $\frac{\partial m}{\partial \epsilon }=0$, using code

ContourPlot[D[min20[β][ϵ], ϵ] == 0, {β, 0, 10}, {ϵ, -1, 1}, Evaluated -> True]

while I cannot get the software to solve the related equation:

Solve[D[min20[β][ϵ], ϵ] == 0, β, MaxExtraConditions -> All]

yields a Solve::inex: message, while

Solve[D[minimizeme[Ω][β][ϵ], ϵ] == 0, β, MaxExtraConditions -> All]

yields Solve::nsmet:. Furthermore,

Table[NSolve[D[min20[β][ϵ], ϵ] == 0, ϵ], {β, 1.*^-10, 10}]

results in a few NSolve::nsmet: until General::stop: kicks in.

Finally, Reduce behaves just like Solve, resulting in either a Reduce::inex: or a Reduce::nsmet: message. Yet I think that since ContourPlot manages to produce a result somehow, Mathematica has a way to find what I need. So, to summarize:

1. How does ContourPlot work "internally"? How does it manage to produce the results I see when even Solve, NSolve and Reduce fail?
2. Can I turn the results of this ContourPlot into one (or more) InterpolatingFunction objects?
3. (bonus question) Are there any alternatives for solving my equation beyond Taylor series expansion? (I have already solved this with Taylor expansion, so I'd like to avoid that at all if possible)

Thanks a lot.

Answer 1

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.