#13267: Proposal of a DifferentialAlgebra package, relying on the C BLAD
libraries
------------------------+---------------------------------------------------
Reporter: boulier | Owner: Francois.Boulier@…
Type: task | Status: new
Priority: major | Milestone: sage-5.3
Component: packages | Keywords: package, differential algebra,
elimination theory
Work issues: | Report Upstream: N/A
Reviewers: | Authors:
Merged in: | Dependencies:
Stopgaps: |
------------------------+---------------------------------------------------
'''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))
}}}
--
Ticket URL: <http://trac.sagemath.org/sage_trac/ticket/13267>
Sage <http://www.sagemath.org>
Sage: Creating a Viable Open Source Alternative to Magma, Maple, Mathematica,
and MATLAB
--
You received this message because you are subscribed to the Google Groups
"sage-trac" group.
To post to this group, send email to [email protected].
To unsubscribe from this group, send email to
[email protected].
For more options, visit this group at
http://groups.google.com/group/sage-trac?hl=en.
