diff --git a/examples/Exact_and_optimal_quadratization_of_nonlinear_finite-dimensional_non-autonomous_dynamical_systems.ipynb b/examples/Exact_and_optimal_quadratization_of_nonlinear_finite-dimensional_non-autonomous_dynamical_systems.ipynb new file mode 100644 index 0000000..de810fd --- /dev/null +++ b/examples/Exact_and_optimal_quadratization_of_nonlinear_finite-dimensional_non-autonomous_dynamical_systems.ipynb @@ -0,0 +1,766 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "fec3a854-55f3-4d0a-8e03-69c75cf8f658", + "metadata": {}, + "source": [ + "## Demo for [\"Exact and optimal quadratization of nonlinear finite-dimensional non-autonomous dynamical systems\"](https://arxiv.org/abs/2303.10285)" + ] + }, + { + "cell_type": "markdown", + "id": "59ec86e9-bf01-4a33-8158-cfeb0886156a", + "metadata": {}, + "source": [ + "This notebook contains several examples of usage of QBee following Sections 5-7 of the paper" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "d37a9871-469b-4b77-a887-4e42e4a23ba4", + "metadata": {}, + "outputs": [], + "source": [ + "# Loading packages\n", + "import sympy as sp\n", + "from qbee import *" + ] + }, + { + "cell_type": "markdown", + "id": "af2a0786-407f-469f-b55f-ec7bee5e5b24", + "metadata": {}, + "source": [ + "## Introductory examples (Section 5)" + ] + }, + { + "cell_type": "markdown", + "id": "77146346-40f9-4a09-81e2-ad03791f54f6", + "metadata": {}, + "source": [ + "### Basic usage of QBee (Section 5.1)\n", + "\n", + "First we show how one can use QBee to quadratize an autonomous finite-dimensional ODE system such as\n", + "$$\n", + "\\begin{cases}\n", + "\\dot{x}_1 = x_1^2 + x_2^2,\\\\\n", + "\\dot{x}_2 = x_1 + x_2\n", + "\\end{cases}\n", + "$$" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "e9b594b0-c422-48a6-85d6-554b94c67933", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + " \r" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Introduced variables:\n", + "w0 = x1**2\n", + "w1 = x2**2\n", + "\n", + "x1' = w0*x1 + w1\n", + "x2' = x1 + x2\n", + "w0' = 2*w0**2 + 2*w1*x1\n", + "w1' = 2*w1 + 2*x1*x2\n" + ] + } + ], + "source": [ + "x1, x2 = functions(\"x1, x2\")\n", + "system = [ # pairs of the form (x, x')\n", + " (x1, x1**3 + x2**2),\n", + " (x2, x1 + x2)\n", + "]\n", + "quadratize(system).print()" + ] + }, + { + "cell_type": "markdown", + "id": "a0262b45-4e17-4fdc-b868-a753c653b8a5", + "metadata": {}, + "source": [ + "We have thus obtained a monomial quadratization of order two which is guaranteed to be optimal among the monomial quadratizations." + ] + }, + { + "cell_type": "markdown", + "id": "cb7364ee-8a97-47a7-98d1-39a9d0918915", + "metadata": {}, + "source": [ + "### Polynomial ODEs with inputs (Section 5.2)\n", + "\n", + "Now we proceed to the new functionality such as systems with inputs. Consider a scalar ODE from Example 5.2\n", + "$$\n", + "\\dot{x} = x^2u\n", + "$$\n", + "It can be quadratized using QBee as follows." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "9695569a-b183-411a-9fca-1ee21923d667", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + " " + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Introduced variables:\n", + "w0 = u*x\n", + "\n", + "x' = w0*x\n", + "w0' = u'*x + w0**2\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "\r" + ] + } + ], + "source": [ + "x, u = functions(\"x, u\")\n", + "system = [ # pairs of the form (x, x')\n", + " (x, x**2 * u)\n", + "]\n", + "quadratize(system).print() # input-free=False by default" + ] + }, + { + "cell_type": "markdown", + "id": "27dc8a19-7260-4b19-a8d5-6008de0e17f5", + "metadata": {}, + "source": [ + "We observe that the resulting quadratic system involves the derivative of the input, that is, $\\dot{u}$. As explained in Example 3.8 of the paper, in this particular case, this was unavoidable. However, in many cases one can find a quadratization which does not involve the derivatives of the inputs (we call such quadratizations *input-free*). Search for input-free quadratizations can be forced by setting the `input_free` keyword argument to be `True`. \n", + "\n", + "As an example, consider a system\n", + "$$\n", + "\\begin{cases}\n", + " \\dot{x}_1 = x_1 + x_1u,\\\\\n", + " \\dot{x}_2 = x_1^2 u.\n", + "\\end{cases}\n", + "$$\n", + "The following code can be used to find an input-free quadratization for it." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "0ce4746e-d670-4a8c-b74a-7059390f1a6b", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + " " + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Introduced variables:\n", + "w0 = x1**2\n", + "\n", + "x1' = u*x1 + x1\n", + "x2' = u*w0\n", + "w0' = 2*u*w0 + 2*w0\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "\r" + ] + } + ], + "source": [ + "x1, x2, u = functions(\"x1, x2, u\")\n", + "system = [ # pairs of the form (x, x')\n", + " (x1, x1 + x1 * u),\n", + " (x2, x1**2 * u)\n", + "]\n", + "quadratize(system, input_free=True).print()" + ] + }, + { + "cell_type": "markdown", + "id": "2333abce-4513-4e94-b6b1-02ee782d58c1", + "metadata": {}, + "source": [ + "And we have obtained an optimal monomial input-free quadratization of order one." + ] + }, + { + "cell_type": "markdown", + "id": "aa140281-9a1f-4b7b-928d-21d248b6dbf9", + "metadata": {}, + "source": [ + "### Dimension-agnostic quadratization (Section 5.3)\n", + "\n", + "We will now show how QBee can be applied to systems of variable dimension occurring, for example, as discretizations of PDEs. We will consider a PDE coming from traffic modeling (see Example 4.1 for further details):\n", + "$$\n", + "\\frac{\\partial \\rho}{\\partial t}(t, \\xi) = \\rho(t, \\xi) + \\rho^2(t, \\xi) \\frac{\\partial \\rho(t, \\xi)}{\\partial \\xi}, \\qquad \\rho(t, 0) = 0, \\quad \\rho(t, 1) = 1.\n", + "$$\n", + "For any integer $n$, we can discretize the spatial domain uniformly $x_{i}^{[n]}(t) = \\rho(t, i/(n + 1))$.\n", + "Then the vector $\\mathbf{x}^{[n]}(t) = [x_{1}^{[n]}(t), \\ldots, x_{n}^{[n]}(t)]^\\top$ satisfies an ODE system of dimension $n$\n", + "$$\n", + "\\dot{\\mathbf{x}}^{[n]} = \\mathbf{x}^{[n]} + (\\mathbf{x}^{[n]})^2 \\odot (\\mathbf{D}\\mathbf{x}^{[n]}),\n", + "$$\n", + "where $\\mathbf{D}$ is the matrix for taking the finite difference.\n", + "\n", + "QBee allows to describe (in finite terms) a set of new variables which will yield a quadratization for any $n$ and any matrix $\\mathbf{D}$ (we call this *dimension-agnostic* quadratization). The system above can be encoded as follows (see Example 5.4)." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "224b18f6-6bea-4147-8b49-35f991e5a023", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + " " + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Every ODE system in the family can be quadratized by adding the following variables\n", + "\t* For each i, we take variables\n", + "\tx_i**2\n", + "\t* For every i and j such that the variables from the j-th block (that is, x_j) appear in the right-hand sides of the i-th block (that is, in x_i'), we take\n", + "\tx_i*x_j\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "\r" + ] + } + ], + "source": [ + "x, Dx = functions(\"x Dx\")\n", + "system = [\n", + " (x, x + x**2 * Dx)\n", + "]\n", + "quadratize_dimension_agnostic(system);" + ] + }, + { + "cell_type": "markdown", + "id": "e34cb2e0-8b9d-4d22-9325-9d9aa822766b", + "metadata": {}, + "source": [ + "And we see that the result of the computation is a recipy how to construct a quadratization depending on $n$ and the positions of the nonzero entries in $\\mathbf{D}$." + ] + }, + { + "cell_type": "markdown", + "id": "ef7e9b8d-bd33-4a16-b363-9156175d54d7", + "metadata": {}, + "source": [ + "### Polynomialization with QBee (Section 5.4)\n", + "\n", + "QBee can also take as input a systems with nonpolynomial right-hand side and lift it to a polynomial one. This can be used as a preprocessing for quadratization but also as a feature on its own. Consider, for example, a scalar nonpolynomial ODE from Example 5.5 \n", + "$$\n", + "\\dot{x} = e^{-x} + e^{-2x}. \n", + "$$\n", + "This ODE can be lifted to a polynomial system using QBee as follows." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "102a5a6e-497f-4538-86fd-7cd5e9641e28", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + " " + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "x' = w_0**2 + w_0\n", + "w_0' = -w_0*(w_0**2 + w_0)\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "\r" + ] + } + ], + "source": [ + "x = functions(\"x\")\n", + "system = [(x, sp.exp(-x) + sp.exp(-2 * x))] # list of pairs (x, x')\n", + "print(polynomialize(system))" + ] + }, + { + "cell_type": "markdown", + "id": "eeb8f772-baa0-4b1b-bcfc-3a95b8204b02", + "metadata": {}, + "source": [ + "Note that we got a cubic system which can be further quadratized. In fact, this can be done by a single call." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "17437783-6fc8-456a-bca9-1b637f7c780b", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + " \r" + ] + }, + { + "data": { + "text/plain": [ + "w_0 = exp(-x)\n", + "w_1 = w_0**2\n", + "\n", + "x' = w_0 + w_1\n", + "w_0' = -w_0*w_1 - w_1\n", + "w_1' = -2*w_0*w_1 - 2*w_1**2" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "polynomialize_and_quadratize(system)" + ] + }, + { + "cell_type": "markdown", + "id": "b75a5fbc-6298-4de1-b153-970fa055b5ba", + "metadata": {}, + "source": [ + "## Examples from the literature" + ] + }, + { + "cell_type": "markdown", + "id": "b982b709-fef6-4567-b5c3-70a2d22260a9", + "metadata": {}, + "source": [ + "### Tubular reactor model (Section 6.1)\n", + "\n", + "We consider two versions of a non-adiabatic tubular reactor model with a single reaction: with polynomial and exponential reaction terms.\n", + "In both cases, the ODE system of interest will be a result of spatial discretization of a PDE system with two vectors of states: species concentration $\\boldsymbol{\\psi}(t)$ and temperature $\\boldsymbol{\\theta}(t)$." + ] + }, + { + "cell_type": "markdown", + "id": "b9a5eff7-f87a-4e9e-980b-0d00ce642c0a", + "metadata": {}, + "source": [ + "The *polynomial* version (Section 6.1.1) will be a $2n$ dimensional system of ODEs with scalar input $u(t)$ and several scalar or vector parameters:\n", + "$$\n", + "\\begin{cases}\n", + " \\boldsymbol{\\dot{\\psi}} = \\mathbf{A}_\\psi \\boldsymbol{\\psi} + \\mathbf{b}_\\psi - \\mathcal{D}\\boldsymbol{\\psi} \\odot (\\mathbf{c}_0 + \\mathbf{c}_1 \\odot\\boldsymbol{\\theta} + \\mathbf{c}_2 \\odot\\boldsymbol{\\theta}^2 + \\mathbf{c}_3 \\odot\\boldsymbol{\\theta}^3), \\\\\n", + " \\boldsymbol{\\dot{\\theta}} = \\mathbf{A}_\\theta \\boldsymbol{\\theta} + \\mathbf{b}_\\theta + \\mathbf{b} u + \\mathcal{B} \\mathcal{D} \\boldsymbol{\\psi} \\odot (\\mathbf{c}_0 + \\mathbf{c}_1 \\odot\\boldsymbol{\\theta} + \\mathbf{c}_2 \\odot\\boldsymbol{\\theta}^2 + \\mathbf{c}_3 \\odot\\boldsymbol{\\theta}^3),\n", + "\\end{cases}\n", + "$$\n", + "We can search for a dimension-agnostic quadratization of this model using QBee as follows." + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "5b5d98e3-ae3e-4a1e-a661-5573c5780fd9", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Found inputs u, they will be considered space-independent\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + " \r" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Every ODE system in the family can be quadratized by adding the following variables\n", + "\t* For each i, we take variables\n", + "\ttheta_i**2, psi_i*theta_i, psi_i*theta_i**2, theta_i**3\n" + ] + } + ], + "source": [ + "# Dpsi and Dtheta below refer to the matrix-vector products A_psi * psi and A_theta * theta\n", + "# the letter D in from of the variable name should be used so that QBee matches the vector and\n", + "# its transforms\n", + "psi, theta, w1, Dpsi, Dtheta = functions(\"psi theta w1 Dpsi Dtheta\")\n", + "u = functions(\"u\")\n", + "B, D, gamma, bpsi, btheta, b = parameters(\"B D gamma bpsi btheta b\")\n", + "c0, c1, c2, c3 = parameters(\"c0 c1 c2 c3\")\n", + "\n", + "poly = c0 + c1 * theta + c2 * theta**2 + c3 * theta**3\n", + "system = [\n", + " (psi, Dpsi + bpsi - D * psi * poly),\n", + " (theta, Dtheta + btheta + b * u + B * D * psi * poly),\n", + "]\n", + "quadratize_dimension_agnostic(system);" + ] + }, + { + "cell_type": "markdown", + "id": "41f1c288-b6a6-4ede-be02-d7f62f0439a9", + "metadata": {}, + "source": [ + "This means that the system can be always quadratized (independently of the shapes of $\\mathbf{A}_\\psi$ and $\\mathbf{A}_\\theta$) with $4n$ extra variables $\\boldsymbol{\\theta}^2,\\; \\boldsymbol{\\psi}\\odot \\boldsymbol{\\theta},\\; \\boldsymbol{\\psi}\\odot\\boldsymbol{\\theta}^2,\\; \\boldsymbol{\\theta}^3$." + ] + }, + { + "cell_type": "markdown", + "id": "243314ae-3e8d-40f6-98e4-129fe86b40e3", + "metadata": {}, + "source": [ + "The *exponential* version of the model uses the same states but involves exponential function.\n", + "$$\n", + "\\begin{cases}\n", + "\\boldsymbol{\\dot{\\psi}} = \\mathbf{A}_{\\psi} \\boldsymbol{\\psi} + \\mathbf{b}_\\psi - \\mathcal{D} \\boldsymbol{\\psi} \\odot e^{\\gamma - \\gamma / \\boldsymbol{\\theta}},\\\\\n", + "\\boldsymbol{\\dot{\\theta}} = \\mathbf{A}_{\\theta} \\boldsymbol{\\theta} + \\mathbf{b}_\\theta + \\mathbf{b} u(t) + \\mathcal{B} \\mathcal{D} \\boldsymbol{\\psi} \\odot e^{\\gamma - \\gamma / \\bm{\\theta}}\n", + "\\end{cases}\n", + "$$\n", + "QBee cannot find dimension-agnostic quadratizations with arbitrary non-polynomial right-hand side but can handle the case of Laurent monomials (i.e. negative powers). In order to reduce to this case, we manually introduce $\\mathbf{w}_1 = e^{\\gamma - \\gamma / \\boldsymbol{\\theta}}$ and obtain\n", + "$$\n", + "\\begin{cases}\n", + "\\boldsymbol{\\dot{\\psi}} = \\mathbf{A}_{\\psi} \\boldsymbol{\\psi} + \\mathbf{b}_\\psi - \\mathcal{D} \\boldsymbol{\\psi} \\odot \\mathbf{w}_1,\\\\\n", + "\\boldsymbol{\\dot{\\theta}} & = \\mathbf{A}_{\\theta} \\boldsymbol{\\theta} + \\mathbf{b}_\\theta + \\mathbf{b} u(t) + \\mathcal{B} \\mathcal{D} \\boldsymbol{\\psi} \\odot \\mathbf{w}_1,\\\\\n", + " \\dot{\\mathbf{w}}_1 & = \\gamma\\boldsymbol{\\dot{\\theta}} \\odot \\frac{1}{\\boldsymbol{\\theta}^2} \\odot \\boldsymbol{w}_1.\n", + "\\end{cases}\n", + "$$\n", + "Now we can, similarly to the polynomial version, obtain a dimension-agnostic quadratization using QBee" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "97a663d8-0ba2-458b-872d-a02a388b4721", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Found inputs u, they will be considered space-independent\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Nodes processed: 383 nodes [00:04, 83.85 nodes/s, Current best order = 36, To return current best press Ctrl+C] " + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "The quadratization algorithm has been interrupted. Returning the current best...\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Nodes processed: 248 nodes [00:03, 77.36 nodes/s, Current best order = 36, To return current best press Ctrl+C]" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "The quadratization algorithm has been interrupted. Returning the current best...\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + " \r" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Every ODE system in the family can be quadratized by adding the following variables\n", + "\t* For each i, we take variables\n", + "\t1/theta_i, theta_i**(-2), u/theta_i, u/theta_i**2, psi_i*w1_i/theta_i**2, psi_i*w1_i/theta_i\n", + "\t* For every i and j such that the variables from the j-th block (that is, w1_j, theta_j, psi_j) appear in the right-hand sides of the i-th block (that is, in w1_i', theta_i', psi_i'), we take\n", + "\ttheta_j/theta_i, psi_j/theta_i, theta_j/theta_i**2, psi_j/theta_i**2\n" + ] + } + ], + "source": [ + "psi, theta, w1, Dpsi, Dtheta = functions(\"psi theta w1 Dpsi Dtheta\")\n", + "u = functions(\"u\")\n", + "\n", + "B, D, gamma, bpsi, btheta, b = parameters(\"B D gamma bpsi btheta b\")\n", + "\n", + "theta_diff = Dtheta + btheta + b * u + B * D * psi * w1\n", + "system = [\n", + " (psi, Dpsi + bpsi - D * psi * w1),\n", + " (theta, theta_diff),\n", + " (w1, gamma * w1 / theta**2 * theta_diff)\n", + "]\n", + "quadratize_dimension_agnostic(system);" + ] + }, + { + "cell_type": "markdown", + "id": "07150d8c-43ac-4554-bc81-fadffdd85ce5", + "metadata": {}, + "source": [ + "The full computation for this example will take some time (~20 minutes) but already after one minute it will display `Current best order = 36` meaning that some quadratization has been found (the exact number 36 can be ignored, it originates from the depths of the algorithm) and the algorithm is trying to refine it. One can interrupt the computation, and the current best option will be displayed (as above). Interestingly, it turns out that this is indeed the best quadratization to be found. Compared to the polynomial version, the shapes of the matrices matter now." + ] + }, + { + "cell_type": "markdown", + "id": "36ac8d95-b7e1-40bd-a6f0-29abaaf8e518", + "metadata": {}, + "source": [ + "### Species’ reaction model for combustion (Section 6.2)\n", + "\n", + "We consider the rocket engine combustion model which is a non-polynomial non-autonomous ODE system in four variables $x_1, x_2, x_3, x_4$ denoting the concentrations of different species and with one external input $u(t)$.\n", + "$$\n", + "\\begin{cases}\n", + " \\dot{x}_1 = -A \\exp{\\Big(-\\frac{E}{R u(t)}\\Big)} x_1^{0.2} x_2^{1.3}\\\\\n", + " \\dot{x}_2 = 2 \\dot{x}_1\\\\\n", + " \\dot{x}_3 = -\\dot{x}_1 \\\\\n", + " \\dot{x}_4 = -2 \\dot{x}_1.\n", + "\\end{cases}\n", + "$$\n", + "The system can be polynomialized and then quadratized using QBee as follows." + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "4bfd24e3-eb62-4784-b8f8-071b7365e6dc", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + " \r" + ] + }, + { + "data": { + "text/plain": [ + "w_0 = x2**1.3\n", + "w_1 = exp(-E/(R*u))\n", + "w_2 = x1**0.2\n", + "w_3 = w_0*w_1\n", + "w_4 = u'/u**2\n", + "w_5 = u**(-2)\n", + "w_6 = u'/u\n", + "w_7 = 1/u\n", + "w_8 = w_0*w_1*w_2/x1\n", + "w_9 = w_0*w_1*w_2/x2\n", + "\n", + "x1' = -A*w_2*w_3\n", + "x2' = -2*A*w_2*w_3\n", + "x3' = A*w_2*w_3\n", + "x4' = 2*A*w_2*w_3\n", + "w_0' = -13*A*w_0*w_9/5\n", + "w_1' = E*w_1*w_4/R\n", + "w_2' = -A*w_2*w_8/5\n", + "w_3' = -13*A*w_3*w_9/5 + E*w_3*w_4/R\n", + "w_4' = u''*w_5 - 2*w_4*w_6\n", + "w_5' = -2*w_5*w_6\n", + "w_6' = u''*w_7 - w_6**2\n", + "w_7' = -w_6*w_7\n", + "w_8' = 4*A*w_8**2/5 - 13*A*w_8*w_9/5 + E*w_4*w_8/R\n", + "w_9' = -A*w_8*w_9/5 - 3*A*w_9**2/5 + E*w_4*w_9/R" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "x1, x2, x3, x4 = functions(\"x1, x2, x3, x4\")\n", + "u = functions(\"u\")\n", + "A, E, R = parameters(\"A, E, R\")\n", + "x1_diff = -A * sp.exp(-E / (R * u)) * x1**0.2 * x2**1.3\n", + "\n", + "system = [\n", + " (x1, x1_diff),\n", + " (x2, 2 * x1_diff),\n", + " (x3, -x1_diff),\n", + " (x4, -2 * x1_diff)\n", + "]\n", + "# we allow the input order go up to two because u' will already appear at the polynomialization step\n", + "polynomialize_and_quadratize(system, input_der_orders={u: 2})" + ] + }, + { + "cell_type": "markdown", + "id": "39abf3b8-c900-4511-aa28-50d3293d8e73", + "metadata": {}, + "source": [ + "As a result, we get a quadratization of order 10 which involves the second derivative of the input (it is possible to prove that it is impossible to have a lower order derivative of the input)." + ] + }, + { + "cell_type": "markdown", + "id": "51655191-587b-4fca-a9bb-bd6624e1d270", + "metadata": {}, + "source": [ + "### Solar wind modeling (Section 7)" + ] + }, + { + "cell_type": "markdown", + "id": "55dd03d9-b4e7-4abd-a721-5d44adf3ac93", + "metadata": {}, + "source": [ + "The solar wind model we consider is a result of semi-discretization of a PDE and it is an $n$-dimensional ODE system in a vector of unknown functions $\\mathbf{v}(r)$ in the $r$ variable:\n", + "$$\n", + "\\frac{\\text{d}\\mathbf{v}(r)}{\\text{d}r} = \\mathbf{D}\\ln \\left[\\mathbf{v}(r)\\right] - \\frac{\\xi}{\\Omega_\\text{rot}(\\hat{\\theta})}\\mathbf{D} \\mathbf{v}(r),\n", + "$$\n", + "where $\\mathbf{D}$ is the matrix encoding the finite difference scheme. In order to use QBee to search for a dimension-agnostic quadratization, we eliminate the logarithmic term by introducing $\\mathbf{w}_0 := \\ln (\\mathbf{v})$. We also set $C_1 := \\frac{\\xi}{\\Omega_{\\text{rot}}(\\hat{\\theta})}$ to simplify the notation. We therefore obtain a $2n$-dimensional ODE model\n", + "$$\n", + "\\begin{cases}\n", + "\\frac{\\text{d} \\mathbf{v}}{\\text{d}r} = \\mathbf{D} \\mathbf{w}_0 - C_1 \\mathbf{D} \\mathbf{v},\\\\\n", + " \\frac{\\text{d} \\mathbf{w}_0}{\\text{d}r} = \\frac{1}{\\mathbf{v}}\\mathbf{D} \\mathbf{w}_0 - \\frac{C_1}{\\mathbf{v}}\\mathbf{D} \\mathbf{v}.\n", + "\\end{cases}\n", + "$$\n", + "Now we can use QBee to find a dimension-agnostic quadratization." + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "a1854811-b864-4f44-8097-63de45b98e6a", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + " \r" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Every ODE system in the family can be quadratized by adding the following variables\n", + "\t* For each i, we take variables\n", + "\t1/v_i, u_i/v_i\n", + "\t* For every i and j such that the variables from the j-th block (that is, v_j, u_j) appear in the right-hand sides of the i-th block (that is, in v_i', u_i'), we take\n", + "\tu_j/v_i, v_j/v_i\n" + ] + } + ], + "source": [ + "u, v, Du, Dv = functions(\"u v Du Dv\")\n", + "C1 = parameters(\"C1\")\n", + "system = [\n", + " (v, Du + C1 * Dv),\n", + " (u, Du / v + C1 * Dv / v)\n", + "]\n", + "quadratize_dimension_agnostic(system);" + ] + }, + { + "cell_type": "markdown", + "id": "cd89ab39-d637-4bc0-8cbf-33d822f033f7", + "metadata": {}, + "source": [ + "As explained in the paper, for the particular matrix $\\mathbf{D}$ of interest, one can reduce the produced quadratization further by hand." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "912206a3-f86b-4cc5-9802-fa2d3b458a84", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.18" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}