Help with DifferentialEquations.jl time bases and dependent functions

New to Julia
1

Hi everyone.
This is my first foray into the Julia world, and I’m a little stuck on the best way to frame a problem I would like to solve with DifferentialEquations.jl.

Consider a set of coupled ODEs:

\begin{align} \dot{C} &= E-\phi C\\ \dot{T} &= log(1+C) -T \\ \dot{E} &= \chi(T) E \end{align}

which are constructed to be non-dimensional (for example T = \hat{T}/\alpha, where \hat{T} has some associated unit removed by a known, constant \alpha), and solved in the time basis t=\hat{t}/\tau.
Here, \chi(T) = \tau[\hat{\chi}(\hat{T})-\mu], where \tau and \mu are constant, but as you can see, this is a conversion to the dimensional basis.

To start with, I’d create a representative function and set up the ODE solver:

function closed_system(du,u,p,t)
 du[1] = u[3]-p[1]*u[1}
 du[2] = log(1+u[1])-u[2]
 du[3] = χ(u[2])*u[3]
end

u0 = [380.0;0.7;160.0]
tspan = (2000.0/τ,22000.0/τ)
p = [1.0]
prob = ODEProblem(closed_system,u0,tspan,p)

with χ not implemented yet, since we need the dimensional value \hat{\chi}(\hat{T}) = \xi(1-\delta\cdot\hat{T}) (where \xi and \delta are known constants).
With a bit of conversion we can arrive at a complete solution:

function closed_system(du,u,p,t)
 χhat(That) = p[5]*(1-p[6]*That)
 χ(T) = p[2]*(χhat(T*p[4])-p[3])
 du[1] = u[3]-p[1]*u[1]
 du[2] = log(1+u[1])-u[2]
 du[3] = χ(u[2])*u[3]
end

u0 = [380.0;0.7;160.0]
tspan = (2000.0/τ,22000.0/τ)
#ϕ, τ, μ, α, ξ, δ
p = [1.0;50.0;0.01;4.328;0.01;0.5]
prob = ODEProblem(closed_system,u0,tspan,p)
soln = solve(prob)

So that all seems to be fine, but here’s where I hit a roadblock:

\begin{align} d\hat{W}/d\hat{t} &= \hat{\chi}(\hat{T})\hat{W} &\quad (1)\\ E &= [\tau\eta\exp(-\mu\hat{t})\hat{W}]/C_{pi} &\quad (2) \end{align}

New values I’ve not introduced yet are known constants, but as you can see I now need to couple this dimensional value of \hat{W} into the system of equations, since it is dependent on the \hat{\chi}(\hat{T}) result.

As its an equation in the dimensional basis though, I can’t just set it up as du[4] = χhat(T*p[4])*u[4], since t in closed_system is non-dimensional.

Question 1. Is it possible to add the ODE (1) into the system of equations if is uses a different time basis like this?

The second problem here is the equation (2) dependency. I’m fairly certain it would be possible to put function calls into the parameters portion for example, but since this dependency is again in connected to the dimensional basis through \hat{t} and \hat{W}, I’m at a loss how one would go about this at all.

Question 2. How could one frame the two equations (1) and (2) such that they could affect the closed_system function appropriately?

Any help would be appreciated. If you need any more information about my problem that I’ve left out for simplicity let me know.

2

I’m having trouble reading your full description, but it sounds like E is an algebraic variable in that last description. In that case, you have a DAE instead of an ODE Differential-algebraic system of equations - Wikipedia . There are two ways to solve this. The first way is you can add a mass matrix to the ODE solver so that it’s Mu' = f(u,p,t). That M can be singular, so if it’s a row of zeros then that corresponds to a constraint equation. The other way is to use the implicit DAE solvers. This is shown in the tutorial:

http://docs.juliadiffeq.org/latest/tutorials/dae_example.html

3

Thanks Chris. I’ll take a look into that and see how I get along.

4

