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

CosmoProf · June 15, 2020, 11:39pm

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:

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?
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.)
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!

ChrisRackauckas · June 15, 2020, 11:43pm

Try:

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

https://diffeq.sciml.ai/latest/solvers/ode_solve/

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.

CosmoProf · June 15, 2020, 11:45pm

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

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

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


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]),
    [tspanForward[2]])
#
t_fineBack = append!(append!([tspanBackward[1]],
    [(tspanForward[1] - (stepN*fineStepSize)) for stepN = 1:fineBackStepsWOends]),
    [tspanBackward[2]])
#
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…

CosmoProf · June 16, 2020, 12:04am

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.
Thanks.

CosmoProf · June 16, 2020, 11:17pm

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")

ChrisRackauckas · June 16, 2020, 11:37pm

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.

CosmoProf · June 17, 2020, 2:27am

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

Thanks for all the advice!