I’m trying to solve a boundary value problem but still have some gaps in my understanding after reading the tutorial. Specifically, I’m trying to understand how boundary constraints should be encoded. I believe I’ve “figured it out” in that it appears to be working, but I won’t say I’m confident in the implications of what I’ve done.
As a MWE, I have defined the following problem
using DifferentialEquations
using LinearAlgebra
function dynamics(u,p,t)
x = u[1:3]
y = u[4:6]
dx = y
dy = cross(x, y)
return [dx;dy]
end
function boundaryconditions!(resid,u,p,t)
xinit = u[begin][1:3]
yfinal = u[end][4:6]
resid[1:3] = xinit - p.x0
resid[4:6] = yfinal - p.yf
return resid
end
tspan = (0.0,1.0)
x0 = Float64[1,0,0]
y0 = Float64[0,1,0]
# with y0 known, this is a simple ODE
truesol = ODEProblem(dynamics, [x0; y0], tspan) |> solve
yf = truesol.u[end][4:6]
# now fix yf and attempt to find the y0 that yields it
bvsol = BVProblem(dynamics, boundaryconditions!, [x0; zero(yf)], tspan, (;x0,yf)) |> solve
The dynamics are just made up for the purpose of this example. I have intended to impose the vector equality constraints that x(0) = \text{x0} and y(1) = \text{yf}. Or at least the solution appears to indicate I have imposed this correctly…
From examining the evaluation of this problem, it appears that resid has the same dimension as the state (size(resid) == (6,), in this case). I don’t exactly understand why this is. I suppose extra constraints are likely to yield infeasible problems and for fewer constraints one could simply set the extra residuals to zero, so that matching the state dimension is suitable for most problems?
It appears that the problem still yields an identical solution if I re-order the values within resid. Is there any significance to the ordering of the entries of resid? Are there implementation concerns such that one ordering might be more efficient? I suppose what I might be looking for here is a reference to how the residual is used in solving the problem.
The tutorial and documentation seems to allege that the main difference in defining a function for a TwoPointBVProblem is that (within the boundaryconditions! function?) u will be a NTuple{2} of the initial and final state, rather than a Vector of states. However, this does not appear to be true and I don’t notice any other difference. What is the practical difference between a general BVProblem and a TwoPointBVProblem? Will the TwoPointBVProblem be more efficient to solve, where applicable?
Thanks for answers to any subset of these questions!