I have the following ODE with complex $\omega$:

$\qquad f(r)\frac{d}{dr}(f(r)\frac{dR}{dr})+(\omega^2-f(r)(\frac{1}{r^2}+\mu^2))R(r)=0,$

and with $f(r)=1-2/r$. I want to solve it starting at $r=2$ outwards for a specific set of $\{\omega,\mu\}$. For that I use NDSolve and the code

f[r_] := 1 - 2/r (*M=1*)
Diff = ((f[r] D[f[r] D[#, r],r] + (ω^2 - f[r] (1/r^2+μ^2)) #)) &;
eq = ((Diff @ R[r]) //. {ω -> ωnr + I ωni, μ->0.01} // Simplify) == 0;
sol = 
    {eq, R[2.01] == 1 + I, R'[2.01] == 1 + I}, 
    R, {r, 2.01, 100}, {ωnr, ωni}, 
    MaxStepSize -> 0.001];
Rsol = R /. sol;
Shouldbezero = 
  Diff@(Rsol[0.5, 0.5][r]) //. {ω -> 0.5 + I 0.5, μ ->0.01} /. {M -> 1} // Simplify;
LogPlot[Shouldbezero // Abs, {r, 2.01, 50}]

As can be seen above, I plugged the solution that NDSolve provides back into the ODE, to check, whether the result is actually a solution. To my surprise, the numerical solution seems to be off, by quite a bit.

Plot of the absolute value of R(r) plugged back into the ODE:

Plot of the absolute value of R(r)

I have played around with AccuracyGoal, etc., but could only achieve minimal improvement. How come Mathematica gives me such a poor result? Or is this due to the nature of the equation?

  • $\begingroup$ I don't think it's a poor result. Do notice the order of magnitude of solution is 14 around r==50. $\endgroup$ – xzczd Mar 22 '19 at 4:19
  • 2
    $\begingroup$ Looks like an issue of relative vs. absolute error. $\endgroup$ – J. M.'s discontentment Mar 22 '19 at 12:06

Plotting the absolute value of the solution,

LogPlot[Abs@Rsol[0.5, 0.5][r], {r, 2.01, 50}]

enter image description here

shows that it too grows exponentially and is about eight orders of magnitude larger than the error curve. So, the solution, Rsol[0.5, 0.5][r] actually is quite accurate.

| improve this answer | |

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.