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

**URL:** <https://discourse.julialang.org/t/best-solver-for-bessel-like-differential-equations-with-singularity-unexpected-interpolation-behavior/41482>\
**Category:** Numerics\
**Created:** [June 15, 2020, 11:39pm UTC](https://discourse.julialang.org/t/best-solver-for-bessel-like-differential-equations-with-singularity-unexpected-interpolation-behavior/41482 "2020-06-15T23:39:00Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![CosmoProf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cosmoprof/32/15516_2.png) [@CosmoProf](https://discourse.julialang.org/u/CosmoProf)\
**Post date:** [June 15, 2020, 11:39pm UTC](https://discourse.julialang.org/t/best-solver-for-bessel-like-differential-equations-with-singularity-unexpected-interpolation-behavior/41482/1 "2020-06-15T23:39:00Z")

</div>

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!

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 15, 2020, 11:43pm UTC](https://discourse.julialang.org/t/best-solver-for-bessel-like-differential-equations-with-singularity-unexpected-interpolation-behavior/41482/2 "2020-06-15T23:43:22Z")

</div>

> [@CosmoProf](#):
>
> oscillatory in general, with a singularity in the function and its derivatives at `t=0` .

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/](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.

---

<div class="post-metadata">

**Author:** ![CosmoProf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cosmoprof/32/15516_2.png) [@CosmoProf](https://discourse.julialang.org/u/CosmoProf)\
**Post date:** [June 15, 2020, 11:45pm UTC](https://discourse.julialang.org/t/best-solver-for-bessel-like-differential-equations-with-singularity-unexpected-interpolation-behavior/41482/3 "2020-06-15T23:45:57Z")

</div>

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

```julia
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…

 ![InterpResids_Plot1](https://global.discourse-cdn.com/julialang/original/3X/3/8/383919d0feb3cd0fe13d9b8359bff2be3f9f0d23.jpeg)

---

<div class="post-metadata">

**Author:** ![CosmoProf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cosmoprof/32/15516_2.png) [@CosmoProf](https://discourse.julialang.org/u/CosmoProf)\
**Post date:** [June 16, 2020, 12:04am UTC](https://discourse.julialang.org/t/best-solver-for-bessel-like-differential-equations-with-singularity-unexpected-interpolation-behavior/41482/4 "2020-06-16T00:04:42Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![CosmoProf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cosmoprof/32/15516_2.png) [@CosmoProf](https://discourse.julialang.org/u/CosmoProf)\
**Post date:** [June 16, 2020, 11:17pm UTC](https://discourse.julialang.org/t/best-solver-for-bessel-like-differential-equations-with-singularity-unexpected-interpolation-behavior/41482/5 "2020-06-16T23:17:25Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 16, 2020, 11:37pm UTC](https://discourse.julialang.org/t/best-solver-for-bessel-like-differential-equations-with-singularity-unexpected-interpolation-behavior/41482/6 "2020-06-16T23:37:38Z")

</div>

> [@CosmoProf](#):
>
> 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` .

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.

---

<div class="post-metadata">

**Author:** ![CosmoProf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cosmoprof/32/15516_2.png) [@CosmoProf](https://discourse.julialang.org/u/CosmoProf)\
**Post date:** [June 17, 2020, 2:27am UTC](https://discourse.julialang.org/t/best-solver-for-bessel-like-differential-equations-with-singularity-unexpected-interpolation-behavior/41482/7 "2020-06-17T02:27:20Z")

</div>

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

Thanks for all the advice!