Hey, Im trying to solve a Boundary Value Problems with Shooting() and DifferentialEquations.jl.
My problem is formulated as this:
This problem is an IVP with known initial values and has also a known analytical solution.
The problem is solved as an ODE:
using DifferentialEquations
const μ = 0.04 # equity premium
const σ = 0.13 # volatility of equity return (std)
const γ = 5.0 # risk aversion parameter
const ρ = 0.03 # discount rate
const r = 0.02 # risk-free rate
t_min = 1.0 # start
t_max = 50.0 # end
tspan = (t_min,t_max) # fund value interval
step_size = 1 # saved solutin every step of size
η = (1/γ)ρ+(1-1/γ)*(r+(1/2)*(μ^2/(γ*σ^2))) # calculate eta
# Defining the IVP and solving it.
x₀ = η^(-γ)*(t_min^(1-γ)/(1-γ)) # initial value for the value function
y₀ = η^(-γ)*t_min^(-γ) # initial value for the first derivative of the value function
u₀ = [x₀, y₀] # initial state vector
function nohabits!(du,u,p,t)
V = u[1]
V_a = u[2]
du[1] = V_a # the first derivative of the value (dx) function equals y
du[2] = (0.5*V_a^2*(μ/σ)^2)/((γ/(1-γ))*V_a^(1-1/γ)-ρ*V+r*t*V_a) # equation 12a from the notes
end
prob = ODEProblem(nohabits!,u₀,tspan,maxiters=10e5) # define the problem
@time sol = solve(prob,reltol=1e-12,abstol=1e-12, saveat=step_size) # solve the problem
The solution time is: 35.300468 seconds (86.01 M allocations: 4.058 GiB, 6.43% gc time) for the first run. And for the 2nd run: 0.000710 seconds (881 allocations: 75.578 KiB).
If the tolerance level is increased to reltol=1e-8, abstol=1e-8, the error becomes present.
And at a tolerance level of reltol=1e-9,abstol=1e-9 it becomes insignificant.
My true problem has the same form as the previous, but the initial conditions are not known.
Hence, I want to solve the previous problem as a BVP, using Shooting() to find the initial conditions before attempting to solve my true problem.
This problem is formulated as this:
# Building and solve the BVP
function nohabits2!(du,u,p,t)
du[1] = u[2]
du[2] = (0.5*u[2]^2*μ^2/σ^2)/((γ/(1-γ))*max(u[2],0)^(1-1/γ)-ρ*u[1]+r*t*u[2])
end
off_set = 0
x_guess = x₀ + off_set # My guess of V_init
V_a_end = sol.u[end][2] # Last value of V_a in the known solution
u₀_2 = [x_guess, y₀] # The initial guess
function bc!(residual,u,p,t)
residual[1] = u[end][2] - V_a_end
end
bvp3 = BVProblem(nohabits2!, bc!, u₀_2, tspan, maxiters=1e5)
@time sol3 = solve(bvp3, Shooting(Tsit5()), isoutofdomain=(u,p,t) -> u[2]<0, saveat=step_size)
And solved by guessing the correct initial conditions.
This problems solves in: 10.299134 seconds (8.04 M allocations: 378.491 MiB, 1.45% gc time)
And, retcode: MaxIters.
By increasing maxiters=1e7
This challenge diapers, and the problem solves quickly but with an high error:
10.142641 seconds (3.72 M allocations: 178.774 MiB, 0.46% gc time)
retcode: Success
When trying to reduce the error by reducing the tolerance level and max iterations: maxiters=1e8 and reltol=1e-5,abstol=1e-5. The solution times become: 199.300782 seconds (3.72 M allocations: 178.776 MiB, 0.02% gc time)
A tolerance level at reltol=1e-5,abstol=1e-5 yields a high error. And from the first problem we know that the tolerance level should be reltol=1e-9,abstol=1e-9.
In the search for a proper solution I have tried several other algorithms, such as: AutoTsit5(Rosenbrock23()), AutoVern7(Rodas5()), Vern7(), Vern9() etc.
with no success.
So the question is:
Why is the BVP so hard to solve when the initial guess is correct?
From DifferentialEquations.jl 3.0 and a Roadmap for 4.0, DifferentialEquations.jl 3.0 and a Roadmap for 4.0 | juliabloggers.com
Flesh out the BVP solvers
Our Shooting methods are very flexible, but they will never do well on problems which are sensitive to the initial condition. We need to instead make our MIRK methods get all of the bells and whistles: continuous extension, adaptivity, etc., and we need to wrap some of the classic Netlib solvers into the same interface so we can start to do some comprehensive benchmarking.
Is this the wrong track to find the initial conditions for my BVP?
If so, is there any methods of recomandation?
Or is it just a blunder in my code? 