How to remove oscillations in a solution to Differential Equation having indeterminate form using DifferentialEquations.jl

Question

I am trying to solve differential equations with indeterminate forms and I am playing with Julia and DifferentialEquations.jl package to do so.

I started with a simple differential equation with indeterminate form that I know of (as a test case) to see if this might work ( My eventual goal is to solve a complex system of differential equations and couple it with Machine Learning. This requires substantial amount of effort on my part to analyze and code, which I would only like to begin once I know that some basic test cases work).

Here is the test case code I started with:

bondi(u,p,t) = -(2*k/t^2)*(1 - (2*a^2*t/k))/(1 - a^2/u)
u0 = 0.01
p = (1e4, 10)
k,a = p
tc = k/(2*a^2)
tspan = (tc/10, tc*5)
prob = ODEProblem(bondi,u0,tspan,p)

The analytical solution to this problem (known as the Bondi flow problem in Astrophysics) is well known (even for the indeterminate cases). I am interested in the indeterminate solutions obtained from the solver.

When I solve using sol = solve(prob) I am getting an oscillating solution quite different from the smooth analytical solution I am aware exists (see the figure in the link below).

Plot u vs t

I was expecting to encounter some 'issues' as t approaches 50 (while simultaneously the y axis variable (representing speed) denoted by u would approach 100) as only then the numerator (and denominator) would vanish together. Any ideas why the solutions starts to oscillate?

I also tried with sol = solve(prob, alg_hints = [:stiff]) and got the following warning:

Warning: Interrupted. Larger maxiters is needed. @ DiffEqBase C:\Users\User.julia\packages\DiffEqBase\1yTcS\src\integrator_interface.jl:329

The solution still oscillates (similar to the solutions obtained without imposing stiffness).

Am I doing something incorrect here? Is there another way to solve such indeterminate equations with the DifferentialEquations.jl package?

Answer 1

That ODE is highly numerically unstable at the singularity. Time stepping based methods are prone to explode when encountering that kind of problem. You need to use a collocation based method, like something in ApproxFun.jl, to avoid directly evaluating at the singularity in order to stabilize the solve.