Solving ODEs arising from spherical and cylindrical systems?

General Usage
1

I want to know what is the correct way of solving an ODE of the following form using DifferentialEquations.

Consider for example the ODE satisfied by a function u(r) (prime indicates a derivative)

u'' + \frac{u'}{r} - \frac{u}{r^2} = 0 \qquad r \in (0,1) \\ u(0) = 0 \\ u'(0) = g

Such ODEs often arise when solving problems under the assumption of spherical or cylindrical symmetry. Solutions are found by assuming a power series form for u. We can verify that the function u(r) = gr satisfies the ODE as well as the initial conditions.

A typical way of solving such a second order ODE numerically is to re-write it as a system of two first order ODEs by introducing a variable v = u'. The system is then:

v' = - \frac{v}{r} + \frac{u}{r^2} \\ u' = v \\ v(0) = g \\ u(0) = 0

Using DifferentialEquations I attempt to solve the above problem as follows:

using DifferentialEquations

function g(du, u, p, t)
    du[1] = u[2]
    du[2] = -u[2]/t + u[1]/t^2
end

u0 = [0.0, -1.0]
tspan = (0.0, 1.0)
prob = ODEProblem(g, u0, tspan)
sol = solve(prob)

This produces the following warnings,

┌ Warning: First function call produced NaNs. Exiting.
└ @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/uJP6R/src/initdt.jl:74
┌ Warning: Automatic dt set the starting dt as NaN, causing instability.
└ @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/uJP6R/src/solve.jl:376
┌ Warning: Instability detected. Aborting
└ @ DiffEqBase ~/.julia/packages/DiffEqBase/EFqMn/src/integrator_interface.jl:162

I understand that this is probably because the differential equation cannot be evaluated at r = 0 as this will produce NaNs.

I can offset tspan to start a little bit away from r = 0 and produce a result. My question, however, is how does one go about solving such an ODE using DifferentialEquations?

Note: For the function u to be analytic at r = 0 we can show by Taylor expansion that u''(0) = 0. I use this condition in my regular Finite Difference solver. Perhaps I can introduce such a constraint somehow into the ODEProblem?

2

Yeah this needs to use some form of collocation since using the derivative at 0 will cause the instability.

3 Likes
3

This is something that ApproxFun.jl might be good at. In fact, we have a WIP implementation of Zernike polynomials which correctly capture the smoothness at the origin of the disk, though for your needs the bog standard Chebyshev method in ApproxFun might suffice.

2 Likes
4

Why not check t==0 and branch to du[2]=0, may be like

du[2] = t==0 ? 0 : (-u[2]/t + u[1]/t^2)
1 Like
5

For your particular problem the solution of @tomaklutfu works quite fine. Integrating away from these kind of “singularities” works in most cases, while it will fail in most cases, when you integrate into the “singularity”, e.g. if you had tspan = (1.0, 0.0).
In that case I haven’t found a better solution than expand the solution around the singularity and integrate in tspan = ( 1.0, e ), e > 0. And pad the expanded solution at the end.

@dlfivefifty: Am I doing something wrong or is it just taking very long ( > 10 minutes )?

using ApproxFun
x = Fun( identity, 0..1 )
u0 = one(x)
N = u -> [ u(0), u'(0)-10, u''+u'/x-u/x^2 ]
u = newton( N, u0 )

I have to admit, that this approach is like black magic to me ^^;

6

This is simple, perfect for my case, and one of those painfully obvious things that just don’t strike you. Thanks!

7

newton is very experimental. This worked for me:

julia> x = Fun(0..1);

julia> D = Derivative();

julia> u = [ivp(); x^2*D^2 + x * D - 1] \ [0,10,0]
Fun(Chebyshev(0..1),[5.0, 5.0, -0.0, -0.0, -0.0, -0.0, -0.0])

julia> u(0.1) ≈ 10*0.1
true
8

That looks interesting. Could you give a little Explanation what is happening here? This approximation theory is quite new for me.

9

The ApproxFun docs are a good place to start. Basically we represent the function in Chebyshev series on 0..1 (which are cosine series in disguise under the map u(2cos(θ)-1)). The differential operator corresponds to a sparse (banded) matrix. \ adaptively increases the resolution until convergence is detected, in this case it only requires 2 coefficients.

1 Like