diff --git a/docs/pages.jl b/docs/pages.jl index f6c49a0de3..829c0d1c08 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -10,6 +10,7 @@ pages = [ "tutorials/stochastic_diffeq.md", "tutorials/discrete_system.md", "tutorials/parameter_identifiability.md", + "tutorials/change_independent_variable.md", "tutorials/bifurcation_diagram_computation.md", "tutorials/SampledData.md", "tutorials/domain_connections.md", diff --git a/docs/src/systems/ODESystem.md b/docs/src/systems/ODESystem.md index 6cc34725c4..4be0f6e263 100644 --- a/docs/src/systems/ODESystem.md +++ b/docs/src/systems/ODESystem.md @@ -30,6 +30,7 @@ ODESystem structural_simplify ode_order_lowering dae_index_lowering +change_independent_variable liouville_transform alias_elimination tearing diff --git a/docs/src/tutorials/change_independent_variable.md b/docs/src/tutorials/change_independent_variable.md new file mode 100644 index 0000000000..4b7ba037fe --- /dev/null +++ b/docs/src/tutorials/change_independent_variable.md @@ -0,0 +1,146 @@ +# Changing the independent variable of ODEs + +Ordinary differential equations describe the rate of change of some dependent variables with respect to one independent variable. +For the modeler it is often most natural to write down the equations with a particular independent variable, say time $t$. +However, in many cases there are good reasons for changing the independent variable: + + 1. One may want $y(x)$ as a function of $x$ instead of $(x(t), y(t))$ as a function of $t$ + 2. Some differential equations vary more nicely (e.g. less stiff) with respect to one independent variable than another. + 3. It can reduce the number of equations that must be solved (e.g. $y(x)$ is one equation, while $(x(t), y(t))$ are two). + +To manually change the independent variable of an ODE, one must rewrite all equations in terms of a new variable and transform differentials with the chain rule. +This is mechanical and error-prone. +ModelingToolkit provides the utility function [`change_independent_variable`](@ref) that automates this process. + +## 1. Get one dependent variable as a function of another + +Consider a projectile shot with some initial velocity in a vertical gravitational field with constant horizontal velocity. + +```@example changeivar +using ModelingToolkit +using ModelingToolkit: t_nounits as t, D_nounits as D +@variables x(t) y(t) +@parameters g=9.81 v # gravitational acceleration and horizontal velocity +eqs = [D(D(y)) ~ -g, D(x) ~ v] +initialization_eqs = [D(x) ~ D(y)] # 45° initial angle +M1 = ODESystem(eqs, t; initialization_eqs, name = :M) +M1s = structural_simplify(M1) +@assert length(equations(M1s)) == 3 # hide +M1s # hide +``` + +This is the standard parametrization that arises naturally from kinematics and Newton's laws. +It expresses the position $(x(t), y(t))$ as a function of time $t$. +But suppose we want to determine whether the projectile hits a target 10 meters away. +There are at least three ways of answering this: + + - Solve the ODE for $(x(t), y(t))$ and use a callback to terminate when $x$ reaches 10 meters, and evaluate $y$ at the final time. + - Solve the ODE for $(x(t), y(t))$ and use root finding to find the time when $x$ reaches 10 meters, and evaluate $y$ at that time. + - Solve the ODE for $y(x)$ and evaluate it at 10 meters. + +We will demonstrate the last method by changing the independent variable from $t$ to $x$. +This transformation is well-defined for any non-zero horizontal velocity $v$, so $x$ and $t$ are one-to-one. + +```@example changeivar +M2 = change_independent_variable(M1, x) +M2s = structural_simplify(M2; allow_symbolic = true) +# a sanity test on the 10 x/y variables that are accessible to the user # hide +@assert allequal([x, M1s.x]) # hide +@assert allequal([M2.x, M2s.x]) # hide +@assert allequal([y, M1s.y]) # hide +@assert allunique([M1.x, M1.y, M2.y, M2s.y]) # hide +@assert length(equations(M2s)) == 2 # hide +M2s # display this # hide +``` + +The derivatives are now with respect to the new independent variable $x$, which can be accessed with `M2.x`. + +!!! warn + + At this point `x`, `M1.x`, `M1s.x`, `M2.x`, `M2s.x` are *three* different variables. + Meanwhile `y`, `M1.y`, `M1s.y`, `M2.y` and `M2s.y` are *four* different variables. + It can be instructive to inspect these yourself to see their subtle differences. + +Notice how the number of equations has also decreased from three to two, as $\mathrm{d}x/\mathrm{d}t$ has been turned into an observed equation. +It is straightforward to evolve the ODE for 10 meters and plot the resulting trajectory $y(x)$: + +```@example changeivar +using OrdinaryDiffEq, Plots +prob = ODEProblem(M2s, [M2s.y => 0.0], [0.0, 10.0], [v => 8.0]) # throw 10 meters with x-velocity 8 m/s +sol = solve(prob, Tsit5()) +plot(sol; idxs = M2.y) # must index by M2.y = y(x); not M1.y = y(t)! +``` + +!!! tip "Usage tips" + + Look up the documentation of [`change_independent_variable`](@ref) for tips on how to use it. + + For example, if you also need $t(x)$, you can tell it to add a differential equation for the old independent variable in terms of the new one using the [inverse function rule](https://en.wikipedia.org/wiki/Inverse_function_rule) (here $\mathrm{d}t/\mathrm{d}x = 1 / (\mathrm{d}x/\mathrm{d}t)$). If you know an analytical expression between the independent variables (here $t = x/v$), you can also pass it directly to the function to avoid the extra differential equation. + +## 2. Alleviating stiffness by changing to logarithmic time + +In cosmology, the [Friedmann equations](https://en.wikipedia.org/wiki/Friedmann_equations) describe the expansion of the universe. +In terms of conformal time $t$, they can be written + +```@example changeivar +@variables a(t) Ω(t) +a = GlobalScope(a) # global var needed by all species +function species(w; kw...) + eqs = [D(Ω) ~ -3(1 + w) * D(a) / a * Ω] + return ODESystem(eqs, t, [Ω], []; kw...) +end +@named r = species(1 // 3) # radiation +@named m = species(0) # matter +@named Λ = species(-1) # dark energy / cosmological constant +eqs = [Ω ~ r.Ω + m.Ω + Λ.Ω, D(a) ~ √(Ω) * a^2] +initialization_eqs = [Λ.Ω + r.Ω + m.Ω ~ 1] +M1 = ODESystem(eqs, t, [Ω, a], []; initialization_eqs, name = :M) +M1 = compose(M1, r, m, Λ) +M1s = structural_simplify(M1) +``` + +Of course, we can solve this ODE as it is: + +```@example changeivar +prob = ODEProblem(M1s, [M1s.a => 1.0, M1s.r.Ω => 5e-5, M1s.m.Ω => 0.3], (0.0, -3.3), []) +sol = solve(prob, Tsit5(); reltol = 1e-5) +@assert Symbol(sol.retcode) == :Unstable # surrounding text assumes this was unstable # hide +plot(sol, idxs = [M1.a, M1.r.Ω / M1.Ω, M1.m.Ω / M1.Ω, M1.Λ.Ω / M1.Ω]) +``` + +But the solver becomes unstable due to stiffness. +Also notice the interesting dynamics taking place towards the end of the integration (in the early universe), which gets compressed into a very small time interval. +These ODEs would benefit from being defined with respect to a logarithmic "time" that better captures the evolution of the universe through *orders of magnitude* of time, rather than linear time. + +It is therefore common to write these ODEs in terms of $b = \ln a$. +To do this, we will change the independent variable in two stages; first from $t$ to $a$, and then from $a$ to $b$. +Notice that $\mathrm{d}a/\mathrm{d}t > 0$ provided that $\Omega > 0$, and $\mathrm{d}b/\mathrm{d}a > 0$, so the transformation is well-defined since $t \leftrightarrow a \leftrightarrow b$ are one-to-one. +First, we transform from $t$ to $a$: + +```@example changeivar +M2 = change_independent_variable(M1, M1.a) +@assert !ModelingToolkit.isautonomous(M2) # hide +M2 # hide +``` + +Unlike the original, notice that this system is *non-autonomous* because the independent variable $a$ appears explicitly in the equations! +This means that to change the independent variable from $a$ to $b$, we must provide not only the rate of change relation $db(a)/da = \exp(-b)$, but *also* the equation $a(b) = \exp(b)$ so $a$ can be eliminated in favor of $b$: + +```@example changeivar +a = M2.a # get independent variable of M2 +Da = Differential(a) +@variables b(a) +M3 = change_independent_variable(M2, b, [Da(b) ~ exp(-b), a ~ exp(b)]) +``` + +We can now solve and plot the ODE in terms of $b$: + +```@example changeivar +M3s = structural_simplify(M3; allow_symbolic = true) +prob = ODEProblem(M3s, [M3s.r.Ω => 5e-5, M3s.m.Ω => 0.3], (0, -15), []) +sol = solve(prob, Tsit5()) +@assert Symbol(sol.retcode) == :Success # surrounding text assumes the solution was successful # hide +plot(sol, idxs = [M3.r.Ω / M3.Ω, M3.m.Ω / M3.Ω, M3.Λ.Ω / M3.Ω]) +``` + +Notice that the variables vary "more nicely" with respect to $b$ than $t$, making the solver happier and avoiding numerical problems. diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl index e6b8130af3..41f9292c73 100644 --- a/src/ModelingToolkit.jl +++ b/src/ModelingToolkit.jl @@ -255,7 +255,8 @@ export isinput, isoutput, getbounds, hasbounds, getguess, hasguess, isdisturbanc tunable_parameters, isirreducible, getdescription, hasdescription, hasunit, getunit, hasconnect, getconnect, hasmisc, getmisc, state_priority -export ode_order_lowering, dae_order_lowering, liouville_transform +export ode_order_lowering, dae_order_lowering, liouville_transform, + change_independent_variable export PDESystem export Differential, expand_derivatives, @derivatives export Equation, ConstrainedEquation diff --git a/src/systems/diffeqs/basic_transformations.jl b/src/systems/diffeqs/basic_transformations.jl index e2be889c5e..a08c83ffb6 100644 --- a/src/systems/diffeqs/basic_transformations.jl +++ b/src/systems/diffeqs/basic_transformations.jl @@ -14,26 +14,20 @@ the final value of `trJ` is the probability of ``u(t)``. Example: ```julia -using ModelingToolkit, OrdinaryDiffEq, Test +using ModelingToolkit, OrdinaryDiffEq @independent_variables t @parameters α β γ δ @variables x(t) y(t) D = Differential(t) +eqs = [D(x) ~ α*x - β*x*y, D(y) ~ -δ*y + γ*x*y] +@named sys = ODESystem(eqs, t) -eqs = [D(x) ~ α*x - β*x*y, - D(y) ~ -δ*y + γ*x*y] - -sys = ODESystem(eqs) sys2 = liouville_transform(sys) -@variables trJ - -u0 = [x => 1.0, - y => 1.0, - trJ => 1.0] - -prob = ODEProblem(complete(sys2),u0,tspan,p) -sol = solve(prob,Tsit5()) +sys2 = complete(sys2) +u0 = [x => 1.0, y => 1.0, sys2.trJ => 1.0] +prob = ODEProblem(sys2, u0, tspan, p) +sol = solve(prob, Tsit5()) ``` Where `sol[3,:]` is the evolution of `trJ` over time. @@ -46,12 +40,159 @@ Optimal Transport Approach Abhishek Halder, Kooktae Lee, and Raktim Bhattacharya https://abhishekhalder.bitbucket.io/F16ACC2013Final.pdf """ -function liouville_transform(sys::AbstractODESystem) +function liouville_transform(sys::AbstractODESystem; kwargs...) t = get_iv(sys) @variables trJ D = ModelingToolkit.Differential(t) neweq = D(trJ) ~ trJ * -tr(calculate_jacobian(sys)) neweqs = [equations(sys); neweq] vars = [unknowns(sys); trJ] - ODESystem(neweqs, t, vars, parameters(sys), checks = false) + ODESystem( + neweqs, t, vars, parameters(sys); + checks = false, name = nameof(sys), kwargs... + ) +end + +""" + change_independent_variable( + sys::AbstractODESystem, iv, eqs = []; + add_old_diff = false, simplify = true, fold = false + ) + +Transform the independent variable (e.g. ``t``) of the ODE system `sys` to a dependent variable `iv` (e.g. ``u(t)``). +The transformation is well-defined when the mapping between the new and old independent variables are one-to-one. +This is satisfied if one is a strictly increasing function of the other (e.g. ``du(t)/dt > 0`` or ``du(t)/dt < 0``). + +Any extra equations `eqs` involving the new and old independent variables will be taken into account in the transformation. + +# Keyword arguments + +- `add_old_diff`: Whether to add a differential equation for the old independent variable in terms of the new one using the inverse function rule ``dt/du = 1/(du/dt)``. +- `simplify`: Whether expanded derivative expressions are simplified. This can give a tidier transformation. +- `fold`: Whether internal substitutions will evaluate numerical expressions. + +# Usage before structural simplification + +The variable change must take place before structural simplification. +In following calls to `structural_simplify`, consider passing `allow_symbolic = true` to avoid undesired constraint equations between between dummy variables. + +# Usage with non-autonomous systems + +If `sys` is non-autonomous (i.e. ``t`` appears explicitly in its equations), consider passing an algebraic equation relating the new and old independent variables (e.g. ``t = f(u(t))``). +Otherwise the transformed system can be underdetermined. +If an algebraic relation is not known, consider using `add_old_diff` instead. + +# Usage with hierarchical systems + +It is recommended that `iv` is a non-namespaced variable in `sys`. +This means it can belong to the top-level system or be a variable in a subsystem declared with `GlobalScope`. + +# Example + +Consider a free fall with constant horizontal velocity. +Physics naturally describes position as a function of time. +By changing the independent variable, it can be reformulated for vertical position as a function of horizontal position: +```julia +julia> @variables x(t) y(t); + +julia> @named M = ODESystem([D(D(y)) ~ -9.81, D(D(x)) ~ 0.0], t); + +julia> M = change_independent_variable(M, x); + +julia> M = structural_simplify(M; allow_symbolic = true); + +julia> unknowns(M) +3-element Vector{SymbolicUtils.BasicSymbolic{Real}}: + xˍt(x) + y(x) + yˍx(x) +``` +""" +function change_independent_variable( + sys::AbstractODESystem, iv, eqs = []; + add_old_diff = false, simplify = true, fold = false +) + iv2_of_iv1 = unwrap(iv) # e.g. u(t) + iv1 = get_iv(sys) # e.g. t + + if is_dde(sys) + error("System $(nameof(sys)) contains delay differential equations (DDEs). This is currently not supported!") + elseif isscheduled(sys) + error("System $(nameof(sys)) is structurally simplified. Change independent variable before structural simplification!") + elseif !iscall(iv2_of_iv1) || !isequal(only(arguments(iv2_of_iv1)), iv1) + error("Variable $iv is not a function of the independent variable $iv1 of the system $(nameof(sys))!") + end + + # Set up intermediate and final variables for the transformation + iv1name = nameof(iv1) # e.g. :t + iv2name = nameof(operation(iv2_of_iv1)) # e.g. :u + iv2, = @independent_variables $iv2name # e.g. u + iv1_of_iv2, = GlobalScope.(@variables $iv1name(iv2)) # inverse, e.g. t(u), global because iv1 has no namespacing in sys + D1 = Differential(iv1) # e.g. d/d(t) + div2_of_iv1 = GlobalScope(default_toterm(D1(iv2_of_iv1))) # e.g. uˍt(t) + div2_of_iv2 = substitute(div2_of_iv1, iv1 => iv2) # e.g. uˍt(u) + div2_of_iv2_of_iv1 = substitute(div2_of_iv2, iv2 => iv2_of_iv1) # e.g. uˍt(u(t)) + + # If requested, add a differential equation for the old independent variable as a function of the old one + if add_old_diff + eqs = [eqs; Differential(iv2)(iv1_of_iv2) ~ 1 / div2_of_iv2] # e.g. dt(u)/du ~ 1 / uˍt(u) (https://en.wikipedia.org/wiki/Inverse_function_rule) + end + @set! sys.eqs = [get_eqs(sys); eqs] # add extra equations we derived + @set! sys.unknowns = [get_unknowns(sys); [iv1, div2_of_iv1]] # add new variables, will be transformed to e.g. t(u) and uˍt(u) + + # Create a utility that performs the chain rule on an expression, followed by insertion of the new independent variable: + # e.g. (d/dt)(f(t)) -> (d/dt)(f(u(t))) -> df(u(t))/du(t) * du(t)/dt -> df(u)/du * uˍt(u) + function transform(ex) + # 1) Replace the argument of every function; e.g. f(t) -> f(u(t)) + for var in vars(ex; op = Nothing) # loop over all variables in expression (op = Nothing prevents interpreting "D(f(t))" as one big variable) + is_function_of_iv1 = iscall(var) && isequal(only(arguments(var)), iv1) # of the form f(t)? + if is_function_of_iv1 && !isequal(var, iv2_of_iv1) # prevent e.g. u(t) -> u(u(t)) + var_of_iv1 = var # e.g. f(t) + var_of_iv2_of_iv1 = substitute(var_of_iv1, iv1 => iv2_of_iv1) # e.g. f(u(t)) + ex = substitute(ex, var_of_iv1 => var_of_iv2_of_iv1; fold) + end + end + # 2) Repeatedly expand chain rule until nothing changes anymore + orgex = nothing + while !isequal(ex, orgex) + orgex = ex # save original + ex = expand_derivatives(ex, simplify) # expand chain rule, e.g. (d/dt)(f(u(t)))) -> df(u(t))/du(t) * du(t)/dt + ex = substitute(ex, D1(iv2_of_iv1) => div2_of_iv2_of_iv1; fold) # e.g. du(t)/dt -> uˍt(u(t)) + end + # 3) Set new independent variable + ex = substitute(ex, iv2_of_iv1 => iv2; fold) # set e.g. u(t) -> u everywhere + ex = substitute(ex, iv1 => iv1_of_iv2; fold) # set e.g. t -> t(u) everywhere + return ex + end + + # Use the utility function to transform everything in the system! + function transform(sys::AbstractODESystem) + eqs = map(transform, get_eqs(sys)) + unknowns = map(transform, get_unknowns(sys)) + unknowns = filter(var -> !isequal(var, iv2), unknowns) # remove e.g. u + ps = map(transform, get_ps(sys)) + ps = filter(!isinitial, ps) # remove Initial(...) # TODO: shouldn't have to touch this + observed = map(transform, get_observed(sys)) + initialization_eqs = map(transform, get_initialization_eqs(sys)) + parameter_dependencies = map(transform, get_parameter_dependencies(sys)) + defaults = Dict(transform(var) => transform(val) + for (var, val) in get_defaults(sys)) + guesses = Dict(transform(var) => transform(val) for (var, val) in get_guesses(sys)) + assertions = Dict(transform(ass) => msg for (ass, msg) in get_assertions(sys)) + systems = get_systems(sys) # save before reconstructing system + wascomplete = iscomplete(sys) # save before reconstructing system + sys = typeof(sys)( # recreate system with transformed fields + eqs, iv2, unknowns, ps; observed, initialization_eqs, + parameter_dependencies, defaults, guesses, + assertions, name = nameof(sys), description = description(sys) + ) + systems = map(transform, systems) # recurse through subsystems + sys = compose(sys, systems) # rebuild hierarchical system + if wascomplete + wasflat = isempty(systems) + sys = complete(sys; flatten = wasflat) # complete output if input was complete + end + return sys + end + return transform(sys) end diff --git a/test/basic_transformations.jl b/test/basic_transformations.jl index d9b59408e6..544a89cd29 100644 --- a/test/basic_transformations.jl +++ b/test/basic_transformations.jl @@ -1,34 +1,217 @@ -using ModelingToolkit, OrdinaryDiffEq, Test +using ModelingToolkit, OrdinaryDiffEq, DataInterpolations, Test @independent_variables t -@parameters α β γ δ -@variables x(t) y(t) D = Differential(t) -eqs = [D(x) ~ α * x - β * x * y, - D(y) ~ -δ * y + γ * x * y] +@testset "Liouville transform" begin + @parameters α β γ δ + @variables x(t) y(t) + eqs = [D(x) ~ α * x - β * x * y, D(y) ~ -δ * y + γ * x * y] + @named sys = ODESystem(eqs, t) + sys = complete(sys) -sys = ODESystem(eqs) + u0 = [x => 1.0, y => 1.0] + p = [α => 1.5, β => 1.0, δ => 3.0, γ => 1.0] + tspan = (0.0, 10.0) + prob = ODEProblem(sys, u0, tspan, p) + sol = solve(prob, Tsit5()) -u0 = [x => 1.0, - y => 1.0] + sys2 = liouville_transform(sys) + sys2 = complete(sys2) -p = [α => 1.5, - β => 1.0, - δ => 3.0, - γ => 1.0] + u0 = [x => 1.0, y => 1.0, sys2.trJ => 1.0] + prob = ODEProblem(sys2, u0, tspan, p, jac = true) + sol = solve(prob, Tsit5()) + @test sol[end, end] ≈ 1.0742818931017244 +end -tspan = (0.0, 10.0) -prob = ODEProblem(sys, u0, tspan, p) -sol = solve(prob, Tsit5()) +@testset "Change independent variable (trivial)" begin + @variables x(t) y(t) + eqs1 = [D(D(x)) ~ D(x) + x, D(y) ~ 1] + M1 = ODESystem(eqs1, t; name = :M) + M2 = change_independent_variable(M1, y) + @variables y x(y) yˍt(y) + Dy = Differential(y) + @test Set(equations(M2)) == Set([ + yˍt^2 * (Dy^2)(x) + yˍt * Dy(yˍt) * Dy(x) ~ x + Dy(x) * yˍt, + yˍt ~ 1 + ]) +end -sys2 = liouville_transform(sys) -@variables trJ +@testset "Change independent variable" begin + @variables x(t) y(t) z(t) s(t) + eqs = [ + D(x) ~ y, + D(D(y)) ~ 2 * x * D(y), + z ~ x + D(y), + D(s) ~ 1 / (2 * s) + ] + initialization_eqs = [x ~ 1.0, y ~ 1.0, D(y) ~ 0.0] + M1 = ODESystem(eqs, t; initialization_eqs, name = :M) + M2 = change_independent_variable(M1, s) -u0 = [x => 1.0, - y => 1.0, - trJ => 1.0] + M1 = structural_simplify(M1; allow_symbolic = true) + M2 = structural_simplify(M2; allow_symbolic = true) + prob1 = ODEProblem(M1, [M1.s => 1.0], (1.0, 4.0), []) + prob2 = ODEProblem(M2, [], (1.0, 2.0), []) + sol1 = solve(prob1, Tsit5(); reltol = 1e-10, abstol = 1e-10) + sol2 = solve(prob2, Tsit5(); reltol = 1e-10, abstol = 1e-10) + ts = range(0.0, 1.0, length = 50) + ss = .√(ts) + @test all(isapprox.(sol1(ts, idxs = M1.x), sol2(ss, idxs = M2.x); atol = 1e-7)) && + all(isapprox.(sol1(ts, idxs = M1.y), sol2(ss, idxs = M2.y); atol = 1e-7)) +end -prob = ODEProblem(sys2, u0, tspan, p, jac = true) -sol = solve(prob, Tsit5()) -@test sol[end, end] ≈ 1.0742818931017244 +@testset "Change independent variable (Friedmann equation)" begin + @independent_variables t + D = Differential(t) + @variables a(t) ȧ(t) Ω(t) ϕ(t) + a, ȧ = GlobalScope.([a, ȧ]) + species(w; kw...) = ODESystem([D(Ω) ~ -3(1 + w) * D(a) / a * Ω], t, [Ω], []; kw...) + @named r = species(1 // 3) + @named m = species(0) + @named Λ = species(-1) + eqs = [ + Ω ~ r.Ω + m.Ω + Λ.Ω, + D(a) ~ ȧ, + ȧ ~ √(Ω) * a^2, + D(D(ϕ)) ~ -3 * D(a) / a * D(ϕ) + ] + M1 = ODESystem(eqs, t, [Ω, a, ȧ, ϕ], []; name = :M) + M1 = compose(M1, r, m, Λ) + + # Apply in two steps, where derivatives are defined at each step: first t -> a, then a -> b + M2 = change_independent_variable(M1, M1.a) + M2c = complete(M2) # just for the following equation comparison (without namespacing) + a, ȧ, Ω, ϕ, aˍt = M2c.a, M2c.ȧ, M2c.Ω, M2c.ϕ, M2c.aˍt + Ωr, Ωm, ΩΛ = M2c.r.Ω, M2c.m.Ω, M2c.Λ.Ω + Da = Differential(a) + @test Set(equations(M2)) == Set([ + aˍt ~ ȧ, # dummy equation + Ω ~ Ωr + Ωm + ΩΛ, + ȧ ~ √(Ω) * a^2, + Da(aˍt) * Da(ϕ) * aˍt + aˍt^2 * (Da^2)(ϕ) ~ -3 * aˍt^2 / a * Da(ϕ), + aˍt * Da(Ωr) ~ -4 * Ωr * aˍt / a, + aˍt * Da(Ωm) ~ -3 * Ωm * aˍt / a, + aˍt * Da(ΩΛ) ~ 0 + ]) + + @variables b(M2.a) + extraeqs = [Differential(M2.a)(b) ~ exp(-b), M2.a ~ exp(b)] + M3 = change_independent_variable(M2, b, extraeqs) + + M1 = structural_simplify(M1) + M2 = structural_simplify(M2; allow_symbolic = true) + M3 = structural_simplify(M3; allow_symbolic = true) + @test length(unknowns(M3)) == length(unknowns(M2)) == length(unknowns(M1)) - 1 +end + +@testset "Change independent variable (simple)" begin + @variables x(t) y1(t) # y(t)[1:1] # TODO: use array variables y(t)[1:2] when fixed: https://github.com/JuliaSymbolics/Symbolics.jl/issues/1383 + Mt = ODESystem([D(x) ~ 2 * x, D(y1) ~ y1], t; name = :M) + Mx = change_independent_variable(Mt, x) + @variables x xˍt(x) xˍtt(x) y1(x) # y(x)[1:1] # TODO: array variables + Dx = Differential(x) + @test Set(equations(Mx)) == Set([xˍt ~ 2 * x, xˍt * Dx(y1) ~ y1]) +end + +@testset "Change independent variable (free fall with 1st order horizontal equation)" begin + @variables x(t) y(t) + @parameters g=9.81 v # gravitational acceleration and constant horizontal velocity + Mt = ODESystem([D(D(y)) ~ -g, D(x) ~ v], t; name = :M) # gives (x, y) as function of t, ... + Mx = change_independent_variable(Mt, x; add_old_diff = true) # ... but we want y as a function of x + Mx = structural_simplify(Mx; allow_symbolic = true) + Dx = Differential(Mx.x) + u0 = [Mx.y => 0.0, Dx(Mx.y) => 1.0, Mx.t => 0.0] + p = [v => 10.0] + prob = ODEProblem(Mx, u0, (0.0, 20.0), p) # 1 = dy/dx = (dy/dt)/(dx/dt) means equal initial horizontal and vertical velocities + sol = solve(prob, Tsit5(); reltol = 1e-5) + @test all(isapprox.(sol[Mx.y], sol[Mx.x - g * (Mx.t)^2 / 2]; atol = 1e-10)) # compare to analytical solution (x(t) = v*t, y(t) = v*t - g*t^2/2) +end + +@testset "Change independent variable (free fall with 2nd order horizontal equation)" begin + @variables x(t) y(t) + @parameters g = 9.81 # gravitational acceleration + Mt = ODESystem([D(D(y)) ~ -g, D(D(x)) ~ 0], t; name = :M) # gives (x, y) as function of t, ... + Mx = change_independent_variable(Mt, x; add_old_diff = true) # ... but we want y as a function of x + Mx = structural_simplify(Mx; allow_symbolic = true) + Dx = Differential(Mx.x) + u0 = [Mx.y => 0.0, Dx(Mx.y) => 1.0, Mx.t => 0.0, Mx.xˍt => 10.0] + prob = ODEProblem(Mx, u0, (0.0, 20.0), []) # 1 = dy/dx = (dy/dt)/(dx/dt) means equal initial horizontal and vertical velocities + sol = solve(prob, Tsit5(); reltol = 1e-5) + @test all(isapprox.(sol[Mx.y], sol[Mx.x - g * (Mx.t)^2 / 2]; atol = 1e-10)) # compare to analytical solution (x(t) = v*t, y(t) = v*t - g*t^2/2) +end + +@testset "Change independent variable (crazy 3rd order nonlinear system)" begin + @independent_variables t + D = Differential(t) + @variables x(t) y(t) + eqs = [ + (D^3)(y) ~ D(x)^2 + (D^2)(y^2) |> expand_derivatives, + D(x)^2 + D(y)^2 ~ x^4 + y^5 + t^6 + ] + M1 = ODESystem(eqs, t; name = :M) + M2 = change_independent_variable(M1, x; add_old_diff = true) + @test_nowarn structural_simplify(M2) + + # Compare to pen-and-paper result + @variables x xˍt(x) xˍt(x) y(x) t(x) + Dx = Differential(x) + areequivalent(eq1, eq2) = isequal(expand(eq1.lhs - eq2.lhs), 0) && + isequal(expand(eq1.rhs - eq2.rhs), 0) + eq1lhs = xˍt^3 * (Dx^3)(y) + xˍt^2 * Dx(y) * (Dx^2)(xˍt) + + xˍt * Dx(y) * (Dx(xˍt))^2 + + 3 * xˍt^2 * (Dx^2)(y) * Dx(xˍt) + eq1rhs = xˍt^2 + 2 * xˍt^2 * Dx(y)^2 + + 2 * xˍt^2 * y * (Dx^2)(y) + + 2 * y * Dx(y) * Dx(xˍt) * xˍt + eq1 = eq1lhs ~ eq1rhs + eq2 = xˍt^2 + xˍt^2 * Dx(y)^2 ~ x^4 + y^5 + t^6 + eq3 = Dx(t) ~ 1 / xˍt + @test areequivalent(equations(M2)[1], eq1) + @test areequivalent(equations(M2)[2], eq2) + @test areequivalent(equations(M2)[3], eq3) +end + +@testset "Change independent variable (registered function / callable parameter)" begin + @independent_variables t + D = Differential(t) + @variables x(t) y(t) + @parameters f::LinearInterpolation (fc::LinearInterpolation)(..) # non-callable and callable + callme(interp::LinearInterpolation, input) = interp(input) + @register_symbolic callme(interp::LinearInterpolation, input) + eqs = [ + D(x) ~ 2t, + D(y) ~ 1fc(t) + 2fc(x) + 3fc(y) + 1callme(f, t) + 2callme(f, x) + 3callme(f, y) + ] + M1 = ODESystem(eqs, t; name = :M) + + # Ensure that interpolations are called with the same variables + M2 = change_independent_variable(M1, x, [t ~ √(x)]) + @variables x xˍt(x) y(x) t(x) + Dx = Differential(x) + @test Set(equations(M2)) == Set([ + t ~ √(x), + xˍt ~ 2t, + xˍt * Dx(y) ~ 1fc(t) + 2fc(x) + 3fc(y) + + 1callme(f, t) + 2callme(f, x) + 3callme(f, y) + ]) + + _f = LinearInterpolation([1.0, 1.0], [-100.0, +100.0]) # constant value 1 + M2s = structural_simplify(M2; allow_symbolic = true) + prob = ODEProblem(M2s, [M2s.y => 0.0], (1.0, 4.0), [fc => _f, f => _f]) + sol = solve(prob, Tsit5(); abstol = 1e-5) + @test isapprox(sol(4.0, idxs = M2.y), 12.0; atol = 1e-5) # Anal solution is D(y) ~ 12 => y(t) ~ 12*t + C => y(x) ~ 12*√(x) + C. With y(x=1)=0 => 12*(√(x)-1), so y(x=4) ~ 12 +end + +@testset "Change independent variable (errors)" begin + @variables x(t) y z(y) w(t) v(t) + M = ODESystem([D(x) ~ 1, v ~ x], t; name = :M) + Ms = structural_simplify(M) + @test_throws "structurally simplified" change_independent_variable(Ms, y) + @test_throws "not a function of" change_independent_variable(M, y) + @test_throws "not a function of" change_independent_variable(M, z) + @variables x(..) # require explicit argument + M = ODESystem([D(x(t)) ~ x(t - 1)], t; name = :M) + @test_throws "DDE" change_independent_variable(M, x(t)) +end diff --git a/test/runtests.jl b/test/runtests.jl index 0ae66d3755..8f1f2a0be4 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -46,6 +46,7 @@ end @safetestset "Model Parsing Test" include("model_parsing.jl") @safetestset "Error Handling" include("error_handling.jl") @safetestset "StructuralTransformations" include("structural_transformation/runtests.jl") + @safetestset "Basic transformations" include("basic_transformations.jl") @safetestset "State Selection Test" include("state_selection.jl") @safetestset "Symbolic Event Test" include("symbolic_events.jl") @safetestset "Stream Connect Test" include("stream_connectors.jl")