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: | |
---|---|---|---|
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: |
Description
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))
Change History (5)
comment:1 follow-up: ↓ 2 Changed 9 years ago by
comment:2 in reply to: ↑ 1 Changed 9 years ago by
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
- Status changed from new to needs_review
comment:4 Changed 9 years ago by
- 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
- Resolution set to duplicate
- Status changed from positive_review to closed
Duplicate of #13268?