Describing a particular second-order equation in Differentialequations.jl

I am using DifferentialEquations.jl . If I wanted to solve a simple second order equation of the form:
a t^2 f’’(t) + b t f’(t) +c = 0 (for a,b,c real constants)
I think the way to describe the problem in differentialequations.jl is to use the following function and ODEProblem (for appropriate values of p):

function eqns(du,u,p,t)
du[1] = u[2]
du[2] = -p[1]*u[2]/t-p[2]*u[1]/t^2.0
end

prob=ODEProblem(eqns,u0,tspan,p)

Clearly, this involves re-arranging the formula for f’’(t), including dividing across by t^2. This causes a difficulty with solving the equation, since, for example if we want to start the solution at t=0, the right hand side is undefined at that point.

Am I doing something wrong or is there another way to describe the problem to avoid this issue? Thanks.

(The full example code is below for reference.)

using DifferentialEquations
using Plots
starttime=0.1
endtime=1.1
steptime=0.1
p=[200.0,100.0] #The parameters
function eqns(du,u,p,t)
du[1] = u[2]
du[2] = -p[1]*u[2]/t-p[2]*u[1]/t^2.0
end
u0=[-5, -5] # The initial conditions
tspan=(starttime,endtime)
prob=ODEProblem(eqns,u0,tspan,p)
sol=solve(prob)
display(plot(sol))

The equation is singular at t=0. Solutions for this equation will exist on intervals that do not include t=0.

It’s a Cauchy-Euler equation, which can be solved in closed form by assuming f(t) = t^m and deriving an equation for m.

1 Like

Thanks very much for looking at this, John. I am just posting this simple equation to describe the problem minimally. The actual equations I have in mind to solve are considerably more complicated and, unfortunately, do not have closed-form solutions.

1 Like

You’re welcome. If your real problem has the same kind of singularity (coefficients of t^n in front of the nth derivative of f), you should not expect to be able to produce a solution starting at t=0, and you should be wary of trying to produce a numerical solution of a problem for which a solution does not exist.

You may be able to use the tstops argument to force the solver to hit points right around the undefined point (look for singularity on this page).

Also, see Best Solver for Bessel-like Differential Equations with Singularity; unexpected Interpolation Behavior - #2 by ChrisRackauckas

1 Like

Thanks for that, John.

Thanks very much for your two suggestions, klaff. The second, especially, I believe could be particularly relevant to my issues…but I’ll have to look closer at both ideas. Thanks.