Quadratization of differential equations in python
Project description
QBee
Online playground
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
- Or, if you want our bleeding edge version, clone
- 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
Citation
If you find this code useful in your research, please consider citing the above paper that works best for you.
Project details
Download files
Download the file for your platform. If you're not sure which to choose, learn more about installing packages.