Best Solver for Bessel-like Differential Equations with Singularity; unexpected Interpolation Behavior

Hi, I’m starting a project solving equations (using the DifferentialEquations.jl package) that will be similar to Bessel equations: oscillatory in general, with a singularity in the function and its derivatives at t=0. (I will ultimately extend the work to modified types of Bessel eq’ns, entangled systems of Inhomogeneous Diff Eq’s, and possibly 4th-order Diff Eq’s via order-reduction.)

Some questions:

  1. What might be the best solver(s) and parameters for me to use with SecondOrderODE, etc., to ensure the best possible accuracy and precision – especially letting me get high accuracy as close as possible to t -> 0. (I am willing to sacrifice computing speed & time.) Does anyone have specific knowledge and experience with this?

  2. When doing plots of residuals (vs. theory) of the discrete values of sol.u, plotted along with interpolations of sol(time) to see how the interpolated solution fits the discrete solution values, the interpolation results are super-weird, with strange oscillatory behavior to much higher error values. (Especially when sampling the interpolations too finely.) Does anyone know why this might be happening? (I might have to rely on those interpolations at some point.)

  3. Is it worth it going to BigFloat and BigInt everywhere, and possibly do setprecision() very tightly?

In my next post, I’ll include a portion of my code so you can see what I mean, and (if possible) a couple pictures of those weird interpolation plots.

Thanks in advance, everyone!


RKO65 - Olver’s 6 stage 5th order method. This method is robust on problems which have a singularity at t=0 .

If that doesn’t work for your problem, you might want to try ApproxFun. Most of the methods in DiffEq are not safe with a singularity at t=0: generally that might need a discretization that’s fully implicit over time.


Hi again, some of my relevant code and images, to illustrate my question above:

BesselJInitAmount = 1.75
BesselYInitAmount = 0.37

FnOrder = [1.27]

tspanForward = (5.0, 20.0)
tspanBackward = (5.0, 0.01)

function OrdBesselDiffEqToSolve!(ddu,du,u,p,t)
    OrderParam = p[1]
    ddu[1] = -(du[1]/t) - (1 - ((OrderParam/t)^2))*u[1]

function OrdBesselInitCondCalc(BesJamt, BesYamt, Ord, tVal)
    StartVal = (BesJamt*besselj(Ord,tVal)) + (BesYamt*bessely(Ord,tVal))
    StartDeriv = ((BesJamt/2)*(besselj((Ord-1),tVal) - besselj((Ord+1),tVal))) +
                 ((BesYamt/2)*(bessely((Ord-1),tVal) - bessely((Ord+1),tVal)))
    return [StartVal, StartDeriv]

function OrdBesselTheorySoln(t, BesJamt, BesYamt, Ord)::BigFloat
    return ((BesJamt*besselj(Ord,t)) + (BesYamt*bessely(Ord,t)))

InitialConditions = OrdBesselInitCondCalc(BesselJInitAmount, BesselYInitAmount,
    FnOrder[1], tspanForward[1])

xInit = [InitialConditions[1]]
vInit = [InitialConditions[2]]

probFor = SecondOrderODEProblem(OrdBesselDiffEqToSolve!,vInit,xInit,tspanForward,FnOrder)
probBack = SecondOrderODEProblem(OrdBesselDiffEqToSolve!,vInit,xInit,tspanBackward,FnOrder)
#@time solFor = solve(probFor, dense=true, reltol=1e-12, abstol=1e-12)
#@time solBack = solve(probBack, dense=true, reltol=1e-12, abstol=1e-12)
#@time solFor = solve(probFor, Tsit5(), dense=true, reltol=1e-12, abstol=1e-12)
#@time solBack = solve(probBack, Tsit5(), dense=true, reltol=1e-12, abstol=1e-12)
@time solFor = solve(probFor, RKO65(), dense=true, reltol=1e-12, abstol=1e-12, dt=0.001)
@time solBack = solve(probBack, RKO65(), dense=true, reltol=1e-12, abstol=1e-12, dt=-0.001)

TvalsFor = getindex.(solFor.t, 1)
FprimeOfTFor = getindex.(solFor.u, 1)
FofTFor = getindex.(solFor.u, 2)
TvalsBack = getindex.(solBack.t, 1)
FprimeOfTBack = getindex.(solBack.u, 1)
FofTBack = getindex.(solBack.u, 2)

u_analyticalFor = @. OrdBesselTheorySoln(TvalsFor,
                        BesselJInitAmount, BesselYInitAmount, FnOrder[1])
errFor = abs.(u_analyticalFor - FofTFor)
@. errFor[errFor == 0 ] = NaN # set zeros to NaN to avoid problems caused by log(0) = -Inf
u_analyticalBack = @. OrdBesselTheorySoln(TvalsBack,
                        BesselJInitAmount, BesselYInitAmount, FnOrder[1])
errBack = abs.(u_analyticalBack - FofTBack)
@. errBack[errBack == 0 ] = NaN # set zeros to NaN to avoid problems caused by log(0) = -Inf
fineStepSize = 0.00299
fineForStepsWOends = BigInt(((tspanForward[2]-tspanForward[1]) ÷ fineStepSize)-1)
fineBackStepsWOends = BigInt(((tspanBackward[1]-tspanBackward[2]) ÷ fineStepSize)-1)
t_fineFor = append!(append!([tspanForward[1]],
    [(tspanForward[1] + (stepN*fineStepSize)) for stepN = 1:fineForStepsWOends]),
t_fineBack = append!(append!([tspanBackward[1]],
    [(tspanForward[1] - (stepN*fineStepSize)) for stepN = 1:fineBackStepsWOends]),
fineForResids = [abs((solFor(time)[2] -
    OrdBesselTheorySoln(time, BesselJInitAmount, BesselYInitAmount,
    FnOrder[1]))) for time in t_fineFor]; @. fineForResids[fineForResids == 0] = NaN
fineBackResids = [abs((solBack(time)[2] -
    OrdBesselTheorySoln(time, BesselJInitAmount, BesselYInitAmount,
    FnOrder[1]))) for time in t_fineBack]; @. fineBackResids[fineBackResids == 0] = NaN
InterpDataPtResidPlots =
    scatter(TvalsFor, errFor, <plot options>); 
    plot!(t_fineFor, fineForResids, <plot options>); 
    scatter!(TvalsBack, errBack, <plot options>); 
    plot!(t_fineBack, fineBackResids, <plot options>)

Hopefully this code snippet (and the attached image) clarifies my question…

Hi Chris, thanks for the quick reply!

I did see that notation in the documentation, and in my code in my follow-up post you can see that I used that solver. But RKO65 is a fixed-timestep method, so I was wary about it. But if you’re sure that’s the best one for t=0 singularities, I could stick with it.

What you said about “a discretization that’s fully implicit over time” is interesting, could you explain what you mean? Also I haven’t heard of ApproxFun, so I guess I’ll have to look into that.

By the way, I have found that (fixed-timestep) method RK065 works about a billion (10^9) times worse for this problem than Tsit5 at t = 0.001.

I am currently also trying solvers AutoVern9(Rodas4()), AutoVern9(Rodas5()), AutoVern9(KenCarp4()), and AutoTsit5(Rosenbrock23()) with SecondOrderODE, but this is creating error messages (as I have raised in another thread):

DimensionMismatch("parent has 4 elements, which is incompatible with length 2")

That’s odd. But yeah… things with singularities at t=0 is a pretty special problem and it’s not something most of these methods are robust for, so really check out ApproxFun.jl with this.

Yeah, I’m trying out everything, looking for the best results.

Thanks for all the advice!