# Two-point boundary value problem in Julia DiffEq: "control hamiltonian"

**URL:** <https://discourse.julialang.org/t/two-point-boundary-value-problem-in-julia-diffeq-control-hamiltonian/36477>\
**Category:** Numerics\
**Tags:** diffeq\
**Created:** [March 25, 2020, 5:37am UTC](https://discourse.julialang.org/t/two-point-boundary-value-problem-in-julia-diffeq-control-hamiltonian/36477 "2020-03-25T05:37:48Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![jacob-roth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jacob-roth/32/1862_2.png) [@jacob-roth](https://discourse.julialang.org/u/jacob-roth)\
**Post date:** [March 25, 2020, 5:37am UTC](https://discourse.julialang.org/t/two-point-boundary-value-problem-in-julia-diffeq-control-hamiltonian/36477/1 "2020-03-25T05:37:48Z")

</div>

I’m interested in solving Hamilton’s equations of motion in DiffEq with two boundary conditions, an initial condition on the state variables and a terminal condition on the momentum variables. First I tried setting up a toy problem (P1) with these boundary conditions but was unable to solve it. Next I tried changing the boundary conditions (P2) to initial _and_ terminal conditions on the state variables but was again unable to solve it (I chose the final state to condition to be a point that was reachable using another solution method).

As an overview:  
Let x \in \mathbb{R}^n denote the state variables, y \in \mathbb{R}^n denote the momentum variables, and H(x,y) denote the Hamiltonian function H:\mathbb{R}^{2n} \to \mathbb{R}. I want to solve Hamilton’s equations of motion over the interval (0,T) by solving

\dot{x} = \nabla\_y H(x,y), \quad x(0) = x\_0 \\ \dot{y} = -\nabla\_x H(x,y), \quad y(T) = y\_T \tag{P1}

where y\_T := g(x\_T) for some function g: \mathbb{R}^n \to \mathbb{R}  
I wasn’t able to solve P1 so I tried solving P2

\dot{x} = \nabla\_y H(x,y), \quad x(0) = x\_0, x(T) = x\_T \\ \dot{y} = -\nabla\_x H(x,y) \tag{P2}.

I used `BVProblem` in both cases but didn’t get satisfactory answers and was wondering if I’m misusing `BVP`. I’ve posted my code below:

```julia
using DifferentialEquations
using ForwardDiff
using LinearAlgebra
using SparseArrays

##
## data
##

D = Dict()

D[:m] = 1.0

D[:A] = [
-1.0 0.0 1.0 1.0 0.0;
 0.0 -1.0 -1.0 0.0 0.0;
-1.0 1.0 0.0 0.0 0.0;
-1.0 0.0 0.0 -1.0 0.0;
 0.0 0.0 0.0 0.0 -1.0;
]

D[:S] = [
1.0 0.0 0.0 0.0 0.0;
0.0 1.0 0.0 0.0 0.0;
0.0 0.0 0.0 0.0 0.0;
0.0 0.0 0.0 1.0 0.0;
0.0 0.0 0.0 0.0 1.0;
]

D[:J] = [
 0.0 0.0 1.0 1.0 0.0;
 0.0 0.0 -1.0 0.0 0.0;
-1.0 1.0 0.0 0.0 0.0;
-1.0 0.0 0.0 0.0 0.0;
 0.0 0.0 0.0 0.0 0.0;
]

D[:B] = [
 20.0 -10.0 -10.0;
-10.0 20.0 -10.0;
-10.0 -10.0 20.0;
]

D[:G] = [
0.0 0.0 0.0;
0.0 0.0 0.0;
0.0 0.0 0.0;
]

D[:Y] = D[:G] + im*D[:B]

D[:P] = [-1.0, -2.0, +3.0]
D[:Q] = [-0.72716, +0.1629, +0.1]
D[:θ] = [0.0]
D[:V] = [1.05, 1.02]
D[:limit] = 5.8

D[:fr] = 2
D[:to] = 3

##
## hamiltonian
##

function helper(x::AbstractArray{T,1}, p) where T
    ww = x[1:2]
    a = x[3:4]; aa = [p[:D][:θ]; a]
    v = x[5]; vv = [p[:D][:V][1]; p[:D][:V][2]; v]
    i1 = 2
    i2 = 3
    i3 = 3
    df = -[
    -p[:D][:m] * ww ;
    vv[i1] * sum(
        vv[k] * ( p[:D][:G][i1,k] * cos(aa[i1] - aa[k]) + p[:D][:B][i1,k] * sin(aa[i1] - aa[k]) )
        for k in 1:3
                ) - p[:D][:P][i1] ;
    vv[i2] * sum(
        vv[k] * ( p[:D][:G][i2,k] * cos(aa[i2] - aa[k]) + p[:D][:B][i2,k] * sin(aa[i2] - aa[k]) )
        for k in 1:3
                ) - p[:D][:P][i2] ;
    1.0 * sum(
        vv[k] * ( p[:D][:G][i3,k] * sin(aa[i3] - aa[k]) - p[:D][:B][i3,k] * cos(aa[i3] - aa[k]) )
        for k in 1:3
                ) - p[:D][:Q][i3]/vv[i3] ;
    ]
    return p[:D][:A] * df
end
function dhelper_dx(x::AbstractArray{T,1}, p) where T
    ## = A * jac(df)
    return ForwardDiff.jacobian(xx -> helper(xx, p), x)
end

function H_x!(H_x, u, p, t)
    x = u[1:5]
    y = u[6:10]
    H_x .= helper(x, p) .+ (p[:D][:S] * y)
    nothing
end
function H_y!(H_y, u, p, t)
    x = u[1:5]
    y = u[6:10]
    H_y .= dhelper_dx(x, p)' * y
    nothing
end

function f!(f, u, p, t)
    x = u[1:5]
    y = u[6:10]
    fx = f[1:5]
    fy = f[6:10]
    H_y!(fx, u, p, t)
    H_x!(fy, u, p, t); fy .*= -1.0
    nothing
end

function F(x::AbstractArray{T,1}, p, t) where T
    a = x[3:4]; aa = [0.0; a]
    v = x[5]; vv = [p[:D][:V][1]; p[:D][:V][2]; v]
    fr = p[:D][:fr]
    to = p[:D][:to]
    return abs(p[:D][:limit] - (p[:D][:G][fr,to]^2 + p[:D][:B][fr,to]^2) * (
           vv[fr]^2 - 2.0 * vv[fr] * vv[to] * cos(aa[fr] - aa[to]) + vv[to]^2)
                               )
end
function g(x::AbstractArray{T,1}, p, t) where T
    ForwardDiff.gradient(xx -> F(xx, p, t), x)
end

function bc_p1!(resid, u, p, t)
    u0 = u[1]
    uT = u[end]
    x0 = u0[1:5]
    xT = uT[1:5]
    y0 = u0[6:10]
    yT = uT[6:10]
    resid[1:5] .= x0 .- p[:x0]
    resid[6:10] .= yT .- p[:λ] * g(xT, p, t)
    nothing
end
function bc_p2!(resid, u, p, t)
    u0 = u[1]
    uT = u[end]
    x0 = u0[1:5]
    xT = uT[1:5]
    y0 = u0[6:10]
    yT = uT[6:10]
    resid[1:5] .= x0 .- p[:x0]
    resid[6:10] .= xT .- p[:xT]
    nothing
end

##
## BVP
##

tspan = (0.0, 1.0)
x0 = [0.0, 0.0, 0.03275487325306676, -0.12652594142990567, 1.019442869573113]
xT = [0.0, 0.0, 0.07112546267753951, -0.1661538337407222, 1.014115273212668]

p = Dict()
p[:D] = D
p[:x0] = x0
p[:xT] = xT
p[:Tmin] = tspan[1]
p[:Tmax] = tspan[2]
p[:λ] = 10.0

y0 = [ones(2); 0.0; ones(2)]
u0 = [x0; y0]

```

Here’s the first problem:

```julia
bvp1 = BVProblem(f!, bc_p1!, u0, tspan, p)
@time sol1 = solve(bvp1, GeneralMIRK4(), dt=0.05)
594.490035 seconds (4.67 G allocations: 209.837 GiB, 4.05% gc time)
retcode: Failure
Interpolation: 1st order linear
t: 20-element Array{Float64,1}:
 0.0
 0.05263157894736842
 0.10526315789473684
 0.15789473684210525
 0.21052631578947367
 0.2631578947368421
 0.3157894736842105
 0.3684210526315789
 0.42105263157894735
 0.47368421052631576
 0.5263157894736842
 0.5789473684210527
 0.631578947368421
 0.6842105263157895
 0.7368421052631579
 0.7894736842105263
 0.8421052631578947
 0.8947368421052632
 0.9473684210526315
 1.0
u: 20-element Array{Array{Float64,1},1}:
 [-4.239414800591688e-9, 3.0325839589094507e-7, -0.4802163263936531, 0.9147598440962671, 0.8478514830618534, 2.002252652421146e7, 2.593265697402368e6, -2.4439318337548118e7, 2.0709904944606456e6, 3.433285100266265e7]
 [-8.450848765489579e-9, 5.97840061842031e-7, -0.45735552623518144, 0.8640305883108456, 0.8569290467849295, 2.002252652285557e7, 2.5932656974023674e6, -2.4439318337548118e7, 2.0709903922141853e6, 3.433285118328228e7]
 [-1.2607438274102086e-8, 8.873558311224565e-7, -0.4344975079942275, 0.8133066058151732, 0.8660057986188268, 2.002252652150169e7, 2.5932656974023674e6, -2.443931833754811e7, 2.0709902900265823e6, 3.433285136393691e7]
 [-1.6682239229703818e-8, 1.1609346110535076e-6, -0.4116435740543389, 0.7625903718963596, 0.8750813630212912, 2.0022526520149466e7, 2.5932656974023664e6, -2.4439318337548107e7, 2.0709901877998041e6, 3.433285154456614e7]
 [-2.064924426577092e-8, 1.4415543125599669e-6, -0.3887949323984538, 0.7118841702963683, 0.8841553867661249, 2.0022526518797923e7, 2.593265697402366e6, -2.4439318337548096e7, 2.0709900856376044e6, 3.43328517251444e7]
 [-2.4480190283580932e-8, 1.7185449869743696e-6, -0.36595265280792527, 0.6611900355213642, 0.8932275575674994, 2.002252651744593e7, 2.5932656974023655e6, -2.443931833754808e7, 2.0709899834948273e6, 3.433285190574506e7]
 [-2.8158793482526193e-8, 1.996807209777585e-6, -0.34311765805819033, 0.6105097083197393, 0.9022976093193785, 2.0022526516092375e7, 2.593265697402364e6, -2.4439318337548066e7, 2.0709898812840465e6, 3.433285208633519e7]
 [-3.1647098982555986e-8, 2.2570888952330476e-6, -0.32029068632040664, 0.5598445876120427, 0.9113653261844182, 2.002252651473649e7, 2.5932656974023636e6, -2.4439318337548047e7, 2.0709897791398256e6, 3.433285226689161e7]
 [-3.49411055772337e-8, 2.4899534409900337e-6, -0.29747227713057517, 0.5091956911861446, 0.9204305513443887, 2.0022526513378013e7, 2.5932656974023627e6, -2.443931833754803e7, 2.0709896769321298e6, 3.433285244748132e7]
 [-3.8017233411041604e-8, 2.700597969727869e-6, -0.27466276917972726, 0.4585636502832126, 0.9294931877679541, 2.002252651201736e7, 2.5932656974023622e6, -2.443931833754801e7, 2.0709895747855399e6, 3.433285262806584e7]
 [-4.085496575541288e-8, 2.8877283630091196e-6, -0.2518622646869529, 0.4079486601496, 0.938553205724572, 2.0022526510655455e7, 2.5932656974023604e6, -2.443931833754799e7, 2.0709894725849093e6, 3.433285280861921e7]
 [-4.343298874288074e-8, 3.075970972444515e-6, -0.22907067207172987, 0.3573505222570041, 0.9476106332970388, 2.0022526509293456e7, 2.59326569740236e6, -2.443931833754797e7, 2.0709893704436808e6, 3.43328529892419e7]
 [-4.5734012394993197e-8, 3.2347288116968842e-6, -0.20628765008927538, 0.30676859440540283, 0.9566655706611404, 2.002252650793238e7, 2.5932656974023585e6, -2.443931833754795e7, 2.0709892682467552e6, 3.433285316979684e7]
 [-4.774799371395633e-8, 3.386284351741266e-6, -0.18351267150499398, 0.25620185295719083, 0.9657181732764821, 2.0022526506572865e7, 2.593265697402358e6, -2.4439318337547936e7, 2.0709891660996801e6, 3.433285335039366e7]
 [-4.946690130702644e-8, 3.502980437935417e-6, -0.1607450039935337, 0.20564889789525265, 0.974768653962489, 2.0022526505215023e7, 2.5932656974023567e6, -2.4439318337547917e7, 2.070989063893104e6, 3.433285353099324e7]
 [-5.088819200456722e-8, 3.6176556537542164e-6, -0.13798372609723702, 0.15510797251160804, 0.9838172850654285, 2.002252650385853e7, 2.593265697402356e6, -2.4439318337547906e7, 2.0709889617288164e6, 3.433285371156347e7]
 [-5.2000766667494496e-8, 3.70649364667698e-6, -0.11522778473899299, 0.10457705385166409, 0.9928643730491353, 2.0022526502502795e7, 2.5932656974023553e6, -2.4439318337547895e7, 2.070988859506491e6, 3.433285389217903e7]
 [-5.279632178699331e-8, 3.788482814224316e-6, -0.09247596300768987, 0.054053830253660305, 1.001910273180109, 2.0022526501147173e7, 2.5932656974023553e6, -2.4439318337547887e7, 2.070988757332218e6, 3.433285407280616e7]
 [-5.3273693068031426e-8, 3.836980112030325e-6, -0.06972697687952568, 0.0035358392543325084, 1.0109553631900086, 2.002252649979109e7, 2.5932656974023543e6, -2.4439318337547883e7, 2.0709886551638795e6, 3.43328542534252e7]
 [-5.3434140036942706e-8, 3.8573638465594695e-6, -0.04697945044250697, -0.0469795402054251, 1.0200000429095584, 2.0022526498434436e7, 2.593265697402355e6, -2.443931833754788e7, 2.070988552905876e6, 3.4332854434119605e7]

```

And for the second problem, I decreased the timespan (and and `dt`), but it still took a while to return “Failure”:

```julia
tspan = (0.0, 0.01)
julia> @time sol2 = solve(bvp2, GeneralMIRK4(), dt=0.001)

100.947457 seconds (773.25 M allocations: 35.162 GiB, 4.00% gc time)
retcode: Failure
Interpolation: 1st order linear
t: 10-element Array{Float64,1}:
 0.0
 0.0011111111111111111
 0.0022222222222222222
 0.0033333333333333335
 0.0044444444444444444
 0.005555555555555556
 0.006666666666666667
 0.0077777777777777776
 0.008888888888888889
 0.01
u: 10-element Array{Array{Float64,1},1}:
 [-3.1719990334149255e-17, -1.8364175759621842e-17, 0.026909781904971872, -0.1111419104467956, 1.0162830927379543, -120233.16769388679, -310399.71017440595, 62807.773764772726, 669887.7669387929, -116763.06774559719]
 [-6.625685415860106e-17, 4.388983021918772e-17, 0.032472089055429615, -0.11896368225301811, 1.0163933103366456, -120233.16769389932, -310399.71017440595, 62807.773764772726, 669887.7651410811, -116763.06447560257]
 [-7.501206766719988e-17, 5.76884033248743e-17, 0.038034396323445345, -0.12678545418065085, 1.0165035279190149, -120233.16769391185, -310399.71017440595, 62807.773764772726, 669887.7633433691, -116763.06120560844]
 [-1.1116918486654339e-16, 7.250701254091926e-17, 0.043596703755792936, -0.1346072262780001, 1.016613745478567, -120233.16769392438, -310399.71017440595, 62807.773764772726, 669887.7615456574, -116763.05793561385]
 [-1.2308312421807804e-16, 1.4890042753488815e-17, 0.04915901138593286, -0.14242899857962282, 1.0167239630106568, -120233.1676939369, -310399.71017440595, 62807.773764772726, 669887.7597479458, -116763.05466561933]
 [-1.2286881881487558e-16, -1.8774318087602835e-17, 0.05472131923130163, -0.15025077110352694, 1.016834180512863, -120233.16769394944, -310399.71017440595, 62807.773764772726, 669887.757950234, -116763.05139562505]
 [-1.1730688355096783e-16, -5.320517827904788e-17, 0.06028362729189938, -0.1580725438497125, 1.0169443979851855, -120233.16769396196, -310399.71017440595, 62807.773764772726, 669887.756152522, -116763.0481256305]
 [-1.217251857208342e-16, -2.880356627782075e-17, 0.06584593555028938, -0.16589431680017158, 1.0170546154300453, -120233.1676939745, -310399.71017440595, 62807.773764772726, 669887.7543548103, -116763.0448556359]
 [-7.034791934430188e-17, -1.0769135627919558e-17, 0.07140824397301115, -0.1737160899203473, 1.0171648328520884, -120233.16769398702, -310399.71017440595, 62807.773764772726, 669887.7525570986, -116763.04158564172]
 [-5.67879223938636e-17, -4.573095647813079e-17, 0.0769705525132911, -0.18153786316193318, 1.0172750502578092, -120233.16769399955, -310399.71017440595, 62807.773764772726, 669887.750759387, -116763.03831564705]

```

Any help/perspective would be great!

---

<div class="post-metadata">

**Author:** ![jacob-roth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jacob-roth/32/1862_2.png) [@jacob-roth](https://discourse.julialang.org/u/jacob-roth)\
**Post date:** [March 25, 2020, 6:04am UTC](https://discourse.julialang.org/t/two-point-boundary-value-problem-in-julia-diffeq-control-hamiltonian/36477/2 "2020-03-25T06:04:45Z")

</div>

And as a follow up, I tried to encode this as a `DAEProblem` with multivariate algebraic function in an older version of Julia (1.1) but am not able to work out the dimensions:

```julia
p = Dict()
p[:D] = D
p[:x0] = x0
p[:xT] = xT
p[:Tmin] = tspan[1]
p[:Tmax] = tspan[2]
y0 = [ones(2); 0.0; ones(2)]
u0 = [x0; y0]
du0 = [zeros(size(u0)); zeros(5)]
dv = [fill(true, 10); fill(false, 5)]
prob = DAEProblem(dae!, du0, u0, tspan, p, differential_vars=dv)
sol = solve(prob, IDA())

julia> ERROR: DimensionMismatch("array could not be broadcast to match destination")
Stacktrace:
 [1] check_broadcast_shape at ./broadcast.jl:456 [inlined]
 [2] check_broadcast_axes at /Users/jakeroth/.julia/packages/DiffEqBase/DqkH4/src/diffeqfastbc.jl:19 [inlined]
 [3] instantiate at ./broadcast.jl:258 [inlined]
 [4] materialize! at ./broadcast.jl:756 [inlined]
 [5] (::getfield(Sundials, Symbol("##60#64")){DAEProblem{Array{Float64,1},Array{Float64,1},Tuple{Float64,Float64},false,Dict{Any,Any},DAEFunction{false,typeof(dae!),Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,Array{Bool,1}}})(::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::Dict{Any,Any}, ::Float64) at /Users/jakeroth/.julia/packages/Sundials/eAJ5O/src/common_interface/solve.jl:687
 [6] #__init#58(::Bool, ::Nothing, ::Float64, ::Bool, ::Bool, ::Nothing, ::Float64, ::Float64, ::Array{Float64,1}, ::Array{Float64,1}, ::Int64, ::Bool, ::Bool, ::Bool, ::Bool, ::Nothing, ::Bool, ::Bool, ::Bool, ::Nothing, ::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::typeof(DiffEqBase.__init), ::DAEProblem{Array{Float64,1},Array{Float64,1},Tuple{Float64,Float64},false,Dict{Any,Any},DAEFunction{false,typeof(dae!),Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,Array{Bool,1}}, ::IDA{:Dense}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}) at /Users/jakeroth/.julia/packages/Sundials/eAJ5O/src/common_interface/solve.jl:814
 [7] __init(::DAEProblem{Array{Float64,1},Array{Float64,1},Tuple{Float64,Float64},false,Dict{Any,Any},DAEFunction{false,typeof(dae!),Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,Array{Bool,1}}, ::IDA{:Dense}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}) at /Users/jakeroth/.julia/packages/Sundials/eAJ5O/src/common_interface/solve.jl:641
 [8] #init#380 at /Users/jakeroth/.julia/packages/DiffEqBase/DqkH4/src/solve.jl:20 [inlined]
 [9] init at /Users/jakeroth/.julia/packages/DiffEqBase/DqkH4/src/solve.jl:8 [inlined]
 [10] #__solve#19(::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::Function, ::DAEProblem{Array{Float64,1},Array{Float64,1},Tuple{Float64,Float64},false,Dict{Any,Any},DAEFunction{false,typeof(dae!),Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,Array{Bool,1}}, ::IDA{:Dense}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at /Users/jakeroth/.julia/packages/Sundials/eAJ5O/src/common_interface/solve.jl:10
 [11] __solve at /Users/jakeroth/.julia/packages/Sundials/eAJ5O/src/common_interface/solve.jl:10 [inlined] (repeats 5 times)
 [12] #solve#381 at /Users/jakeroth/.julia/packages/DiffEqBase/DqkH4/src/solve.jl:39 [inlined]
 [13] solve(::DAEProblem{Array{Float64,1},Array{Float64,1},Tuple{Float64,Float64},false,Dict{Any,Any},DAEFunction{false,typeof(dae!),Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,Array{Bool,1}}, ::IDA{:Dense}) at /Users/jakeroth/.julia/packages/DiffEqBase/DqkH4/src/solve.jl:27
 [14] top-level scope at none:0

```

---

<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:** [March 25, 2020, 10:16am UTC](https://discourse.julialang.org/t/two-point-boundary-value-problem-in-julia-diffeq-control-hamiltonian/36477/3 "2020-03-25T10:16:34Z")

</div>

Have you tried a Shooting method instead? Not sure that the returned solution makes sense.

```julia
julia> @time sol1 = solve(bvp1, Shooting(Tsit5()))
  0.389356 seconds (1.82 M allocations: 95.004 MiB, 2.54% gc time)
retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 6-element Array{Float64,1}:
 0.0                  
 9.999999999999999e-5 
 0.0010999999999999998
 0.011099999999999997 
 0.11109999999999996  
 1.0                  
u: 6-element Array{Array{Float64,1},1}:
 [0.0, 0.0, -0.04688552948990853, -0.04688553868683625, 1.0199999998607172, 1.0, 1.0, 0.0, 1.0, 1.0]
 [0.0, 0.0, -0.04688552948990853, -0.04688553868683625, 1.0199999998607172, 1.0, 1.0, 0.0, 1.0, 1.0]
 [0.0, 0.0, -0.04688552948990853, -0.04688553868683625, 1.0199999998607172, 1.0, 1.0, 0.0, 1.0, 1.0]
 [0.0, 0.0, -0.04688552948990853, -0.04688553868683625, 1.0199999998607172, 1.0, 1.0, 0.0, 1.0, 1.0]
 [0.0, 0.0, -0.04688552948990853, -0.04688553868683625, 1.0199999998607172, 1.0, 1.0, 0.0, 1.0, 1.0]
 [0.0, 0.0, -0.04688552948990853, -0.04688553868683625, 1.0199999998607172, 1.0, 1.0, 0.0, 1.0, 1.0]

```

---

<div class="post-metadata">

**Author:** ![mschauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mschauer/32/13946_2.png) [@mschauer](https://discourse.julialang.org/u/mschauer)\
**Post date:** [March 25, 2020, 12:01pm UTC](https://discourse.julialang.org/t/two-point-boundary-value-problem-in-julia-diffeq-control-hamiltonian/36477/4 "2020-03-25T12:01:46Z")

</div>

Would it make sense to consider a stochastic relaxation of the problem with additive white noise (with Brownian motions W^{(1)} and W^{(2)}, or only one BM on the generalised momentum variable)?

\dot{x} = \nabla\_y H(x,y) + d W^{(1)}, \quad x(0) = x\_0 \\ \dot{y} = -\nabla\_x H(x,y) + dW^{(2)}, \quad y(T) = y\_T \tag{P1}

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [March 25, 2020, 12:02pm UTC](https://discourse.julialang.org/t/two-point-boundary-value-problem-in-julia-diffeq-control-hamiltonian/36477/5 "2020-03-25T12:02:33Z")

</div>

Generally I would recommend spectral methods for BVPs, especially with control. I find

_Boyd, J. P. (2001). Chebyshev and fourier spectral methods._

a very nice reference. You can find it online.

---

<div class="post-metadata">

**Author:** ![jacob-roth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jacob-roth/32/1862_2.png) [@jacob-roth](https://discourse.julialang.org/u/jacob-roth)\
**Post date:** [March 25, 2020, 2:59pm UTC](https://discourse.julialang.org/t/two-point-boundary-value-problem-in-julia-diffeq-control-hamiltonian/36477/6 "2020-03-25T14:59:22Z")

</div>

> [@rveltz](#):
>
> Shooting(Tsit5())

Thanks @rveltz, I was able to increase the tspan to `tspan = (0.0, 16.0)` and decrease the step size to `dt=1e-4` (parameters that work in another method I have used), but I return a similar solution to yours:

```julia
julia> tspan = (0.0, 16.0)
(0.0, 16.0)

julia> bvp2 = BVProblem(f!, bc_p2!, u0, tspan, p)
BVProblem with uType Array{Float64,1} and tType Float64. In-place: true
timespan: (0.0, 16.0)
u0: [0.0, 0.0, 0.03275487325306676, -0.12652594142990567, 1.019442869573113, 1.0, 1.0, 0.0, 1.0, 1.0]

julia> @time sol2 = solve(bvp2, Shooting(Tsit5()), dt=1e-4)
  0.384874 seconds (1.36 M allocations: 74.673 MiB, 8.70% gc time)
retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 8-element Array{Float64,1}:
  0.0
  0.0001
  0.0011
  0.0111
  0.11109999999999999
  1.1111
 11.111099999999999
 16.0
u: 8-element Array{Array{Float64,1},1}:
 [0.0, 0.0, 0.05194016796530314, -0.14633988758531394, 1.0167790713928904, 1.0, 1.0, 0.0, 1.0, 1.0]
 [0.0, 0.0, 0.05194016796530314, -0.14633988758531394, 1.0167790713928904, 1.0, 1.0, 0.0, 1.0, 1.0]
 [0.0, 0.0, 0.05194016796530314, -0.14633988758531394, 1.0167790713928904, 1.0, 1.0, 0.0, 1.0, 1.0]
 [0.0, 0.0, 0.05194016796530314, -0.14633988758531394, 1.0167790713928904, 1.0, 1.0, 0.0, 1.0, 1.0]
 [0.0, 0.0, 0.05194016796530314, -0.14633988758531394, 1.0167790713928904, 1.0, 1.0, 0.0, 1.0, 1.0]
 [0.0, 0.0, 0.05194016796530314, -0.14633988758531394, 1.0167790713928904, 1.0, 1.0, 0.0, 1.0, 1.0]
 [0.0, 0.0, 0.05194016796530314, -0.14633988758531394, 1.0167790713928904, 1.0, 1.0, 0.0, 1.0, 1.0]
 [0.0, 0.0, 0.05194016796530314, -0.14633988758531394, 1.0167790713928904, 1.0, 1.0, 0.0, 1.0, 1.0]

```

and neither the initial nor final conditions are satisfied.

* * *

Thanks @mschauer, I think you’re right; the problems P1 and P2 are actually control versions of an underlying SDE. I guess that I could use the `SDE` package to try and simulate this, and then use Monte Carlo to get an estimate of the discrete \Delta W^{(1)} increments at each time step, but I was hoping to solve this deterministically.

* * *

@Tamas_Papp: thanks for the reference. I was leaning toward collocation schemes for enforcing the dynamics and will take a look through the book!

---

<div class="post-metadata">

**Author:** ![jacob-roth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jacob-roth/32/1862_2.png) [@jacob-roth](https://discourse.julialang.org/u/jacob-roth)\
**Post date:** [March 25, 2020, 3:04pm UTC](https://discourse.julialang.org/t/two-point-boundary-value-problem-in-julia-diffeq-control-hamiltonian/36477/7 "2020-03-25T15:04:47Z")

</div>

[moved previous replies into one post (per the Discourse suggestion)]

---

<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:** [March 25, 2020, 3:15pm UTC](https://discourse.julialang.org/t/two-point-boundary-value-problem-in-julia-diffeq-control-hamiltonian/36477/8 "2020-03-25T15:15:48Z")

</div>

> [@jacob-roth](#):
>
> bc\_p1

I am not sure your `bc_p1!` is right. It should be ` resid[6:10] .= yT .- g(xT, p, t)`, no?

Also, as a very minor point, I would advise you to use named tuple for the parameters : D = (m = 1.0, fr = 2, to = 3, etc

---

<div class="post-metadata">

**Author:** ![mschauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mschauer/32/13946_2.png) [@mschauer](https://discourse.julialang.org/u/mschauer)\
**Post date:** [March 25, 2020, 4:04pm UTC](https://discourse.julialang.org/t/two-point-boundary-value-problem-in-julia-diffeq-control-hamiltonian/36477/9 "2020-03-25T16:04:30Z")

</div>

In that case you problem seems to be accessible by our methods in

> **[Diffusion bridges for stochastic Hamiltonian systems and shape evolutions](https://arxiv.org/abs/2002.00885)**
>
> Stochastically evolving geometric systems are studied in shape analysis and computational anatomy for modelling random evolutions of human organ shapes. The notion of geodesic paths between shapes is central to shape analysis and has a natural...

Diffusion bridge = 2-sided boundary problems for SDEs

---

<div class="post-metadata">

**Author:** ![jacob-roth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jacob-roth/32/1862_2.png) [@jacob-roth](https://discourse.julialang.org/u/jacob-roth)\
**Post date:** [March 25, 2020, 4:44pm UTC](https://discourse.julialang.org/t/two-point-boundary-value-problem-in-julia-diffeq-control-hamiltonian/36477/10 "2020-03-25T16:44:18Z")

</div>

Thanks, yes, I mistakenly forgot to add `yT` to the BCs in P1! Have updated the original question. I can get it to return a “Success” code now for P1 (using Shooting), but the solution is in some sense trivial (so perhaps I’ve formulated it incorrectly…).

Also, why use a named tuple vs a dictionary?

---

<div class="post-metadata">

**Author:** ![jacob-roth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jacob-roth/32/1862_2.png) [@jacob-roth](https://discourse.julialang.org/u/jacob-roth)\
**Post date:** [March 25, 2020, 4:53pm UTC](https://discourse.julialang.org/t/two-point-boundary-value-problem-in-julia-diffeq-control-hamiltonian/36477/11 "2020-03-25T16:53:30Z")

</div>

Oh, very cool, I wasn’t aware of your approach!