CODEWITHSUNDEEP
juliaodedifferential-equationsdifferentialequations.jl
Ekaterina Ivanushkina
Ekaterina Ivanushkina

Reputation: 1

Julia - Second-order ODE gives wrong results

I am trying to solve a 2nd order ODE (r+delta)V=a* exp(bz)+ mu * V' +1/2 sigma^2*V''

This is the code that I wrote

using DifferentialEquations, Plots

p = -0.3, 0.04, 0.05, 0.1, 0.66, 3.0
μ, σ2, r, δ, α, w = p
a = (1-α)*(α/w)^(α/(1-α))
b = 1/(1-α)

function value_firm(dv, v, p, z)
    μ, σ2, r, δ, α, w = p
    a = (1-p[5])*(p[5]/p[6])^(p[5]/(1-p[5]))
    b = 1/(1-p[5])
    2/σ2*((p[3]+p[4])*v-a*exp(b*z)-μ*dv)   
end


nu  = r+δ-b*μ-b^2*σ2/2
e = μ/σ2+sqrt((μ/σ2)^2+r*2/σ2)
z_bar = log(w/(r+δ)*e/(e+b)/a*nu)/b


zspan = (z_bar, 2.0)
v_init = w/(r+δ)
dv_init = 0.0

an_sol(z) = a/nu*exp(b*z) + v_init*b/(b+e)*exp(-e*(z-z_bar))

prob = SecondOrderODEProblem(value_firm, dv_init, v_init, zspan, p)
sol = solve(prob, Vern7(), reltol=1e-8, abstol=1e-8)

v_sol = [u[2] for u in sol.u]
plot(sol.t, v_sol, linewidth=2.0, label="numerical")
plot!(sol.t, an_sol.(sol.t), color="black", linewidth=2.0, label="analytical")

However, for some reason, the solution that I get is very different from the (correct) analytical solution an_sol. I am not that experienced with differential equations in Julia, so any help will be appreciated.

Upvotes: 0

Views: 50

Answers (0)

Related Questions