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.