# Fast solution of a two-point boundary value problem with Shooting

**URL:** <https://discourse.julialang.org/t/fast-solution-of-a-two-point-boundary-value-problem-with-shooting/110531>\
**Category:** Performance\
**Tags:** diffeq, differentialequation\
**Created:** [February 21, 2024, 5:23pm UTC](https://discourse.julialang.org/t/fast-solution-of-a-two-point-boundary-value-problem-with-shooting/110531 "2024-02-21T17:23:26Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![flibe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/flibe/32/8382_2.png) [@flibe](https://discourse.julialang.org/u/flibe)\
**Post date:** [February 21, 2024, 5:23pm UTC](https://discourse.julialang.org/t/fast-solution-of-a-two-point-boundary-value-problem-with-shooting/110531/1 "2024-02-21T17:23:26Z")

</div>

Dear All,  
I am trying to find a way to solve as fast as possible a two point boundary value problem (TPBVP) related to the optimal control of a storage system (it arises from the associated optimality necessary conditions - Pontryagin’s minimum principle).

The differential equation is (state and costate equations - the state is the level of energy in the storage):

```julia
function area_ODE!(du, u, p, t)
    a, _ = p
    @inbounds @fastmath du[1] = min( a[3](t)::Float64 , max( a[4](t) , -1/a[2]::Float64*u[2] - a[5](t)::Float64 ) )
    @inbounds @fastmath du[2] = -a[1]::Float64*u[1]
    nothing
end

```

`a` contains parameters and functions of time.

The residual function is:

```julia
function bc2!(residual, u, p, t)
    _, b = p
    @inbounds @fastmath residual[1] = u[1][1] - b[1]::Float64
    @inbounds @fastmath residual[2] = u[end][2] - b[2]::Float64*u[end][1]
end

```

The first residual is the initial state condition (initial level of energy equal to zero), while the second one is the transversality condition, defined at at the final time.

I define the TPBVP as:

```julia
bvp2 = TwoPointBVProblem(area_ODE!, bc2!, u0, tspan, p)

```

Now, when I solve it with command:

```julia
sol = @time solve(bvp2, GeneralMIRK4(), dt = 0.01)

```

The solution takes very long, but seems correct (residuals very close to zero and control is as expected):

```julia
julia> include("solve_problem.jl")
┌ Warning: Backwards compatability support of the new return codes to Symbols will be deprecated with the Julia v1.9 release. Please see https://docs.sciml.ai/SciMLBase/stable/interfaces/Solutions/#retcodes for more information
└ @ SciMLBase ~/.julia/packages/SciMLBase/szsYq/src/retcodes.jl:354
1023.778347 seconds (24.47 G allocations: 1.718 TiB, 7.64% gc time, 0.42% compilation time: <1% of which was recompilation)

```

If instead I use Shooting:

```julia
sol = @time solve(bvp2, Shooting(AutoTsit5(Rosenbrock23())), dt = 0.000001, abstol = 1e-12, reltol = 1e-12)

```

a solution is found very fast:

```julia
julia> include("solve_problem.jl")
┌ Warning: Backwards compatability support of the new return codes to Symbols will be deprecated with the Julia v1.9 release. Please see https://docs.sciml.ai/SciMLBase/stable/interfaces/Solutions/#retcodes for more information
└ @ SciMLBase ~/.julia/packages/SciMLBase/szsYq/src/retcodes.jl:354
  2.176208 seconds (26.81 M allocations: 2.303 GiB, 7.71% gc time)

```

but the final (second) residual is very far from zero.

**My question is, is there a way to get Shooting return fast a correct solution?** Or to speed up `solve(bvp2, GeneralMIRK4(), dt = 0.01)`.

I have seen here [https://docs.sciml.ai/DiffEqDocs/dev/solvers/bvp\_solve/#Recommended-Methods](https://docs.sciml.ai/DiffEqDocs/dev/solvers/bvp_solve/#Recommended-Methods) that `MultipleShooting(5, FBDF())` could be better than single shooting, but when I try:

```julia
sol = @time solve(bvp2, MultipleShooting(5, FBDF()), dt = 0.000001, abstol = 1e-12, reltol = 1e-12)

```

I get:

```julia
julia> include("solve_problem.jl")
ERROR: LoadError: UndefVarError: `MultipleShooting` not defined

```

Finally, my ODE seems to be a DynamicalODE. I am wandering if it is necessary or not to write it as explained in [https://docs.sciml.ai/DiffEqDocs/dev/solvers/dynamical\_solve/#Dynamical,-Hamiltonian,-and-2nd-Order-ODE-Solvers](https://docs.sciml.ai/DiffEqDocs/dev/solvers/dynamical_solve/#Dynamical,-Hamiltonian,-and-2nd-Order-ODE-Solvers).

Any help would be appreciated, sorry for the long post.  
Thank you!  
Best  
Francesco

---

<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:** [February 23, 2024, 3:01pm UTC](https://discourse.julialang.org/t/fast-solution-of-a-two-point-boundary-value-problem-with-shooting/110531/2 "2024-02-23T15:01:05Z")

</div>

The MIRK methods should get a major speed improvement in about a month or two, there’s a lot that hasn’t been merged. If you open an issue on the repo we can track it there.