QBee

Tutorial in Jupyter (Colab)

QBee is a Python library for transforming systems of differential equations into a systems with quadratic right-rand side.

Installation

PyPI

pip install qbee

Manual

Clone repository: https://github.com/AndreyBychkov/QBee.git
- Or, if you want our bleeding edge version, clone https://github.com/AndreyBychkov/QBee/tree/dev
Change directory: cd QBee
Install package: pip install .

If you use poetry you can alternately install if with poetry install

What is quadratization?

The problem of quadratization is, given a system of ODEs with polynomial right-hand side, reduce the system to a system with quadratic right-hand side by introducing as few new variables as possible. We will explain it using toy example. Consider the system

$$ \begin{cases} x_1' = x_1 x_2 \ x_2' = -x_1 x_2^3 \end{cases} $$

An example of quadratization of this system will be a singular vector of new variables

$$ [y' = x_2 y - 2y^2] $$

leading to the following ODE

$$ y' = x_2 y - 2y^2 $$

Thus, we attained the system with quadratic right-hand side

$$ \begin{cases} x_1' = x_1 x_2 \ x_2 ' = -x^2 y \ y' = x_2 y - 2y^2 \end{cases} $$

We used only one new variable, so we achieved an optimal quadratization.

Qbee usage

QBee implements algorithms that take system of ODEs with elementary functions right-hand side and return optimal monomial quadratization - optimal quadratization constructed from monomial substitutions.

We will demonstrate usage of QBee on the example below. Other interactive examples you can find in examples section.

1. Importing QBee

QBee relies on Sympy for a high-level API.

import sympy as sp
from qbee import *

2. System definition

For example, we will take the A1 system from Swischuk et al.'2020

$$ \begin{cases} c_1' = -A \exp(-E_a / (R_u T)) c_1 ^{0.2} c_2^{1.3} \ c_2 ' = 2c_1' \ c_3' = -c_1' \ c_4' = -2 c_1' \end{cases} $$

The parameters in the system are A, Ea and Ru, and the others are either state variables or inputs. So, let's define them with the system in code:

A, Ea, Ru = parameters("A, Ea, Ru")
c1, c2, c3, c4, T = functions("c1, c2, c3, c4, T")  

eq1 = -A * sp.exp(-Ea / (Ru * T)) * c1 ** 0.2 * c2 ** 1.3
system = [
    (c1, eq1),
    (c2, 2 * eq1),
    (c3, -eq1),
    (c4, -2 * eq1)
]

3. Polynomialization and Quadratization

When we work with ODEs with the right-hand side being a general continuous function, we utilize the following pipeline:

Input system -> Polynomial System -> Quadratic System

and the transformations are called polynomialization and quadratization accordingly.

The example system is not polynomial, so we use the most general method for achieving optimal monomial quadratization.

# {T: 2} means than T can have a derivative of order at most two => T''
quadr_system = polynomialize_and_quadratize(system, input_der_orders={T: 2})
if quadr_system:
    quadr_system.print()

Sample output:

Introduced variables:
w_0 = exp(-Ea/(Ru*T))
w_1 = c1**0.2
w_2 = c2**1.3
w_3 = w_0*w_1
w_4 = T'/T**2
w_5 = T**(-2)
w_6 = T'/T
w_7 = 1/T
w_8 = w_0*w_1*w_2/c1
w_9 = w_0*w_1*w_2/c2

c1' = -A*w_2*w_3
c2' = -2*A*w_2*w_3
c3' = A*w_2*w_3
c4' = 2*A*w_2*w_3
w_0' = Ea*w_0*w_4/Ru
w_1' = -A*w_1*w_8/5
w_2' = -13*A*w_2*w_9/5
w_3' = -A*w_3*w_8/5 + Ea*w_3*w_4/Ru
w_4' = T''*w_5 - 2*w_4*w_6
w_5' = -2*w_5*w_6
w_6' = T''*w_7 - w_6**2
w_7' = -w_6*w_7
w_8' = 4*A*w_8**2/5 - 13*A*w_8*w_9/5 + Ea*w_4*w_8/Ru
w_9' = -A*w_8*w_9/5 - 3*A*w_9**2/5 + Ea*w_4*w_9/Ru

Introduced variables are the optimal monomial quadratization.

Papers

Optimal Monomial Quadratization for ODE systems: arxiv, Springer

Citation

If you find this code useful in your research, please consider citing the above paper that works best for you.

AndreyBychkov / QBee

readme