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 ẇarnings]:=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)ˆ (nk)*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);

G := [( −x^{2} + 1 ) Dx + 2 k + n x − n,
( −2 n^{2} − 4 n + 4 n k − 2 + 4 k − 2 k^{2} ) Sn + n^{2} x + n^{2} + 2 n x + 2 n + x + 1,
( −k^{2} x − k^{2} − 2 k x − 2 k − x − 1 ) Sk + n^{2} x − n^{2} − 2 n k x + 2 n k + k^{2} x − k^{2}]

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 Sn n x + 3 Sn x − n − 1 − Sn^{2} n − 2 Sn^{2}, −n x + Sn n − Dx x^{2} + 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

GL_{ϕ} := [ 2 x y Dx − y^{2} Dx + y − Dx, 2 x y Dy + x − y^{2} 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);

y ϕ( x, y) + ( −y^{2} + 2 x y − 1 ) (\frac{∂}{∂x} ϕ( x, y))

> pdesolve(”,phi(x,y));

ϕ( x, y) = \frac{F1( y)}{\sqrt{−y^{2} + 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 Ḟ1(y).

> numer(normal(applyopr(GL[FG][2],op(2,”),A)));
GL_{FG}_{2} F1( y)

Since the previous expression is identically zero, Ḟ1(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);

ϕ := \frac{1}{\sqrt{1 − 2 x y + y^{2}}}

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);

h := \frac{1}{\sqrt{1 − 2 x u + u^{2}} \sqrt{1 − 2 x v + v^{2}}}

> G:=hypergeomtoholon(h,A);

G := [( −1 + 2 x v − v^{2} + 2 x u − 4 x^{2} u v + 2 x u v^{2} − u^{2} + 2 u^{2} x v − u^{2} v^{2} ) Dx + u
− 4 u x v + u v^{2} + v + v u^{2}, ( −1 + 2 x u − u^{2} ) Du + x − u,
( −1 + 2 x v − v^{2} ) 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 Du u Dv − 2 Du u Dv v^{2} − 2 Du u v + 2 Dv v Du + 2 Dv v Du u^{2} + 2 Dv v u
+ Du + Du u^{2} + u − Dv − Dv v^{2} − v, ( v − 1 ) ( v + 1 ) ( u − 1 ) ( u + 1 ) Dx + u v^{2}
− 2 Dv v^{2} + 2 Dv u^{2} v^{2} + 2 Du u^{2} v^{2} − v + v u^{2} − u − 2 Du u^{2}]

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);
C_{0} := u v^{2} − 2 Dv v^{2} + 2 Dv u^{2} v^{2} + 2 Du u^{2} v^{2} − v + v u^{2} − u − 2 Du u^{2}

> 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));

( u v^{2} − 2 Dv v^{2} + 2 Dv u^{2} v^{2} + 2 Du u^{2} v^{2} − v + v u^{2} − u − 2 Du u^{2} )

∫_{−1}^{1} \frac{1}{\sqrt{1 − 2 x u + u^{2}} \sqrt{1 − 2 x v + v^{2}}} dx =

−( 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);

G := [u − v + 2 ( u v − 1 ) ( u − v) Du Dv + ( −v^{2} + 2 u v − 1 ) Dv + ( −2 u v + 1 + u^{2} ) Du,

−v ( u + v)^{2} − 2 v^{2} ( u − 1 ) ( u + 1 ) ( u + v) Du Dv − 2 v^{2} ( u^{2} + 2 u v + 1 ) Dv
+ ( 6 u v − 4 u^{2} v^{2} + 3 u^{2} + v^{2} − v u^{3} − 5 u v^{3} ) Du
− 2 u^{2} ( v − 1 ) ( v + 1 ) ( u + v) Du^{2}, −u ( u + v)^{2}
− 2 u^{2} ( v − 1 ) ( v + 1 ) ( u + v) Du Dv
+ ( 3 v^{2} − 4 u^{2} v^{2} + 6 u v − u v^{3} − 5 v u^{3} + u^{2} ) Dv − 2 u^{2} ( 1 + v^{2} + 2 u v) Du
− 2 v^{2} ( u − 1 ) ( u + 1 ) ( u + v) Dv^{2}]

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);

GL_{ψ} := [−u v ( u − v) ( u + v) − v^{2} ( −v − 3 v u^{2} + 3 u + u^{3} ) Dv
− u ( −5 u v^{3} + 4 v^{2} + 3 u^{2} v^{2} + 4 v u^{3} − 3 u v − 3 u^{2} ) Du
− 2 u^{2} ( u − v) ( u v − 1 ) ( u + v) Du^{2},
−u + v − 2 ( u v − 1 ) ( u − v) Du Dv + ( v^{2} + 1 − 2 u v) Dv + ( −u^{2} − 1 + 2 u v) Du,
u v ( u − v) ( u + v) + v ( 5 v u^{3} − 3 u^{2} v^{2} + 3 u v − 4 u v^{3} − 4 u^{2} + 3 v^{2} ) Dv
+ u^{2} ( −3 v + 3 u v^{2} − v^{3} + u) Du + 2 v^{2} ( u − v) ( u v − 1 ) ( u + v) Dv^{2}]

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 Ḋu:=collect(dependency(psi(u,v),u,3,GL,T,BR),Du,factor);

