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

2b5de62a-05ae-46d3-a99b-c6405c33a11e commented 12 years ago

This ticket proposes the inclusion in SAGE of a new package, dedicated to differential algebra.

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).

Upstream: BLAD BMI

Installation

After building Sage, install the experimental packages blad and bmi, then run make or make build again, i.e:

make build
SAGE_ROOT=$(pwd) ./local/bin/sage -i blad
SAGE_ROOT=$(pwd) ./local/bin/sage -i bmi
make build

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 Rosenfeld-Groebner algorithm 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))

CC: @sagetrac-sage-combinat @sagetrac-boulier @burcin @kcrisman @BrentBaccala

Component: packages: optional

Keywords: package, differential algebra, elimination theory

Author: Nicolas M. Thiéry, François Boulier, Charles Bouillaguet

Branch/Commit: public/optional_spkg/differential_algebra-13268 @ 81aef63

Reviewer: Charles Bouillaguet, Karl-Dieter Crisman

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

9d1e8f33-9c5a-487d-bc38-8944bef1a4ae commented 2 years ago

comment:39

In the above example there is supposed to be a DifferentialAlgebra.pyx in this branch of sage:

https://github.com/sagemath/sagetrac-mirror/tree/public/optional_spkg/differential_algebra-13268/src/sage

see the examples at the very beginning which has the line:

from sage.libs.blad.DifferentialAlgebra import DifferentialRing, RegularDifferentialChain, BaseFieldExtension

Where is this now?

BrentBaccala commented 2 years ago

comment:40

Replying to @tdupu:

In the above example there is supposed to be a DifferentialAlgebra.pyx in this branch of sage:

https://github.com/sagemath/sagetrac-mirror/tree/public/optional_spkg/differential_algebra-13268/src/sage

see the examples at the very beginning which has the line:
from sage.libs.blad.DifferentialAlgebra import DifferentialRing, RegularDifferentialChain, BaseFieldExtension
Where is this now?

DifferentialRing is defined in src/sage/calculus/DifferentialAlgebra.pyx

You also need to install the BLAD and BMI libraries.

If you checkout the public/optional_spkg/differential_algebra-13268 branch and run git diff 9.0.beta4 --stat, that should give you a good idea of where things are in this branch.

sagemath / sage

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