# Two point boundary value problem solution not obeying boundary condition

Hello,

I am trying to replicate the solution of an optimal control problem (Redirecting).

The canonical system is given by

\begin{aligned} & \dot{x}(t)=x(t)\left(1-\frac{x(t)}{a}\right)-\frac{1}{(c+\tau)} \frac{x(t)^2}{1+x(t)^2}-u(t) x(t) \\ & \dot{\lambda}(t)=\lambda(t)\left(r-1+\frac{2 x(t)}{a}+\frac{2 x(t)}{(c+\tau)\left(1+x(t)^2\right)^2}\right)-u(t)(p-\lambda(t)) \end{aligned}

where x is the state variable, \lambda the costate, and u the control.

The maximizer is given by u=\frac{x}{2}(p-\lambda).

The boundary conditions are the initial condition x(0) =x_0 and the transversality condition \lim_{t\rightarrow \infty} e^{-rt}\lambda(t)=0.

This system has a hyperbolic saddle point \hat y and I am trying to solve for the stable path from the initial state to the steady state.

Let

\dot y(t) = \begin{pmatrix} \dot x(t)\\ \dot \lambda(t) \end{pmatrix} = g(y(t))

and denote by J(\hat y) the Jacobian matrix at the the steady state.

Let E^s(\hat y) be the eigenspace associated with the negative eigenvalues, which is tangent to the stable manifold at \hat y.

Then, for T large enough, the transversality condition can be approximated by F^\prime(y(T)-\hat y) =0, where F^\prime is a matrix spanning the orthogonal complement to the stable eigenspace.

My try of implementing this so far:

using DifferentialEquations
using Optim
using ForwardDiff
using LinearAlgebra
using Plots

r = 0.1
τ = 0.25
c = 0.375
a = 7.5
pr = 0.1

par = [r, τ, c, a, pr]

T = 500

G(X) = X * (1 - X / a)
P(X) = (1 / (c + τ)) * (X^2 / (1 + X^2))


Getting the Jacobian matrix and F

function f(u)

fish = u[1] / 2 * (pr - u[2])
du1 = G(u[1]) - P(u[1]) - u[1] * fish
du2 = u[2] * (r - 1 + 2 * u[1] / a + 2 * u[1] / ((c + τ) * (1 + u[1]^2)^2)) - fish * (pr - u[2])
return [du1, du2]
end
Df(u) = ForwardDiff.jacobian(f, u)

r, τ, c, a, pr = par
fish(x, barλ) = x * (pr - barλ) / 2
eqx(x, λ) = G(x) - P(x) - x * fish(x, λ)
eqλ(x, λ) = λ * (r - 1 + 2 * x / a + 2 * x / ((c + τ) * (1 + x^2)^2)) - fish(x, λ) * (pr - λ)
eqs(x, λ) = [eqx(x, λ), eqλ(x, λ)]
res = Optim.optimize(x -> sum(abs.(eqs(x[1], x[2]))), [3, 0.5], Optim.Options(g_tol=1e-15))
x = res.minimizer[1]
λ = res.minimizer[2]
return x, λ
end

J = Df([barx, barλ])
ev = eigvals(J)
stable_evecs = eigvecs(J)[:, ev.<0]
F = nullspace(stable_evecs')


Boundary conditions

x0 = barx + 0.1

function bc1!(residuala, sola, p)
par, x0, barx, barλ = p
residuala[1] = sola[1][1] - x0
end

function bc2!(residualb, solb, p)
par, x0, barx, barλ = p
residualb[1] = (F'*(solb[end].-[barx, barλ]))[1]
end



Here I also ran into a problem when trying to call sola(0.0)[1] as in the tutorial (Boundary Value Problems · DifferentialEquations.jl) it returns an error

ERROR: MethodError: objects of type Vector{Float64} are not callable
Use square brackets [] for indexing an Array.
Stacktrace:
[1] bc1!(residuala::SubArray{…}, sola::Vector{…}, p::Tuple{…})


Solving


strt = [barx, barλ]
prob = TwoPointBVProblem(f!, (bc1!, bc2!), strt, (0.0, 500.0), (par, x0, barx, barλ); bcresid_prototype=(zeros(1), zeros(1)))
@time sol = solve(prob, MIRK4(), dt=0.01)


Which yields

where for some reason the solution diverges at the end, violating the boundary condition.

Help would be appreciated!

I misunderstood the solution object in the boundary conditions in the two point boundary value problem vs the general boundary value problem. sola and solb only contain the value of the solution at the start and end points respectively, so solb[end] was a scalar.

The correct boundary conditions are

function bc1!(residuala, sola, p)
par, x0, barx, barλ = p
residuala[1] = sola[1] - x0
end

function bc2!(residualb, solb, p)
par, x0, barx, barλ = p
residualb[1] = (F'*(solb.-[barx, barλ]))[1]
end


which yield the expected solution.