I would like to solve the following boundary value problem for $y(x)$ for a fixed value of $k$ between $0 < k <1$:

$$y'' + \frac{3}{x} y' - y + \frac{3}{2}y^2 - \frac{k}{2}y^3=0 \\ y'(0) = 0,\qquad y(\infty)\rightarrow0$$

To make this work, I truncate the domain of the problem to $x_0 < x < x_\text{Max}$ to avoid the singularity at $x=0$ and so that the second boundary condition is brought to a finite $y(x_\text{Max})=0$. I will be working with $x_0 = 0.01$ and $x_\text{Max} = 10.$

I want to write a fast Mathematica program to solve this by using undershoot/overshoot on the related initial-value problem $y(x_0) = y_0$ to achieve the second boundary condition $y(x_\text{Max}) = 0$.

So far, I wrote a nice little code that solves the initial-value problem for a given $k$ and $y_0$:

sol = With[{x0 = .01, xMax = 10},
          {y''[x] + 3/x y'[x] - y[x] + 3/2 y[x]^2 - k/2 y[x]^3 == 0 (*diff.eq.*),
           y[x0] == y0 + 1/8 x0^2 (y0 - 3/2 y0^2 + k/2 y0^3)        (*init.cond.1*),
           y'[x0] == 1/4 x0 (y0 - 3/2 y0^2 + k/2 y0^3)              (*init.cond.2*),
           WhenEvent[y[x] == 0, {Print["Overshot"], "StopIntegration"}],
           WhenEvent[y'[x] == 0, {Print["Undershot"], "StopIntegration"}]
          y, {x, x0, xMax}, {k, y0}]]

The solutions are obtained by calling sol[k,y0][x] and can be plotted by giving numbers for k and y0.

Example: Take $k=0.4$

Then, $y_0 = 6.0$ (blue) leads to an "overshoot", $y_0 = 3.0$ (red) leads to an "undershoot", and $y_0 = 4.51114$ (black) is the "sweetspot" which gives the solution to the boundary value problem.

enter image description here

Question: What is the best way (in terms of performance) to automate getting the solution to the boundary value problem in Mathematica?

I definitely know that the correct $y_0$ lies somewhere between $\frac{3-\sqrt{9-8k}}{2k} < y_0 < \frac{3+\sqrt{9-8k}}{2k}$.

Appendix: the following code generates the graphic above.

plot1 = Plot[sol[0.4, 6][x], {x, .01, sol[0.4, 6]["Domain"][[1, 2]]}, PlotStyle -> Blue];
plot2 = Plot[sol[0.4, 3][x], {x, .01, sol[0.4, 3]["Domain"][[1, 2]]}, PlotStyle -> Red];
plot3 = Plot[sol[0.4, 4.51114][x], {x, .01, sol[0.4, 4.51114]["Domain"][[1, 2]]}, PlotStyle -> {Black, Thick}];
Show[plot1, plot2, plot3, PlotRange -> {{0, 10}, {-1, 6}}]
  • $\begingroup$ @ QuantumDot, how do you write the (init.cond.1 and 2), that is, y[x0] and y'[x0] as a polynomial in your code? Thanks! $\endgroup$ – Enter Jan 8 '16 at 4:25
  • $\begingroup$ @can It's the small x behavior of the solution at x=x0. I can't start integration right at x=0 because that is a singular point. $\endgroup$ – QuantumDot Jan 8 '16 at 14:11
  • $\begingroup$ @ QuantumDot, I understand why you do this, but I cannot figure out how you do this. Could you give a hint, say, how $y(x_0)$ approximated as $y_0 + \frac{1}{8} x_0^2 (y_0 - \frac{3}{2} y_0^2 + \frac{k}{2} y_0^3)$. Sorry if it is a silly question. $\endgroup$ – Enter Jan 9 '16 at 3:22

A mundane but effective approach is to use the Shooting Method built into NDSolve.

x0 = .01; xMax = 11.5; k = .4;
c = N[D[(BesselK[1, x]/x), x]/(BesselK[1, x]/x) /. x -> xMax];
sol = NDSolveValue[{y''[x] + 3/x y'[x] - y[x] + 3/2 y[x]^2 - 
     k/2 y[x]^3 == 0, y'[xMax] == c y[xMax], y'[x0] == 0}, y, {x, x0, xMax}, 
  Method -> {"Shooting", "StartingInitialConditions" -> {y[x0] == 4.5, y'[x0] == 0}}]

Note that the boundary condition y'[xMax] == c y[xMax] is employed to match the exponentially decreasing asymptotic solution to the ODE, BesselK[1, x]/x. Also, xMax is increased to 11.5.

The k = .4 solution, compared to that given in the question is,

LogPlot[{sol[x], solt[x]}, {x, x0, xMax}, AxesLabel -> {x, y}]

enter image description here

Two features are evident. First, the solution already is exponentially decreasing at x = 3. Second, both approaches become unstable, the method in the Question at about x = 9.5 and that given here at about x = 11. This presumably occurs due to unavoidable round-off errors that couple into the second, exponentially growing, asymptotic solution to the ODE.

The value of y near the origin agrees quite closely with that in the Question.

(* 4.51095 *)

The next plot depicts y[x0] (points) as a function of k, along with the bounds on y[x0] (solid curves) given in the Question. (These bounds correspond to the values of y[x0] which satisfy the ODE for all x, apart from the boundary condition at infinity.)

data = {{0., 5.78018}, {0.1, 5.49104}, {0.2, 5.18459}, {0.3, 4.85871}, {0.4, 4.51095}, 
        {0.5, 4.13862}, {0.6, 3.73916}, {0.7, 3.31151}, {0.75, 3.08818}, 
        {0.8, 2.86076}, {0.85, 2.63322}, {0.88, 2.49923}, {0.9, 2.41187}};
Show[ListPlot[data, AxesLabel -> {"k", y[0]}, PlotStyle -> Black], 
  Plot[{(3 - Sqrt[9 - 8 k])/(2 k), (3 + Sqrt[9 - 8 k])/(2 k)}, {k, 0, 1},
  PlotStyle -> Black], PlotRange -> {{0, 1}, {0, 6}}]

enter image description here

Solutions for k < .8 are obtained with moderate ease by first running the code above with xMax = 5, which is fairly forgiving of inaccurate initial guesses for y[x0], and then using the resulting actual y[x0] with xMax = 10 to obtain y[x0] accurate to five significant figures or better. The process becomes rather more delicate at larger k, because y[0] lies very close to its upper bound. I typically would need to guess about five times (moving the point at which NDSolve declared the equations stiff toward xMax) and also increase xMax somewhat to obtain a solution that was both stable and accurate. Any particular calculation takes less than a second of computer time.

A comparison among y for k of 0(blue), .6(brown), .8(green), and .9 (red) follows. That the k = 0.9 solution is essentially constant out to about x = 10 explains why the calculation at larger k is so delicate and also why xMax = 15 is necessary for accuracy.

enter image description here

Note that I have not experimented with NDSolve Methods for handling stiffness, which may improve the robustness of this solution. Nonetheless, I was able to obtain the curves above without excessive difficulty.

  • $\begingroup$ I have mathematica 9, and I can't get your solution to work for any value of $k$ except for 0.4. I always get "stiff system" errors. $\endgroup$ – QuantumDot May 25 '15 at 19:43
  • $\begingroup$ @QuantumDot I have added information to my answer to assist you in obtaining solutions for various k and also corrected a minor (and numerically inconsequential) error in the upper boundary condition. $\endgroup$ – bbgodfrey May 27 '15 at 21:36
  • $\begingroup$ @bbgodfrey, a digression: I am interested in how to estimate or obtain the asymptotic solution to the ODE. If I use DSolve[3/x y'[x] - y[x] == 0, y[x], x], it gives (*y[x] -> E^(x^2/6) C[1]*) $\endgroup$ – Enter Jan 7 '16 at 10:25
  • $\begingroup$ @can Because we are seeking the asymptotic solution that goes to zero for large x, drop the nonlinear terms in the ODE and solve with DSolve. It returns {y[x] -> -((BesselJ[1, I x] C[1])/x) + (BesselY[1, -I x] C[2])/x}, and an appropriate choice of the two constants yields BesselK[1, x]/x, which satisfies the y[Infinity] == 0 boundary condition. $\endgroup$ – bbgodfrey Jan 7 '16 at 11:10
  • $\begingroup$ @bbgodfrey, thanks again for your explanation. You appear to know all the maths about numerical and analytical one! $\endgroup$ – Enter Jan 7 '16 at 15:39

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.