How to read a system of equations?

Specific Domains Modelling & Simulations
3

One way to define the dynamics is as follows:

const c1 = 1.0
const c2 = 2.0
const c3 = 3.0
const c4 = 4.0
const c5 = 5.0
const c6 = 6.0

function dynamics(Y)
    x, y = Y
    [c1 + c2 * x + c3 * y,
     c4 + c5 * x + c6 * y]
end

So that you can evaluate the dynamics as:

julia> Y = [2.0, 3.0];

julia> dynamics(Y)
2-element Array{Float64,1}:
 14.0
 32.0

Of course you can find the constant offset by simply evaluating the function at Y = 0:

julia> A = dynamics(zeros(2))
2-element Array{Float64,1}:
 1.0
 4.0

Since this system is linear, B = \frac{\partial Y(t)}{\partial Y(t - 1)}. So you could use (automatic) differentiation to extract B. One option is to use ForwardDiff.jl (install using Pkg.add("ForwardDiff")) to extract B:

julia> using ForwardDiff

julia> B = ForwardDiff.jacobian(dynamics, Y)
2×2 Array{Float64,2}:
 2.0  3.0
 5.0  6.0

Here, you’re essentially implicitly telling ForwardDiff that Y contains the variables by having Y be the argument of the dynamics function while the c’s are global constants. Note that ForwardDiff.jacobian expects as function that maps a vector to a vector. If you don’t want the c’s to be global constants, you can use a ‘closure’:

function dynamics2(Y, C)
    x, y = Y
    [C[1] + C[2] * x + C[3] * y,
     C[4] + C[5] * x + C[6] * y]
end

A = ForwardDiff.jacobian(Y -> dynamics2(Y, [1.0, 2.0, 3.0, 4.0, 5.0, 6.0]), Y)

By the way, Discourse allows you to write math equations in LaTeX form using e.g. $x = \frac{y}{z}$, and see this post for how to get nice code formatting.

4 Likes