# Solving a Petersen Matrix with Julia

**URL:** <https://discourse.julialang.org/t/solving-a-petersen-matrix-with-julia/67674>\
**Category:** New to Julia\
**Tags:** diffeq\
**Created:** [September 4, 2021, 2:35am UTC](https://discourse.julialang.org/t/solving-a-petersen-matrix-with-julia/67674 "2021-09-04T02:35:43Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![Kester](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kester/32/28870_2.png) [@Kester](https://discourse.julialang.org/u/Kester)\
**Post date:** [September 4, 2021, 2:35am UTC](https://discourse.julialang.org/t/solving-a-petersen-matrix-with-julia/67674/1 "2021-09-04T02:35:43Z")

</div>

Julia has been great so far for my work in environmental engineering and modeling. I have only used it for simple systems thus far, but it is very performant.

Maybe I am missing something simple, but I have not found a good succinct way to input system rate equations in the form of a Petersen Matrix:  
[https://en.wikipedia.org/wiki/Petersen\_matrix](https://en.wikipedia.org/wiki/Petersen_matrix)

As in the example from Wikipedia, I would normally write out the system of equations and solve

As the system gets more complex I would like to solve in matrix format:  
 ![image](https://global.discourse-cdn.com/julialang/original/3X/9/d/9d3ff8944136022d3ec11190e5a903d04b76f912.png)

But I am not sure how to do the syntax. Thanks!

---

<div class="post-metadata">

**Author:** ![Erasmo98](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/erasmo98/32/26098_2.png) [@Erasmo98](https://discourse.julialang.org/u/Erasmo98)\
**Post date:** [September 4, 2021, 3:55am UTC](https://discourse.julialang.org/t/solving-a-petersen-matrix-with-julia/67674/2 "2021-09-04T03:55:47Z")

</div>

Welcome Kester, I hope you have a great time in discourse.

Well to solve the equation you would need to use some package like [`DifferentialEquations.jl`](https://diffeq.sciml.ai/stable/) which I haven’t use myself.

But, lets say that you a Julia function named `solve` which solves the equation

\frac{\text{d}}{\text{d}t}\vec{x} = f(\vec{x}),

then in your case you have that

\vec{x} = ([A], [B], [S], [E], [ES], [P])^T.

Enough maths, lets get into Julia. The function in the right hand of the ODE, which I call “f”, has a matrix form in your case. You can write a Julia function `f` which takes an array `x` as above but also an array of parameters `k`given by

\vec{k} = (k\_1, k\_r, k\_f, k\_\mathrm{cat})^T.

Then your Julia code would look something like

```julia
f(x, k) = [-1 0 0 0;
           -2 0 0 0;
           1 -1 1 0;
           0 -1 1 1;
           0 1 -1 -1;
           0 0 0 1] * [k[1]*x[1]*x[2]^2,
                         k[2]*x[4]*x[3],
                         k[3]*x[5],
                         k[4]*x[6]]

```

I’m explaining myself? I hope that this is useful…

Edit: syntax correction

---

<div class="post-metadata">

**Author:** ![Kester](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kester/32/28870_2.png) [@Kester](https://discourse.julialang.org/u/Kester)\
**Post date:** [September 4, 2021, 4:25am UTC](https://discourse.julialang.org/t/solving-a-petersen-matrix-with-julia/67674/3 "2021-09-04T04:25:59Z")

</div>

Yes that is very helpful, thank you. I have experience with simpler equations in DifferentialEquations.jl, I don’t know why I was not understanding how to do this… your example makes perfect sense though.

Part of my error was with syntax. I thought matrix format as like:

```julia
A = [1. 0 0 -5
      4 -2 4 -3
     -4 0 0 1
      5 -2 2 3]

```

But I see you are using an array? Or an array of vectors? (Sorry I am still very new with data types in Julia).

EDIT: I tried the format in your answer and it gave me a dimension mismatch when attempting to use DiffEq.jl. However this format worked:

```julia
f(x, k) = [-1 0 0 0;
              -2 0 0 0; 
              1 -1 1 0; 
              0 -1 1 1; 
              0 1 -1 -1; 
              0 0 0 1] * [k[1]*x[1]*x[2]^2,
                            k[2]*x[4]*x[3],
                            k[3]*x[5],
                            k[4]*x[6]]

```

Still think there is more I need to learn about data types and structures in Julia. Thanks again for the help. If anyone has any more pointers it would be great!

---

<div class="post-metadata">

**Author:** ![Erasmo98](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/erasmo98/32/26098_2.png) [@Erasmo98](https://discourse.julialang.org/u/Erasmo98)\
**Post date:** [September 4, 2021, 5:14am UTC](https://discourse.julialang.org/t/solving-a-petersen-matrix-with-julia/67674/4 "2021-09-04T05:14:53Z")

</div>

Sorry, for the wrong syntax; I didn’t run the example. 😅

---

<div class="post-metadata">

**Author:** ![yewalenikhil65](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yewalenikhil65/32/26873_2.png) [@yewalenikhil65](https://discourse.julialang.org/u/yewalenikhil65)\
**Post date:** [September 4, 2021, 5:48am UTC](https://discourse.julialang.org/t/solving-a-petersen-matrix-with-julia/67674/5 "2021-09-04T05:48:41Z")

</div>

@Kester  
You can check out how to directly write reactions as shown in Wikipedia page you have cited, using `Catalyst.jl` package [Home · Catalyst.jl](https://catalyst.sciml.ai/dev/) and, convert it into `ODESystem` automatically.

---

<div class="post-metadata">

**Author:** ![Kester](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kester/32/28870_2.png) [@Kester](https://discourse.julialang.org/u/Kester)\
**Post date:** [September 4, 2021, 6:34pm UTC](https://discourse.julialang.org/t/solving-a-petersen-matrix-with-julia/67674/6 "2021-09-04T18:34:12Z")

</div>

Youre example was very helpful though, thank you!

Do you know what the difference is between the formats? I think one is a matrix, while one is a vector of vectors?

---

<div class="post-metadata">

**Author:** ![Kester](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kester/32/28870_2.png) [@Kester](https://discourse.julialang.org/u/Kester)\
**Post date:** [September 4, 2021, 6:39pm UTC](https://discourse.julialang.org/t/solving-a-petersen-matrix-with-julia/67674/7 "2021-09-04T18:39:30Z")

</div>

Thanks for the suggestion! I looked at this package, as well as ModelingToolkit.jl. I see how these would function, but I don’t quite see the benefit compared to just explicitly writing the ODEs myself and solving them. Is the goal with these packages simplicity when you have very complex systems?

With Catalyst the benefit may be based on how the problems are formulated? Since my field of research I am normally constructing a Petersen Matrix, I only need to transpose to get the system rate equation.

---

<div class="post-metadata">

**Author:** ![Erasmo98](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/erasmo98/32/26098_2.png) [@Erasmo98](https://discourse.julialang.org/u/Erasmo98)\
**Post date:** [September 4, 2021, 8:35pm UTC](https://discourse.julialang.org/t/solving-a-petersen-matrix-with-julia/67674/8 "2021-09-04T20:35:46Z")

</div>

Yes you are right. You wrote a 6x4 matrix while I wrote a 6 element array of 1x4 matrices.

---

<div class="post-metadata">

**Author:** ![yewalenikhil65](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yewalenikhil65/32/26873_2.png) [@yewalenikhil65](https://discourse.julialang.org/u/yewalenikhil65)\
**Post date:** [September 5, 2021, 4:48am UTC](https://discourse.julialang.org/t/solving-a-petersen-matrix-with-julia/67674/9 "2021-09-05T04:48:34Z")

</div>

`Catalyst.jl` benefits for automatic conversion of reaction system to ODE, SDEs etc. If you know explicitly the matrix form it’s nice(as in your case)

But if you do not know explicitly ODEs from chemical reaction (say for large networks) , you might find `Catalyst.jl` beneficial as you do not have to do tedious task of checking whether the entries in matrix are correct or not. It gives you a way to explicitly find this matrix and rate expressions. Do take a look at its documentation.

---

<div class="post-metadata">

**Author:** ![genkuroki](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/genkuroki/32/18030_2.png) [@genkuroki](https://discourse.julialang.org/u/genkuroki)\
**Post date:** [September 8, 2021, 3:55am UTC](https://discourse.julialang.org/t/solving-a-petersen-matrix-with-julia/67674/10 "2021-09-08T03:55:47Z")

</div>

**Simple working example**

Load packages:

```julia
using DifferentialEquations
using Parameters
using StaticArrays
using Plots

```

Describe the diff. eq. with StaticArrays (Assume that `u` is a `SVector`of size (6,), `M` a `SMatrix` of size (6, 4), and `p` a named tuple (or struct) of the all model parameters):

```julia
function f(u, p, t)
    @unpack M, k1, kf, kr, kcat = p
    A, B, S, E, ES, P = u
    v = SVector(k1*A*B^2, kf*E*S, kr*ES, kcat*ES)
    M*v
end

```

The reason to use StaticArrays (SVector, SMatrix, etc.) is to improve performance.

cf. [https://diffeq.sciml.ai/stable/tutorials/ode\_example/#ode\_other\_types](https://diffeq.sciml.ai/stable/tutorials/ode_example/#ode_other_types)

Parameters:

```julia
M = [
    -1 0 0 0
    -2 0 0 0
    +1 -1 +1 0
     0 -1 +1 +1
     0 +1 -1 -1
     0 0 0 +1
]

p = (
    M = SMatrix{6, 4, Float64}(M),
    k1 = 1.0, 
    kf = 1.0, 
    kr = 1.0, 
    kcat = 1.0, 
)

u0 = SVector(
    #= A =# 1.5, 
    #= B =# 1.2, 
    #= S =# 0.0, 
    #= E =# 1.0, 
    #= ES =# 0.0, 
    #= P =# 0.0,
)

tspan = (0.0, 10.0)

```

Solve the problem and plot the result:

```julia
prob = ODEProblem(f, u0, tspan, p)
sol = solve(prob)
label = ["A" "B" "S" "E" "ES" "P"]
linestyle = [:dash :dash :dashdot :dashdot :solid :solid]
plot(sol; label, linestyle, lw=1.5)

```

 ![2021-09-08](https://global.discourse-cdn.com/julialang/original/3X/a/2/a2329721aa14a7d149a68d7deb7feff9a6876098.png)

You can obtain the same result with

```julia
function g(u, p, t)
    @unpack M, k1, kf, kr, kcat = p
    A, B, S, E, ES, P = u
    dA = -k1*A*B^2
    dB = -2k1*A*B^2
    dS = k1*A*B^2 - kf*E*S + kr*ES
    dE = -kf*E*S + kr*ES + kcat*ES
    dES = kf*E*S - kr*ES - kcat*ES
    dP = kcat*ES
    SVector(dA, dB, dS, dE, dES, dP)
end

prob = ODEProblem(g, u0, tspan, p)
sol = solve(prob)
label = ["A" "B" "S" "E" "ES" "P"]
linestyle = [:dash :dash :dashdot :dashdot :solid :solid]
plot(sol; label, linestyle, lw=1.5)

```

Jupyter notebook: [https://github.com/genkuroki/public/blob/main/0019/Petersen%20matrix.ipynb](https://github.com/genkuroki/public/blob/main/0019/Petersen%20matrix.ipynb)

---

<div class="post-metadata">

**Author:** ![Kester](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kester/32/28870_2.png) [@Kester](https://discourse.julialang.org/u/Kester)\
**Post date:** [September 11, 2021, 4:29pm UTC](https://discourse.julialang.org/t/solving-a-petersen-matrix-with-julia/67674/11 "2021-09-11T16:29:54Z")

</div>

Thank you very much for the detailed examples, and for the info about StaticArrays. I learned a lot looking through your code here!