#13268: Proposal of a DifferentialAlgebra package, relying on the C BLAD
libraries
---------------------------------------------------------------------+------
Reporter: boulier |
Owner: boulier
Type: enhancement |
Status: new
Priority: major |
Milestone: sage-5.6
Component: optional packages |
Resolution:
Keywords: package, differential algebra, elimination theory | Work
issues:
Report Upstream: N/A |
Reviewers: Charles Bouillaguet, Karl-Dieter Crisman
Authors: Nicolas M. Thiéry, François Boulier |
Merged in:
Dependencies: |
Stopgaps:
---------------------------------------------------------------------+------
Changes (by Bouillaguet):
* type: task => enhancement
Old description:
> 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).
>
> '''Release manager : '''
>
> * Make [[http://www.lifl.fr/~bouillaguet/download/blad-3.9.spkg]] an
> experimental/optional package
> * Apply [attachment:13268_differential_algebra.patch]
>
> '''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))
> }}}
New description:
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).
'''Release manager : '''
* Make [[http://www.lifl.fr/~bouillaguet/download/blad-3.9.spkg]] an
experimental/optional package
* Apply [attachment:13268_differential_algebra.patch]
* Apply [attachment:13268_differential_algebra_support_files.patch]
* Apply [attachment:13268_module_list.patch]
'''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))
}}}
--
--
Ticket URL: <http://trac.sagemath.org/sage_trac/ticket/13268#comment:10>
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.
