# BVP Shooting method with free parameters

Hey, Im trying to solve a Boundary Value Problems with DifferentialEquations.jl.

My system of equations is composed by 4 equations of which I know 3 initial conditions. One initial condition u_0 and 2 free parameters (omega and sigma) need to be determined. The 4th boundary conditions are known.

As I understood from the minimal examples in the documentation, BVProblem can only find the initial conditions that satisfy the boundaries.

A possible solution seems to be the one proposed here: ODE-BVP’s with Unknown Parameter using an additional “dumb” equation with no dynamics for each unknown parameters.

Here’s my tentative with an encountered error:

function ShapeIntegrator!(dy,y,p,t)
omega = y[5]
sigma = y[6]

dy[1] = omega * y[2]

a1 = -y[2] / y[4] * cos(y[1])
a2 = (sin(y[1]) * cos(y[1])) / y[4]^2
a3 = (sin(y[1]) * y[3]) / (2π * y[4])
dy[2] = omega * (a1 + a2 + a3)

b1 = π *(y[2]^2 - sin(y[1])^2 / y[4]^2)
b2 = 2π * sigma
dy[3] = omega * (b1 + b2)

dy[4] = omega * cos(y[1])
du[5] = 0 #omega
du[6] = 0 #sigma
end

function bc1!(residual, y, p, t)
residual[1] = y[end][1] - psistar
residual[2] = y[end][2] - ustar
residual[4] = y[end][4] - xstar
end

# boundary conditions
psistar = pi
ustar = 1
xstar = 0.25

# free params
omega0 = 3
sigma0 = 0.0
u0 = 1.0

# init conditions
init = [0.035,u0,0.004, 0.035,omega0,sigma0]

s_span = [0,1.0]

bvp2 = BVProblem(ShapeIntegrator!, bc1!, init, s_span)

At this point BVProblem return an error:

MethodError: objects of type Vector{Float64} are not callable
Use square brackets [] for indexing an Array.

Is this the right direction?
If not, is there any methods of recomandation?

1 Like

I see a few typos in the code you posted, but even before fixing those I don’t get the error you mentioned. Here is a version that is able to construct the problem, but then crashes with an error when attempting to solve. Are you sure your problem definition is correct? Specifically, it is worrying that your initial guess produces NaN values.

using DifferentialEquations

function ShapeIntegrator!(dy, y, p, t)
psistar, ustar, xstar = p
omega = y[5]
sigma = y[6]

dy[1] = omega * y[2]

a1 = -y[2] / y[4] * cos(y[1])
a2 = (sin(y[1]) * cos(y[1])) / y[4]^2
a3 = (sin(y[1]) * y[3]) / (2π * y[4])
dy[2] = omega * (a1 + a2 + a3)

b1 = π * (y[2]^2 - sin(y[1])^2 / y[4]^2)
b2 = 2π * sigma
dy[3] = omega * (b1 + b2)

dy[4] = omega * cos(y[1])
dy[5] = 0 #omega
dy[6] = 0 #sigma
end

function bc1!(residual, y, p, t)
psistar, ustar, xstar = p
residual[1] = y[end][1] - psistar
residual[2] = y[end][2] - ustar
residual[4] = y[end][4] - xstar
end

function test()
# boundary conditions
psistar = pi
ustar = 1
xstar = 0.25

p = [psistar, ustar, xstar]

# free params
omega0 = 3
sigma0 = 0.0
u0 = 1.0

# init conditions
init = [0.035, u0, 0.004, 0.035, omega0, sigma0]

s_span = [0, 1.0]

