# How to read a system of equations?

**URL:** <https://discourse.julialang.org/t/how-to-read-a-system-of-equations/10207>\
**Category:** Modelling & Simulations\
**Tags:** question\
**Created:** [April 6, 2018, 9:33pm UTC](https://discourse.julialang.org/t/how-to-read-a-system-of-equations/10207 "2018-04-06T21:33:46Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![tomj](https://avatars.discourse-cdn.com/v4/letter/t/dfb087/32.png) [@tomj](https://discourse.julialang.org/u/tomj)\
**Post date:** [April 6, 2018, 9:33pm UTC](https://discourse.julialang.org/t/how-to-read-a-system-of-equations/10207/1 "2018-04-06T21:33:46Z")

</div>

Hello Julia users!

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

Imagine a system of equations like

```julia
 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

---

<div class="post-metadata">

**Author:** ![Nosferican](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nosferican/32/9275_2.png) [@Nosferican](https://discourse.julialang.org/u/Nosferican)\
**Post date:** [April 7, 2018, 2:38am UTC](https://discourse.julialang.org/t/how-to-read-a-system-of-equations/10207/2 "2018-04-07T02:38:36Z")

</div>

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,

```julia
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 `Array`s 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).

---

<div class="post-metadata">

**Author:** ![tkoolen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tkoolen/32/1603_2.png) [@tkoolen](https://discourse.julialang.org/u/tkoolen)\
**Post date:** [April 7, 2018, 3:17am UTC](https://discourse.julialang.org/t/how-to-read-a-system-of-equations/10207/3 "2018-04-07T03:17:41Z")

</div>

One way to define the dynamics is as follows:

```julia
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
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
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](https://github.com/JuliaDiff/ForwardDiff.jl) (install using `Pkg.add("ForwardDiff")`) to extract `B`:

```julia
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`](http://www.juliadiff.org/ForwardDiff.jl/stable/user/api.html#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’:

```julia
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](https://discourse.julialang.org/t/psa-how-to-quote-code-with-backticks/7530) for how to get nice code formatting.

---

<div class="post-metadata">

**Author:** ![tomj](https://avatars.discourse-cdn.com/v4/letter/t/dfb087/32.png) [@tomj](https://discourse.julialang.org/u/tomj)\
**Post date:** [April 7, 2018, 12:19pm UTC](https://discourse.julialang.org/t/how-to-read-a-system-of-equations/10207/4 "2018-04-07T12:19:01Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [April 7, 2018, 12:41pm UTC](https://discourse.julialang.org/t/how-to-read-a-system-of-equations/10207/5 "2018-04-07T12:41:14Z")

</div>

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

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

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [April 7, 2018, 1:13pm UTC](https://discourse.julialang.org/t/how-to-read-a-system-of-equations/10207/6 "2018-04-07T13:13:13Z")

</div>

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](http://www.dynare.org/), 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.