share/contrib/meijer_gv2.mac

   1 /* Copyright 2007 by Edmond Orignac.
   2  * This file is released under the terms of the GNU General Public License, version 2.
   3  */
   4
   5 /* MacRobert E-function  as defined in Bateman Manuscript Project.
   6 a=[a_1,...a_n]
   7 b=[b_1,...b_m]
   8 */
   9
  10     /* Quoting Edmond Orignac (the author) from an email 2007-03-26:
  11
  12     "The second script [i.e. this one] attempts to reduce the E-function of MacRobert
  13     and the G-function of Meijer to simpler formulas with generalized
  14     hypergeometric functions or Bessel/Trig/exponentials/powers.
  15     The E-function and G-function generalize hypergeometric functions
  16     and are described in the Bateman Manuscript Project (vol.1).
  17     The simplification  is done by application of some simplification
  18     formulas valid for some special cases and
  19     by application of Slater's theorem when no two indices differ by
  20     an integer. I have done some simple tests
  21     comparing the output of the program with the formulas from
  22     functions.wolfram.com and the result was satisfactory.
  23     However, my script remains incomplete."
  24
  25     "First, I have not succeeded to implement some of the formulas
  26     quoted by CS Meijer in the Compositio Mathematicae (available
  27     at http://www.numdam.com) article mentioned in the comments
  28     of my program. The problem is simply that I could not devise
  29     a simple way to find the indices of two elements in a list of 4
  30     coefficients whose sum would be zero.
  31     I was thinking of using nested loops with a while condition
  32     as I would do in a C program but I was unsure how to exit
  33     from the nested loops while preserving the values of the
  34     two loop indices. Maybe there is a more natural way of doing
  35     this in Maxima, that does not use nested loops at all, but
  36     I did not find it."
  37
  38     "Second, there exists an article by Kelly Roach (presented at
  39     Issac in 1997) and available from citeseer at Penn State
  40     which is discussing another method based on contiguity relations
  41     to simplify Meijer G-functions. I have not understood how exactly
  42     the contiguity relations are used in that case, but I think that
  43     quite a few extra cases could be treated by this method."
  44
  45      */
  46
  47 macrobert_ered(a,b,x):=block([aa,bb,c,d,y],
  48 c:[1],
  49 d:[],
  50 aa:b,
  51 bb:a,
  52 y:meijer_gred(x,bb,c,aa,d),
  53 return(y))$
  54
  55
  56
  57 /* Meijer G-function as defined in Bateman Manuscript Project.
  58 a=[a_1,...a_n]
  59 b=[b_1,...b_m]
  60 c=[a_{n+1},...,a_p]
  61 d=[b_{m+1},...,b_q]
  62
  63 Some formulas have been taken from http://functions.wolfram.com,
  64 page on Meijer G-function.
  65
  66
  67 More formulas have to be implemented.
  68 */
  69 meijer_gred(x,b,a,c,d):=block([m,n,p,q,sa,sb,sc,sd,s1,s2,i,j,k,k1,k2,nu],
  70  if (not(listp(a) and listp(b) and listp(c) and
  71 listp(d))) then return(false),
  72 s1:[],
  73 s2:[],
  74
  75 /* We put all the elements common to lists a and d in s1 */
  76 for i in a do if member(i,d) then s1:cons(i,s1),
  77 /* We put all  the elements common to lists b and c in s2 */
  78 for i in b do if member(i,c) then s2:cons(i,s2),
  79
  80
  81 /* We delete the common elements in lists a and d */
  82 for i in s1 do if (member(i,a) and member(i,d)) then block([],
  83 a:delete(i,a,1),
  84 d:delete(i,d,1)
  85 ),
  86
  87 /* We delete the common elements in lists b and c */
  88 for i in s2 do if (member(i,b) and member(i,c)) then block([],
  89 b:delete(i,b,1),
  90 c:delete(i,c,1)
  91 ),
  92
  93 /* In case a check of the contents of the lists would be necessary
  94 print (b),
  95 print (a),
  96 print (c),
  97 print (d),
  98 */
  99
 100
 101 m:length(b),
 102 n:length(a),
 103 p:n+length(c),
 104 q:m+length(d),
 105
 106
 107
 108 /* The case of m=0 */
 109 if (m>0) then go(m1),
 110
 111 if (n=0 and p=0 and q=0) then return(delta(x-1)),
 112
 113 /* The case of m=0, n=1 */
 114
 115 if (n>1) then go (m1),
 116 if (p=1 and q=0) then return(x^(first(a)-1)*exp(-1/x)),
 117 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)),
 118 if (p=2 and q=0) then return (x^((first(a)+first(c))/2-1)*bessel_j(first(c)-first(a),2/sqrt(x))),
 119 if (p=2 and q=1) then return (hgfred([1-first(a)+first(d)],[1-first(a)+first(c)],1/x)),
 120 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)),
 121 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)),
 122 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))),
 123
 124
 125 /* The case of m=1 */
 126
 127 m1,
 128 if (m>1) then go (m2),
 129 if (p=0 and q=1) then return(x^first(b)*exp(-x)),
 130 if (p=0 and q=2) then return (x^((first(b)+first(d))/2)*bessel_j(first(b)-first(d),2*sqrt(x))),
 131 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)),
 132
 133 m2,
 134 /* formula from L. T. Wille J. Math. Phys. 29 p. 599 (1987) */
 135 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))),
 136
 137 /* Formulas from C. S. Meijer Compositio Mathematica 8, p. 49 (1951) */
 138 if (m#2 or p#0 or q#4) then go(bsl2),
 139 if (expand((b[1]-b[2])^2)#1/4 or expand((d[1]-d[2]))^2#1/4) then go(bsl2),
 140 if (expand(b[1]+b[2]+d[1]+d[2])#1) then go(bsl1a),
 141 nu:expand(b[1]+b[2]-1/2),
 142 return(bessel_j(2*nu,4*x^(1/4))),
 143 bsl1a,
 144 if (expand(b[1]+b[2]+d[1]+d[2])#0) then go(bsl2),
 145 nu:expand(b[1]+b[2]-1/2),
 146 return(bessel_j(2*nu+1,4*x^(1/4))/x^(1/4)),
 147
 148 bsl2,
 149
 150
 151
 152
 153 /* 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. */
 154
 155 if ((p<=q) and not(test_int_diff(b))) then return (block([y,z,k],
 156
 157     z:0,
 158     for k:1 thru m do z:z+block([],
 159       w:x^(b[k]),
 160       for j:1 thru length(c) do w:w/gamma(c[j]-b[k]),
 161       for j:1 thru length(a) do w:w*gamma(1+b[k]-a[j]),
 162       for j:1 thru (k-1) do w:w*gamma(b[j]-b[k]),
 163       for j:(k+1) thru m do w:w*gamma(b[j]-b[k]),
 164       for j:1 thru length(d) do w:w/gamma(1+b[k]-d[j]),
 165       sa:[],
 166       sb:[],
 167       for j:1 thru length(a) do sa:cons(1+b[k]-a[j],sa),
 168       for j:1 thru length(c) do sa:cons(1+b[k]-c[j],sa),
 169       for j:1 thru (k-1) do sb:cons(1+b[k]-b[j],sb),
 170       for j:(k+1) thru length(b) do sb:cons(1+b[k]-b[j],sb),
 171       for j:1 thru length(d) do sb:cons(1+b[k]-d[j],sb),
 172       return (w*hgfred(sa,sb,(-1)^(m+n-p)*x))),
 173     return(z)
 174       )),
 175
 176
 177 /* We try to apply Slater's theorem after analytic continuation when p>q  */
 178
 179 if ((p>q) and not(test_int_diff(a))) then return (block([y,z,k],
 180
 181     z:0,
 182     for k:1 thru n do z:z+block([j,w],
 183       w:x^(a[k]-1),
 184       for j:1 thru length(d) do w:w/gamma(a[k]-d[j]),
 185       for j:1 thru length(b) do w:w*gamma(1+b[j]-a[k]),
 186       for j:1 thru (k-1) do w:w*gamma(a[k]-a[j]),
 187       for j:(k+1) thru n do w:w*gamma(a[k]-a[j]),
 188       for j:1 thru length(c) do w:w/gamma(1+c[j]-a[k]),
 189       sa:[],
 190       sb:[],
 191       for j:1 thru length(b) do sa:cons(1+b[j]-a[k],sa),
 192       for j:1 thru length(d) do sa:cons(1+d[j]-a[k],sa),
 193       for j:1 thru (k-1) do sb:cons(1+a[j]-a[k],sb),
 194       for j:(k+1) thru length(a) do sb:cons(1+a[j]-a[k],sb),
 195       for j:1 thru length(c) do sb:cons(1+c[j]-a[k],sb),
 196       return (w*hgfred(sa,sb,(-1)^(m+n-q)/x))),
 197     return(z)
 198       )),
 199
 200
 201 /* Catch all case */
 202
 203
 204
 205 return (meijer_G(x,b,a,c,d)))$
 206
 207
 208 test_int_diff(a):=block([i,j,n,z],
 209   if (not(listp(a))) then return(false),
 210   n:length(a),
 211   if (n=1) then return(false),
 212   z:false,
 213   for i:1 thru n do for j:i+1 thru n do z:(z or integerp(a[j]-a[i])),
 214   return(z))$