Numerical solution of Bessel-like equation using NDSolve

Question

I need to calculate solution of Bessel-like equation having general form: $\frac{d^2F}{dr^2}+\frac{1}{r}\frac{dF}{dr}+Q(r)F(r)=0$. Problems come from the points near $r=0$ leading to numeric errors.

For example, simple Bessel equation can be solved using DSolve:

DSolve[{y''[x] + 1/x y'[x] + y[x] == 0, y[0] == 1, y'[0] == 0}, y[x], x]

with output:

{{y[x] -> BesselJ[0, x]}}

However same equation using NDSolve produces errors:

sol = NDSolve[{y''[x] + 1/x y'[x] + y[x] == 0, y[0] == 1, y'[0] == 0}, y[x], {x, 0, 10}]

Power::infy: Infinite expression 1/0. encountered. >>
Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered. >>
NDSolve::ndnum: Encountered non-numerical value for a derivative at x == 0.`. >>

The question is: how to overcome erros in NDSolve routine and get numeric solution?

Answer 1

There is a problem with the initial conditions. At $x=0$ there is no solution. The analytical solution given by Mathematica seems to have used the standard solution for standard Bessel ODE http://mathworld.wolfram.com/BesselDifferentialEquation.html, which is a Bessel function. But I am not sure for $n=0$ one can just use Bessel(0,x), I think BesselY(0,x) is also needed. This is what I get.

But to answer you, I think you need to use series solution for this. Or avoid $x=0$ as was suggested in the comment for a numerical solution.

\begin{align*} y^{\prime\prime}\left( x\right) +\frac{1}{x}y^{\prime}\left( x\right) +y\left( x\right) & =0\\ y\left( 0\right) & =1\\ y^{\prime}\left( 0\right) & =0 \end{align*}

Let $x=e^{z}$ or $\ln\left( x\right) =z$, hence $\frac{dy}{dx}=\frac{dy} {dz}\frac{dz}{dx}=\frac{dy}{dz}\frac{1}{x}$ and $\frac{d^{2}y}{dx^{2}} =\frac{d^{2}y}{dz^{2}}\frac{dz}{dx}\frac{1}{x}+\frac{dy}{dz}\left( \frac {-1}{x^{2}}\right) =\frac{d^{2}y}{dz^{2}}\frac{1}{x^{2}}-\frac{1}{x^{2}} \frac{dy}{dz}$

Substituting all these back into the original ODE gives

\begin{align*} \frac{d^{2}y}{dz^{2}}\frac{1}{x^{2}}-\frac{1}{x^{2}}\frac{dy}{dz}+\frac{1} {x}\frac{dy}{dz}\frac{1}{x}+y\left( z\right) & =0\\ \frac{d^{2}y}{dz^{2}}+x^{2}y\left( z\right) & =0\\ \frac{d^{2}y}{dz^{2}}+e^{2z}y\left( z\right) & =0 \end{align*}

The solution to the above is

$$ y\left( z\right) = BesselJ_{0}\left( \sqrt{e^{2z}}\right) c_{1}+2 BesselY_{0}\left( \sqrt{e^{2z}}\right) c_{2} $$

Replacing back

\begin{equation} y\left( x\right) = BesselJ_{0}\left( x\right) c_{1} +2 BesselY_{0}\left( x\right) c_{2}\tag{1} \end{equation}

You see that there is a BesselY[0,x] function in the solution which do not show up the solution given by Mathematica for some reason.

Now to find the constants $c_{1},c_{2}$ we use initial conditions. At $x=0,y=1\,\ $,hence

$$ 1=c_{1}-\infty $$

For the derivative at $x=0$,

$$ y^{\prime}\left( x\right) =- BesselJ_{1}\left( x\right) c_{1}-2 BesselY_{1}\left( x\right) c_{2} $$

At $x=0$

$$ 0=0-\infty $$

sol = DSolve[{y''[z] + Exp[2 z] y[z] == 0}, y[z], z]

Btw, Maple 17 gives a numerical solution for this, without the 1/0 issue. But I am not sure how it avoided the point $x=0$ now. It uses RK45 standard numerical method. Here it is:

Answer 2

This can be done with "EquationSimplification" -> "Residual", provided the denominators are cleared from the equation.

ndsol = NDSolve[{x y''[x] + y'[x] + x y[x] == 0,
   y[0] == 1, y'[0] == 0}, y[x], {x, 0, 10},
  Method -> {"EquationSimplification" -> "Residual"}];

Check with DSolve:

dsol = DSolve[{x y''[x] + y'[x] + x y[x] == 0, y[0] == 1, y'[0] == 0},
    y[x], {x, 0, 10}];

Plot[y[x] /. Join[dsol, ndsol] // Evaluate, {x, 0, 10}, 
 PlotStyle -> {AbsoluteThickness[5], Automatic}]

Answer 3

In this case MMa won't do it directly, but we can use L'Hospital's rule to find the value of the undefined term at x = 0.

The undefined term is y'[x]]/x and taking the derivative of numerator and denominator gives us

y''[x]/1

and the PDE at x = 0 becomes

y''[0] + y''[0]/1 + y[0] == 0 /. y[0] -> 1

Solve[%, y''[0]] // Flatten
(*{(y''[0] -> -(1/2)}*)

which tells us that the limit of y'[x]/x as x -> 0 is -1/2, so we can form a Piecewise expression for the PDE.

pde = y''[x] + Piecewise[{{-1/2, x == 0}, {y'[x]/x, x > 0}}] + y[x] == 0

Now NDSolve has no problem at x = 0

sol = NDSolve[{pde, y[0] == 1, y'[0] == 0}, y[x], {x, 0, 10}] // Flatten;

And we get a match to the J0 solution.

Plot[{BesselJ[0, x], y[x] /. sol}, {x, 0, 10}]