Extract linear dimension reduction in a matrix form

[mleprovost]

Hello,

First, this is a very cool package! This is not to report a bug, but to ask a user question.

Is there a routine to return the matrix of coefficients from the old coordinate system to the new coordinate system for a given reduction?

For instance,

using ExactODEReduction

odes = @ODEsystem(
  x1'(t) = x1^2 + 2x1*x2,
  x2'(t) = x2^2 + x3 + x4,
  x3'(t) = x2 + x4,
  x4'(t) = x1 + x3
)
reductions = find_reductions(odes)

red1 = reductions[1]

new_system(red1)
## Prints:
y1'(t) = y1(t)^2 + y2(t)
y2'(t) = y1(t) + y2(t)

new_vars(red1)
## Prints:
Dict{Nemo.fmpq_mpoly, Nemo.fmpq_mpoly} with 2 entries:
  y2 => x3 + x4
  y1 => x1 + x2

Is there a simple way to get the matrix L = [1, 1, 0, 0; 0, 0, 1, 1] from reductions[1] ?

Can we also get the linear combinations of the old variables which are conserved?

Thank you for your help

[sumiya11]

Hi @mleprovost , thank you for your kind words.

I don't think there is an existing method for getting the matrix.

As a quick fix, you can use this branch #31, like this:

using ExactODEReduction

odes = @ODEsystem(
  x1'(t) = x1^2 + 2x1*x2,
  x2'(t) = x2^2 + x3 + x4,
  x3'(t) = x2 + x4,
  x4'(t) = x1 + x3
)
reductions = find_reductions(odes)

red1 = reductions[1]

red1.lumping_matrix
# prints
2×4 Matrix{Any}:
 1  1  0  0
 0  0  1  1

It can fail for some cases, we will perhaps fix it over the weekend :^).

Can we also get the linear combinations of the old variables which are conserved?

If I understand your question correctly, you can take the rows of the matrix.

Feel free to let me know if you have other questions.

[mleprovost]

Awesome, thank you.
Let me rephrase my second question:

I am interested in finding the linear conserved quantities of a system of ODEs.

Can we use ExactODEReduction.jl to find linear combinations of the original variables which have zero dynamics?

In the example above:
We complete the rows of L with the matrix Lcomplement = [1 -1 0 0; 0 0 -1 1] to form a basis of R^4

Can we say that the new variables y_invariant = Lcomplement times x have zero dynamics i.e. d(Lcomplement*x)/dt = 0?

[sumiya11]

Ah, I see, thanks for the clarification.

I have merged #31, together with a very simple function for finding linear conservation laws:

using ExactODEReduction

ode = @ODEsystem(
         x1'(t) = x1 + x2*x3,
         x2'(t) = x2*x3 + x3,
         x3'(t) = x1 - x3,
       )

ExactODEReduction.find_conservation_laws(ode)
# returns
1-element Vector{Any}:
 (basis = Nemo.fmpq_mpoly[x1, x2, x3], matrix = [-1 1 1], law = Nemo.fmpq_mpoly[-x1 + x2 + x3])

Let me know if that's what you were asking about.

In the example above:
We complete the rows of L with the matrix Lcomplement = [1 -1 0 0; 0 0 -1 1] to form a basis of R^4

I don't think a complement will have zero dynamics in general. In fact, in the matrix L itself, some of the rows may correspond to quantities with zero dynamics in the standard basis, but this is not guaranteed.

[mleprovost]

That's awesome, thank you so much for your help! What is the best way to cite this package and your work on this topic?

[sumiya11]

You can cite this paper :-)

@article{DDP24,
title = {Exact hierarchical reductions of dynamical models via linear transformations},
journal = {Communications in Nonlinear Science and Numerical Simulation},
volume = {131},
pages = {107816},
year = {2024},
issn = {1007-5704},
doi = {https://doi.org/10.1016/j.cnsns.2024.107816},
url = {https://www.sciencedirect.com/science/article/pii/S1007570424000029},
author = {Alexander Demin and Elizaveta Demitraki and Gleb Pogudin},
keywords = {Ordinary differential equations, Exact reduction, Lumping, Dimensionality reduction, Matrix algebras},
abstract = {Dynamical models described by ordinary differential equations (ODEs) are a fundamental tool in the sciences and engineering. Exact reduction aims at producing a lower-dimensional model in which each macro-variable can be directly related to the original variables, and it is thus a natural step towards the model’s formal analysis and mechanistic understanding. We present an algorithm which, given a polynomial ODE model, computes a longest possible chain of exact linear reductions of the model such that each reduction refines the previous one, thus giving a user control of the level of detail preserved by the reduction. This significantly generalizes over the existing approaches which compute only the reduction of the lowest dimension subject to an approach-specific constraint. The algorithm reduces finding exact linear reductions to a question about representations of finite-dimensional algebras. We provide an implementation of the algorithm, demonstrate its performance on a set of benchmarks, and illustrate the applicability via case studies. Our implementation is freely available at https://github.com/x3042/ExactODEReduction.jl.}
}

[sumiya11]

Out of curiosity, do you use these linear reductions in some application ?

[mleprovost]

AbsoIutely, I am looking at the conservation of linear constraints in data assimilation problems. I wanted to experiment my ideas on biological and chemical ODES which have many invariants arising from the conservation of elements.

I have notice that many of these biological ODEs (from BIOMEd) have wide disparity in the scaling of the initial conditions.
Let's imagine that I would like to rescale my variables to be of order O(1).

Consider a dynamical model $\frac{dx}{dt} = f(x, t)$

We consider the change of variable $\tilde{x} = D x$ with D a diagonal matrix of scaling numbers, then $\frac{d\tilde{x}}{dt} = D \frac{dx}{dt} = D f(x,t) = D f(D^{-1} \tilde{x}, t)$

Can we do define this rescaled system of ODEs with ExactODEReduction.jl?

[pogudingleb]

For just making a substitution, you can use some general-purpose tool for manipulating differential equations like ModelingToolkit, here is an example of code (I am not an expert in MTK though)

using ModelingToolkit

# define the starting system
@variables t, x1(t), x2(t)
D = Differential(t)
eqs = [D(x1) ~ x1 + x2, D(x2) ~ x1^2 + 3 * x2]

# this is the diagonal of your matrix
scaling = [2, 0.5]

# this is the map x -> D^-1 \tilde{x}
scaling_map = (x1, x2) .=> scaling.^(-1) .* (x1, x2)

# and here are the new equations
scaled_eqs = [e.lhs ~ s * substitute(e.rhs, scaling_map) for (s, e) in zip(scaling, eqs)]

Note that you can use the result of such substitution with ExactODEReduction.

[mleprovost]

Thank you all for your help !!