In many fields it is common to encounter systems of differential equations with boundary conditions that occur at limits (rather than at finite points).

For example, the Neoclassical growth model boils down to a system of two boundary value ODEs:

\dot{C}(t)= \frac{ 1 }{ \gamma } \left( \alpha z (K(t))^{\alpha -1} - \delta - \rho \right) \times C(t)

\dot{K}(t) =z (K(t))^{\alpha} - \delta K(t) - C(t)

K(0)=K_{0} given (Initial Condition)

\lim_{T\to \infty} e^{-\rho T} \lambda(T) K(T) =
e^{-\rho T} C(T)^{-\gamma} K(T)
\to 0 (Terminal Condition)

This system has two variables, one state variable K(t), one control variable C(t).

The initial condition depends only on the state variable.

The terminal condition depends on both variables.

This is standard in optimal control problems.

AFAIK currently SciML does not allow boundary conditions at limits.

This post will attempt to “hack” a solution.

First, let’s look at the closed form solution (for a subset of parameters where we know `C0`

):

```
using DifferentialEquations, Plots;
function ngm!(dX, X, θ, t)
C, K = X # variables: Consumption, Capital
ρ, γ, z, α, δ = θ # parameters
KN = abs(K)
#
dC = ((α*z*(KN^(α-1.0)) -δ -ρ)/γ)*C # Euler Eqn
dK = z*KN^(α) - δ*KN - C # Law of Motion
dX .= (dC, dK) # Time derivatives
end
α=0.20; γ=6.0; δ=0.25; z=1.0; ρ = δ*(α*γ-1.0)
K0 = 1.0; C0 = (1-(1/γ))*z*(K0)^α;
KSS= ((α)/(ρ+δ))^(1/(1-α)); CSS = (KSS)^α - δ*KSS;
X0 = [C0, K0] # initial conditions
T = 100.0 # time horizon
tspan = (0.0, T) # time span
θ = [ρ, γ, z, α, δ] # parameters
#
TC(t,c,k;γ=γ,ρ=ρ) = exp(-ρ*t) * (c^(-γ)) * k # Terminal condition
#
prob = ODEProblem(ngm!, [C0, K0], tspan, θ); # problem
sol = solve(prob,abstol=1e-10,reltol=1e-10)
#
sol_C(t) = sol(t, idxs=1)
sol_K(t) = sol(t, idxs=2)
#
tt=0.0:0.01:T
plot()
plot!(tt, sol_C, label="C(t)", color=1)
plot!(tt, sol_K, label="K(t)", color=2)
plot!(tt, t -> TC(t, sol_C(t), sol_K(t)), label="TC(t)", color=3)
plot!([CSS KSS 0.0], seriestype = :hline, lab=["CSS" "KSS" ""], color=[1 2 3], l=:dash)
```

Observe:

1: the control variable C(t) goes to its steady state $C_{SS}

2: the state variable K(t) goes to its steady state K_{SS}

3: the terminal condition holds TC(t)\to 0

In most cases we don’t have a closed form solution, so we don’t know `C0`

.

Consider a shooting method, to find the appropriate `C0`

such that the terminal condition holds TC(t)\to 0.

First consider a grid of `C0`

```
grid_c0 = .75*K0 :0.01 : 1.25*K0;
grid_TC =[]
for c0 in grid_c0
prob = ODEProblem(ngm!, [c0, K0], tspan, θ); # problem
#sol = solve(prob, Tsit5(), abstol=1e-10,reltol=1e-10, maxiters = 10_000_000)
sol = solve(prob, Tsit5())
if sol.retcode == :Success
_tc = TC(T, sol(T, idxs=1), sol(T, idxs=2))
else
_tc = -100
end
push!(grid_TC, _tc)
println(c0, " ", _tc)
end
cl(x) = clamp(x,-100,100)
plot(grid_c0, cl.(grid_TC) )
```

It appears the appropriate `C0`

(which satisfies the TC) is between 0.8 and 0.9.

This is good bc we know the true `C0=0.833333`

Another approach is to use a nonlinear solver (or optimized) to find the unique value of `C0`

at which the terminal condition is satisfied.

```
using Optim;
function f(x; K0=K0,T=100.0,tspan=(0.0, T),θ=θ)
prob = ODEProblem(ngm!, [x[1], K0], tspan, θ); # problem
#sol = solve(prob,Tsit5(),reltol=1e-10,maxiters=100_000,abstol=1e-8)
sol = solve(prob,Tsit5())
if sol.retcode == :Success #C, K
cc(x) = sol(x, idxs=1)
KK(x) = sol(x, idxs=2)
_tc1 = TC(T, abs(cc(T)), abs(KK(T)))
else
_tc1 = -10_000.0
end
return abs(_tc1)
end
f([0.75]) #
result = optimize(f, [0.75])
xx=result.minimizer
result.minimum
```

`Optim.jl`

returns `C0=0.8333338128872135`

which is close to the true `C0=.83333333333`

.

- Does anyone know of a good way to do this in the Julia ecosystem or otherwise
- I’ve tried in on larger problems (2 state/2 control) variables w/o success
- I arbitrarily choose
`T=100.0`

, is there a better way for finding the limit of`TC(t)`

maybe when some tolerance is reached?