Opened 9 years ago

Closed 9 years ago

#13267 closed task (duplicate)

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

Reported by: boulier Owned by: Francois.Boulier@…
Priority: major Milestone: sage-duplicate/invalid/wontfix
Component: packages: standard Keywords: package, differential algebra, elimination theory
Cc: Merged in:
Authors: Reviewers:
Report Upstream: N/A Work issues:
Branch: Commit:
Dependencies: Stopgaps:

Status badges



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.


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:

        E + S  -----------> ES
        ES     -----------> E + S
        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

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

    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

    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

    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

    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

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

Change History (5)

comment:1 follow-up: Changed 9 years ago by jhpalmieri

Duplicate of #13268?

comment:2 in reply to: ↑ 1 Changed 9 years ago by boulier

Replying to jhpalmieri:

Duplicate of #13268?

Yes, sorry for my clumsy manipulations. I am quite new as a Sage user and a Sage contributor.

comment:3 Changed 9 years ago by jhpalmieri

  • Status changed from new to needs_review

comment:4 Changed 9 years ago by jhpalmieri

  • Milestone changed from sage-5.3 to sage-duplicate/invalid/wontfix
  • Status changed from needs_review to positive_review

comment:5 Changed 9 years ago by jdemeyer

  • Resolution set to duplicate
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.