Domain Errors and Log Transforms

[azev77]

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]

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?