Introduction to the Mgfun package, a package
for Multivariate Generating FUnctions.
The Mgfun package is used to perform calculations on multivariate functions and
sequences. In particular, it is well suited for calculations on so-called special functions,
including hypergeometric functions and sequences, their q-analogues and holonomic
functions. In the philosophy of Mgfun, these functions and sequences are described by
systems of linear homogeneous differential or difference equations, or by q-analogues of
them. A natural tool that plays a central role in these calculations is a suitable
generalization of Groebner bases.
When using the Mgfun package, any function or sequence is represented by a set of linear
homogeneous equations it satisfies. When such systems are known for functions f(x,y)
and g(x,y), Mgfun can compute a system for the sum f(x,y)+g(x,y), the product
f(x,y)g(x,y) or definite integrals of f(x,y) with respect to y. In the case of sequences
u[n,m] and v[n,m], it can compute equations for the sum u[n,m]+v[n,m], the product
u[n,m]v[n,m] or definite summations of u[n,m] over all m. Mgfun can also deal with
mixed cases, such as f[n,m](x,y), or with other types of linear equations the user might
want to define.
Using these basic operations, one can:
- find closed forms (in terms of hypergeometric functions and sequences and their
q-analogues);
- compute (ordinary, exponential, Bessel...) generating functions,
- extract coefficients of a series (for any orthogonal basis of the algebra of series in a
general class);
- check identities and prove new ones;
- compute diagonals of series.
This worksheet gives an overview of Mgfun through some applications.
Once you have installed Mgfun, the following call should run with no error.
> with(Mgfun);
[algtoholon, annihilators, applyopr, commalg, dependency, fglm, gbasis,
hdefint, hdefqsum, hdefsum, hdiag, hilbertdim, hpowerrs, hprod, hprodrs,
hsum, hsumrs, hypergeomtoholon, leadcoeff, leading_term, leadmon, makeopr,
oppower, opprod, orealg, qbinomial, qfactorial, qpochhammer, randopr,
reduce, reducelist, reducescale, setparam, shiftalg, skewelim, skewgcdex,
skewpdiv, skewprem, spoly, termorder, testorder, weylalg]
We want Mgfun to be verbose.
> infolevel[Mgfun_warnings]:=1:
The Legendre polynomials.
As mentioned above, definite summations can be computed using the Mgfun package. As
an example, we show how to calculate equations that describe the Legendre polynomials
P[n](x). The n-th polynomial of this sequence is defined as the sum for all k of
binomial(n,k)^2*(x-1)^k*(x+1)^(n-k)*2^(-n). Let us now compute this sum.
We first declare an algebra of operators to enter equations that describe the summand.
We need differential equations with respect to x and recurrence equations with respect to
n and k. It will also be possible to enter mixed differential-recurrence equations.
> A:=orealg(diff=[Dx,x],shift=[Sn,n],shift=[Sk,k]):
Most calculations by Mgfun on multivariate functions make intensive use of Groebner
bases computations. Here, we declare an elimination term order to prepare an elimination
of k by Groebner basis computation, since k is the index of the summation we want to
compute.
> T:=termorder(A,lexdeg=[[k],[Dx,Sn,Sk]]):
We now enter equations satisfied by the term to be summed.
> G:=hypergeomtoholon(binomial(n,k)^2*(x-1)^k*(x+1)^(n-k)*2^(-n),A);
2
G := [(- x + 1) Dx + 2 k + n x - n,
2 2 2 2
(- 2 n - 4 n + 4 n k - 2 + 4 k - 2 k ) Sn + n x + n + 2 n x + 2 n + x + 1,
2 2 2 2 2 2
(- k x - k - 2 k x - 2 k - x - 1) Sk + n x - n - 2 n k x + 2 n k + k x - k
]
And we compute equations satisfied by the Legendre polynomials P[n](x).
> GL[P]:=hdefsum(G,k=0..infinity,T);
Mgfun/Holonomy/definitesum:
Assume the summand and its shifts vanish in 0 and infinity
GL[P] :=
2 2 2
[2 Sn n x + 3 Sn x - n - 1 - Sn n - 2 Sn , - n x + Sn n - Dx x + Dx + Sn - x]
A particular use of summation is to compute generating function. As an example, we
give here the generating function of Legendre polynomials. Recall that it is defined as
phi(x,y)=sum(P[n](x)y^n,n=0..infinity).
We follow the same method as above, beginning by declaring another algebra and another
term order well-suited for elimination of the index of summation n.
> A:=orealg(diff=[Dx,x],diff=[Dy,y],shift=[Sn,n]):\
T:=termorder(A,lexdeg=[[n],[Dx,Dy,Sn]]):\
BR:=ratpoly(rational,[x,y]):
From a description of P[n](x), we derive a description of P[n](x)y^n.
> G:=[op(numer(subs(Sn=Sn/y,GL[P]))),y*Dy-n]:
We next derive a description for their generating function phi(x,y).
> GL[phi]:=hdefsum(G,n=0..infinity,T);
Mgfun/Holonomy/definitesum:
Assume the summand and its shifts vanish in 0 and infinity
2 2
GL[phi] := [2 x y Dx - y Dx + y - Dx, 2 x y Dy + x - y Dy - y - Dy]
Solving these differential equations by separation of variables yields a closed form for the
generating function.
> applyopr(GL[phi][1],phi(x,y),A);
2 / d \
y phi(x, y) + (- y + 2 x y - 1) |---- phi(x, y)|
\ dx /
> pdesolve(",phi(x,y));
_F1(y)
phi(x, y) = ---------------------
2 1/2
(- y + 2 x y - 1)
Any solution must be of the previous form. The equation in derivatives with respect to y
makes it possible to compute _F1(y).
> numer(normal(applyopr(GL[FG][2],op(2,"),A)));
GL[FG][2] _F1(y)
Since the previous expression is identically zero, _F1(y) is a constant (at least in a
neighbourhood of 0) and phi(x,0)=P[0](x)=1 yields the following generating function.
> phi:=1/sqrt(1-2*x*y+y^2);
1
phi := -------------------
2 1/2
(1 - 2 x y + y )
Computing integrals using Mgfun is very similar to computing sums. As an example, we
prove that Legendre polynomials form an orthogonal family for the scalar product
<f(x),g(x)>=int(f(x)g(x),x=-1..1).
We prove this by integration of the bivariate generating function h of the products
J[n](x)J[m](x) between -1 and +1, so as to find the generating function of the scalar
products <J[n](x),J[m](x)>. So we first declare another algebra and define h.
> A:=weylalg([Dx,x],[Du,u],[Dv,v]):\
h:=subs(y=u,phi)*subs(y=v,phi);
1
h := ---------------------------------------
2 1/2 2 1/2
(1 - 2 x u + u ) (1 - 2 x v + v )
> G:=hypergeomtoholon(h,A);
G := [
2 2 2 2 2 2 2
(- 1 + 2 x v - v + 2 x u - 4 x u v + 2 x u v - u + 2 u x v - u v ) Dx + u
2 2
- 4 u x v + u v + v + v u ,
2 2
(- 1 + 2 x u - u ) Du + x - u, (- 1 + 2 x v - v ) Dv + x - v]
Once again, we will perform an elimination by calculation of a Groebner basis. However,
the existing function hdefint of Mgfun is not suited for this computation, since h does not
vanish at -1 and +1. Here we need to eliminate the variable of integration, x, by
lower-level tools.
> T:=termorder(A,lexdeg=[[x],[Dx,Du,Dv]]):
The Groebner basis calculation is performed by a call to gbasis. This function is Mgfun's
own routine, which performs more efficiently than grobner[gbasis] in the case of
commutative calculations. We keep only those polynomials that are free of x.
> EL:=collect(remove(has,gbasis(G,T,ratpoly(rational,[u,v])),x),Dx,factor);
EL := [
2 2
- 2 Du u Dv - 2 Du u Dv v - 2 Du u v + 2 Dv v Du + 2 Dv v Du u + 2 Dv v u + Du
2 2
+ Du u + u - Dv - Dv v - v,
2 2 2 2 2 2
(v - 1) (v + 1) (u - 1) (u + 1) Dx + u v - 2 Dv v + 2 Dv u v + 2 Du u v
2 2
- v + v u - u - 2 Du u ]
The previous operations are known as Zeilberger's method of creative telescoping. They
are usually followed by the substitution Dx=0. This would correspond to cases where
h(+1)=h(-1)=0. Here, we need to take h(+1)-h(-1) into account. Let us give the
following notation.
> C[0]:=coeff(EL[2],Dx,0);
2 2 2 2 2 2 2 2
C[0] := u v - 2 Dv v + 2 Dv u v + 2 Du u v - v + v u - u - 2 Du u
> C[1]:=coeff(EL[2],Dx,1);
C[1] := (v - 1) (v + 1) (u - 1) (u + 1)
By integration, EL[2] yields the following non-homogeneous equation.
> C[0]*Int(h,x=-1..1)=-C[1]*('h'(x,1)-'h'(x,-1));
2 2 2 2 2 2 2 2
(u v - 2 Dv v + 2 Dv u v + 2 Du u v - v + v u - u - 2 Du u )
1
/
| 1
| --------------------------------------- dx =
| 2 1/2 2 1/2
/ (1 - 2 x u + u ) (1 - 2 x v + v )
-1
- (v - 1) (v + 1) (u - 1) (u + 1) (h(x, 1) - h(x, -1))
Let us compute C[1]*(h(x,1)-h(x,-1)).
> A:=weylalg([Du,u],[Dv,v]):\
simplify(C[1]*normal(subs(x=1,h)-subs(x=-1,h)),symbolic);
2 u + 2 v
Generally, this difference has to be holonomic (i.e., it is always possible to get equations
that describe it using the theory of holonomy). In the present case, it is a polynomial. We
thus get simple linear operators that cancel on them.
> hypergeomtoholon(",A);
[(- u - v) Du + 1, (- u - v) Dv + 1]
Applying these operators to the right-hand side of the integral equation above yields zero.
Therefore, it is the same when we apply them to the left-hand side. Composition of
operators translates into product.
> G:=map(collect,[EL[1],op(map(opprod,",C[0],A))],[Du,Dv],distributed,factor);
2
G := [u - v + 2 (u v - 1) (u - v) Du Dv + (- v + 2 u v - 1) Dv
2
+ (- 2 u v + 1 + u ) Du,
2 2 2 2
- v (u + v) - 2 v (u - 1) (u + 1) (u + v) Du Dv - 2 v (u + 2 u v + 1) Dv
2 2 2 2 3 3
+ (6 u v - 4 u v + 3 u + v - v u - 5 u v ) Du
2 2
- 2 u (v - 1) (v + 1) (u + v) Du ,
2 2
- u (u + v) - 2 u (v - 1) (v + 1) (u + v) Du Dv
2 2 2 3 3 2
+ (3 v - 4 u v + 6 u v - u v - 5 v u + u ) Dv
2 2 2 2
- 2 u (1 + v + 2 u v) Du - 2 v (u - 1) (u + 1) (u + v) Dv ]
Recall that these operators cancel on the integral psi(u,v) we want to determine. Now,
we compute a Groebner basis for another purpose, namely to get a normal form for this
set of equations.
> T:=termorder(A,tdeg=[Du,Dv],max):\
BR:=ratpoly(rational,[u,v]):\
GL[psi]:=collect(gbasis(G,T,BR),[Du,Dv],distributed,factor);
2 2 3
GL[psi] := [- u v (u - v) (u + v) - v (- v - 3 v u + 3 u + u ) Dv
3 2 2 2 3 2
- u (- 5 u v + 4 v + 3 u v + 4 v u - 3 u v - 3 u ) Du
2 2
- 2 u (u - v) (u v - 1) (u + v) Du ,
2
- u + v - 2 (u v - 1) (u - v) Du Dv + (v + 1 - 2 u v) Dv
2
+ (- u - 1 + 2 u v) Du,
3 2 2 3 2 2
u v (u - v) (u + v) + v (5 v u - 3 u v + 3 u v - 4 u v - 4 u + 3 v ) Dv
2 2 3 2 2
+ u (- 3 v + 3 u v - v + u) Du + 2 v (u - v) (u v - 1) (u + v) Dv
]
To solve this system in psi(u,v), we separate variables. To compute an equation that
involves derivatives with respect to u only, we use the function dependency, that searches
for a dependency between successive derivatives with respect to u.
> eq_Du:=collect(dependency(psi(u,v),u,3,GL,T,BR),Du,factor);
2 2 3 3 5
eq_Du := - 4 u (- v - 3 v u + 3 u + u ) (u - v) (u v - 1) Du - 4 u (5 u v
4 2 4 3 3 3 2 2 2 3
- 21 u v - 3 u + 12 u v + 34 v u - 29 u v - 15 u + 6 u v + 15 u v
2 2 5 4 2 4 3 3 3 2 2
- 4 v ) Du + (- 17 u v + 78 u v + 3 u - 21 u v - 129 v u + 89 u v
2 3 2
+ 45 u - 27 u v - 30 u v + 9 v ) Du
4 3 2 2 2 2
- v (u - 6 v u - 3 u v + 15 u - 10 u v + 3 v )
For symmetry reasons, a similar equation could be obtained for derivatives with respect to
v.
At this stage, we have to solve to get a solution with parameter v. Maple's dsolve does
not return any solution. So we look for a function of the form sigma(uv). We change
variables by z=uv.
> collect(numer(subs([u=z/v,Du=v*Dz],eq_Du)),Dz,factor);
2 4 2 2 2 3 2 3 5
- 4 z v (v + 3 z v - 3 z v - z ) (- z + v ) (z - 1) Dz - 4 z v (5 z
4 2 4 3 4 3 2 2 4 2 2 6
- 21 z v - 3 z + 12 z v + 34 z v - 29 z v - 15 z v + 6 z v
4 6 2 5 4 2 4 3 4 3 2
+ 15 z v - 4 v ) Dz - (17 z - 78 z v - 3 z + 21 z v + 129 z v
2 4 2 2 6 4 6
- 89 z v - 45 z v + 27 z v + 30 z v - 9 v ) v Dz
4 3 2 2 4 2 2 4 6
- v (z - 6 z v - 3 z v + 15 z v - 10 z v + 3 v )
Now, we are working in the (z,v) plane, so we set v=1 to get an equation satisfied by any
solutions of the form sigma(uv).
> collect(primpart(subs(v=1,"),Dz),Dz,factor);
2 2 3 2 2
- 4 z (z - 1) Dz - 4 z (5 z - 4) (z - 1) Dz + (- 17 z + 30 z - 9) Dz - z
+ 3
Maple is now able to solve.
> A:=weylalg([Dz,z]):\
dsol:=dsolve({applyopr("",sigma(z),A)},sigma(z));
1/2
_C1 _C2 arctanh(z ) _C3 (- ln(z) + 2 ln(z - 1))
dsol := sigma(z) = ---- + ----------------- + ---------------------------
1/2 1/2 1/2
z z z
Let us compute an asymptotic series at 0 to be able to use initial conditions in 0, so as to
determine the integration constants _C1, _C2 and _C3.
> series(op(2,dsol),z,2);
_C1 + _C3 (ln(1/z) - 2 I csgn(I (z - 1)) Pi) 1/2
-------------------------------------------- + _C2 - 2 _C3 z + 1/3 _C2 z
1/2
z
3/2
+ O(z )
Since <P[0](x),P[0](x)>=2, i.e., sigma(z)=2, we find the generating function of the scalar
products.
> sigma:=2*arctanh(sqrt(z))/sqrt(z):\
subs(z=u*v,sigma);
1/2
arctanh((u v) )
2 -----------------
1/2
(u v)
This proves that the family of Legendre polynomials is orthogonal. It is now possible to
compute the squares of the modules ||P[n](x)||^2=<P[n](x),P[n](x)>. We first reduce the
order of the equation that defines GF. To do so, we use the method of indetermined
coefficients to get a second order differential equation.
> solve({coeffs(numer(normal(a*sigma+b*diff(sigma,z)+c*diff(sigma,z,z))),arctanh(sqrt(z)))},{a,b,c});
2
{b = 5 a z - 3 a, c = 2 a z - 2 a z, a = a}
> DE_sigma:=normal(primpart(subs(",a+b*Dz+c*Dz^2),Dz));
2 2 2
DE_sigma := 1 + 5 Dz z - 3 Dz + 2 Dz z - 2 Dz z
We easily check this equation.
> normal(applyopr(DE_sigma,sigma,A));
0
The following series expansion makes us conjecture that there is a closed form for the
coefficients, namely 2/(2*n+1). We now proceed to prove this conjecture.
> series(sigma,z,10);
2 3 4 5 6 7 8
2 + 2/3 z + 2/5 z + 2/7 z + 2/9 z + 2/11 z + 2/13 z + 2/15 z + 2/17 z
9 19/2
+ 2/19 z + O(z )
The differential equation that defines sigma can easily be translated in terms of
coefficients. To do so, we use Salvy and Zimmermann's gfun package to derive a
recurrence equation satisfied by the coefficients of sigma.
> with(share):
See ?share and ?share,contents for information about the share library
> readshare(gfun,analysis):
> with(gfun):
> diffeqtorec(applyopr(DE_sigma,f(z),A),f(z),u(n));
(2 n + 1) u(n) + (- 2 n - 3) u(n + 1)
And Maple is able to solve this recurrence equation.
> rsolve({",u(0)=2},u(n));
2
-------
2 n + 1
A particular case of calculation of integrals is taking the coefficients of a series by
Cauchy's formula. We exemplify it by computing again the squares of the modules
||P[n](x)||^2=<P[n](x),P[n](x)> in this more direct way.
We first divide the generating function GF of the squares of the modules by x^(n+1).
> A:=orealg(diff=[Dz,z],shift=[Sn,n]):\
numer(normal(applyopr(collect(DE_sigma,Dz,factor),z^(n+1)*f(z),A)/z^(n+1)));
2 / d \ 2 / d \ 2
6 f(z) z + 2 z f(z) n + 7 z f(z) n + 4 |---- f(z)| z n + 9 |---- f(z)| z
\ dz / \ dz /
/ 2 \
| d | 3 2 / d \
+ 2 |----- f(z)| z - 2 f(z) n - 5 f(z) n - 4 |---- f(z)| z n
| 2 | \ dz /
\ dz /
/ 2 \
/ d \ | d | 2
- 7 |---- f(z)| z - 2 |----- f(z)| z - 3 f(z)
\ dz / | 2 |
\ dz /
> collect(subs([diff(f(z),z,z)=Dz^2,diff(f(z),z)=Dz,f(z)=1],"),Dz,factor);
2 2
2 z (z - 1) Dz + z (9 z + 4 z n - 7 - 4 n) Dz + (2 n + 3) (z n - n + 2 z - 1)
We eliminate z since the integration is with respect to this variable. Note the trick of
integrating from 1 to 1! This means integration on any contour from and to 1 on which
the integrand is well defined.
> T:=termorder(A,lexdeg=[[z],[Dz,Sn]],max):\
hdefint(["",z*Sn-1],z=1..1,T);
Mgfun/Holonomy/definiteintegral:
Assume the integrand and its derivatives vanish in 1 and 1
[2 n - 2 Sn n - 3 Sn + 1]
And Maple is able to solve the corresponding recurrence.
> rsolve({applyopr(op("),u(n),A),u(0)=2},u(n));
2
-------
2 n + 1