Describing a particular second-order equation in DifferentialEquations.jl

Specific Domains Numerics
,
#1

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

#2

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.

2 Likes
#3

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
#4

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.

#5

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

2 Likes
#6

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

3 Likes
#7

Thanks for that, John.

#8

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.

#9

Apologies for the delay in closing this out. Both the “tstops” and the “RSO65()” suggestions by klaff helped greatly with my problem. I will try to indicate them as solutions.

1 Like