How to read a system of equations?

Specific Domains Modelling & Simulations
1

Hello Julia users!

I’m a beginner in Julia, and am currently blocked with the following problem.

Imagine a system of equations like

 x(t) = c1 + c2*x(t-1) + c3*y(t-1)
 y(t) = c4 + c5*y(t-1) + c6*x(t-1)

I would like to write those equations in a script for Julia, and next code an operation that would transform this system of equation into matrix form:
Y(t) = A + B*Y(t-1), where Y = [x y]' with A = [c1 c4]' and B = [c2 c3;c5 c6]

My questions:

  1. Is there a way to tell Julia that x and y are variables, and all the c’s coefficients?
  2. How can I make Julia read the system, identify variables and coefficients in order to construct A and B matrices?

Thanks a lot in advance!

Thomas

2

Please backtick code to make it easier to read. Without knowing more about your needs is hard to assist. I made the following assumptions: linear system of equations. Here is a toy example,

macro parse(ex::Expr)
    findvals(ex)
end
function findvals(ex::Expr)
    args = ex.args
    if args[2] isa Number
        return float.(push!(args[2], args[1]))
    else
        return float.(push!(args[2].args[2:end], args[1]))
    end
end
function update(A::AbstractMatrix{<:Number},
                b::AbstractVector{<:Number},
        xy::AbstractVector{<:Number})
    return vcat(A, xy[2:end]'), push!(copy(b), xy[1])
end
A = Matrix{Float64}(0, 3)
b = Vector{Float64}(0)
A, b = update(A, b, @parse(3 = 1 + 2 + -4))
A, b = update(A, b, @parse(2 = 0 + -2 + 3))

I assume the equation are ordered, but could potentially could be parse with various variables through value * variable. You could create the rows one by one or if you know ahead of time allocate the Arrays for efficiency. The toy model doesn’t have full safeguards (e.g., notice that 2 = 1 + 2 - 4 will not work while 2 = 1 + 2 + -4 will).

1 Like
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
4

Hello!

Thanks for your quick answers! I guess I was a little too vague about what I wanted to do, here are more details (I don’t have any proper code yet, only pieces like what follows):

I start from the system of 2 equations, given as 2 ::Expr
eq1 = “x = alpha + beta * x(-1) + beta/gamma * y(-1)”
eq2 = “y = sigma + rho * y(-1) + gamma*x(-1)”

where x(-1) indicates a lag for the variable x. For the moment there is not value assigned to coefficients, they could be Symbol for instance.

In a first step, I would like Julia to construct coefficient matrices such that
A = [alpha; sigma]
B = [beta beta/gamma; gamma rho]
Here I don’t know how to do this. “parse” and “.args” help to separate the terms, but how to tell Julia that beta is associated with x(-1), beta/gamma is the coefficient of y(-1) in the equation for x, etc… ?

Next, in a second step, I would assign values to the symbols beta, gamma etc… in order to evaluate the matrices. This step is straightforward.

Best,

Thomas

5

In general, I wouldn’t recommend inputting equations as strings that you have to parse. Instead just pass them as a function argument.

1 Like
6

I think that @tkoolen gave a really useful answer to your questions — it addresses everything, you may want to study it a bit.

For some toolkits, particularly Matlab-based ones, it is common to provide an interface where the user enters strings which are then parsed as you suggest. Eg Dynare, a (mainly) Matlab toolkit for DSGE models has a very similar syntax.

In Julia, you may be better off coding this directly as a function, as suggested above. You will be able to test and optimize it directly (for more complex cases, of course), and it will be just compiled as is as a first-class object without further steps.

2 Likes