19
$\begingroup$

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.

Improve this question
$\endgroup$

2 Answers 2

Reset to default
19
$\begingroup$

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

contours

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.

Improve this answer
$\endgroup$
2
  • 1
    $\begingroup$ A nice solution and still relevant and helpful, +1. $\endgroup$ Commented Nov 24, 2014 at 19:52
  • $\begingroup$ @Lou thanks - I saw there was a duplicate question just now, and marked it as such... $\endgroup$ Commented Nov 24, 2014 at 20:19
10
$\begingroup$

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:

Mathematica graphics

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

Mathematica graphics

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

Mathematica graphics

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:

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.