Proposal of a DifferentialAlgebra package, relying on the C BLAD libraries

Introduction

The DifferentialAlgebra Sage package is an analogue of the MAPLE 14 DifferentialAlgebra package. The underlying theory is the differential algebra of Ritt and Kolchin. Its main tool is a simplifier for systems of polynomial differential equations, ordinary or with partial derivatives, called RosenfeldGroebner. It is related to the differential elimination theory. This simplifier decomposes the radical differential ideal I generated by an input system, as an intersection of radical differential ideals presented by regular differential chains (a slight generalization of Ritt characteristic sets). The output permits to test membership in the differential ideal I.

Further developments

The package should be developed to enhance numerical solvers of DAE (computation of the underlying ODE and of the equations which give the constraints on the initial values).
The package provides an implementation of differential fields, which could be very important for many developments (Ore algebra, differential modules, Kähler differentials, ...).
The package involves an implementation of Ritt's Low Power Theorem, for analysing the solutions of a single polynomial differential equation (general, particular, singular solution).
The development of the package was undertaken as first step, towards a control theory package.

Software

The package is written in Cython. The computations are performed by the BLAD libraries (C libraries, 60000 lines, LGPL license). The interface between Sage and BLAD is handled by the BMI library (C library, 10000 lines, LGPL license).

Getting started

The attached rebuild file is a shell command file which should help to build the whole stuff. This file was tested on Linux architectures.

An example