bvp2 = BVProblem(ShapeIntegrator!, bc1!, init, s_span, p)
end
\$ julia
_
_       _ _(_)_     |  Documentation: https://docs.julialang.org
(_)     | (_) (_)    |
_ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
| | | | | | |/ _` |  |
| | |_| | | | (_| |  |  Version 1.10.4 (2024-06-04)
_/ |\__'_|_|_|\__'_|  |  Official https://julialang.org/ release
|__/                   |

(@v1.10) pkg> activate --temp
Activating new project at `/tmp/jl_EgfZA0`

# --- Snip installation output --

julia> include("bvp.jl")
test (generic function with 1 method)

julia> bvp2 = test()
BVProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 1.0)
u0: 6-element Vector{Float64}:
0.035
1.0
0.004
0.035
3.0
0.0

julia> soln = solve(bvp2)

Which results in

Warnings and error message
┌ Warning: First function call produced NaNs. Exiting. Double check that none of the
│ initial conditions, parameters, or timespan values are NaN.
└ @ OrdinaryDiffEq /home/user/.julia/packages/OrdinaryDiffEq/s27pa/src/initdt.jl:131
┌ Warning: Potential Rank Deficient Matrix Detected. Attempting to solve using
│ Pivoted QR Factorization.
└ @ NonlinearSolve /home/user/.julia/packages/NonlinearSolve/82B4C/src/internal/linear_solve.jl:162
┌ Warning: At t=0.8699345110024782, dt was forced below floating point epsilon
│ 1.1102230246251565e-16, and step error estimate = 2.529307595032869. Aborting.
│ There is either an error in your model specification or the true solution is
│ unstable (or the true solution can not be represented in the precision of
│ Float64).
└ @ SciMLBase /home/user/.julia/packages/SciMLBase/vhP5T/src/integrator_interface.jl:600
┌ Warning: Interrupted. Larger maxiters is needed. If you are using an integrator
│ for non-stiff ODEs or an automatic switching algorithm (the default), you may
│ want to consider using a method for stiff equations. See the solver pages for
│ more details (e.g.
│ https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/#Stiff-Problems).
└ @ SciMLBase /home/user/.julia/packages/SciMLBase/vhP5T/src/integrator_interface.jl:566
┌ Warning: Potential Rank Deficient Matrix Detected. Attempting to solve using
│ Pivoted QR Factorization.
└ @ NonlinearSolve /home/user/.julia/packages/NonlinearSolve/82B4C/src/internal/linear_solve.jl:162
┌ Warning: At t=0.8699345110024786, dt was forced below floating point epsilon
│ 1.1102230246251565e-16, and step error estimate = 2.5293053452384693. Aborting.
│ There is either an error in your model specification or the true solution is
│ unstable (or the true solution can not be represented in the precision of
│ Float64).
└ @ SciMLBase /home/user/.julia/packages/SciMLBase/vhP5T/src/integrator_interface.jl:600

# Above repeats many times with slightly variying details ...

┌ Warning: Potential Rank Deficient Matrix Detected. Attempting to solve using
│ Pivoted QR Factorization.
└ @ NonlinearSolve /home/user/.julia/packages/NonlinearSolve/82B4C/src/internal/linear_solve.jl:162
ERROR: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{SciMLBase.JacobianWrapper{…}, Float64}, Float64, 1})

Closest candidates are:
(::Type{T})(::Real, ::RoundingMode) where T<:AbstractFloat
@ Base rounding.jl:207
(::Type{T})(::T) where T<:Number
@ Core boot.jl:792
Float64(::IrrationalConstants.Sqrt2π)
@ IrrationalConstants ~/.julia/packages/IrrationalConstants/vp5v4/src/macro.jl:112
...

Stacktrace:
[1] convert(::Type{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{SciMLBase.JacobianWrapper{…}, Float64}, Float64, 1})
@ Base ./number.jl:7
[2] setindex!(A::Vector{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{…}, Float64, 1}, i1::Int64)
@ Base ./array.jl:1021
[3] _unsafe_copyto!(dest::Vector{Float64}, doffs::Int64, src::Vector{ForwardDiff.Dual{…}}, soffs::Int64, n::Int64)
@ Base ./array.jl:299
[4] unsafe_copyto!
@ ./array.jl:353 [inlined]
[5] _copyto_impl!
@ ./array.jl:376 [inlined]
[6] copyto!
@ ./array.jl:363 [inlined]
[7] copyto!
@ ./array.jl:385 [inlined]
[8] recursivecopy!
@ ~/.julia/packages/RecursiveArrayTools/K1bCr/src/utils.jl:62 [inlined]
[9] reinit!(integrator::OrdinaryDiffEq.ODEIntegrator{…}, u0::Vector{…}; t0::Float64, tf::Float64, erase_sol::Bool, tstops::Tuple{}, saveat::Tuple{}, d_discontinuities::Tuple{}, reset_dt::Bool, reinit_callbacks::Bool, initialize_save::Bool, reinit_cache::Bool, reinit_retcode::Bool)
@ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/s27pa/src/integrators/integrator_interface.jl:393
[10] reinit!
@ ~/.julia/packages/OrdinaryDiffEq/s27pa/src/integrators/integrator_interface.jl:380 [inlined]
[11] __single_shooting_loss!(resid_::Vector{…}, u0_::Vector{…}, p::Vector{…}, cache::OrdinaryDiffEq.ODEIntegrator{…}, bc::typeof(bc1!), u0_size::Tuple{…}, pt::SciMLBase.StandardBVProblem, resid_size::Tuple{…})
@ BoundaryValueDiffEq ~/.julia/packages/BoundaryValueDiffEq/KejJ4/src/solve/single_shooting.jl:106
[12] #51
@ ~/.julia/packages/BoundaryValueDiffEq/KejJ4/src/solve/single_shooting.jl:30 [inlined]
[13] NonlinearFunction
@ ~/.julia/packages/SciMLBase/vhP5T/src/scimlfunctions.jl:2300 [inlined]
[14] JacobianWrapper
@ ~/.julia/packages/SciMLBase/vhP5T/src/function_wrappers.jl:108 [inlined]
[15] auto_jacvec!(dy::Vector{…}, f::SciMLBase.JacobianWrapper{…}, x::Vector{…}, v::Vector{…}, cache1::Vector{…}, cache2::Vector{…})
@ SparseDiffTools ~/.julia/packages/SparseDiffTools/hj4gQ/src/differentiation/jaches_products.jl:17
[16] #90
@ ~/.julia/packages/NonlinearSolve/82B4C/src/internal/operators.jl:112 [inlined]
[17] JacobianOperator
@ ~/.julia/packages/NonlinearSolve/82B4C/src/internal/operators.jl:196 [inlined]
[18] (::NonlinearSolve.var"#157#169"{…})(du::Vector{…}, u::Vector{…}, fu::Vector{…}, p::Vector{…})
@ NonlinearSolve ~/.julia/packages/NonlinearSolve/82B4C/src/globalization/line_search.jl:163
[19] (::NonlinearSolve.var"#161#173"{…})(f::Function, p::Vector{…}, u::Vector{…}, du::Vector{…}, α::Float64, u_cache::Vector{…}, fu_cache::Vector{…}, deriv_op::NonlinearSolve.var"#157#169"{…})
@ NonlinearSolve ~/.julia/packages/NonlinearSolve/82B4C/src/globalization/line_search.jl:191
[20] #177
@ ~/.julia/packages/NonlinearSolve/82B4C/src/globalization/line_search.jl:205 [inlined]
[21] #__internal_solve!#174
@ ~/.julia/packages/NonlinearSolve/82B4C/src/globalization/line_search.jl:212 [inlined]
[22] __internal_solve!(cache::NonlinearSolve.LineSearchesJLCache{…}, u::Vector{…}, du::Vector{…})
@ NonlinearSolve ~/.julia/packages/NonlinearSolve/82B4C/src/globalization/line_search.jl:200
[23] __step!(cache::NonlinearSolve.GeneralizedFirstOrderAlgorithmCache{…}; recompute_jacobian::Nothing, kwargs::@Kwargs{})
@ NonlinearSolve ~/.julia/packages/NonlinearSolve/82B4C/src/core/generalized_first_order.jl:265
[24] __step!
@ ~/.julia/packages/NonlinearSolve/82B4C/src/core/generalized_first_order.jl:215 [inlined]
[25] #step!#218
@ ~/.julia/packages/NonlinearSolve/82B4C/src/core/generic.jl:50 [inlined]
[26] step!
@ ~/.julia/packages/NonlinearSolve/82B4C/src/core/generic.jl:45 [inlined]
[27] solve!(cache::NonlinearSolve.GeneralizedFirstOrderAlgorithmCache{…})
@ NonlinearSolve ~/.julia/packages/NonlinearSolve/82B4C/src/core/generic.jl:13
[28] #__solve#217
@ ~/.julia/packages/NonlinearSolve/82B4C/src/core/generic.jl:4 [inlined]
[29] __solve
@ ~/.julia/packages/NonlinearSolve/82B4C/src/core/generic.jl:1 [inlined]
[30] macro expansion
@ ~/.julia/packages/NonlinearSolve/82B4C/src/default.jl:282 [inlined]
[31] __solve(::NonlinearProblem{…}, ::NonlinearSolvePolyAlgorithm{…}; stats::SciMLBase.NLStats, alias_u0::Bool, verbose::Bool, kwargs::@Kwargs{…})
@ NonlinearSolve ~/.julia/packages/NonlinearSolve/82B4C/src/default.jl:255
[32] __solve(prob::BVProblem{…}, alg_::Shooting{…}; odesolve_kwargs::@NamedTuple{}, nlsolve_kwargs::@NamedTuple{}, verbose::Bool, kwargs::@Kwargs{…})
@ BoundaryValueDiffEq ~/.julia/packages/BoundaryValueDiffEq/KejJ4/src/solve/single_shooting.jl:77
[33] __solve
@ ~/.julia/packages/BoundaryValueDiffEq/KejJ4/src/solve/single_shooting.jl:1 [inlined]
[34] solve_call(_prob::BVProblem{…}, args::Shooting{…}; merge_callbacks::Bool, kwargshandle::Nothing, kwargs::@Kwargs{…})
@ DiffEqBase ~/.julia/packages/DiffEqBase/V6SCE/src/solve.jl:612
[35] solve_call
@ ~/.julia/packages/DiffEqBase/V6SCE/src/solve.jl:569 [inlined]
[36] #solve_up#53
@ ~/.julia/packages/DiffEqBase/V6SCE/src/solve.jl:1080 [inlined]
[37] solve_up
@ ~/.julia/packages/DiffEqBase/V6SCE/src/solve.jl:1066 [inlined]
[38] solve(prob::BVProblem{…}, args::Shooting{…}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{…}, kwargs::@Kwargs{…})
@ DiffEqBase ~/.julia/packages/DiffEqBase/V6SCE/src/solve.jl:1003
[39] __solve(::BVProblem{…}, ::Nothing; default_set::Bool, kwargs::@Kwargs{…})
@ DifferentialEquations ~/.julia/packages/DifferentialEquations/zEcqZ/src/default_solve.jl:14
[40] __solve
@ ~/.julia/packages/DifferentialEquations/zEcqZ/src/default_solve.jl:1 [inlined]
[41] #__solve#72
@ ~/.julia/packages/DiffEqBase/V6SCE/src/solve.jl:1394 [inlined]
[42] __solve
@ ~/.julia/packages/DiffEqBase/V6SCE/src/solve.jl:1386 [inlined]
[43] solve_call(::BVProblem{…}; merge_callbacks::Bool, kwargshandle::Nothing, kwargs::@Kwargs{})
@ DiffEqBase ~/.julia/packages/DiffEqBase/V6SCE/src/solve.jl:612
[44] solve_call(::BVProblem{…})
@ DiffEqBase ~/.julia/packages/DiffEqBase/V6SCE/src/solve.jl:569
[45] #solve_up#53
@ ~/.julia/packages/DiffEqBase/V6SCE/src/solve.jl:1072 [inlined]
[46] solve_up
@ ~/.julia/packages/DiffEqBase/V6SCE/src/solve.jl:1066 [inlined]
[47] solve(::BVProblem{…}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{…}, kwargs::@Kwargs{})
@ DiffEqBase ~/.julia/packages/DiffEqBase/V6SCE/src/solve.jl:1003
[48] solve(::BVProblem{…})
@ DiffEqBase ~/.julia/packages/DiffEqBase/V6SCE/src/solve.jl:993
[49] top-level scope
@ REPL[5]:1
Some type information was truncated. Use `show(err)` to see complete types.

You need to make sure you add the right amount of conditions so you don’t end up with a singular Jacobian. This is something that the least squares formalism would simplify, but it doesn’t exist yet.

@contradict The issue with this problem lies in some divergences with the initial values, which require an expansion to resolve. I initially left out the expansion for simplicity, but I’ll add it now.

Regarding the MethodError, this was caused by using an older version of Julia. After updating to 1.10, the error disappeared.

@ChrisRackauckas I’ve been using the least squares method, but it’s highly sensitive to the initial conditions of the free parameters, often leading to different local minima. This is one of the reasons I’m exploring an alternative approach using BVProblem!.

The current issue with the BVProblem approach is that I’m encountering a rank deficiency matrix error. The BVP seems to interpret the problem as having 3 boundary conditions but 6 unknowns (4 initial conditions + 2 free parameters). In reality, only one initial condition, u0, is unknown, making the problem well-posed with 3 boundary conditions and 3 unknowns (1 initial condition and 2 free parameters).

My questions are:

• How can I specify to BVProblem that not all initial conditions are unknown, and that three of them should remain constant (or be evaluated at each integration call, as I did with the least squares method in the residual function)?
• If this isn’t possible, would it be more convenient to stick with the least squares method? If so, do you have any suggestions for managing the instability of the parameters found?

Using Least squares method:

using DifferentialEquations
using Optim
using Interpolations

function ShapeIntegrator!(dy, y, p, t)
omega,sigma = p
k = 1.0
# dy[1] = DPsi_DS(omega, u)
dy[1] = omega * y[2]

a1 = -y[2] / y[4] * cos(y[1])
a2 = (sin(y[1]) * cos(y[1])) / y[4]^2
a3 = (sin(y[1]) * y[3]) / (2π * y[4] * k)
dy[2] = omega * (a1 + a2 + a3)

b1 = π * k * (y[2]^2 - sin(y[1])^2 / y[4]^2)
b2 = 2π * sigma * k
dy[3] = omega * (b1 + b2)

dy[4] = omega * cos(y[1])
end

# South Pole expansion
function South_X(s, omega, u0)
x1 = omega
x3 = -omega^3 * u0^2
return x1 * s + 1/6 * x3 * s^3
end

function South_Psi(s, omega, u0, sigma)
psi1 = omega * u0
psi3 = omega^3 * (3 * sigma * u0 - 2 * u0^3)
return psi1 * s + 1/6 * psi3 * s^3
end

function South_U(s, omega, u0, sigma)
return u0 + 1/2 * omega^2 * 0.5 * (3 * u0 * sigma - 2 * u0^3) * s^2
end

function South_Gamma(s, omega, u0, sigma, k)
gamma1 = omega * (2π * sigma * k)
gamma3 = 4/3 * π * k * u0 * omega^3 * (3 * sigma * u0 - 2 * u0^3)
return gamma1 * s + 1/6 * gamma3 * s^3
end

function InitialArcLength(p)
# Find the minimimum length s_init to not have a divergence in x(s).
omega = p[1]
u0 = p[2]
threshold = 0.035  # Changed from 0.01 for better stability
n = 1
delta_s = 0.0001
while true
if South_X(n * delta_s, omega, u0) > threshold
break
end
n += 1
end
return n * delta_s
end

function InitialValues(s_init, p)
omega = p[1]
u0 = p[2]
sigma = p[3]
k = p[4]

return [South_Psi(s_init, omega, u0, sigma),
South_U(s_init, omega, u0, sigma),
South_Gamma(s_init, omega, u0, sigma, k),
South_X(s_init, omega, u0),
]
end

function Residuales(parameters, boundary_conditions)
k = 1
omega, sigma, u0 = parameters

# Calculate initial arc length
s_init = InitialArcLength((omega, u0))

# Calculate initial values
z_init = InitialValues(s_init, (omega, u0, sigma, k))

# Evaluation points for the solution
s = range(s_init, stop=1, length=1000)

# Solve the IVP using Radau method
f = ODEFunction(ShapeIntegrator!,jac=ShapeJacobian!)
prob = ODEProblem(f, z_init, (s_init, 1), (omega, sigma))
sol = solve(prob, saveat=s)

if sol.retcode != :Success
println("integration failed")
return [1e3, 1e3, 1e3]
end

z_fina_num = sol.u[end]
psif, uf, xf = boundary_conditions

psi = z_fina_num[1] - psif
u = z_fina_num[2] - uf
x = z_fina_num[4] - xf
return [psi, u, x]
end

function main()
# Constitutive relations
Rparticle = 3
Rvesicle = 30
rpa = Rparticle / Rvesicle

# Wrapping angle
deg = 60

# Free parameters initial values
omega = 3.0
sigma = 0.019
u0 = 1.0

# Boundary conditions
psistar = π + phi
ustar = 1 / rpa
xstar = rpa * sin(phi)

# Check on xstar value
if xstar < 0.035
println("xstar: \$xstar")
error("xstar too small, can lead to divergences")
else
println("xstar: \$xstar")
end

boundary_conditions = [psistar, ustar, xstar]
free_params_extended = [omega, sigma, u0]

# Shooting algorithm and solver
result = optimize(x -> sum(abs2, Residuales(x, boundary_conditions)), free_params_extended, NelderMead())

@printf "Err: %.5f\n" result.minimum

best_parameters = result.minimizer
print(best_parameters)

Using BVProblem approach:

function ShapeIntegrator!(dy, y, p, t)
# psistar, ustar, xstar = p
omega = y[5]
sigma = y[6]

dy[1] = omega * y[2]

a1 = -y[2] / y[4] * cos(y[1])
a2 = (sin(y[1]) * cos(y[1])) / y[4]^2
a3 = (sin(y[1]) * y[3]) / (2π * y[4])
dy[2] = omega * (a1 + a2 + a3)

b1 = π * (y[2]^2 - sin(y[1])^2 / y[4]^2)
b2 = 2π * sigma
dy[3] = omega * (b1 + b2)

dy[4] = omega * cos(y[1])
dy[5] = 0 #omega
dy[6] = 0 #sigma
end

function bc1!(residual, y, p, t)
psistar, ustar, xstar = p
residual[1] = y[end][1] - psistar
residual[2] = y[end][2] - ustar
residual[4] = y[end][4] - xstar
end

# South Pole expansion
function South_X(s, omega, u0)
x1 = omega
x3 = -omega^3 * u0^2
return x1 * s + 1/6 * x3 * s^3
end

function South_Psi(s, omega, u0, sigma)
psi1 = omega * u0
psi3 = omega^3 * (3 * sigma * u0 - 2 * u0^3)
return psi1 * s + 1/6 * psi3 * s^3
end

function South_U(s, omega, u0, sigma)
return u0 + 1/2 * omega^2 * 0.5 * (3 * u0 * sigma - 2 * u0^3) * s^2
end

function South_Gamma(s, omega, u0, sigma, k)
gamma1 = omega * (2π * sigma * k)
gamma3 = 4/3 * π * k * u0 * omega^3 * (3 * sigma * u0 - 2 * u0^3)
return gamma1 * s + 1/6 * gamma3 * s^3
end

function InitialArcLength(p)
# Find the minimimum length s_init to not have a divergence in x(s).
omega = p[1]
u0 = p[2]
threshold = 0.035  # Changed from 0.01 for better stability
n = 1
delta_s = 0.0001
while true
if South_X(n * delta_s, omega, u0) > threshold
break
end
n += 1
end
return n * delta_s
end

function InitialValues(s_init, p)
omega = p[1]
u0 = p[2]
sigma = p[3]
k = p[4]

return [South_Psi(s_init, omega, u0, sigma),
South_U(s_init, omega, u0, sigma),
South_Gamma(s_init, omega, u0, sigma, k),
South_X(s_init, omega, u0),
]
end

function test()

# constitutive relations
Rparticle = 3
Rvesicle = 30
rpa = Rparticle / Rvesicle
# wrapping angle
phi = pi/6

# Boundary conditions
psistar = pi + phi
ustar = 1/rpa
xstar = rpa*sin(phi)
p = [psistar, ustar, xstar]

# check on xstar value
if xstar < 0.035
print("xstar:",xstar)
throw(DomainError())
else
print("xstar:",xstar)
end

# free params
omega0 = 3.0
sigma0 = 0.019
u0 = 1.0

k = 1
s_init = InitialArcLength((omega0, u0))
s_span = (s_init, 1.0)

# init conditions
init = InitialValues(s_init, (omega0, u0, sigma0, k))
init = [init; [omega0, sigma0]]

bvp2 = BVProblem(ShapeIntegrator!, bc1!, init, s_span, p)
end

bvp2 = test()
soln = solve(bvp2)

Note that BVPs are one of our biggest areas of work right now, so some things are moving fast. In particular, there is a (currently undocumented?) argument nlls

which can be used to switch the internal nonlinear solve into a NonlinearLeastSquaresProblem which should remove the square Jacobian requirement and thus the rank issue.

Could you expand a little bit more? How can I fix some of the initial conditions using this nlls argument?

And secondly, from the source code (and docs) it seems that I can use a function to set the initial condition, is this function called at every integration step?

See this test:

You say it’s an nlls and and give the residual vector definition (so that it knows the size of the constraints).

I have tried using the nlls keyword but I still got errors.

I checked using valid initial conditions and valid free params I can perform a successfull integration using ODEProblem without any NaN alert.

using DifferentialEquations
using Plots
using Interpolations

function ShapeIntegrator!(dy, y, p, t)
# omega = y[5]
# sigma = y[6]

omega, u0, sigma, k, psistar, ustar, xstar = p

dy[1] = omega * y[2]

a1 = -y[2] / y[4] * cos(y[1])
a2 = (sin(y[1]) * cos(y[1])) / y[4]^2
a3 = (sin(y[1]) * y[3]) / (2π * y[4])
dy[2] = omega * (a1 + a2 + a3)

b1 = π * (y[2]^2 - sin(y[1])^2 / y[4]^2)
b2 = 2π * sigma
dy[3] = omega * (b1 + b2)

dy[4] = omega * cos(y[1])
# dy[5] = 0 #omega
# dy[6] = 0 #sigma
end

function ShapeIntegratorExtended!(dy, y, p, t)

# omega, u0, sigma, k, psistar, ustar, xstar = p
omega = y[5]
sigma = y[6]
# print(omega)

dy[1] = omega * y[2]

a1 = -y[2] / y[4] * cos(y[1])
a2 = (sin(y[1]) * cos(y[1])) / y[4]^2
a3 = (sin(y[1]) * y[3]) / (2π * y[4])
dy[2] = omega * (a1 + a2 + a3)

b1 = π * (y[2]^2 - sin(y[1])^2 / y[4]^2)
b2 = 2π * sigma
dy[3] = omega * (b1 + b2)

dy[4] = omega * cos(y[1])
dy[5] = 0 #omega
dy[6] = 0 #sigma
end

function bc1!(residual, y, p, t)
# psistar, ustar, xstar = p
omega, u0, sigma, k, psistar, ustar, xstar = p

residual[1] = y[end][1] - psistar
residual[2] = y[end][2] - ustar
residual[3] = y[end][4] - xstar
end

# South Pole expansion
function South_X(s, omega, u0)
x1 = omega
x3 = -omega^3 * u0^2
return x1 * s + 1/6 * x3 * s^3
end

function South_Psi(s, omega, u0, sigma)
psi1 = omega * u0
psi3 = omega^3 * (3 * sigma * u0 - 2 * u0^3)
return psi1 * s + 1/6 * psi3 * s^3
end

function South_U(s, omega, u0, sigma)
return u0 + 1/2 * omega^2 * 0.5 * (3 * u0 * sigma - 2 * u0^3) * s^2
end

function South_Gamma(s, omega, u0, sigma, k)
gamma1 = omega * (2π * sigma * k)
gamma3 = 4/3 * π * k * u0 * omega^3 * (3 * sigma * u0 - 2 * u0^3)
return gamma1 * s + 1/6 * gamma3 * s^3
end

function InitialArcLength(p)
# Find the minimimum length s_init to not have a divergence in x(s).
omega, u0, sigma, k, psistar, ustar, xstar = p
threshold = 0.035  # Changed from 0.01 for better stability
n = 1
delta_s = 0.0001
while true
if South_X(n * delta_s, omega, u0) > threshold
break
end
n += 1
end
return n * delta_s
end

function InitialValues(p,s_init)
omega, u0, sigma, k, psistar, ustar, xstar = p
return [South_Psi(s_init, omega, u0, sigma),
South_U(s_init, omega, u0, sigma),
South_Gamma(s_init, omega, u0, sigma, k),
South_X(s_init, omega, u0),
]
end

function test()

# constitutive relations
Rparticle = 3
Rvesicle = 30
rpa = Rparticle / Rvesicle
# wrapping angle
phi = pi/3

# Boundary conditions
psistar = pi + phi
ustar = 1/rpa
xstar = rpa*sin(phi)

# check on xstar value
if xstar < 0.035
print("xstar:",xstar,"\n")
throw(DomainError())
else
print("xstar:",xstar,"\n")

end

# free params
omega0 = 3.419653288280937
sigma0 = 0.03890024246105399
u0 = 0.8710355880503611

# init conditions
k = 1
p = (omega0, u0, sigma0, k, psistar, ustar, xstar)
init = [0.030671103390489153, 0.8706571792130107, 0.008576555548218276, 0.035216903281918656]

s_init = InitialArcLength(p)
s_init = 0.0103
s_span = (s_init, 1.0)

# successful integration without any NaN values
prob = ODEProblem(ShapeIntegrator!, init, s_span,p)
sol = solve(prob)

# BVP problem, using 2 equations to represent the free params.
init = [0.030671103390489153, u0, 0.008576555548218276, 0.035216903281918656,omega0,sigma0]
bvp = BVProblem(BVPFunction{true}(ShapeIntegratorExtended!, bc1!; bcresid_prototype = zeros(3)),init, s_span, p; nlls=Val(true))
solve(bvp,MIRK6(),dt=1e-4)

end

bvp2 = test()

With this code it simply gets stuck (?) I am not sure if during the compilation or during the integration procedure

I tried your example with a simple numerical shooting approach to solve the boundary value problem directly (three equations, three unknowns); see below for the code. The main change is in bc1!, though I also changed your parameters to be a NamedTuple to make mistakes in ordering less likely.

The solver doesn’t converge very well, and I have to use a high tolerance for the ODE solver, which suggests that your problem is ill conditioned in some way.

using OrdinaryDiffEq
using NonlinearSolve

function ShapeIntegrator!(dy, y, p, t)
(; omega, sigma) = p  # use a named tuple to pass parameters

dy[1] = omega * y[2]

a1 = -y[2] / y[4] * cos(y[1])
a2 = (sin(y[1]) * cos(y[1])) / y[4]^2
a3 = (sin(y[1]) * y[3]) / (2π * y[4])
dy[2] = omega * (a1 + a2 + a3)

b1 = π * (y[2]^2 - sin(y[1])^2 / y[4]^2)
b2 = 2π * sigma
dy[3] = omega * (b1 + b2)

dy[4] = omega * cos(y[1])
end

# South Pole expansion
function South_X(s, omega, u0)
x1 = omega
x3 = -omega^3 * u0^2
return x1 * s + 1/6 * x3 * s^3
end

function InitialArcLength(p)
# Find the minimimum length s_init to not have a divergence in x(s).
(; omega, u0) = p
threshold = 0.035  # Changed from 0.01 for better stability
n = 1
delta_s = 0.0001
while true
if South_X(n * delta_s, omega, u0) > threshold
break
end
n += 1
end
return n * delta_s
end

function bc1!(residual, prob, u, p)
# psistar, ustar, xstar = p
(; psistar, ustar, xstar) = p

init = [0.030671103390489153, u[1], 0.008576555548218276, 0.035216903281918656]  # update the initial conditions
p′ = merge(p, (omega = u[2], sigma = u[3]))  # update the parameters

sol = solve(prob; u0 = init, p = p′, save_everystep = false, abstol=1e-10, reltol=1e-10)

residual[1] = sol[end][1] - psistar
residual[2] = sol[end][2] - ustar
residual[3] = sol[end][4] - xstar

return residual
end

function test()

# constitutive relations
Rparticle = 3
Rvesicle = 30
rpa = Rparticle / Rvesicle
# wrapping angle
phi = pi/3

# Boundary conditions
psistar = pi + phi
ustar = 1/rpa
xstar = rpa*sin(phi)

# check on xstar value
if xstar < 0.035
print("xstar:",xstar,"\n")
throw(DomainError())
else
print("xstar:",xstar,"\n")

end

# free params
omega0 = 3.419653288280937
sigma0 = 0.03890024246105399
u0 = 0.8710355880503611

# init conditions
k = 1
p = (; omega = omega0, u0, sigma = sigma0, k, psistar, ustar, xstar)
init = [0.030671103390489153, 0.8706571792130107, 0.008576555548218276, 0.035216903281918656]

s_init = InitialArcLength(p)
s_init = 0.0103
s_span = (s_init, 1.0)

# successful integration without any NaN values
ode_prob = ODEProblem(ShapeIntegrator!, init, s_span,p)

# BVP problem using numerical shooting
nl_prob = NonlinearProblem((res, u, p) -> bc1!(res, ode_prob, u, p), [init[2], omega0, sigma0], p)

return solve(nl_prob, TrustRegion(; autodiff = AutoFiniteDiff()); show_trace = Val(true))
end

bvp2 = test()
1 Like

In BifurcationKit, I tested a lot the shooting paradigm. My “general” conclusion is

• you need multiple shooting (unless the orbit is super attracting)
• you need an implicit solver with relatively high tolerances, I often use Rodas5P

You can see that in the tutorials. This is especially true when you do something fancy like computing curves of bifurcations of periodic orbits.

It would be nice to try orthogonal collocation on this, I will see if I have some time.

If it’s not well-conditioned then you really want to make use of collocation methods. Did you try any MIRK tableaus?

I tried MIRK methods in the last script I posted but it gets stuck for hours

@dawbarton Thank you for spending some time helping me!
This approach is similar to what I implemented (and posted above) using the least square optimizer to minimize the residual. Your approach completes the solve but the final result is not the one expected: Omega should be close to pi.

I think this is related to the fact that the initial conditions for the ODE depends from omega,sigma,u0 values. I will spend some time trying to use your approach but including an initial values update for each optimization step.

Try with a fixed dt, i.e. adaptive=false. And try something like coldae as well, it should just swap in as an algorithm choice. Any of those working?