[maxima program] Meijer's G function and MacRobert's E function
Subject: [maxima program] Meijer's G function and MacRobert's E function
From: Edmond Orignac
Date: Sat, 24 Mar 2007 14:59:37 +0100
Dear Maxima Users,
Below find as attachment a program to transform
the Meijer G function and the MacRobert E function
into simpler expressions by application of the
theorem of L. J. Slater and some miscellaneous
identities. Definitions for E and G follow
the Bateman Manuscript project.
The theorem of Slater allows the E and G functions
to be reduced to sums of generalized hypergeometric
functions that are themselves simplified by the hgfred
function of Maxima. I also used two miscellaneous
identities from CS Meijer to simplify some expressions
in terms of Bessel functions.
To use these functions, one has to enter:
meijer_gred(x,[a1...an],[b1...bm],[bm+1...bq],[an+1...ap])
or
macrobert_ered([a1...ap],[r1...rq],x)
The program is called meijer_gv2.mac
Best,
Edmond Orignac
-------------- next part --------------
/* MacRobert E-function as defined in Bateman Manuscript Project.
a=[a_1,...a_n]
b=[b_1,...b_m]
*/
macrobert_ered(a,b,x):=block([aa,bb,c,d,y],
c:[1],
d:[],
aa:b,
bb:a,
y:meijer_gred(x,bb,c,aa,d),
return(y))$
/* Meijer G-function as defined in Bateman Manuscript Project.
a=[a_1,...a_n]
b=[b_1,...b_m]
c=[a_{n+1},...,a_p]
d=[b_{m+1},...,b_q]
Some formulas have been taken from http://functions.wolfram.com,
page on Meijer G-function.
More formulas have to be implemented.
*/
meijer_gred(x,b,a,c,d):=block([m,n,p,q,sa,sb,sc,sd,s1,s2,i,j,k,k1,k2,nu],
if (not(listp(a) and listp(b) and listp(c) and
listp(d))) then return(false),
s1:[],
s2:[],
/* We put all the elements common to lists a and d in s1 */
for i in a do if member(i,d) then s1:cons(i,s1),
/* We put all the elements common to lists b and c in s2 */
for i in b do if member(i,c) then s2:cons(i,s2),
/* We delete the common elements in lists a and d */
for i in s1 do if (member(i,a) and member(i,d)) then block([],
a:delete(i,a,1),
d:delete(i,d,1)
),
/* We delete the common elements in lists b and c */
for i in s2 do if (member(i,b) and member(i,c)) then block([],
b:delete(i,b,1),
c:delete(i,c,1)
),
/* In case a check of the contents of the lists would be necessary
print (b),
print (a),
print (c),
print (d),
*/
m:length(b),
n:length(a),
p:n+length(c),
q:m+length(d),
/* The case of m=0 */
if (m>0) then go(m1),
if (n=0 and p=0 and q=0) then return(delta(x-1)),
/* The case of m=0, n=1 */
if (n>1) then go (m1),
if (p=1 and q=0) then return(x^(first(a)-1)*exp(-1/x)),
if (p=1 and q=1) then return (x^(first(a)-1)*(1-1/x)^(first(a)-first(d) -1)/gamma(first(a)-first(d))*theta(abs(x)-1)),
if (p=2 and q=0) then return (x^((first(a)+first(c))/2-1)*bessel_j(first(c)-first(a),2/sqrt(x))),
if (p=2 and q=1) then return (hgfred([1-first(a)+first(d)],[1-first(a)+first(c)],1/x)),
if (p=2 and q=2) then return (hgfred([1-first(a)+first(d),1-first(a)+last(d)],[1-first(a)+first(c)],-1/x)),
if (p=3 and q=0) then return (x^(first(a)-1)*hgfred([],[1-first(a)+first(c)],[1-first(a)+last(c)],-1/x)),
if (p=3 and q=1) then return (x^(first(a)-1)*hgfred([1-first(a)+first(b)],[1-first(a)+first(c),1-first(a)+last(c)],1/x)/gamma(first(a)-first(b))),
/* The case of m=1 */
m1,
if (m>1) then go (m2),
if (p=0 and q=1) then return(x^first(b)*exp(-x)),
if (p=0 and q=2) then return (x^((first(b)+first(d))/2)*bessel_j(first(b)-first(d),2*sqrt(x))),
if (m=1 and n=1 and p=1 and q=1) then return(gamma(1-first(a)+first(b))*x^(first(b))*(1+x)^(first(a)-first(b)-1)),
m2,
/* formula from L. T. Wille J. Math. Phys. 29 p. 599 (1987) */
if (m=2 and n=0 and p=1 and q=2 and d[1]=1 and ((b[1]=0 and b[2]=1/2) or (b[1]=1/2 and b[2]=0))) then return(sqrt(%pi)*(1-erf(x))),
/* Formulas from C. S. Meijer Compositio Mathematica 8, p. 49 (1951) */
if (m#2 or p#0 or q#4) then go(bsl2),
if (expand((b[1]-b[2])^2)#1/4 or expand((d[1]-d[2]))^2#1/4) then go(bsl2),
if (expand(b[1]+b[2]+d[1]+d[2])#1) then go(bsl1a),
nu:expand(b[1]+b[2]-1/2),
return(bessel_j(2*nu,4*x^(1/4))),
bsl1a,
if (expand(b[1]+b[2]+d[1]+d[2])#0) then go(bsl2),
nu:expand(b[1]+b[2]-1/2),
return(bessel_j(2*nu+1,4*x^(1/4))/x^(1/4)),
bsl2,
/* If Slater's theorem is applicable the Meijer G function becomes a sum of power of x times pFq generalized hypergeometric function of x. Note that for p=q we should really apply 07.34.26.0006.01 from Wolfram functions site. Here we are implicitly assuming abs(x)<1. */
if ((p<=q) and not(test_int_diff(b))) then return (block([y,z,k],
z:0,
for k:1 thru m do z:z+block([],
w:x^(b[k]),
for j:1 thru length(c) do w:w/gamma(c[j]-b[k]),
for j:1 thru length(a) do w:w*gamma(1+b[k]-a[j]),
for j:1 thru (k-1) do w:w*gamma(b[j]-b[k]),
for j:(k+1) thru m do w:w*gamma(b[j]-b[k]),
for j:1 thru length(d) do w:w/gamma(1+b[k]-d[j]),
sa:[],
sb:[],
for j:1 thru length(a) do sa:cons(1+b[k]-a[j],sa),
for j:1 thru length(c) do sa:cons(1+b[k]-c[j],sa),
for j:1 thru (k-1) do sb:cons(1+b[k]-b[j],sb),
for j:(k+1) thru length(b) do sb:cons(1+b[k]-b[j],sb),
for j:1 thru length(d) do sb:cons(1+b[k]-d[j],sb),
return (w*hgfred(sa,sb,(-1)^(m+n-p)*x))),
return(z)
)),
/* We try to apply Slater's theorem after analytic continuation when p>q */
if ((p>q) and not(test_int_diff(a))) then return (block([y,z,k],
z:0,
for k:1 thru n do z:z+block([j,w],
w:x^(a[k]-1),
for j:1 thru length(d) do w:w/gamma(a[k]-d[j]),
for j:1 thru length(b) do w:w*gamma(1+b[j]-a[k]),
for j:1 thru (k-1) do w:w*gamma(a[k]-a[j]),
for j:(k+1) thru n do w:w*gamma(a[k]-a[j]),
for j:1 thru length(c) do w:w/gamma(1+c[j]-a[k]),
sa:[],
sb:[],
for j:1 thru length(b) do sa:cons(1+b[j]-a[k],sa),
for j:1 thru length(d) do sa:cons(1+d[j]-a[k],sa),
for j:1 thru (k-1) do sb:cons(1+a[j]-a[k],sb),
for j:(k+1) thru length(a) do sb:cons(1+a[j]-a[k],sb),
for j:1 thru length(c) do sb:cons(1+c[j]-a[k],sb),
return (w*hgfred(sa,sb,(-1)^(m+n-q)/x))),
return(z)
)),
/* Catch all case */
return (meijer_G(x,b,a,c,d)))$
test_int_diff(a):=block([i,j,n,z],
if (not(listp(a))) then return(false),
n:length(a),
if (n=1) then return(false),
z:false,
for i:1 thru n do for j:i+1 thru n do z:(z or integerp(a[j]-a[i])),
return(z))$