# Representation of ODEs with mix of Matrix and Vector Components

**URL:** <https://discourse.julialang.org/t/representation-of-odes-with-mix-of-matrix-and-vector-components/48687>\
**Category:** New to Julia\
**Tags:** diffeq\
**Created:** [October 20, 2020, 2:16pm UTC](https://discourse.julialang.org/t/representation-of-odes-with-mix-of-matrix-and-vector-components/48687 "2020-10-20T14:16:02Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![Sam\_Power](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sam_power/32/18972_2.png) [@Sam\_Power](https://discourse.julialang.org/u/Sam_Power)\
**Post date:** [October 20, 2020, 2:16pm UTC](https://discourse.julialang.org/t/representation-of-odes-with-mix-of-matrix-and-vector-components/48687/1 "2020-10-20T14:16:02Z")

</div>

Hi,

I’m currently using DifferentialEquations.jl to solve some ODEs arising in Chemical Reaction Networks, in particular the Linear Noise Approximation to the Chemical Langevin Equation.

In my example, the ODE involves a pair of 2-dimensional vectors (z, m), and a 2x2 matrix V, as follows

```julia
# dz/dt = \alpha(z)
# dm/dt = F(z) m 
# dV/dt = V F(z)' + F(z) V + \beta(z) \beta(z)'

```

where `(\alpha, \beta, F)` are defined below, and ’ denotes taking the transpose.

I wanted to ask about the syntax for ODEs which have ‘mixed form’ in this sense of containing both matrix and vector components, in case there are any relevant worked examples which I can look at, or if there are any extra details of which I should be wary.

To be concrete: given expressions for `(\alpha, \beta, F)` as below, how should I best go about aggregating them, in order to represent and solve the above ODE jointly for `(z, m , F)`

My current code is as follows:

```julia
using DifferentialEquations

function parameterised_α_LNA!(dx, x, c, t)
    dx[1] = c[1] * x[1] - c[2] * x[1] * x[2]
    dx[2] = c[2] * x[1] * x[2] - c[3] * x[2]
end

function parameterised_β_LNA!(σx, x, c, t)
    σx[1, 1] = c[1] * x[1] + c[2] * x[1] * x[2]
    σx[1, 2] = -c[2] * x[1] * x[2]
    σx[2, 1] = -c[2] * x[1] * x[2]
    σx[2, 2] = c[2] * x[1] * x[2] + c[3] * x[2]
end

function parameterised_F_LNA!(jx, x, c, t)
    jx[1, 1] = c[1] - c[2] * x[2]
    jx[1, 2] = -c[2] * x[1]
    jx[2, 1] = c[2] * x[2]
    jx[2, 2] = c[2] * x[1] - c[3]
end

x0 = [70.0; 80,0]
m0 = [0.0; 0.0]
V0 = [0.0 0.0; 0.0 0.0]
tspan = (0.0, 50.0)
c = [1.0, 0.005, 0.6]

```

---

<div class="post-metadata">

**Author:** ![mschauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mschauer/32/13946_2.png) [@mschauer](https://discourse.julialang.org/u/mschauer)\
**Post date:** [October 20, 2020, 2:21pm UTC](https://discourse.julialang.org/t/representation-of-odes-with-mix-of-matrix-and-vector-components/48687/2 "2020-10-20T14:21:55Z")

</div>

Hehe, I asked that yesterday on Zulip as well. @ChrisRackauckas suggested `ArrayPartition` from [https://github.com/SciML/RecursiveArrayTools.jl](https://github.com/SciML/RecursiveArrayTools.jl) to represent the different Matrices/Vector valued coordinate as long vector you can unpack into its constituents or access in a flattened way.

---

<div class="post-metadata">

**Author:** ![Sam\_Power](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sam_power/32/18972_2.png) [@Sam\_Power](https://discourse.julialang.org/u/Sam_Power)\
**Post date:** [October 20, 2020, 3:09pm UTC](https://discourse.julialang.org/t/representation-of-odes-with-mix-of-matrix-and-vector-components/48687/3 "2020-10-20T15:09:17Z")

</div>

Thanks, this looks useful - trying to follow along with some of the examples in the test folder.

Currently I’m trying the following:

```julia
using DifferentialEquations

z0 = [70.0; 80.0]
m0 = [0.0, 0.0]
V0 = [0.0 0.0; 0.0 0.0]
u0 = ArrayPartition(z0, m0, V0)

tspan = (0.0, 50.0)

function α(x, c)
    dx[1] = c[1] * x[1] - c[2] * x[1] * x[2]
    dx[2] = c[2] * x[1] * x[2] - c[3] * x[2]
    dx
end

function β(x, c)
    σx[1, 1] = c[1] * x[1] + c[2] * x[1] * x[2]
    σx[1, 2] = -c[2] * x[1] * x[2]
    σx[2, 1] = -c[2] * x[1] * x[2]
    σx[2, 2] = c[2] * x[1] * x[2] + c[3] * x[2]
    σx
end

function F(x, c)
    jx[1, 1] = c[1] - c[2] * x[2]
    jx[1, 2] = -c[2] * x[1]
    jx[2, 1] = c[2] * x[2]
    jx[2, 2] = c[2] * x[1] - c[3]
    jx
end

function drift(du, u, c, t)
    z = u[1,:]
    m = u[2,:]
    V = u[3,:]
    du[1] = α(z, c)
    du[2] = F(z, c) * m
    du[3] = V * transpose(F(z, c)) + F(z, c) * V + β(z, c) * transpose(β(z, c))
end

c = [1.0, 0.005, 0.6]
prob = ODEProblem(drift, u0, tspan, c)
sol = solve(prob)

```

and this runs into problems at the final stage; upon running `solve(prob)`, I get an error of `UndefVarError: dx not defined`, which points back to the definition of α.

I’m guessing that I need to either

- structure (α, β, F) differently,
- pass them to `drift` in a different format,
- apply `ArrayPartition` to `du` somehow

or something along these lines?

---

<div class="post-metadata">

**Author:** ![mschauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mschauer/32/13946_2.png) [@mschauer](https://discourse.julialang.org/u/mschauer)\
**Post date:** [October 20, 2020, 3:14pm UTC](https://discourse.julialang.org/t/representation-of-odes-with-mix-of-matrix-and-vector-components/48687/4 "2020-10-20T15:14:04Z")

</div>

`dx` is not defined in `α`, try

```julia
function α(x, c)
     dx = [c[1]*x[1] - c[2]*x[1]*x[2], c[2]*x[1]*x[2] - c[3]*x[2], ?]
     return dx
end

```

(not sure what `?` should be here)

---

<div class="post-metadata">

**Author:** ![Sam\_Power](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sam_power/32/18972_2.png) [@Sam\_Power](https://discourse.julialang.org/u/Sam_Power)\
**Post date:** [October 22, 2020, 11:41am UTC](https://discourse.julialang.org/t/representation-of-odes-with-mix-of-matrix-and-vector-components/48687/5 "2020-10-22T11:41:30Z")

</div>

Thanks, I’ve amended as suggested (also for the other functions similarly). `α` is just a vector of length 2, so I think it should be okay to leave `?` empty, unless I’ve misunderstood.

I’ve resolved as below. It seems that I had been a bit careless in how I recalled the matrix part of the ArrayPartition; I had been calling it as `u[3, :]`, but it needs to be called as `u[3, :, :]`, which makes sense in hindsight.

```julia
using DifferentialEquations

z0 = [70.0, 80.0]
m0 = [0.0, 0.0]
V0 = [0.0 0.0; 0.0 0.0]
u0 = ArrayPartition(z0, m0, V0)

tspan = (0.0, 50.0)

function α(x, c)
    dx = [c[1] * x[1] - c[2] * x[1] * x[2], c[2] * x[1] * x[2] - c[3] * x[2]]
    return dx
end

function β_squared(x, c)
    σx = Array{Float64}(undef, 2, 2)
    σx[1, 1] = c[1] * x[1] + c[2] * x[1] * x[2]
    σx[1, 2] = -c[2] * x[1] * x[2]
    σx[2, 1] = -c[2] * x[1] * x[2]
    σx[2, 2] = c[2] * x[1] * x[2] + c[3] * x[2]
    return σx
end

function F(x, c)
    jx = Array{Float64}(undef, 2, 2)
    jx[1, 1] = c[1] - c[2] * x[2]
    jx[1, 2] = -c[2] * x[1]
    jx[2, 1] = c[2] * x[2]
    jx[2, 2] = c[2] * x[1] - c[3]
    return jx
end

function drift(du, u, c, t)
    z = u[1,:]
    m = u[2,:]
    V = u[3,:, :]
    du[1,:] = α(z, c)
    du[2, :] = F(z, c) * m
    du[3, :, :] = V * transpose(F(z, c)) + F(z, c) * V + β_squared(z, c)
end

c = [1.0, 0.005, 0.6]

prob = ODEProblem(drift, u0, tspan, c)
sol = solve(prob)

```

I think I’m sorted for now; I’ll return to this thread if there are any other issues related to the mixed format. Thanks!

---

<div class="post-metadata">

**Author:** ![Sam\_Power](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sam_power/32/18972_2.png) [@Sam\_Power](https://discourse.julialang.org/u/Sam_Power)\
**Post date:** [October 26, 2020, 3:21pm UTC](https://discourse.julialang.org/t/representation-of-odes-with-mix-of-matrix-and-vector-components/48687/6 "2020-10-26T15:21:26Z")

</div>

Hi again,

I’ve streamlined the code slightly as follows (using some model-specific facts about the coefficients), and am now working with

```julia
using DifferentialEquations, LinearAlgebra, ForwardDiff

S = [1 -1 0; 0 1 -1]

function h(x, c)
    return h = [c[1] * x[1]; c[2] * x[1] * x[2]; c[3] * x[2]]
end

function α(x, c)
    return dx = S * h(x, c)
end

function β(x, c)
    return σx = S * Diagonal(h(x, c)) * transpose(S)
end

function F(x, c)
    return jx = ForwardDiff.jacobian(x -> α(x, c), x)
end

# LNA ODEs
# dz = α(z) dt
# dm = F(z) m dt
# dV = ( V F(z)' + F(z) V + β(z) ) dt

function lotvol_LNA_drift(du, u, c, t)
    z = u[1,:]
    m = u[2,:]
    V = u[3,:, :]
    du[1, :] = α(z, c)
    du[2, :] = F(z, c) * m
    du[3, :, :] = V * transpose(F(z, c)) + F(z, c) * V + β(z, c)
end

z0 = [70.0, 80.0]
m0 = [0.0, 0.0]
V0 = [0.0 0.0; 0.0 0.0]
u0 = ArrayPartition(z0, m0, V0)

tspan = (0.0, 50.0)

c = [1.0, 0.005, 0.6]

lotvol_LNA_ODE = ODEProblem(lotvol_LNA_drift, u0, tspan, c)
LNA_sol = solve(lotvol_LNA_ODE)

```

and this is all working great.

The remaining question I have concerns extracting the separate `(z, m, V)` coordinates of the ODE from `LNA_sol`. At the moment I’m doing

```julia
z_LNA = LNA_sol[1:2, :]
m_LNA = LNA_sol[3:4, :]
V_LNA = LNA_sol[5:8, :]

```

which _essentially_ gets the job done, but

1. extracts the `V` component as a 4-vector, rather than as an array, and
2. doesn’t make use of the `ArrayPartition` setup of the problem.

As such, I’m guessing / hoping that there’s a more elegant solution here, where I can pull the `V` component directly as a matrix, and treat it as such.

---

<div class="post-metadata">

**Author:** ![mschauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mschauer/32/13946_2.png) [@mschauer](https://discourse.julialang.org/u/mschauer)\
**Post date:** [October 26, 2020, 3:45pm UTC](https://discourse.julialang.org/t/representation-of-odes-with-mix-of-matrix-and-vector-components/48687/7 "2020-10-26T15:45:41Z")

</div>

If you feel bold you can call `dump(sol)` and see what you can figure out from that.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [October 26, 2020, 4:50pm UTC](https://discourse.julialang.org/t/representation-of-odes-with-mix-of-matrix-and-vector-components/48687/8 "2020-10-26T16:50:46Z")

</div>

[GitHub - SciML/LabelledArrays.jl: Arrays which also have a label for each element for easy scientific machine learning (SciML)](https://github.com/SciML/LabelledArrays.jl) might be a nice way to get a naming scheme on there to make the code look a little nicer. Or ArrayPartition. Basically, whatever array library you want.

Our plan is to automate the generation of the LNA with Catalyst.jl/ModelingToolkit.jl soon though. When that’s done, I’d recommend just using that.

---

<div class="post-metadata">

**Author:** ![Sam\_Power](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sam_power/32/18972_2.png) [@Sam\_Power](https://discourse.julialang.org/u/Sam_Power)\
**Post date:** [October 26, 2020, 4:52pm UTC](https://discourse.julialang.org/t/representation-of-odes-with-mix-of-matrix-and-vector-components/48687/9 "2020-10-26T16:52:51Z")

</div>

Thanks, I’ll take a look at these and see if there’s a neat solution.

Re: LNA, that’s great to hear - the project I’m working on at the moment uses LNA pretty heavily, so that would be super useful 🙂.