Borrowed from DifferentialAlgebra.pyx, to motivate (hopefully) reviewers.

        sage: from sage.libs.blad.DifferentialAlgebra import DifferentialRing, RegularDifferentialChain, BaseFieldExtension
        sage: leader,order,rank = var ('leader,order,rank')
        sage: derivative = function ('derivative')

    This example shows how to build the Henri Michaelis Menten formula by differential elimination. One considers a chemical reaction system describing the enzymatic reaction:

                   k(1)
        E + S  -----------> ES
                   k(-1)
        ES     -----------> E + S
                   k(2)
        ES     -----------> E + P

    A substrate S is transformed into a product P, in the presence of an enzyme E. An intermediate complex ES is formed.

        sage: t = var('t')
        sage: k,F_1,E,S,ES,P = function('k,F_1,E,S,ES,P')
        sage: params = [k(-1),k(1),k(2)]
        sage: params
        [k(-1), k(1), k(2)]

    The main assumption is that k(1), k(-1) >> k(2) i.e. that the revertible reaction is much faster than the last one. One performs a quasi-steady state approximation by considering the following differential-algebraic system (it comes from the mass-action law kinetics, replacing the contribution of the fast reactions by an unknown function F_1(t), on the algebraic variety where the fast reaction would equilibrate if they were alone).

        sage: syst = [diff(E(t),t) == - F_1(t) + k(2)*ES(t), diff(S(t),t) == - F_1(t), diff (ES(t),t) == - k(2)*ES(t) + F_1(t), diff (P(t),t) == k(2)*ES(t), 0 == k(-1)*E(t)*S(t) - k(1)*ES(t) ]
        sage: syst
        [D[0](E)(t) == k(2)*ES(t) - F_1(t), D[0](S)(t) == -F_1(t), D[0](ES)(t) == -k(2)*ES(t) + F_1(t), D[0](P)(t) == k(2)*ES(t), 0 == k(-1)*E(t)*S(t) - k(1)*ES(t)]

    Differential elimination permits to simplify this DAE. To avoid discussing the possible vanishing of ``params``, one moves them to the base field of the equations.

        sage: Field = BaseFieldExtension (generators = params)
        sage: Field
        differential_field

        sage: R = DifferentialRing (derivations = [t], blocks = [F_1, [E,ES,P,S], params], parameters = params)
        sage: R
        differential_ring

    The RosenfeldGroebner considers three cases. The two last ones are degenerate cases.

        sage: ideal = R.RosenfeldGroebner (syst, basefield = Field)
        sage: ideal
        [regular_differential_chain, regular_differential_chain, regular_differential_chain]
        sage: [ C.equations (solved = true) for C in ideal ]
        [[E(t) == k(1)*ES(t)/(k(-1)*S(t)), D[0](S)(t) == -(k(-1)*k(2)*S(t)^2*ES(t) + k(1)*k(2)*S(t)*ES(t))/(k(-1)*S(t)^2 + k(1)*S(t) + k(1)*ES(t)), D[0](P)(t) == k(2)*ES(t), D[0](ES)(t) == -k(1)*k(2)*ES(t)^2/(k(-1)*S(t)^2 + k(1)*S(t) + k(1)*ES(t)), F_1(t) == (k(-1)*k(2)*S(t)^2*ES(t) + k(1)*k(2)*S(t)*ES(t))/(k(-1)*S(t)^2 + k(1)*S(t) + k(1)*ES(t))], [S(t) == -k(1)/k(-1), ES(t) == 0, E(t) == 0, D[0](P)(t) == 0, F_1(t) == 0], [S(t) == 0, ES(t) == 0, D[0](P)(t) == 0, D[0](E)(t) == 0, F_1(t) == 0]]

    The sought equation, below, is not yet the Henri-Michaelis-Menten formula. This is expected, since some minor hypotheses have not yet been taken into account

        sage: ideal [0].equations (solved = true, selection = leader == derivative (S(t)))
        [D[0](S)(t) == -(k(-1)*k(2)*S(t)^2*ES(t) + k(1)*k(2)*S(t)*ES(t))/(k(-1)*S(t)^2 + k(1)*S(t) + k(1)*ES(t))]

    Let us take them into account. First create two new constants. Put them among ``params``, together with initial values.

        sage: K,V_max = var ('K,V_max')
        sage: params = [k(-1),k(1),k(2),E(0),ES(0),P(0),S(0),K,V_max]
        sage: params
        [k(-1), k(1), k(2), E(0), ES(0), P(0), S(0), K, V_max]

        sage: R = DifferentialRing (blocks = [F_1, [ES,E,P,S], params], parameters = params, derivations = [t])
        sage: R
        differential_ring

    There are relations among the parameters: initial values supposed to be zero, and equations meant to rename constants.

        sage: relations_among_params = RegularDifferentialChain ([P(0) == 0, ES(0) == 0, K == k(1)/k(-1), V_max == k(2)*E(0)], R)
        sage: relations_among_params
        regular_differential_chain

    Coming computations will be performed over a base field defined by generators and relations

        sage: Field = BaseFieldExtension (generators = params, relations = relations_among_params)
        sage: Field
        differential_field

    Extend the DAE with linear conservation laws. They could have been computed from the stoichimetry matrix of the chemical system.

        sage: newsyst = syst
        sage: newsyst.append (E(t) + ES(t) == E(0) + ES(0))
        sage: newsyst.append (S(t) + ES(t) + P(t) == S(0) + ES(0) + P(0))
        sage: newsyst
        [D[0](E)(t) == k(2)*ES(t) - F_1(t), D[0](S)(t) == -F_1(t), D[0](ES)(t) == -k(2)*ES(t) + F_1(t), D[0](P)(t) == k(2)*ES(t), 0 == k(-1)*E(t)*S(t) - k(1)*ES(t), E(t) + ES(t) == E(0) + ES(0), S(t) + ES(t) + P(t) == S(0) + ES(0) + P(0)]

    Simplify again. Only one case is left.

        sage: ideal = R.RosenfeldGroebner (newsyst, basefield = Field)
        sage: ideal
        [regular_differential_chain]

    To get the traditional Henri-Michaelis-Menten formula, one still needs to neglect the term K*E(0)

        sage: ideal[0].equations (solved = true, selection = leader == derivative (S(t)))
        [D[0](S)(t) == -(K*V_max*S(t) + V_max*S(t)^2)/(K^2 + K*E(0) + 2*K*S(t) + S(t)^2)]

    One can also get it by computing the right hand side of the equation which gives the evolution of the product P

        sage: ideal[0].normal_form (diff(P(t),t))
        V_max*S(t)/(K + S(t))

Component: packages: standard

Keywords: package, differential algebra, elimination theory

Issue created by migration from https://trac.sagemath.org/ticket/13267

sagemath / sage

Proposal of a DifferentialAlgebra package, relying on the C BLAD libraries #13267