# \[DifferentialEquations.jl\] Solving Berman's equation

**URL:** <https://discourse.julialang.org/t/differentialequations-jl-solving-bermans-equation/63223>\
**Category:** General Usage\
**Tags:** diffeq\
**Created:** [June 20, 2021, 11:48am UTC](https://discourse.julialang.org/t/differentialequations-jl-solving-bermans-equation/63223 "2021-06-20T11:48:30Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![bmxam](https://avatars.discourse-cdn.com/v4/letter/b/b38774/32.png) [@bmxam](https://discourse.julialang.org/u/bmxam)\
**Post date:** [June 20, 2021, 11:48am UTC](https://discourse.julialang.org/t/differentialequations-jl-solving-bermans-equation/63223/1 "2021-06-20T11:48:31Z")

</div>

Hi everyone,

I’m trying to solve the axisymetric Berman’s equation (see [here](https://en.wikipedia.org/wiki/Berman_flow)) with DifferentialEquations.jl. This is a fourth order BVP :

\left[\left( \frac{F}{r} \right) \left(\frac{F'}{r} \right)' - \left( \frac{F'}{r} \right)^2 \right]' - \frac{1}{R} \left[\frac{1}{r} \left( r \left( \frac{F'}{r}\right)' \right)' \right]' = 0,

with the associated boundary conditions:

\lim\_{r \rightarrow 0} \left(\frac{F'}{r} \right)' = \lim\_{r \rightarrow 0} \frac{F}{r} = F'(1) = 0,~ F(1) = 1.

R is a constant parameter corresponding to the flow Reynolds number. The solution is very close to F(r) = sin(\pi/2 r^2), especially when R \rightarrow +\infty.

I used to solve this equation using MATLAB bvp4c function. I am know trying DifferentialEquations.jl. I can recover the correct solution but as soon as I ask for a “fine” mesh, i.e \Delta r = 0.001, the computation time is very high (more than 10 minutes on my laptop) whereas it was 60s maximum for MATLAB bvp4c.

I think a clever initialization could help, but I didn’t find out how to specify the \sin(\pi/2r^2) solution as an initial solution.

Could anybody look at my implementation and give me some advice please? Especially:

- what should look like the init. vector?
- how to deal with the jacobian singularity in r=0?

Here is my code:

```julia
function f!(du,u,p,r)
    F0 = u[1]
    F1 = u[2]
    F2 = u[3]
    F3 = u[4]

    if (r < p.myeps)
        du .= [
            F1
            F2
            F3
            0.
        ]
    else
        du .= [
            F1
            F2
            F3
            1/r*( 2*F3-3*F2/r + 3*F1/(r*r) ) - p.Rey/r*( -F1*F2+F0*F3 -3*F0*F2/r + F1*F1/r +3*F0*F1/(r*r) )
        ]
    end
end

function bc!(residual, u, p, r)
    uleft = u[1]
    uright = u[end]
    residual .= [
        uright[1] - 1.
        uright[2]
        uleft[1]
        uleft[2]
    ]
end

rspan = (0.,1.)
p = (Rey = Rey, myeps = 1e-20)
u0 = [0., 0., 0., 0.] # I don't know how to use the init function
bvp = TwoPointBVProblem(f!, bc!, u0, rspan, p)
sol = solve(bvp, MIRK4(), dt=dt) # we need to use the MIRK4 solver for TwoPointBVProblem

```

---

<div class="post-metadata">

**Author:** ![Salmon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/salmon/32/22968_2.png) [@Salmon](https://discourse.julialang.org/u/Salmon)\
**Post date:** [June 20, 2021, 1:26pm UTC](https://discourse.julialang.org/t/differentialequations-jl-solving-bermans-equation/63223/2 "2021-06-20T13:26:07Z")

</div>

I am not an expert on these particular kinds of problem, so maybe there is a bigger issue but I do see some unnecessary allocations in your derivative funciton for example:

> [@bmxam](#):
>
> ```julia
> du .= [
> F1
> F2
> F3
> 0.
> ]
> 
> ```

Unless you have a particular reason not to do that, I would suggest you replace all these assignments with

```julia
        du .= F1,F2,F3,0.

```

which does not allocate anything. If I benchmark your f! function this way, I reduced the time by a factor of 2.

Taking this one step further, maybe it is worth to consider using StaticArrays since your state vector only has four elements. They are stack-allocated and are likely faster in your case as suggested by the DifferentialEquations.jl library:  
[https://diffeq.sciml.ai/stable/basics/faq/#My-ODE-is-solving-really-slow](https://diffeq.sciml.ai/stable/basics/faq/#My-ODE-is-solving-really-slow)  
[https://diffeq.sciml.ai/stable/tutorials/ode\_example/#ode\_other\_types](https://diffeq.sciml.ai/stable/tutorials/ode_example/#ode_other_types)

Ultimately, there may of course still be some issue that has to do with the BV problem in particular, I hope someone else who knows more about these chimes in in that case. I’d probably still at least try out the first option I’ve mentioned since its easily done.

---

<div class="post-metadata">

**Author:** ![bmxam](https://avatars.discourse-cdn.com/v4/letter/b/b38774/32.png) [@bmxam](https://discourse.julialang.org/u/bmxam)\
**Post date:** [June 20, 2021, 6:48pm UTC](https://discourse.julialang.org/t/differentialequations-jl-solving-bermans-equation/63223/3 "2021-06-20T18:48:22Z")

</div>

Thank you. I’m bit ashamed that I did not present a code with SA, I was only focused on the package and not and my function. I will try that.

However I still believe that initializing properly the problem would greatly improve the computation time.

---

<div class="post-metadata">

**Author:** ![bmxam](https://avatars.discourse-cdn.com/v4/letter/b/b38774/32.png) [@bmxam](https://discourse.julialang.org/u/bmxam)\
**Post date:** [June 20, 2021, 7:09pm UTC](https://discourse.julialang.org/t/differentialequations-jl-solving-bermans-equation/63223/4 "2021-06-20T19:09:33Z")

</div>

I have followed your advice. Here is the updated code:

```julia
function f!(du,u,p,r)
    F0 = u[1]
    F1 = u[2]
    F2 = u[3]
    F3 = u[4]

    if (r < p.myeps)
        du .= F1, F2, F3, 0.
    else
        du .= F1, F2, F3, 1/r*( 2*F3-3*F2/r + 3*F1/(r*r) ) - p.Rey/r*( -F1*F2+F0*F3 -3*F0*F2/r + F1*F1/r +3*F0*F1/(r*r) )
    end
end

function bc!(residual, u, p, r)
    uleft = u[1]
    uright = u[end]
    residual .= uright[1] - 1., uright[2], uleft[1], uleft[2]
end

rspan = (0.,1.)
p = (Rey = Rey, myeps = 1e-20)
u0 = [0., 0., 0., 0.] # I don't know how to use the init function
bvp = TwoPointBVProblem(f!, bc!, u0, rspan, p)
sol = solve(bvp, MIRK4(), dt=dt) 

```

Using BenchmarkTools, I get the following results for the `solve` with `dt=0.01` and `Rey=2000`.  
Code from my first post: 2.291 s (39462640 allocations: 4.03 GiB)  
Updated code : 2.229 s (35868878 allocations: 3.65 GiB)

So even if we’ve improved the performance of `f`, the `solve` performance is not significantly improved…

---

<div class="post-metadata">

**Author:** ![stillyslalom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stillyslalom/32/45687_2.png) [@stillyslalom](https://discourse.julialang.org/u/stillyslalom)\
**Post date:** [June 20, 2021, 11:52pm UTC](https://discourse.julialang.org/t/differentialequations-jl-solving-bermans-equation/63223/5 "2021-06-20T23:52:31Z")

</div>

Using [ProfileView.jl](https://github.com/timholy/ProfileView.jl), we can find out where all the runtime is coming from:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/9/0/901278f92722be6018073f6cf259081806af55a8.png)

All the runtime cost appears downstream of `FiniteDiff.finite_diff_jacobian!`, which is allocating tons of temporary arrays despite the availability of a cache for the Jacobian. @YingboMa might have some idea of what’s going wrong here.

---

<div class="post-metadata">

**Author:** ![YingboMa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yingboma/32/2181_2.png) [@YingboMa](https://discourse.julialang.org/u/YingboMa)\
**Post date:** [June 21, 2021, 4:29am UTC](https://discourse.julialang.org/t/differentialequations-jl-solving-bermans-equation/63223/6 "2021-06-21T04:29:34Z")

</div>

Yeah, unfortunately, BoundaryValueDiffEq.jl doesn’t use optimized approximate Jacobian evaluations yet.

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [June 21, 2021, 6:53am UTC](https://discourse.julialang.org/t/differentialequations-jl-solving-bermans-equation/63223/7 "2021-06-21T06:53:53Z")

</div>

Have you tried a shooting method instead?

---

<div class="post-metadata">

**Author:** ![bmxam](https://avatars.discourse-cdn.com/v4/letter/b/b38774/32.png) [@bmxam](https://discourse.julialang.org/u/bmxam)\
**Post date:** [June 21, 2021, 7:17am UTC](https://discourse.julialang.org/t/differentialequations-jl-solving-bermans-equation/63223/8 "2021-06-21T07:17:46Z")

</div>

Is the shooting ok for TwoPointBVP? What algorithm do you recommend in association to the `Shooting` object?

And what should look like the init vector? The `Shooting` should be efficient only if I can provide an initial guess.

---

<div class="post-metadata">

**Author:** ![bmxam](https://avatars.discourse-cdn.com/v4/letter/b/b38774/32.png) [@bmxam](https://discourse.julialang.org/u/bmxam)\
**Post date:** [June 21, 2021, 7:56am UTC](https://discourse.julialang.org/t/differentialequations-jl-solving-bermans-equation/63223/9 "2021-06-21T07:56:11Z")

</div>

Looking at the sources, it seems that `BVPSystem` doesn’t handle user-provided jacobian function. Am I right?

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [June 21, 2021, 10:01am UTC](https://discourse.julialang.org/t/differentialequations-jl-solving-bermans-equation/63223/10 "2021-06-21T10:01:33Z")

</div>

Also, is not `bv4pc` adaptive?

---

<div class="post-metadata">

**Author:** ![YingboMa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yingboma/32/2181_2.png) [@YingboMa](https://discourse.julialang.org/u/YingboMa)\
**Post date:** [June 21, 2021, 6:33pm UTC](https://discourse.julialang.org/t/differentialequations-jl-solving-bermans-equation/63223/11 "2021-06-21T18:33:25Z")

</div>

Yeah, that’s right. Contributions are welcome. I wrote the code when I was in high school, so it’s very limited comparing to other software in SciML.