Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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?

share|improve this question
Is this the whole story? I get errors when I run your DSolve, though I also get the answer you show. – bill s Oct 20 '13 at 17:09
You can substitute the 0 making trouble with a small number : sol[x_] = First[y[x] /. NDSolve[{y''[x] + 1/x y'[x] + y[x] == 0, y[$MachineEpsilon] == 1, y'[$MachineEpsilon] == 0}, y[x], {x, $MachineEpsilon, 10}]]. – b.gatessucks Oct 20 '13 at 17:18

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

Mathematica graphics

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:

Mathematica graphics

share|improve this answer
Thanks so much! Brilliant idea to make change of variables, it helped a lot. – denkorw Oct 21 '13 at 12:41

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

Mathematica graphics

share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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