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

Specific Domains Numerics
1

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:

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:

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”:

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!

2

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:

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
3

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

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]
1 Like
4

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}
1 Like
5

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.

2 Likes
6

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> 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!

7

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

8

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

1 Like
9

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

Diffusion bridge = 2-sided boundary problems for SDEs

2 Likes
10

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?

11

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