I'm trying to solve the following forth-order ODE with the shooting method:

$$\frac{1}{5}(y-2xy^\prime)=\frac{1}{x}\left\{\frac{xy^\prime}{y}+xy^3 \left[\frac{(xy^\prime)^\prime}{x} \right]^\prime \right\}^\prime$$ where $\prime$ denotes differentiation. It is clear that this ODE need 4 boundary conditions (BCs): two of which are given $y^\prime(0)=y^{\prime\prime\prime}(0)=0$, however, the other two values of $y(0)$ and $y^{\prime\prime}(0)$ are determined by shooting for BCs at infinity ($x_\text{Max}$).

It has been found that $y=Ax^{1/2}$ is a far-field asymptotic solution to the ODE by neglecting the nonlinear terms. Now, let me refer to the solutions with $y\sim x^{1/2}$ asymptotic behavior as the $x^{1/2}$ solutions. Here, I want to search just for this $x^{1/2}$ solutions by shooting method in which the ODE is integrated with a 4th-order R-K scheme from the origin to a certain end-value of $x_{max}$. In practice, I truncate the domain to $x_0<x<x_\text{Max}$ to avoid the singularity at $x=0$, and I would like to work with $x_0=10^{-4}$ and $x_\text{Max}=10$.

The shooting parameters $y(0)$ and $y^{\prime\prime}(0)$ will be adjusted so that $y-2xy^\prime=0$ or $y+4x^2y^{\prime\prime}=0$ at $x=x_\text{Max}$. The two BCs were chosen to be independent to involve low-order derivatives and to require $y\propto x^{1/2}$ at the end of the integration interval.

My questions are:

(1) How can I use a Taylor series (say, five-terms) to start the integration at $x_0=10^{-4}$. I am thinking the series should be located in the position of initial conditions, but I can't figure it out;

(2) How do I implement this searching method: firstly, find the solutions that satisfying the 1st BCs $y-2xy^\prime=0$ at $x_\text{Max}$, next, find the solutions that satisfying the 2nd BC $y+4x^2y^{\prime\prime}=0$ at $x_\text{Max}$. Then ParametricPlot the shooting parameters $(y(0), y^{\prime\prime}(0))$ that corresponding to solutions satisfying the above-mentioned BCs respectively, thus the intersections of the two sets of curves correspond to the desired $x^{1/2}$ solutions.

(3) Is it possible to use Do loop and FindRoot to find the parameters, as shown in this answer?



1 Answer 1


Interesting equation. It seems to be necessary to use the asymptotic solution as the boundary at $x=x_0$ if you want to solve the equation correctly.

Thanks to this answer, we can easily get the series solution at $x=0$ with seriesDSolve:

eq1 = 1/5 (y[x] - 2 x y'[x]) == D[(x y'[x])/y[x] + x y[x]^3 D[D[x y'[x], x]/x, x], x]/x; 
bc1 = y[x] - 2 x y'[x] == 0; 
bc2 = y[x] + 4 x^2 y''[x] == 0; 

solSeries = seriesDSolve[eq1, y, {x, 0, 5}, 
                         {y[0] -> a, y'[0] -> 0, y''[0] -> b, y'''[0] -> 0}]    

enter image description here

With the asymptotic solution, we obtain a list of more accurate boundaries at $x=x_0$:

newbclist = Thread[(Derivative[#][y][x0] == 
        (D[Normal@solSeries, {x, #}] /. x -> x0) & ) /@ Range[0, 3]]

Finally, solve the equation and plug it into the right boundary and plot:

x0 = 1/10^4; xMax = 10;
sol = ParametricNDSolveValue[{eq1, newbclist}, y, {x, x0, xMax}, {a, b}]

ContourPlot[{bc1, bc2} /. x -> xMax /. y -> sol[a, b] // Evaluate, 
            {a, 0.4, 0.8}, {b, -0.05, 0.25}]

enter image description here

By contrast, if one simply approximate the left boundary as $y'(x_0)=y'''(x_0)=0$, the resulting plot will be undesired:

bclist = Thread[(Derivative[#][y][x0] & /@ Range[0, 3]) == {a, 0, b, 0}]

solfake = ParametricNDSolveValue[{eq1, bclist}, y, {x, x0, xMax}, {a, b}]

ContourPlot[{bc1, bc2} /. x -> xMax /. y -> solfake[a, b] // Evaluate, 
            {a, 0.4, 0.8}, {b, -0.05, 0.25}]

enter image description here

  • $\begingroup$ @W.Robin This isn't the first time I met a equation that using the asymptotic solution as the boundary can make the solving easier. (But I never encountered one that can't be solved successfully without the asymptotic solution, actually I've struggled for more than half an hour before I noticed this!) $\endgroup$
    – xzczd
    Commented Jan 9, 2016 at 8:03
  • $\begingroup$ @ xzczd, 1. your answer employed NDSolve in which the 4 ICs have been estimated by power series expanding near singularity point in a consistent way. Thus, it essentially is a numerical power series method instead of the shooting method. Am I right? 2. How can I figure out the Method used by the NDSolve. I found that if I explicitly use Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 4}, there are some zigzag curves, so the default method is better. Thanks! $\endgroup$
    – W. Robin
    Commented Jan 9, 2016 at 10:41
  • $\begingroup$ @W.Robin 1. I think this method still belongs to shooting method, the power series is just to find a more accurate approximation of $y′(x_0)$ and $y‴(x_0)$. 2. This post may be helpful, BTW, if you simply set Method -> "ExplicitRungeKutta" the resulting plot will be as good as the default. $\endgroup$
    – xzczd
    Commented Jan 10, 2016 at 2:45
  • $\begingroup$ @W.Robin Check this post. $\endgroup$
    – xzczd
    Commented Jan 11, 2016 at 5:45
  • 2
    $\begingroup$ @W.Robin 1. FindRoots2D is essentially a function that uses points on the plot of equation as starting points to find the accurate roots of equation with FindRoot so it should be as effective as "FindRoot with a bunch of good starting points", if you still feel worried, just plug the roots back to the equation to check if they're correct . 2. Check this post. BTW, you may also give a try to other answers under the FindAllCrossings2D question. $\endgroup$
    – xzczd
    Commented Jan 11, 2016 at 15:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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