SciML / DifferentialEquations.jl

Multi-language suite for high-performance solvers of differential equations and scientific machine learning (SciML) components. Ordinary differential equations (ODEs), stochastic differential equations (SDEs), delay differential equations (DDEs), differential-algebraic equations (DAEs), and more in Julia.
https://docs.sciml.ai/DiffEqDocs/stable/
Other
2.85k stars 226 forks source link

Domain Errors and Log Transforms #886

Open azev77 opened 2 years ago

azev77 commented 2 years ago

It works for the following set of parameters:

using DifferentialEquations, Plots;
function ngm!(dX, X, θ, t)
    C, K          = X
    ρ, γ, z, α, δ = θ
    dC = ((α*z*(K^(α-1.0)) -δ - ρ)/(γ)) * (C)
    dK = z*K^(α) - δ*K - C
    dX .= (dC, dK)
end

α=0.20; γ=6.0; δ=0.25; z=1.0; ρ = δ*(α*γ-1.0)
K0 = 1.0; C0 = (1-(1/γ))*z*(K0)^α;

X0 = [C0, K0]
tspan = (0.0, 2.0)
θ = [ρ, γ, z, α, δ]
prob = ODEProblem(ngm!, X0, tspan, θ)
sol = solve(prob)
plot(sol)

Problem 1: if I widen tspan it reports success but the solution is no longer accurate at all for large t. Incorrectly reporting a success is risky. Unless it's somehow within the tolerance...

tspan = (0.0, 100.0)                   # time span
prob = ODEProblem(ngm!, X0, tspan, θ) # problem
sol = solve(prob)
plot(sol, label=["C(t)" "K(t)"])
plot!([CSS KSS 0.0],  seriestype = :hline, lab=["CSS" "KSS" ""], color="grey", l=:dash)

Partial sol: adjusting the tolerance improves it (slower to solve), but it's still not good

tspan = (0.0, 100.0)                   # time span
prob = ODEProblem(ngm!, X0, tspan, θ) # problem
sol = solve(prob,reltol=1e-12)
plot(sol, label=["C(t)" "K(t)"])
plot!([CSS KSS 0.0],  seriestype = :hline, lab=["CSS" "KSS" ""], color="grey", l=:dash)

Main problem, when I increase α from 0.2 to 0.5 I get a DomainError

tspan = (0.0, 30.0) 
α=0.50; γ=6.0; δ=0.25; z=1.0; ρ = δ*(α*γ-1.0)
θ = [ρ, γ, z, α, δ]
prob = ODEProblem(ngm!, X0, tspan, θ) # problem

julia> sol = solve(prob)
ERROR: DomainError with -0.3526915679295306:
Exponentiation yielding a complex result requires a complex argument.
Replace x^y with (x+0im)^y, Complex(x)^y, or similar.
Stacktrace:
  [1] throw_exp_domainerror(x::Float64)
    @ Base.Math .\math.jl:37
  [2] ^
    @ .\math.jl:901 [inlined]
  [3] ngm!(dX::Vector{Float64}, X::Vector{Float64}, θ::Vector{Float64}, t::Float64)
    @ Main .\Untitled-1:6
  [4] ODEFunction
    @ C:\Users\azevelev\.julia\packages\SciMLBase\YasGG\src\scimlfunctions.jl:1624 [inlined]

The problem is the terms K^(α-1.0) and K^(α) in the DE. I believe this can be solved somehow w/ the transform exp(K) and then computing log of the solution to get back the original variable. Economists encounter these types of domain errors all the time. What are the best practices for doing this?

PS: I can eventually solve the problem w/ the new parameters in this simple problem, my point is I wanna know how to handle these issues in general.

ChrisRackauckas commented 2 years ago

Log transforms can make it numerically easier to stay positive. It's not automated so you'd have to do this by hand, for now, though there's work in MTK to do this.

https://diffeq.sciml.ai/stable/basics/faq/#My-ODE-goes-negative-but-should-stay-positive,-what-tools-can-help?