add the liouville transformation

Author: ChrisRackauckas

docs/src/systems/ODESystem.md

@@ -17,6 +17,7 @@ ODESystem
 
 ```@docs
 ode_order_lowering
+liouville_transform
 ```
 
 ## Applicable Calculation and Generation Functions

src/ModelingToolkit.jl

@@ -176,6 +176,7 @@ include("systems/diffeqs/abstractodesystem.jl")
 include("systems/diffeqs/first_order_transform.jl")
 include("systems/diffeqs/modelingtoolkitize.jl")
 include("systems/diffeqs/validation.jl")
+include("systems/diffeqs/basic_transformations.jl")
 
 include("systems/jumps/jumpsystem.jl")
 

src/systems/diffeqs/basic_transformations.jl

@@ -0,0 +1,17 @@
+"""
+$(TYPEDSIGNATURES)
+
+Generates the Liouville transformed set of ODEs, which is the original
+ODE system with a new variable `trJ` appended, corresponding to the
+-tr(Jacobian). This variable is used for properties like uncertainty
+propagation and should be used with a zero initial condition.
+"""
+function liouville_transform(sys)
+      t = independent_variable(sys)
+      @variables trJ
+      D = ModelingToolkit.Differential(t)
+      neweq = D(trJ) ~ -tr(calculate_jacobian(sys))
+      neweqs = [equations(sys);neweq]
+      vars = [states(sys);trJ]
+      ODESystem(neweqs,t,vars,parameters(sys))
+end

test/basic_transformations.jl

@@ -0,0 +1,38 @@
+using ModelingToolkit, OrdinaryDiffEq, Test
+
+@parameters t σ ρ β
+@variables x(t) y(t) z(t)
+@derivatives D'~t
+
+eqs = [D(D(x)) ~ σ*(y-x),
+       D(y) ~ x*(ρ-z)-y,
+       D(z) ~ x*y - β*z]
+
+sys = ODESystem(eqs)
+sys = ode_order_lowering(sys)
+
+u0 = [D(x) => 2.0,
+      x => 1.0,
+      y => 0.0,
+      z => 0.0]
+
+p  = [σ => 28.0,
+      ρ => 10.0,
+      β => 8/3]
+
+tspan = (0.0,100.0)
+prob = ODEProblem(sys,u0,tspan,p,jac=true)
+sol = solve(prob,Tsit5())
+
+sys2 = louiville_transform(sys)
+@variables trJ
+
+u0 = [D(x) => 2.0,
+      x => 1.0,
+      y => 0.0,
+      z => 0.0,
+      trJ => 0.0]
+
+prob = ODEProblem(sys2,u0,tspan,p,jac=true)
+sol = solve(prob,Tsit5())
+@test sol[end,end] ≈ 366.66666666