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