Hmm, so I’m still uncertain about this system. If I were to define the template in the simplest of terms it would be:
\begin{align} F(\dot{C},\dot{T},\dot{E},t) &= 0\\ E(\hat{T}, \hat{W}, \hat{t}) &= 0 \end{align}
Where t=\hat{t}/\tau and T=\hat{T}/\alpha are simple conversions via normalising constants and the expression for E is tightly coupled to the results in F.
If we follow the form of the Roberts model
\begin{align} du &= f(u,p,t)\\ 0 &= g(u,p,t), \end{align}
E could indeed be considered as a representation of g on the surface, but in my system it’s dependent on T and as such I can’t seem to re-frame the problem to align with the examples here. Based on your live stream Chris, your expertise in these matters towers over mine. I’ll give you my full system and perhaps there’s a pathway you can see.

\begin{align} \dot{C}(t) &= E(t)-\phi C(t)\\ \dot{T}(t) &= \log[1+C(t)] -T(t) \\ \dot{E}(t) &= \chi[T(t)] E(t) \end{align}
where
\begin{align} t &= \hat{t}/\tau, &C(t) &= \hat{C(\hat{t})}/C_{pi} - 1, &T(t) &= \hat{T(\hat{t})}/\alpha,\\ \alpha &= \Delta T_x / log(2), &E(t) &= \frac{\tau\eta\exp(-\mu\hat{t})\hat{W(\hat{t})}}{C_{pi}}, &\phi &= \tau/\tau_C\\ \chi[T(t)] &= \tau\{\hat{\chi}[\hat{T}(\hat{t})]-\mu\}, &\hat{\chi}[\hat{T}(\hat{t})] &= \xi[1-\delta\cdot\hat{T}(\hat{t})], &\dot{\hat{W}}(\hat{t}) &= \hat{\chi}[\hat{T}(\hat{t})]\hat{W}(\hat{t}) \end{align}

I’ve tried to be explicit here with time dependent variables, with either ‘real’ time \hat{t} or dimensionless time t. All other values are constants.

If I attempt a DAE implementation, I get as far as

function f(out,du,u,p,t)
 out[1] = u[3]-p[1]*u[1] - du[1] #dC
 out[2] = log(1+u[1])-u[2] - du[2] #dT
 out[3] = p[2]*(p[5]*(1-p[6]*u[2]*p[4])-p[3])*u[3] - du[3] #dE
 out[4] = p[5]*(1-p[6]*u[2]*p[4])*u[4] - du[4]; #\hat{W}, but this should be against \hat{t}
 out[5] = (p[2]*p[6]*exp(-p[3]*t*p[2])*u[4])/p[7] #E, but out[5] should define u[3] really
end

where the \chi functions have been expanded and put in explicitly here. The two comments there are now my roadblocks. In the Roberts model, the algebraic property is ‘separate’ in the sense that there’s a constraint on the system through 1 = y_1+y_2+y_3, but here I seem to want to define define E (u[3]) rather than constrain it. Additionally, \dot{\hat{W}} (out[4]) is in the t basis rather than \hat{t}, which comes back to my additional question on if changing bases are possible.

Perhaps a simpler form factor is possible with an ODE implementation:

g = @ode_def WTest begin
  dW = ξ*(1-δ*T)*W
  dT = log(C/Cₚᵢ) - T/α
  dC = ((τ*η*exp(μ*t)*W)/Cₚᵢ)-ϕ*((C/Cₚᵢ)-1)
  dE = τ*((ξ*(1-δ*T))-μ)*((τ*η*exp(μ*t)*W)/Cₚᵢ)
    end ξ δ Cₚᵢ α τ η μ ϕ

Here, I think I’m only missing the t = \hat{t}/\tau conversion for dT, dC, dE.

5

Yes, written out in full it looks like an ODE since there isn’t anything defined implicitly. To switch between time you can just define that = t*tau, or the other way around if you let your input time be that. So if you add the *tau in the right spots it looks like you’re good?

(Adding a separate algebraic line that = t*tau in the macro version isn’t something that can be done yet it needs to be in the differential expressions, but it’s one of the things that will come when we do a change in the infrastructure there. If you use the function form then you can add extra lines of course)

6

Hi again Chris.
I think your intuition was indeed correct and this is a DAE. I’ve done away with the normalised unit case to simplify things and kind of have an implementation:

