Representation of ODEs with mix of Matrix and Vector Components

New to Julia
1

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

# 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:

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 ]
2

Hehe, I asked that yesterday on Zulip as well. @ChrisRackauckas suggested ArrayPartition from 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.

3

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

Currently I’m trying the following:

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?

4

dx is not defined in α, try

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)

5

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.

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!

6

Hi again,

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

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

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.

7

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

8

GitHub - SciML/LabelledArrays.jl: Arrays which also have a label for each element for easy scientific machine learning (SciML) 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.

1 Like
9

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 :slightly_smiling_face:.