# Two point boundary value problem solution not obeying boundary condition

**URL:** <https://discourse.julialang.org/t/two-point-boundary-value-problem-solution-not-obeying-boundary-condition/109633>\
**Category:** General Usage\
**Tags:** diffeq, sciml, differentialequation\
**Created:** [February 2, 2024, 5:06pm UTC](https://discourse.julialang.org/t/two-point-boundary-value-problem-solution-not-obeying-boundary-condition/109633 "2024-02-02T17:06:01Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![bast](https://avatars.discourse-cdn.com/v4/letter/b/3e96dc/32.png) [@bast](https://discourse.julialang.org/u/bast)\
**Post date:** [February 2, 2024, 5:06pm UTC](https://discourse.julialang.org/t/two-point-boundary-value-problem-solution-not-obeying-boundary-condition/109633/1 "2024-02-02T17:06:01Z")

</div>

Hello,

I am trying to replicate the solution of an optimal control problem ([Redirecting](https://doi.org/10.1016/j.jedc.2012.04.006)).

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:

```julia
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

```julia
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)

function get_steady_state(par)
    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

barx, barλ = get_steady_state(par)

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

```

Boundary conditions

```julia
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](https://docs.sciml.ai/DiffEqDocs/latest/tutorials/bvp_example/)) it returns an error

```julia
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

```julia

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  
 ![fishery_bvp](https://global.discourse-cdn.com/julialang/original/3X/6/7/670c366f99432b92265dd2ad17edfc5deb747cf6.png)  
where for some reason the solution diverges at the end, violating the boundary condition.

Help would be appreciated!

---

<div class="post-metadata">

**Author:** ![bast](https://avatars.discourse-cdn.com/v4/letter/b/3e96dc/32.png) [@bast](https://discourse.julialang.org/u/bast)\
**Post date:** [February 4, 2024, 10:47am UTC](https://discourse.julialang.org/t/two-point-boundary-value-problem-solution-not-obeying-boundary-condition/109633/2 "2024-02-04T10:47:36Z")

</div>

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

```julia
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.