# Solving ODEs arising from spherical and cylindrical systems?

**URL:** <https://discourse.julialang.org/t/solving-odes-arising-from-spherical-and-cylindrical-systems/23462>\
**Category:** General Usage\
**Tags:** diffeq\
**Created:** [April 24, 2019, 12:30am UTC](https://discourse.julialang.org/t/solving-odes-arising-from-spherical-and-cylindrical-systems/23462 "2019-04-24T00:30:31Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![ArjunNarayanan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/arjunnarayanan/32/6475_2.png) [@ArjunNarayanan](https://discourse.julialang.org/u/ArjunNarayanan)\
**Post date:** [April 24, 2019, 12:30am UTC](https://discourse.julialang.org/t/solving-odes-arising-from-spherical-and-cylindrical-systems/23462/1 "2019-04-24T00:30:31Z")

</div>

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:

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

```julia
┌ 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 `NaN`s.

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`?

---

<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:** [April 24, 2019, 2:38am UTC](https://discourse.julialang.org/t/solving-odes-arising-from-spherical-and-cylindrical-systems/23462/2 "2019-04-24T02:38:57Z")

</div>

> [@ArjunNarayanan](#):
>
> **Note** : For the function uu to be analytic at r=0r = 0 we can show by Taylor expansion that u′′(0)=0u’'(0) = 0. I use this condition in my regular Finite Difference solver. Perhaps I can introduce such a constraint somehow into the `ODEProblem` ?

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

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [April 24, 2019, 2:54am UTC](https://discourse.julialang.org/t/solving-odes-arising-from-spherical-and-cylindrical-systems/23462/3 "2019-04-24T02:54:54Z")

</div>

This is something that ApproxFun.jl might be good at. In fact, we have a WIP [implementation](https://github.com/JuliaApproximation/MultivariateOrthogonalPolynomials.jl/blob/master/src/Disk/Disk.jl) 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.

---

<div class="post-metadata">

**Author:** ![tomaklutfu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tomaklutfu/32/2411_2.png) [@tomaklutfu](https://discourse.julialang.org/u/tomaklutfu)\
**Post date:** [April 24, 2019, 5:49am UTC](https://discourse.julialang.org/t/solving-odes-arising-from-spherical-and-cylindrical-systems/23462/4 "2019-04-24T05:49:49Z")

</div>

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

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

```

---

<div class="post-metadata">

**Author:** ![BeastyBlacksmith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/beastyblacksmith/32/4741_2.png) [@BeastyBlacksmith](https://discourse.julialang.org/u/BeastyBlacksmith)\
**Post date:** [April 24, 2019, 9:48am UTC](https://discourse.julialang.org/t/solving-odes-arising-from-spherical-and-cylindrical-systems/23462/5 "2019-04-24T09:48:33Z")

</div>

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

```julia
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 ^^;

---

<div class="post-metadata">

**Author:** ![ArjunNarayanan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/arjunnarayanan/32/6475_2.png) [@ArjunNarayanan](https://discourse.julialang.org/u/ArjunNarayanan)\
**Post date:** [April 24, 2019, 5:41pm UTC](https://discourse.julialang.org/t/solving-odes-arising-from-spherical-and-cylindrical-systems/23462/6 "2019-04-24T17:41:22Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [April 24, 2019, 6:28pm UTC](https://discourse.julialang.org/t/solving-odes-arising-from-spherical-and-cylindrical-systems/23462/7 "2019-04-24T18:28:31Z")

</div>

`newton` is very experimental. This worked for me:

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

```

---

<div class="post-metadata">

**Author:** ![BeastyBlacksmith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/beastyblacksmith/32/4741_2.png) [@BeastyBlacksmith](https://discourse.julialang.org/u/BeastyBlacksmith)\
**Post date:** [April 25, 2019, 10:05am UTC](https://discourse.julialang.org/t/solving-odes-arising-from-spherical-and-cylindrical-systems/23462/8 "2019-04-25T10:05:14Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [April 25, 2019, 10:13am UTC](https://discourse.julialang.org/t/solving-odes-arising-from-spherical-and-cylindrical-systems/23462/9 "2019-04-25T10:13:30Z")

</div>

The [ApproxFun docs](http://juliaapproximation.github.io/ApproxFun.jl/stable/) 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.