\begin{align} \dot{\hat{C}} &= \hat{E} - \frac{1}{\tau_C}(\hat{C}-C_{pi})\\ \dot{\hat{T}} &= \frac{1}{\tau}\left[\alpha \log\left(\frac{\hat{C}}{C_{pi}}\right)-\hat{T}\right]\\ \dot{\hat{W}} &= \hat{\chi}(\hat{T})\hat{W} = \xi(1-\delta\cdot\hat{T})\hat{W}\\ \dot{\hat{E}} &= [\hat{\chi}(\hat{T})-\mu]\hat{E}\\ \hat{E} &= \eta \exp(-\mu \hat{t})\hat{W}\\ \end{align}

Where the last equation there (for \hat{E}) is the constraint.

This should all transform fairly directly, I’ve included my constants here too:

function closed_system(out,du,u,p,t)
 χ = p[1]*(1-p[2]*u[2])
 out[1] = u[4] - (1/p[8])*(u[1]-p[7]) - du[1] #Chat
 out[2] = (1/p[5])*(p[6]*log(u[1]/p[7])-u[2]) - du[2] #That
 out[3] = χ*u[3] - du[3] #What 
 out[4] = (χ-p[4])*u[4] - du[4] # Ehat
 out[5] = p[3]*exp(-p[4]*t)*u[3] - u[4] # Ehat = η*exp(-μ*that)*What
end

δ = 0.5; #0.5 is equilibrium
τ = 50.0; #yr, characteristic climate time scale
τC = 50.0; #yr, characteristic carbon time scale (consistent with fixed airborne fraction)
η = 0.025; #ppmv/$trillion, initial carbon intensity (consistent with current CO₂ rise of 2ppmv/yr)
Cₚᵢ = 280.0; #ppmv, pre-industrial atmospheric carbon level
ΔT₂ₓ = 3.0; #K, equilibrium climate sensitivity (central estimate)
α = ΔT₂ₓ/log(2.0);

μ = 0.01; #1%/yr, decarbonization rate, based on historical data.
ξ = 0.04; #1% low growth, 4%, high growth
p = [ξ, δ, η, μ, τ, α, Cₚᵢ, τC]

#Initial CO₂ concentration (ppmv), Warming rate (K), wealth ($ trillion), Emissions (GtC)
u₀ = [380.0, 0.7, 160.0, η*160.0*2.14, 0.0]
du₀ = [0.0, 0.0, 0.0, 0.0, 0.0]
tspan = (2000.0,2200.0)
differential_vars = [true,true,true,true,false]
prob = DAEProblem(closed_system,du₀,u₀,tspan,p,differential_vars=differential_vars)
sol = solve(prob,IDA())

Running this throws an error that doesn’t seem to tell me much:

[IDAS ERROR] IDACalcIC
The residual routine or the linear setup or solve routine had a recoverable error, but IDACalcIC was unable to recover.

Is there something glaringly obvious here that I could alter?

7

Hey,
Yeah, Sundials was throwing error code -14 which has that very uninformative output. What’s going on is that IDA cannot handle overdetermined systems. Alan discusses that here:

http://sundials.2283335.n4.nabble.com/IDA-variable-choice-for-equation-td4653781.html#a4653791

One way you can get around that is by manually doing an index reduction routine. While we will be implementing these kinds of routines into ModelingToolkit.jl, it’s not ready. So for now you’d have to apply this algorithm by hand, which shouldn’t be too bad:

https://epubs.siam.org/doi/abs/10.1137/0914043

If you don’t need “full accuracy” and don’t need interpolations, you can make use of projections. You have 4 ODEs and a constraint added on. What you can do is solve the 4 ODE system and then utilize the ManifoldProjection callback.

http://docs.juliadiffeq.org/latest/features/callback_library.html#Manifold-Conservation-and-Projection-1

Only the steps (not the interpolation) are projected, but if you don’t need a continuous solution that should be fine. I would try using this in conjunction with Tsit5() or Rodas5() and see how it goes.

(P.S. Mathematica has some more discussions on overdetermined systems, and recommends pretty much the same thing: Numerical Solution of Differential-Algebraic Equations—Wolfram Language Documentation)

8

Thanks Chris!

The ManifoldProjection solution works enough for my requirements—this was only ever a poking around project.
I think in fact that I need to revisit my derivation, and in the end the solution may be simpler than what’s currently listed here. But even so, this has been a great introduction project on how Julia operates with ODEs for me.

Your work is truly inspiring & I hope we cross paths again in the future.