The following two-point boundary value problem is difficult to solve because the stationary state has a saddle structure, with one positive and one negative eigenvalue: only one path leads to the stationary state — all other paths are unstable. Can the BVProblem options be tweaked to achieve convergence?
using DifferentialEquations
using Plots
"""Parameters (including the stationary state values)"""
p = Dict(:σ => 0.5, :K0 => 0.2, :C0 => 0.1, :Kss => 1, :Css => 0.3)
"""The Production function"""
F(x) = 0.4 * x^(0.5) - 0.1 * x
"""The Ramsey-Cass-Koopmans Dynamic system"""
function sys!(du,u,p::Dict,t)
σ = p[:σ]
K, C = u
du[1] = 0.4 * K^(0.5) - 0.1 * K - C
du[2] = (0.5 * 0.4 * K^(-0.5) - 0.1 - 0.1) * σ * C
end
"""boundary condition
The state variable K has given initial value K0.
The control variable C has known stationary-state value Css."""
function bc!(res,sol,p::Dict,t)
K0, Css = p[:K0], p[:Css]
res[1] = sol[1][1] - K0
res[2] = sol[end][1] - Css
end
# event handler to prevent variables going out of bounds
function set_callback(p::Dict)
Kss, Css = p[:Kss], p[:Css]
c1(u, t, integrator) = u[1] # stay in positive range
c2(u, t, integrator) = u[2] # stay in positive range
c3(u, t, integrator) = 2*Kss - u[1] # don't stray too far from stationary state
c4(u, t, integrator) = 2*Css - u[2] # don't stray too far from stationary state
function terminate_affect!(integrator)
terminate!(integrator)
end
# set up affect upon callback
cb1 = ContinuousCallback(c1, terminate_affect!)
cb2 = ContinuousCallback(c2, terminate_affect!)
cb3 = ContinuousCallback(c3, terminate_affect!)
cb4 = ContinuousCallback(c4, terminate_affect!)
cb = CallbackSet(cb1, cb2, cb3, cb4)
return cb
end
# SOLVE AS TWO-POINT BVP
const u0 = [p[:K0]; p[:C0]]
const tspan=(0, 5)
bvp = BVProblem(sys!, bc!, u0, tspan, p)
cb = set_callback(p)
# attempt to solve bvp with shooting method
sol = solve(bvp, callback=cb, Shooting(Vern7()))
# PHASE PLOT
plt = plot(xlims=(0,2*p[:Kss]), ylims=(0,2*p[:Css]), xlabel="K(t)", ylabel="C(t)",
extra_kwargs=Dict(:subplot=>Dict(:legend_hfactor=>1.5)))
# plot the nullclines
plot!(plt, F, 0, 2*p[:Kss], color=:black, linewidth=2, label="")
vline!(plt, [p[:Kss]], color=:black, linewidth=2, label="")
# mark the position of the stationary state
scatter!(plt, (p[:Kss],p[:Css]), label=["Stationary State" "."], color=:red)
# mark the position of the initial guess u0
scatter!(plt, (u0...,), label=["Initial Guess" "."], color=:green)
# add trajectory to plot
plot!(plt, sol, vars=(1,2), xlims=(0,2*p[:Kss]), ylims=(0,2*p[:Css]), color=:blue, linewidth=3, label="converging path", thickness_scaling=1)
display(plt)
This is a follow-up on solving-boundary-value-differential-equation-problems-for-economics. I decided to start a new topic since the problem was solved there via other methods.
Can the BVP solver find the stable path? This is what the path that converges to the saddle looks like.
And this is what my attempts to use BVProblem yielded:
