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 ofTC(t)maybe when some tolerance is reached?