eq Du := −4 u^{2} ( −v − 3 v u^{2} + 3 u + u^{3} ) ( u − v) ( u v − 1 ) Du^{3} − 4u(5 u^{5} v − 21 u^{4} v^{2}
− 3 u^{4} + 12 u^{3} v^{3} + 34 v u^{3} − 29 u^{2} v^{2} − 15 u^{2} + 6 u v^{3} + 15 u v − 4 v^{2})Du^{2} + (
−17 u^{5} v + 78 u^{4} v^{2} + 3 u^{4} − 21 u^{3} v^{3} − 129 v u^{3} + 89 u^{2} v^{2} + 45 u^{2} − 27 u v^{3}
− 30 u v + 9 v^{2})Du − v ( u^{4} − 6 v u^{3} − 3 u^{2} v^{2} + 15 u^{2} − 10 u v + 3 v^{2} )

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 Ḋu)),Dz,factor);

−4 z^{2} v ( v^{4} + 3 z^{2} v^{2} − 3 z v^{2} − z^{3} ) ( −z + v^{2} ) ( z − 1 ) Dz^{3} − 4z v(5 z^{5} − 21 z^{4} v^{2} − 3 z^{4}
+ 12 z^{3} v^{4} + 34 z^{3} v^{2} − 29 z^{2} v^{4} − 15 z^{2} v^{2} + 6 z v^{6} + 15 z v^{4} − 4 v^{6})Dz^{2} − (17 z^{5}
− 78 z^{4} v^{2} − 3 z^{4} + 21 z^{3} v^{4} + 129 z^{3} v^{2} − 89 z^{2} v^{4} − 45 z^{2} v^{2} + 27 z v^{6} + 30 z v^{4}
− 9 v^{6})v Dz − v ( z^{4} − 6 z^{3} v^{2} − 3 z^{2} v^{4} + 15 z^{2} v^{2} − 10 z v^{4} + 3 v^{6} )

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);
−4 z^{2} ( z − 1 )^{2} Dz^{3} − 4 z ( 5 z − 4 ) ( z − 1 ) Dz^{2} + ( −17 z^{2} + 30 z − 9 ) Dz − z + 3

Maple is now able to solve.

> A:=weylalg([Dz,z]):
> dsol:=dsolve(–applyopr(””,sigma(z),A)˝ ,sigma(z));
dsol := σ( z) = \frac{C1}{\sqrt{z}} + C2 arctanh (\sqrt{z}) ; \sqrt{z} + \frac{C3 ( −ln( z) + 2 ln( z − 1 ) )}{\sqrt{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 Ċ1, Ċ2 and Ċ3.

> series(op(2,dsol),z,2);

\frac{C1 + C3 (ln (\frac{1}{z}) − 2 I csgn( I ( z − 1 ) ) π)}{\sqrt{z}} + C2 − 2 C3 \sqrt{z} + \frac{1}{3} C2 z + O( z^{3/2} )

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);

2 \begin{matrix}arctanh & (\sqrt{u v}) \\ √ & \overline{u v} \end{matrix}

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))),ar
> ctanh(sqrt(z)))˝ ,–a,b,c˝ );
{ b = 5 a z − 3 a, c = 2 a z^{2} − 2 a z, a = a }

> DE ṡigma:=normal(primpart(subs(”,a+b*Dz+c*Dzˆ 2),Dz));
DE sigma := 1 + 5 Dz z − 3 Dz + 2 Dz^{2} z^{2} − 2 Dz^{2} z

We easily check this equation.

> normal(applyopr(DE ṡigma,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 + \frac{2}{3} z + \frac{2}{5} z^{2} + \frac{2}{7} z^{3} + \frac{2}{9} z^{4} + \frac{2}{11} z^{5} + \frac{2}{13} z^{6} + \frac{2}{15} z^{7} + \frac{2}{17} z^{8} + \frac{2}{19} z^{9} + O( z^{19/2} )

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 l

ibrary

> readshare(gfun,analysis):
> with(gfun):
> diffeqtorec(applyopr(DE ṡigma,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 \frac{1}{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 ṡigma,Dz,factor),zˆ (n+1)*f(z),A)/zˆ (n+1)
> ));

6 f( z) z + 2 z f( z) n^{2} + 7 z f( z) n + 4 (\frac{∂}{∂z} f( z)) z^{2} n + 9 (\frac{∂}{∂z} f( z)) z^{2} + 2 (\frac{∂^{2}}{∂z^{2}} f( z)) z^{3}

− 2 f( z) n^{2} − 5 f( z) n − 4 (\frac{∂}{∂z} f( z)) z n − 7 (\frac{∂}{∂z} f( z)) z − 2 (\frac{∂^{2}}{∂z^{2}} f( z)) z^{2} − 3 f( z)

> collect(subs([diff(f(z),z,z)=Dzˆ 2,diff(f(z),z)=Dz,f(z)=1],”),Dz,factor);

2 z^{2} ( z − 1 ) Dz^{2} + 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 \frac{1}{2 n + 1}