4
$\begingroup$

I have an equation of $x$ and $y$, which reads $$ -\alpha^2+\frac{2}{3} \alpha^3+y^2-\frac{2}{3}y^3+x^3f(\frac{y}{x})=0, $$

where $\alpha$ is positive parameter and $f(x)=\int (1-\tanh x) x^2\mathrm{d}x$, with the arbitrary constant of integration defined by $f(\infty)=0$. With Mathematica, $f(x)$ is just what Mathematica gives you.

Now I want to obtain the relation between $x$ and $y$ from that equation, with $x\in(0,1/\ln 2)$, the parameter $\alpha$ should be tested for $\alpha\in(0,5)$, say, especially the values around $2.4$

My Attempt

My attempt is to use FindRoot to do that job, however no matter how I increase the WorkingPrecision, there are complains saying:

The line search decreased the step size to within tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find a sufficient decrease in the merit function. You may need more than 100 digits of working precision to meet these tolerances.

The curve obtained is just wearied, which I can't trust. For your convenience, I will put my attempt code:

ClearAll["Global`*"];
α = 15/10;
stepSize = 1/100;
c0 = Log[2];
eps = 1/100000;
f[x_] = Integrate[(1 - Tanh[x]) x^2, x]; (*f(∞)=0*)
list = 
  Table[
    {c0 x, 
     y /. FindRoot[
            -α^2 + 2/3*α^3 + y^2 - 2/3 y^3 + x^3 f[y/x] == 0, 
            {y, 0.5}, 
            WorkingPrecision -> 100]}, 
    {x, eps, 1/c0, stepSize}];
ListLinePlot[list]

The above code gives a wearied curve as shown below and a lot of complains as mentioned above.

enter image description here

I've also tried the fix $y$ to solve for $x$ instead. However, the code doesn't works well (for some parameters), either.

$\endgroup$
  • 1
    $\begingroup$ There is more than one y for each x, and depending on the initial guess, you will jump between those different roots. Here is a contour plot of your function for the choice of alpha shown in your post, with the zero-contour only shown. $\endgroup$ – march Oct 30 '15 at 5:04
5
$\begingroup$

Depending on what you want, here's how I would go about this. First, using your choice of α, define

α = 15/10;
f[x_] = Integrate[(1 - Tanh[x]) x^2, x]
(* -x^2 Log[1+E^(-2 x)]+x PolyLog[2,-E^(-2 x)]+1/2 PolyLog[3,-E^(-2 x)] *)
g[x_, y_] = -α^2 + 2/3*α^3 + y^2 - 2/3 y^3 + x^3 f[y/x];

Generate the zero-contours of your function for your choice of α using ContourPlot:

p1 = ContourPlot[g[x, y] /. α -> 15/10
  , {x, 0, 1/Log[2]}, {y, -5, 5}
  , Contours -> {0}, ContourShading -> None]

enter image description here

This might be all you want! If you want the points, extract them:

pts = First@Cases[Normal@p1, Line[a_] :> a, Infinity];
ListLinePlot@pts

enter image description here

These aren't perfect, so let's use these points as initial guesses for FindRoot to refine the values:

newPts = {#1, y /. FindRoot[g[#1, y] /. α -> 15/10, {y, #2}]} & @@@ pts;
ListLinePlot@newPts

enter image description here

$\endgroup$
  • $\begingroup$ I need to carefully read this, your results seems very good! In the mean time, can you explain what zero contour means? $\endgroup$ – an offer can't refuse Oct 30 '15 at 5:27
  • $\begingroup$ I just mean that it's the level-curve of the the function that you're setting to zero (which I'm calling g) that corresponds to the value 0. In other words, it's the graph of the solution set of g[x, y] == 0. $\endgroup$ – march Oct 30 '15 at 5:42
  • $\begingroup$ pts = First@Cases[Normal@p1, Line[a_] :> a, Infinity]; In this code, you select the line (points) in the objects Normal@p1, what other things other than line does this Normal@p1 have? $\endgroup$ – an offer can't refuse Oct 30 '15 at 8:54
  • 1
    $\begingroup$ Everything in Mathematica is an expression, including Graphics. Evaluate p1 and then evaluate p1 // FullForm, and you will see all the stuff associated with that Graphics object. Doing Normal@p1 first rearranges the innards of the expression to more directly reflect what's being drawn in the figure; this is necessary because Graphics expressions typically are written in terms of GraphicsComplexes, which (I think) is for the purposes of compression. Anyway, I recommend playing! Make a simple like Plot[{x^2, 1+x^2}, {x, 0, 1}] //Normal//FullForm and examine the results! $\endgroup$ – march Oct 30 '15 at 17:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.