DifferentialAlgebra Examples
Examples
Study of a differential-algebraic system
The following example is from [HW96]. It features a differential-algebraic system of differential index 2, describing a steering wheel. The dependent variables x⁡t and y⁡t are space coordinates. The dependent variable z⁡t is a command and the system is given by
ⅆⅆtx⁡t=7⁢y⁡t10+sin⁡5⁢z⁡t2,ⅆⅆty⁡t=7⁢x⁡t5+cos⁡5⁢z⁡t2,x⁡t2+y⁡t2=1
The RosenfeldGroebner function permits this system to be transformed into an equivalent system of differential index 0. In particular, it permits the following:
To compute a complete set of relations, defining the algebraic variety, over which the integral curves of the system must remain. These relations are essential for choosing consistent initial values.
To compute the missing ODE, which expresses the evolution of z⁡t.
First, the above system is transformed into a polynomial differential form by introducing two new dependent variables s⁡t and c⁡t which encode the sine and the cosine.
sys1≔diff⁡x⁡t,t=710⁢y⁡t+s⁡t,diff⁡y⁡t,t=1410⁢x⁡t+c⁡t,x⁡t2+y⁡t2=1,diff⁡s⁡t,t=2510⁢diff⁡z⁡t,t⁢c⁡t,diff⁡c⁡t,t=−2510⁢diff⁡z⁡t,t⁢s⁡t,s⁡t2+c⁡t2=1
sys1≔ⅆⅆtx⁡t=7⁢y⁡t10+s⁡t,ⅆⅆty⁡t=7⁢x⁡t5+c⁡t,x⁡t2+y⁡t2=1,ⅆⅆts⁡t=5⁢ⅆⅆtz⁡t⁢c⁡t2,ⅆⅆtc⁡t=−5⁢ⅆⅆtz⁡t⁢s⁡t2,s⁡t2+c⁡t2=1
Then, a differential polynomial ring is defined, endowed with an orderly ranking, and the system is simplified.
with⁡DifferentialAlgebra:
R≔DifferentialRing⁡derivations=t,blocks=z,y,x,s,c
R≔differential_ring
ideal≔RosenfeldGroebner⁡sys1,R
ideal≔regular_differential_chain
Because the ranking is orderly, the order 0 differential polynomials of the returned regular differential chain provide the algebraic variety, over which, the integral curves of the system must remain. These differential polynomials form a non-differential squarefree regular chain (see RegularChains), that can be solved from bottom-up. They are important for determining consistent initial values. Indeed, once z⁡0 is chosen, s⁡0 and c⁡0 are fixed, and, x⁡0 cannot be chosen freely.
Equations⁡ideal1,order=0,notation=jet,solved
y=−100⁢c2⁢s⁢x−210⁢c⁢s⁢x2−441⁢s⁢x3+210⁢c⁢s+341⁢s⁢x100⁢c3−100⁢c,x4=−2021⁢x3⁢c+341441⁢x2+2021⁢x⁢c+100441⁢c2,s2=−c2+1
The missing ODE, which expresses the evolution of z⁡t, is the element of the regular differential chain, whose leading derivative is a derivative of z.
Equations⁡ideal1,leader=derivative⁡z⁡t,notation=jet,solved
zt=−55566000⁢c5⁢x2+58344300⁢c4⁢x3−37044000⁢c5−43306200⁢c4⁢x−56831670⁢c3⁢x2−54519507⁢c2⁢x3+36335670⁢c3+40951407⁢c2⁢x+13461000⁢c⁢x2+10892700⁢x3−7987000⁢c−8422700⁢x11025000⁢c5−11025000⁢c3+2500000⁢c
Chemistry Example: The Henri-Michaelis-Menten Reduction
This example features a model reduction for chemical reaction systems, based on a quasi-steady state approximation [BLLM07]. The process is illustrated over a famous case, but has more general applications.
The following system is derived from a basic enzymatic reaction system, by application of the mass-action law. Four chemical species are involved: an enzyme E, a substrate S, an intermediate complex C, and a product P. There are three chemical reactions. The dependent variables represent the concentrations of the four species.
k1
k−1
k2
E+S
--->
C,
C
---->
E+S,
E+P.
sys2≔diff⁡E⁡t,t=−k1⁢E⁡t⁢S⁡t+k−1⁢C⁡t+k2⁢C⁡t,diff⁡S⁡t,t=−k1⁢E⁡t⁢S⁡t+k−1⁢C⁡t,diff⁡C⁡t,t=k1⁢E⁡t⁢S⁡t−k−1⁢C⁡t−k2⁢C⁡t,diff⁡P⁡t,t=k2⁢C⁡t
sys2≔ⅆⅆtE⁡t=−k1⁢E⁡t⁢S⁡t+k−1⁢C⁡t+k2⁢C⁡t,ⅆⅆtS⁡t=−k1⁢E⁡t⁢S⁡t+k−1⁢C⁡t,ⅆⅆtC⁡t=k1⁢E⁡t⁢S⁡t−k−1⁢C⁡t−k2⁢C⁡t,ⅆⅆtP⁡t=k2⁢C⁡t
Under some assumptions, this system can be approximated by the so-called Henri-Michaelis-Menten formula. We are going to show how this formula can be automatically recovered. The main assumption is: the revertible reaction is fast, and, the last one is slow. In other words, k1,k−1 >> k2.
Taking into account this main assumption actually amounts to perform a quasi-steady state approximation of the above system. This can be done as follows.
First, in the system, replace the contribution k1⁢E⁡t⁢S⁡t−k−1⁢C⁡t of the fast reactions, which occurs in the first three ODEs, by a new dependent variable F1. This variable represents the unknown contribution, which cannot be assumed to be 0, of the fast reactions, after the transient step. Somehow, this variable simply encodes that fast reactions are reactions.
Second, add the algebraic equation k1⁢E⁡t⁢S⁡t−k−1⁢C⁡t=0 to the system. This equation forces the dynamics of the system to stays over the algebraic variety that it defines. Somehow, it encodes the speed.
Third, eliminate the new dependent variable F1 from the differential-algebraic system. This elimination, performed using RosenfeldGroebner, computes the simplest ODE which describes the evolution of S⁡t, and which is a consequence of the DAE.
The two first steps yield the following differential-algebraic system. The rate constants are considered to be arbitrary parameters of the problem.
DAE≔diff⁡E⁡t,t=−F1⁡t+k2⁢C⁡t,diff⁡S⁡t,t=−F1⁡t,diff⁡C⁡t,t=F1⁡t−k2⁢C⁡t,diff⁡P⁡t,t=k2⁢C⁡t,k1⁢E⁡t⁢S⁡t−k−1⁢C⁡t=0
DAE≔ⅆⅆtE⁡t=−F1⁡t+k2⁢C⁡t,ⅆⅆtS⁡t=−F1⁡t,ⅆⅆtC⁡t=F1⁡t−k2⁢C⁡t,ⅆⅆtP⁡t=k2⁢C⁡t,k1⁢E⁡t⁢S⁡t−k−1⁢C⁡t=0
species≔C,E,P,S
rate_constants≔k1,k−1,k2
R≔DifferentialRing⁡blocks=F1,species,derivations=t,arbitrary=rate_constants
The third step can now be performed. An inequation, C⁡t≠0, is added to the DAE system passed to RosenfeldGroebner to avoid a spurious discussion on the case C⁡t=0.
ideal≔RosenfeldGroebner⁡op⁡DAE,C⁡t≠0,R
The returned regular differential chain involves a differential polynomial expressing the evolution of S⁡t; however, this relation is not yet the Henri-Michaelis-Menten formula. Indeed, the minor hypotheses have not yet been taken into account. This formula is just the simplest formula that can be deduced from the DAE, by using the main assumption only.
Equations⁡ideal1,leader=diff⁡S⁡t,t,solved
ⅆⅆtS⁡t=−E⁡t⁢S⁡t2⁢k12⁢k2+E⁡t⁢S⁡t⁢k1⁢k−1⁢k2E⁡t⁢k1⁢k−1+S⁡t⁢k1⁢k−1+k−12
A few more parameters, C0,E0,P0 and S0, are introduced and the linear conservation laws of the system involving them:
conservation_laws≔E⁡t+C⁡t=E0+C0,S⁡t+C⁡t+P⁡t=S0+C0+P0
A minor assumption is: P⁡0 = C⁡0 = 0.
hypotheses≔subs⁡P0=0,C0=0,conservation_laws
hypotheses≔E⁡t+C⁡t=E0,S⁡t+C⁡t+P⁡t=S0
At this point, we could either add E0,C0,S0 and P0 to the constants declared as arbitrary when defining the differential ring, or we can call RosenfeldGroebner directly because by default all non declared objects will be taken as arbitrary anyway
One could actually proceed from the formula established above, but, it is simpler to perform the whole simplification anew.
ideal≔RosenfeldGroebner⁡op⁡DAE,op⁡hypotheses,R
formula≔Equations⁡ideal1,leader=diff⁡S⁡t,t,solved
formula≔ⅆⅆtS⁡t=−S⁡t2⁢k12⁢k2⁢E0+S⁡t⁢k1⁢k−1⁢k2⁢E0S⁡t2⁢k12+2⁢S⁡t⁢k1⁢k−1+k1⁢k−1⁢E0+k−12
Not yet: one still needs to perform some renaming on parameters, and to neglect the term K⁢E0, assuming S⁡0 >> E⁡0.
formula≔eval⁡formula,k−1=K⁢k1:
formula≔normal⁡formula:
formula≔algsubs⁡k2⁢E0=Vmax,formula:
Here is the Henri-Michaelis-Menten formula.
formula≔normal⁡eval⁡formula,K⁢E0=0
formula≔ⅆⅆtS⁡t=−Vmax⁢S⁡tK+S⁡t
The pendulum
Consider the set of differential polynomials describing the motion of a pendulum, that is, a point mass m suspended from a massless rod of length l under the influence of gravity g, in Cartesian coordinates x,y. The Lagrangian formulation leads to two second order differential equations with an algebraic constraint. T is the Lagrangian multiplier.
params≔m,l,g:
P≔m⁢diff⁡y⁡t,t,t+T⁡t⁢y⁡t+g,m⁢diff⁡x⁡t,t,t+T⁡t⁢x⁡t,x⁡t2+y⁡t2−l2
P≔m⁢ⅆ2ⅆt2y⁡t+T⁡t⁢y⁡t+g,m⁢ⅆ2ⅆt2x⁡t+T⁡t⁢x⁡t,x⁡t2+y⁡t2−l2
To find the lowest order differential polynomials vanishing on the zeros of the system modeling the pendulum and in particular all the algebraic constraints (that is, the differential polynomials of order 0), use an orderly ranking. The second component in the output of RosenfeldGroebner corresponds to the equilibria of the pendulum. From the first component of the output, it can be seen that the actual motion is described by two first order differential equations with a constraint: the solution depends on only two arbitrary constants.
R≔DifferentialRing⁡derivations=t,blocks=T,x,y,arbitrary=params
ideal≔RosenfeldGroebner⁡P,R
ideal≔regular_differential_chain,regular_differential_chain
Equations⁡ideal1,notation=jet,solved
Tt=−3⁢yt⁢gl2,yt2=−−T⁢l4+T⁢l2⁢y2−g⁢l2⁢y+g⁢y3m⁢l2,x2=l2−y2
Equations⁡ideal2,notation=jet,solved
T=−y⁢gl2,x=0,y2=l2
To obtain the differential system satisfied by x and y alone, eliminate T. Note here that the motion is given by a second order differential equation in y. Then the Lagrange multiplier is given explicitly in terms of y and its first derivative.
Equations⁡ideal,notation=jet,leader<T⁡t,solved
yt,t=−−l2⁢m⁢y⁢yt2−g⁢l4+2⁢g⁢l2⁢y2−g⁢y4−l4⁢m+l2⁢m⁢y2,x2=l2−y2,x=0,y2=l2
General and singular solutions of a single differential equation
In general, when RosenfeldGroebner is applied on a single differential polynomial F, it returns at least one component involving F, and, possibly, some other components. The first component describes the general solution of F. The other ones describe the singular solutions.
R≔DifferentialRing⁡derivations=t,blocks=y
F≔diff⁡y⁡t,t3−4⁢t⁢y⁡t⁢diff⁡y⁡t,t+8⁢y⁡t2
F≔ⅆⅆty⁡t3−4⁢t⁢y⁡t⁢ⅆⅆty⁡t+8⁢y⁡t2
Over this example, RosenfeldGroebner returns the general component plus two singular ones. As we shall see, one of the singular components is essential, while the other one is not [K73,H99].
ideal≔RosenfeldGroebner⁡F,R
ideal≔regular_differential_chain,regular_differential_chain,regular_differential_chain
Equations⁡ideal
ⅆⅆty⁡t3−4⁢t⁢y⁡t⁢ⅆⅆty⁡t+8⁢y⁡t2,27⁢y⁡t−4⁢t3,y⁡t
The three components can be solved explicitly. The solutions of the general component writes y⁡t=c⁢t−c2, where c denotes an arbitrary constant. The solution of the first singular component writes y⁡t=4⁢t327. It is not a particular case of the general solution. This first singular component is said to be essential. The solution of the second singular component writes y⁡t=0. It is a particular case of the general solution, since it is obtained for c=0. This second singular component is said to be non-essential.
The function RosenfeldGroebner with the option singsol = essential only returns the essential components of F:
ideal≔RosenfeldGroebner⁡F,R,singsol=essential
ⅆⅆty⁡t3−4⁢t⁢y⁡t⁢ⅆⅆty⁡t+8⁢y⁡t2,27⁢y⁡t−4⁢t3
Power Series Solutions
The regular differential chains returned by RosenfeldGroebner, permit to compute integral formal power series solutions of the input system. By integral, it is meant that the exponents are integers. By formal, it is meant that convergence issues are not addressed. The following system provides a good example of what a general PDE polynomial system may look like, Its power series solutions turn out to be polynomials.
R≔DifferentialRing⁡derivations=x,y,blocks=v,u
sys3≔ux2−4⁢u,ux,y⁢vy−u+1,vx,x−ux
ideal≔RosenfeldGroebner⁡sys3,R
Equations⁡ideal1,solved
vx,x=ux,vy=−−u⁢ux⁢uy+ux⁢uy4⁢u,ux2=4⁢u,uy2=2⁢u
Looking at the leading derivatives of the returned regular differential chain, one sees that the system solutions depend on three constants, or initial values, which can be chosen almost arbitrary (see [DL84] for related undecidability results). Indeed, they must be solutions of the polynomial system of equations, returned by PowerSeriesSolution with the option conditions:
PowerSeriesSolution⁡conditions,ideal1
D1⁡u⁡0,02−4⁢u⁡0,0=0,D2⁡u⁡0,02−2⁢u⁡0,0=0,u⁡0,0≠0,D1⁡u⁡0,0≠0,D2⁡u⁡0,0≠0
Let us choose the following initial values, and, check that they do satisfy the above polynomial system.
iv≔u=c02,uy=sqrt⁡2⁢c0,ux=2⁢c0,v=c1,vx=c2
iv≔u=c02,uy=2⁢c0,ux=2⁢c0,v=c1,vx=c2
PowerSeriesSolution⁡conditions,ideal1,iv
0=0,0=0,c0≠0
The PowerSeriesSolution function can now compute the integral formal power series solutions of sys3, for these initial values.
sols≔PowerSeriesSolution⁡ideal1,3,iv
sols≔v⁡x,y=c1+c02⁢22−22⁢y+c2⁢x+c0⁢y22+2⁢c0⁢x⁢y+c0⁢x2+2⁢y312+x⁢y22+2⁢x2⁢y2+x33,u⁡x,y=c02+2⁢c0⁢y+2⁢c0⁢x+y22+2⁢x⁢y+x2
equations≔Equations⁡sys3,R,notation=diff
equations≔∂2∂x2v⁡x,y−∂∂xu⁡x,y,∂2∂x∂yu⁡x,y⁢∂∂yv⁡x,y−u⁡x,y+1,∂∂xu⁡x,y2−4⁢u⁡x,y
expand⁡eval⁡equations,sols
0,0,0
Download Help Document