Boundary conditions for `DifferentialEquations.BVProblem`

Specific Domains
1

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!

2

The residual function represents the constraints of the given system in the boundary, as long as we set it correctly, the order doesn’t affect the solving process.

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?

I think yes, if we randomly make up some additional constraints for the given BVP, there might not be numerical solutions.

The detailed usage of the residual function can be seen here for the shooting method and here for the MIRK method, the basics here is to construct a nonlinear problem for the loss function, and use the root finding algorithms to find the zero points of the residual function.

When we only set constraints in the initial state and final state(the start and end point of the time span), the TwoPointBVProblem is the same as BVProblem, so practically there are no differences.