From cee94d52e5012fc8c8fb49644464e72f5c12a0c5 Mon Sep 17 00:00:00 2001 From: Steffen Plunder Date: Wed, 23 Jun 2021 03:48:11 +0200 Subject: [PATCH 1/5] work in progress --- src/ModelingToolkit.jl | 2 +- src/systems/diffeqs/basic_transformations.jl | 82 ++++++++++++++++++ test/changeofvariables.jl | 88 ++++++++++++++++++++ 3 files changed, 171 insertions(+), 1 deletion(-) create mode 100644 test/changeofvariables.jl diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl index 261890d40c..a2f447685b 100644 --- a/src/ModelingToolkit.jl +++ b/src/ModelingToolkit.jl @@ -153,7 +153,7 @@ export JumpProblem, DiscreteProblem export NonlinearSystem, OptimizationSystem export ControlSystem export alias_elimination, flatten, connect, @connector -export ode_order_lowering, liouville_transform +export ode_order_lowering, liouville_transform, changeofvariables export runge_kutta_discretize export PDESystem export Reaction, ReactionSystem, ismassaction, oderatelaw, jumpratelaw diff --git a/src/systems/diffeqs/basic_transformations.jl b/src/systems/diffeqs/basic_transformations.jl index 1478a968b8..8af3ddc041 100644 --- a/src/systems/diffeqs/basic_transformations.jl +++ b/src/systems/diffeqs/basic_transformations.jl @@ -54,3 +54,85 @@ function liouville_transform(sys) vars = [states(sys);trJ] ODESystem(neweqs,t,vars,parameters(sys)) end + + + + + + + + +""" +$(TYPEDSIGNATURES) + +Generates the set of ODEs after change of variables. + + +Example: + +```julia +using ModelingToolkit, OrdinaryDiffEq, Test + +@parameters t α +@variables x(t) +D = Differential(t) +eqs = [D(x) ~ α*x] + +tspan = (0., 1.) +u0 = [x => 1.0] +p = [α => -0.5] + +sys = ODESystem(eqs; defaults=u0) +prob = ODEProblem(sys, [], tspan, p) +sol = solve(prob, Tsit5()) + +@variables z(t) +forward_subs = [exp(x) => z] +backward_subs = [x => log(z)] +new_sys = changeofvariables(sys, forward_subs, backward_subs) +@test equations(new_sys)[1] == (D(z) ~ α*z*log(z)) + +new_prob = ODEProblem(new_sys, [], tspan, p) +new_sol = solve(new_prob, Tsit5()) + +@test isapprox(new_sol[x][end], sol[x][end], atol=1e-4) +``` + +""" +function changeofvariables(sys, forward_subs, backward_subs; simplify=false, t0=missing) + t = independent_variable(sys) + + old_vars = first.(backward_subs) + new_vars = last.(forward_subs) + kept_vars = setdiff(states(sys), old_vars) + rhs = [eq.rhs for eq in equations(sys)] + + # use: dz/dt = ∂z/∂x dx/dt + ∂z/∂t + dzdt = Symbolics.derivative( first.(forward_subs), t ) + new_eqs = [] + for (new_var, ex) in zip(new_vars, dzdt) + for ode_eq in equations(sys) + ex = substitute(ex, ode_eq.lhs => ode_eq.rhs) + end + ex = substitute(ex, Dict(forward_subs)) + ex = substitute(ex, Dict(backward_subs)) + if simplify + ex = Symbolics.simplify(ex, expand=true) + end + push!(new_eqs, Differential(t)(new_var) ~ ex) + end + + defs = get_defaults(sys) + new_defs = Dict() + for f_sub in forward_subs + #TODO call value(...)? + ex = substitute(first(f_sub), defs) + if !ismissing(t0) + ex = substitute(ex, t => t0) + end + new_defs[last(f_sub)] = ex + end + return ODESystem(new_eqs; + defaults=new_defs, + observed=first.(backward_subs) .~ last.(backward_subs)) +end diff --git a/test/changeofvariables.jl b/test/changeofvariables.jl new file mode 100644 index 0000000000..2cb3e747f6 --- /dev/null +++ b/test/changeofvariables.jl @@ -0,0 +1,88 @@ +using ModelingToolkit, OrdinaryDiffEq +using Test, LinearAlgebra + + +# Change of variables: z = exp(x) + +@parameters t α +@variables x(t) +D = Differential(t) +eqs = [D(x) ~ α*x] + +tspan = (0., 1.) +u0 = [x => 1.0] +p = [α => -0.5] + +sys = ODESystem(eqs; defaults=u0) +prob = ODEProblem(sys, [], tspan, p) +sol = solve(prob, Tsit5()) + +@variables z(t) +forward_subs = [exp(x) => z] +backward_subs = [x => log(z)] +new_sys = changeofvariables(sys, forward_subs, backward_subs) +@test equations(new_sys)[1] == (D(z) ~ α*z*log(z)) + +new_prob = ODEProblem(new_sys, [], tspan, p) +new_sol = solve(new_prob, Tsit5()) + +@test isapprox(new_sol[x][end], sol[x][end], atol=1e-4) + + + +# Riccatti equation +@parameters t α +@variables x(t) +D = Differential(t) +eqs = [D(x) ~ t^2 + α - x^2] +sys = ODESystem(eqs, defaults=[x=>1.]) + +@variables z(t) +forward_subs = [t + α/(x+t) => z ] +backward_subs = [ x => α/(z-t) - t] + +new_sys = changeofvariables(sys, forward_subs, backward_subs; simplify=true, t0=0.) +# output should be equivalent to +# t^2 + α - z^2 + 2 (but this simplification is not found automatically) + +tspan = (0., 1.) +p = [α => 1.] +prob = ODEProblem(sys,[],tspan,p) +new_prob = ODEProblem(new_sys,[],tspan,p) + +sol = solve(prob, Tsit5()) +new_sol = solve(new_prob, Tsit5()) + +@test isapprox(sol[x][end], new_sol[x][end], rtol=1e-4) + + +# Linear transformation to diagonal system +@parameters t +@variables x[1:3](t) +D = Differential(t) +A = [0. -1. 0.; -0.5 0.5 0.; 0. 0. -1.] +eqs = D.(x) .~ A*x + +tspan = (0., 10.) +u0 = x .=> [1.0, 2.0, -1.0] + +sys = ODESystem(eqs; defaults=u0) +prob = ODEProblem(sys,[],tspan) +sol = solve(prob, Tsit5()) + +T = eigen(A).vectors + +@variables z[1:3](t) +forward_subs = T \ x .=> z +backward_subs = x .=> T*z + +new_sys = changeofvariables(sys, forward_subs, backward_subs; simplify=true) + +new_prob = ODEProblem(new_sys, [], tspan, p) +new_sol = solve(new_prob, Tsit5()) + +# test RHS +new_rhs = [eq.rhs for eq in equations(new_sys)] +new_A = Symbolics.value.(Symbolics.jacobian(new_rhs, z)) +@test isapprox(diagm(eigen(A).values), new_A, rtol = 1e-10) +@test isapprox( new_sol[x[1],end], sol[x[1],end], rtol=1e-4) From c3037b6dc79ad8da76620f7206f1368b7e35c1ac Mon Sep 17 00:00:00 2001 From: Steffen Plunder Date: Thu, 24 Jun 2021 17:56:52 +0200 Subject: [PATCH 2/5] added correct substitution for text example --- test/changeofvariables.jl | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/test/changeofvariables.jl b/test/changeofvariables.jl index 2cb3e747f6..75fba44a6f 100644 --- a/test/changeofvariables.jl +++ b/test/changeofvariables.jl @@ -2,7 +2,8 @@ using ModelingToolkit, OrdinaryDiffEq using Test, LinearAlgebra -# Change of variables: z = exp(x) +# Change of variables: z = log(x) +# (this implies that x = exp(z) is automatically non-negative) @parameters t α @variables x(t) @@ -18,10 +19,10 @@ prob = ODEProblem(sys, [], tspan, p) sol = solve(prob, Tsit5()) @variables z(t) -forward_subs = [exp(x) => z] -backward_subs = [x => log(z)] +forward_subs = [log(x) => z] +backward_subs = [x => exp(z)] new_sys = changeofvariables(sys, forward_subs, backward_subs) -@test equations(new_sys)[1] == (D(z) ~ α*z*log(z)) +@test equations(new_sys)[1] == (D(z) ~ α) new_prob = ODEProblem(new_sys, [], tspan, p) new_sol = solve(new_prob, Tsit5()) @@ -30,7 +31,7 @@ new_sol = solve(new_prob, Tsit5()) -# Riccatti equation +# Riccati equation @parameters t α @variables x(t) D = Differential(t) From 72a8b3a0d242d3ddadc23177f4137e63f42d441e Mon Sep 17 00:00:00 2001 From: Steffen Plunder Date: Thu, 24 Jun 2021 17:58:47 +0200 Subject: [PATCH 3/5] keep observed variables from original system --- src/systems/diffeqs/basic_transformations.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/systems/diffeqs/basic_transformations.jl b/src/systems/diffeqs/basic_transformations.jl index 8af3ddc041..defa43a524 100644 --- a/src/systems/diffeqs/basic_transformations.jl +++ b/src/systems/diffeqs/basic_transformations.jl @@ -99,7 +99,7 @@ new_sol = solve(new_prob, Tsit5()) ``` """ -function changeofvariables(sys, forward_subs, backward_subs; simplify=false, t0=missing) +function changeofvariables(sys::ODESystem, forward_subs, backward_subs; simplify=false, t0=missing) t = independent_variable(sys) old_vars = first.(backward_subs) @@ -109,7 +109,7 @@ function changeofvariables(sys, forward_subs, backward_subs; simplify=false, t0= # use: dz/dt = ∂z/∂x dx/dt + ∂z/∂t dzdt = Symbolics.derivative( first.(forward_subs), t ) - new_eqs = [] + new_eqs = Equation[] for (new_var, ex) in zip(new_vars, dzdt) for ode_eq in equations(sys) ex = substitute(ex, ode_eq.lhs => ode_eq.rhs) @@ -134,5 +134,6 @@ function changeofvariables(sys, forward_subs, backward_subs; simplify=false, t0= end return ODESystem(new_eqs; defaults=new_defs, - observed=first.(backward_subs) .~ last.(backward_subs)) + observed=vcat(observed(sys),first.(backward_subs) .~ last.(backward_subs)) + ) end From 4a724125ed2c31eb39339419a8d2afa8dc64e34a Mon Sep 17 00:00:00 2001 From: Steffen Plunder Date: Fri, 3 Sep 2021 18:41:44 +0200 Subject: [PATCH 4/5] fixed doc --- src/systems/diffeqs/basic_transformations.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/systems/diffeqs/basic_transformations.jl b/src/systems/diffeqs/basic_transformations.jl index 9da98e185f..5d476f2eb1 100644 --- a/src/systems/diffeqs/basic_transformations.jl +++ b/src/systems/diffeqs/basic_transformations.jl @@ -90,7 +90,7 @@ sol = solve(prob, Tsit5()) forward_subs = [exp(x) => z] backward_subs = [x => log(z)] new_sys = changeofvariables(sys, forward_subs, backward_subs) -@test equations(new_sys)[1] == (D(z) ~ α*z*log(z)) +@test equations(new_sys)[1] == (D(z) ~ α) new_prob = ODEProblem(new_sys, [], tspan, p) new_sol = solve(new_prob, Tsit5()) From 4947f49fb4baedc422825b949e723b48178077bc Mon Sep 17 00:00:00 2001 From: Steffen Plunder Date: Fri, 3 Sep 2021 18:43:46 +0200 Subject: [PATCH 5/5] fix doc --- src/systems/diffeqs/basic_transformations.jl | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/src/systems/diffeqs/basic_transformations.jl b/src/systems/diffeqs/basic_transformations.jl index 5d476f2eb1..d2e1eb8a0a 100644 --- a/src/systems/diffeqs/basic_transformations.jl +++ b/src/systems/diffeqs/basic_transformations.jl @@ -73,6 +73,9 @@ Example: ```julia using ModelingToolkit, OrdinaryDiffEq, Test +# Change of variables: z = log(x) +# (this implies that x = exp(z) is automatically non-negative) + @parameters t α @variables x(t) D = Differential(t) @@ -82,14 +85,15 @@ tspan = (0., 1.) u0 = [x => 1.0] p = [α => -0.5] -sys = ODESystem(eqs; defaults=u0) +@named sys = ODESystem(eqs; defaults=u0) prob = ODEProblem(sys, [], tspan, p) sol = solve(prob, Tsit5()) @variables z(t) -forward_subs = [exp(x) => z] -backward_subs = [x => log(z)] -new_sys = changeofvariables(sys, forward_subs, backward_subs) +forward_subs = [log(x) => z] +backward_subs = [x => exp(z)] + +@named new_sys = changeofvariables(sys, forward_subs, backward_subs) @test equations(new_sys)[1] == (D(z) ~ α) new_prob = ODEProblem(new_sys, [], tspan, p)