Working with a single equation, split up into several

General Usage
1

I’m trying to get into ModelingToolkit.jl but having some trouble understanding the documentation, which doesn’t seem to include basic examples. Solving nonlinear ODEs and plotting solutions in 3D is very impressive, but at the moment the framework still feels very hard to get into, especially if one doesn’t need differential equations.

Let’s say I have a single non-linear equation for a quantity z:

using ModelingToolkit
using NonlinearSolve

@variables u v w x y z
@parameters a b

eqn = z ~ a * u^2 + b * (v + 2w + 3x + 4y)

Question 1: How can I best plug in numbers for variables and get a number as a result? (e.g. all parameters and variables except z)

Note: you might object that ModelingToolkit.jl is not responsible for such trivial things. However, in my opinion, when working with models of any kind, it is essential to be able to do simple things, such as plugging in numbers, plotting equations as functions, etc.

The best solution I came up with after quite a bit of trying seems a bit cumbersome and is not part of ModelingToolkit:

using Symbolics
substitute(eqn, Dict([u=>1., v=>2., w=>3., x=>4., y=>5., a=>6., b=>7.]))
# Result: `z ~ 286.0`. Still need to extract the number, too.

My actual equation is quite long. Publications split it up into several terms and give them names. Some of the terms appear multiple times as sub-expressions of other terms. I would like to be able to present the equations in this form before manipulating them, so readers can verify they are copied from the paper. For example, the same equation above could be written as four equations using helper variables h_1, h_2,h_3 as:

@variables h_1, h_2, h_3

eqns = [
  z ~ h_1 + b * (h_2 + h_3)
h_1 ~ a * u^2
h_2 ~ v + 2w
h_3 ~ 3x + 4y
]

Question 2: How can I get back the single equation by substituting them into one another? (Or, if that’s not possible, at least verify automatically that the equation above is the same as this system?)

I tried turning the equations into a NonlinearSystem and using structural_simplify, but it doesn’t work (“ExtraVariablesSystemException: The system is unbalanced…”). The functions (e.g. simplify or solve_for) in the Symbolics package also didn’t work.

Possibly, I could achieve a similar effect by using functions instead of variables or Equations? What is the right way to approach this?

2

I kind-of solved the problem with ugly hacks. If anyone is interested, here is my haphazard solution, which just calls Symbolics.substitute until all the helper variables are gone.

using Symbolics

"""
Substitutes expressions of `target_eqn` as specified by equations `source_eqns`.
For each equation in `source_eqns`, substitutes its `lhs` by its `rhs` until none of the `lhs` are left in `target_eqn`.

TODO this is just a crappy prototype!
"""
function substitute_all(target_eqn, source_eqns,)
    all_eqns = [target_eqn, source_eqns...]

    eliminated_symbols = Set(map(x->x.lhs, source_eqns))
    length(eliminated_symbols) == length(source_eqns) || error("Redundant equations!")

    for i in 2:length(all_eqns)
        all_eqns = substitute.(all_eqns, (Dict(all_eqns[i].lhs => all_eqns[i].rhs),))
    end

    return all_eqns[1]
end


# test the function using above example

@variables u v w x y z a b

eqn_to_recover = z ~ a * u^2 + b * (v + 2w + 3x + 4y)

# define split-up version
helper_variables = @variables h_1, h_2, h_3

target_eqn = z ~ h_1 + b * (h_2 + h_3)

source_eqns = [
    h_1 ~ a * u^2
    h_2 ~ v + 2w
    h_3 ~ 3x + 4y
]

substitute_all(target_eqn, source_eqns,) == eqn_to_recover || error("Failure")