I am trying to solve the equation below with DSolve. The equation is that of a wave, expected to fall off exponentially as r approaches infinity. The solution is a combination of Spherical Bessel functions, which is great! However, I am unable to obtain a plot of the solution after multiple tries with various ranges of r. Any suggestions will be appreciated.

sol = DSolve[{-w^2 f[r] + f''[r] + 2/r f'[r] - 2/r^2 f[r] == 0, 
        f[10^5] == 0, f'[10^5] == 1}, f, {r, 10^3, 10^5}]
test = Plot[{Evaluate[f[r] /. sol /. w -> 0.05 ]}, {r, 1*10^3, 
       1*10^5}, PlotRange -> {-10, 10}]
  • $\begingroup$ The solution appears to be zero to high order for the parameters chosen. $\endgroup$ – bbgodfrey Jan 9 at 20:13
  • $\begingroup$ Try a scaling regarding $r$ $\endgroup$ – Cesareo Jan 9 at 21:15
  • 1
    $\begingroup$ You can see the exponential behavior by using FunctionExpand, which rewrites the spherical Bessel functions: DSolve[{-w^2 f[r] + f''[r] + 2/r f'[r] - 2/r^2 f[r] == 0, f[10^5] == 0, f'[10^5] == 1}, f[r], r] // FunctionExpand // TrigToExp // FullSimplify. More generally, DSolve[-w^2 f[r] + f''[r] + 2/r f'[r] - 2/r^2 f[r] == 0, f[r], r] // FunctionExpand gives you the solution for arbitrary boundary conditions. $\endgroup$ – Roman Jan 9 at 22:31
  • $\begingroup$ @Roman this certainly helps me see the solution in a form that makes more physical sense. However, the issue of not being able to plot it still remains, even for arbitrary boundary conditions. $\endgroup$ – mv1996 Jan 10 at 18:03

The numerical alternative shows a plot

F = ParametricNDSolveValue[{-w^2 f[r] + f''[r] + 2/r f'[r] -2/r^2 f[r] == 0, f[10^5] == 0, f'[10^5] == 1},f, {r, 10^3, 10^5}, w]
Plot[{Evaluate[F[0.05][r] ]}, {r, 1*10^3, 1*10^5},PlotRange -> {-10, 10} ] 

enter image description here

The line has only real points near the right boundary.

  • $\begingroup$ I used {r, 10^5, 10^35} in the above equation and obtained a plot. Thank you! Interestingly, the same range in DSolve did not give me a plot. $\endgroup$ – mv1996 Jan 15 at 10:16

The function values are very large and a bit tricky to use. They are so large that plotting the function itself seems to overwhelm Mathematica's plot autoscaling. The following works to plot the logarithm of the function, to get an idea of the size of the numbers:

F[r_] = f[r] /. sol[[1]] /. w -> 1/20 // FullSimplify
ListLinePlot[Table[{r, Log[Abs[F[r]]]}, {r, 10^3, 10^5, 10^3}]]

If you just evaluate the function numerically as N[F[1000]] // Re, you get zero (this is wrong). If you use high precision as N[F[1000], 10^3] // Re you get more meaningful answers.


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.