TransWikia.com

NDSolve solution plugged back into ODE quickly diverges

Mathematica Asked on January 29, 2021

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 = 
  ParametricNDSolve[
    {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?

One Answer

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.

Answered by bbgodfrey on January 29, 2021

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP