tests/rtest14.mac

   1 jn(3,4);
   2 jn(3,4); /*0.1320341839226408 ;*/
   3 j0(1);
   4 j0(1); /*0.7651976865579665 ;*/
   5 bessel_j(0,1.0);
   6 0.7651976865579665 ;
   7
   8 /* BESSEL is gone. This makes sure it's gone. */
   9 bessel(2,3);
  10 bessel(2,3)$
  11
  12 bessel(2,3.0);
  13 bessel(2, 3.0);
  14
  15 bessel_j(3,2.0);
  16 0.1289432494744021;
  17
  18 bessel_j(-5,x);
  19 -bessel_j(5,x);
  20 bessel_j(-n,x);
  21 bessel_j(-n,x);
  22
  23 diff(bessel_j(0,x),x);
  24 -bessel_j(1,x);
  25 diff(bessel_j(1,x),x);
  26 (bessel_j(0,x)-bessel_j(2,x))/2;
  27 diff(bessel_y(0,x),x);
  28 -bessel_y(1,x);
  29 diff(bessel_y(1,x),x);
  30 (bessel_y(0,x)-bessel_y(2,x))/2;
  31 diff(bessel_i(0,x),x);
  32 bessel_i(1,x);
  33 diff(bessel_i(1,x),x);
  34 (bessel_i(2,x)+bessel_i(0,x))/2;
  35 diff(bessel_k(0,x),x);
  36 -bessel_k(1,x);
  37 diff(bessel_k(1,x),x);
  38 -((bessel_k(2,x)+bessel_k(0,x))/2);
  39
  40 besselexpand:true;
  41 true$
  42
  43 bessel_j(1/2,x);
  44 sqrt(2/%pi)*sin(x)/sqrt(x);
  45 bessel_j(-1/2,x);
  46 sqrt(2/%pi)*cos(x)/sqrt(x);
  47
  48 bessel_j(3/2,x);
  49 (sqrt(2)*sqrt(x)*(sin(x)/x^2-cos(x)/x))/sqrt(%pi);
  50
  51 bessel_j(5/2,x);
  52 (sqrt(2)*sqrt(x)*((3/x^3-1/x)*sin(x)-(3*cos(x))/x^2))/sqrt(%pi)$
  53
  54 ratsimp(bessel_j(-3/2,x) - (2*(-1/2)/x*bessel_j(-1/2,x)-bessel_j(1/2,x)));
  55 0$
  56
  57 ratsimp(bessel_j(-5/2,x) - (2*(-3/2)/x*bessel_j(-3/2,x)-bessel_j(-1/2,x)));
  58 0$
  59
  60 ratsimp(bessel_y(1/2,x) + sqrt(2/%pi)*cos(x)/sqrt(x));
  61 0$
  62
  63 ratsimp(bessel_y(-1/2,x) - sqrt(2/%pi)*sin(x)/sqrt(x));
  64 0$
  65
  66 ratsimp(bessel_y(3/2,x) - (2*(1/2)/x*bessel_y(1/2,x)-bessel_y(-1/2,x)));
  67 0$
  68
  69 ratsimp(bessel_y(5/2,x) - (2*(3/2)/x*bessel_y(3/2,x)-bessel_y(1/2,x)));
  70 0$
  71
  72 ratsimp(bessel_y(-3/2,x) - (2*(-1/2)/x*bessel_y(-1/2,x)-bessel_y(1/2,x)));
  73 0$
  74
  75 ratsimp(bessel_y(-5/2,x) - (2*(-3/2)/x*bessel_y(-3/2,x)-bessel_y(-1/2,x)));
  76 0$
  77
  78 bessel_i(1/2,x);
  79 sqrt(2*x/%pi)*sinh(x)/x;
  80
  81 bessel_i(-1/2,x);
  82 sqrt(2*x/%pi)*cosh(x)/x;
  83
  84 ratsimp(bessel_i(3/2,x) - (-2*(1/2)/x*bessel_i(1/2,x)+bessel_i(-1/2,x)));
  85 0$
  86
  87 ratsimp(bessel_i(5/2,x) -(-2*(3/2)/x*bessel_i(3/2,x)+bessel_i(1/2,x)));
  88 0$
  89
  90 ratsimp(bessel_i(-3/2,x) -(2*(-1/2)/x*bessel_i(-1/2,x)+bessel_i(1/2,x)));
  91 0$
  92
  93 ratsimp(bessel_i(-5/2,x) - (2*(-3/2)/x*bessel_i(-3/2,x)+bessel_i(-1/2,x)));
  94 0$
  95
  96 ratsimp(bessel_k(1/2,x) - sqrt(%pi/2)*%e^(-x)/sqrt(x));
  97 0$
  98
  99 ratsimp(bessel_k(-1/2,x)- sqrt(%pi/2)*%e^(-x)/sqrt(x));
 100 0$
 101
 102 ratsimp(bessel_k(3/2,x) - (2*(1/2)/x*bessel_k(1/2,x)+bessel_k(-1/2,x)));
 103 0$
 104
 105 ratsimp(bessel_k(5/2,x) -(2*(3/2)/x*bessel_k(3/2,x)+bessel_k(1/2,x)));
 106 0$
 107
 108 ratsimp(bessel_k(-3/2,x) -(-2*(-1/2)/x*bessel_k(-1/2,x)+bessel_k(1/2,x)));
 109 0$
 110
 111 ratsimp(bessel_k(-5/2,x) -(-2*(-3/2)/x*bessel_k(-3/2,x)+bessel_k(-1/2,x)));
 112 0$
 113
 114 (assume(p>0),true);
 115 true$
 116 (assume(4*p+a>0),true);
 117 true$
 118 besselexpand:false;
 119 false$
 120
 121 specint(t^(1/2)*%e^(-a*t/4)*%e^(-p*t),t);
 122 sqrt(%pi)/(2*(p+a/4)^(3/2));
 123
 124 prefer_whittaker:true;
 125 true$
 126
 127 /*
 128  * Reference:  Table of Integral Transforms.
 129  */
 130
 131 /*
 132  * f24p146:
 133  *
 134  * t^(v-1)*exp(-t^2/(8*a))
 135  *     -> gamma(v)*2^v*a^(v/2)*exp(a*p^2)*D[-v](2*p*sqrt(a))
 136  *
 137  * We have v = 7/4, a = b/4 so the result should be
 138  *
 139  *   gamma(7/4)*2^(7/4)*(b/4)^(7/8)*exp(b*p^2/4)*D[-7/4](p*sqrt(b))
 140  *     = 3/4*gamma(3/4)*b^(7/8)*exp(b*p^2/4)*D[-7/4](p*sqrt(b))
 141  *
 142  * But
 143  *
 144  *   D[v](z) = 2^(v/2+1/4)*z^(-1/2)*%w[v/2+1/4,1/4](z^2/2)
 145  *
 146  * and
 147  *
 148  *   %w[k,u](z) = gamma(-2*u)/gamma(1/2-u-k)*%m[k,u](z)
 149  *                 + gamma(2*u)/gamma(1/2+u-k)*%m[k,-u](z)
 150  *
 151  * Thus,
 152  *
 153  *   D[-7/4](p*sqrt(b)) = %w[-5/8,1/4](b*p^2/2)/(2^(5/8)*b^(1/4)*sqrt(p))
 154  *
 155  * and
 156  *
 157  *   %w[-5/8,1/4](b*p^2)
 158  *     = 8/3/gamma(3/8)*%m[-5/8,-1/4](b*p^2/2)
 159  *         - 2*sqrt(%pi)/gamma(7/8)*%m[-5/8,1/4](b*p^2/2)
 160  *
 161  * And finally the transform is
 162  *
 163  * 3*gamma(3/4)*b^(5/8)*exp(b*p^2/4)/(4*2^(5/8)*sqrt(p))
 164  *   *(8*sqrt(%pi)/3/gamma(3/8)*%m[-5/8,-1/4](b*p^2/2)
 165  *      - 2*sqrt(%pi)/gamma(7/8)*%m[-5/8,1/4](b*p^2/2))
 166  */
 167 (assume(b>0),true);
 168 true$
 169
 170 specint(t^(3/4)*%e^(-t^2/2/b)*%e^(-p*t),t);
 171
 172 /*
 173 -sqrt(%pi)*b^(5/8)
 174           *(3*gamma(3/8)*gamma(3/4)*%e^(b*p^2/4)*%m[-5/8,1/4](b*p^2/2)
 175            -4*gamma(3/4)*gamma(7/8)*%e^(b*p^2/4)
 176              *%m[-5/8,-1/4](b*p^2/2))
 177  /(2*2^(5/8)*gamma(3/8)*gamma(7/8)*sqrt(p))$
 178  */
 179
 180 3*gamma(3/4)*b^(7/8)
 181  *%e^(b*p^2/4)
 182  *(2^(19/8)*sqrt(%pi)*%m[-5/8,-1/4](b*p^2/2)/(3*gamma(3/8)*b^(1/4)*sqrt(p))
 183   -2^(3/8)*sqrt(%pi)*%m[-5/8,1/4](b*p^2/2)/(gamma(7/8)*b^(1/4)*sqrt(p)))
 184  /4;
 185
 186 /*
 187  * Sec. 4.5, formula (33):
 188  *
 189  * t^(-1/2)*exp(-2*sqrt(a)*sqrt(t)) ->
 190  *    sqrt(%pi)/sqrt(p)*exp(a/p)*erfc(sqrt(a)/sqrt(p))
 191  */
 192 ratsimp(specint(t^(-1/2)*%e^(-2*a^(1/2)*t^(1/2))*%e^(-p*t),t));
 193 -sqrt(%pi)*(erf(sqrt(a)/sqrt(p))-1)*%e^(a/p)/sqrt(p)$
 194
 195 /*
 196  * The Laplace transform of sin(a*t)*cosh(b*t^2) can be derived by
 197  * expressing sin and cosh in terms of exponential functions.  We end
 198  * up with terms of the form:
 199  *
 200  *   exp(+/-%i*a*t)*exp(+/-b*t^2)
 201  *
 202  * All of these can be computed using formula 24, p. 146 of Tables of
 203  * Integral Transforms, which handle functions of the form
 204  * t^(v-1)*exp(-t^2/8/a).
 205  *
 206  * But, we have terms of the form exp(b*t^2-p*t+%i*a*t).  I don't
 207  * think this converges, so the Laplace transform doesn't exist if b >
 208  * 0.
 209  *
 210  */
 211 /*
 212 radcan(specint(sin(a*t)*cosh(b*t^2)*%e^(-p*t),t));
 213 -%e^-((p^2+2*%i*a*p+a^2)/(4*b))*(sqrt(%pi)*%e^((2*%i*a*p+a^2)/(2*b))*erf((%i
 214  *p+a)/(2*sqrt(b)))-sqrt(%pi)*%e^(a^2/(2*b))*erf((%i*p-a)/(2*sqrt(b)))+sqrt
 215  (%pi)*%i*%e^((p^2+2*%i*a*p)/(2*b))*erf((p+%i*a)/(2*sqrt(b)))-sqrt(%pi)*%i*%e
 216  ^(p^2/(2*b))*erf((p-%i*a)/(2*sqrt(b)))+(sqrt(%pi)*%i-sqrt(%pi)*%i*%e^(%i*a
 217  *p/b))*%e^(p^2/(2*b))-sqrt(%pi)*%e^((2*%i*a*p+a^2)/(2*b))+sqrt(%pi)*%e^(a^2/
 218  (2*b)))/(8*sqrt(b)) $
 219 */
 220
 221 /*
 222  * Sec 4.14, formula (27):
 223  *
 224  * t^(1/2)*bessel_j(1,2*a^(1/2)*t^(1/2)) ->
 225  *    sqrt(a)/p^2*exp(-a/p)
 226  */
 227
 228 specint(t^(1/2)*bessel_j(1,2*a^(1/2)*t^(1/2))*%e^(-p*t),t);
 229 sqrt(a)*%e^-(a/p)/p^2$
 230
 231 /*
 232  * Sec 4.14, formula (3):
 233  *
 234  * t^2*bessel_j(v,a*t) ->
 235  *    ((v^2-1)/r^3 + 3*p*(p+v*r)/r^5)*(a/R)^v
 236  *
 237  * where r = sqrt(p^2+a^2), R = p + r;
 238  *
 239  * (Maxima can't currently compute this transform for general v due to a bug
 240  * in hyp.lisp.)
 241  * This bug is no longer present after correction of legf24 in hyp.lisp.
 242  */
 243 factor(ratsimp(specint(t^2*bessel_j(1,a*t)*%e^(-p*t),t)));
 244 3*a*p/(p^2+a^2)^(5/2) $
 245
 246 (/* This is the Laplace transform of the Struve H function, see
 247   http://dlmf.nist.gov/Draft/ST/about_ST.8.13.html */
 248  2/(%pi*p)-2*p*log(p/(sqrt(p^2+1)-1))/(%pi*sqrt(p^2+1)),
 249  /* And this should be the same as the specint of the next test below */
 250  -diff(%%,p),
 251  ev(fullratsimp(%%),logexpand:all));
 252 -((sqrt(p^2+1)*(2*p^2*log(sqrt(p^2+1)-1)-2*p^2*log(p))-2*p^2-2) /(%pi*p^6+2*%pi*p^4+%pi*p^2))$
 253
 254 (ev(fullratsimp(specint(t*struve_h(1,t)*%e^(-p*t),t)),logexpand:all),
 255  ratsimp(%%/%));
 256 1$
 257
 258 /*
 259  *
 260  * From the comments for hstf in hypgeo.lisp:
 261  *
 262  * struve_h(1,t) = 2/sqrt(%pi)*(t/2)^2/gamma(1+3/2)*%f[1,2]([1],[3/2,5/2],-t^2/4)
 263  *
 264  * So
 265  *
 266  * struve_h(1,sqrt(t)) = 2/(3*%pi)*t*%f[1,2]([1],[3/2,5/2],-t/4)
 267  *
 268  * and the integrand is
 269  *
 270  * 2/(3*%pi)*t^(5/2)*exp(-p*t)*%f[1,2]([1],[3/2,5/2],-t/4).
 271  *
 272  * From the f19p220, the Laplace transform of this, with s = 7/2,
 273  * c=-1/4, k = 1, is
 274  *
 275  * 2/(3*%pi)*gamma(7/2)/p^(7/2)*%f[2,2]([1,7/2],[3/2,5/2],-1/4/p)
 276  *
 277  * From the derivation of SPLITPFQ, we can simplify this
 278  * hypergeometric function.
 279  *
 280  * %f[2,2]([1,7/2],[3/2,5/2],z) =
 281  *
 282  *      1
 283  *     sum z^k/poch(5/2,k)*binomial(1,k) *diff(%f[2,2]([1,5/2],[3/2,5/2],z,k)
 284  *     k=0
 285  *
 286  * But %f[2,2]([1,5/2],[3/2,5/2],z) = %f[1,1]([1],[3/2],z)
 287  * and Maxima knows how to compute this.
 288  */
 289 ratsimp(specint(t^(3/2)*struve_h(1,t^(1/2))*%e^(-p*t),t));
 290 -%e^-(1/(4*p))*(sqrt(%pi)*sqrt(p)
 291                                 *(8*%i*erf(%i/(2*sqrt(p)))*p
 292                                  -%i*erf(%i/(2*sqrt(p))))
 293                       -2*p*%e^(1/(4*p)))
 294         /(8*sqrt(%pi)*p^(9/2)) $
 295
 296 /* Trivial result because %ibes is not bessel_i
 297    After specializing the pattern match of arbpow1 we get a more
 298    correct noun form. The result is adjusted. DK.
 299  */
 300 specint(t*%ibes[0](a*t/2)*%ibes[1](a*t/2)*%e^(-p*t),t);
 301 /* %ibes[0](a*t/2)*%ibes[1](a*t/2)/p^2 $ */
 302 specint(t*%ibes[0](a*t/2)*%ibes[1](a*t/2)*%e^(-p*t),t);
 303
 304 /*
 305  * t^(3/4)*bessel_j(1/2,t)*bessel_j(1/4,t)
 306  *
 307  * Luke gives
 308  *
 309  * bessel_j(u,t)*bessel_j(v,t)
 310  *    = (z/2)^(u+v)/gamma(u+1)/gamma(v+1)
 311  *        * %f[2,3]([(u+v+1)/2,(u+v+2)/2],[u+1,v+1,u+v+1],-z^2)
 312  *
 313  * So the integrand is
 314  *
 315  * 8/2^(3/4)/sqrt(%pi)/gamma(1/4)*t^(3/2)*%f[2,3]([7/8,11/8],[3/2,5/4,7/4],-t^2)
 316  *
 317  * f19p220 gives
 318  *
 319  * t^(3/2)*%f[2,3]([7/8,11/8],[3/2,5/4,7/4],-t^2)
 320  *
 321  *    -> gamma(5/2)*p^(-5/2)
 322  *         *%f[4,3]([7/8,11/8,5/2/2,(5/2+1)/2],[3/2,5/4,7/4],-4/p^2)
 323  *    =  gamma(5/2)*p^(-5/2)
 324  *         $%f[2,1]([7/8,11/8],[3/2],-4/p^2)
 325  *
 326  * And we know %f[2,1]([7/8,11/8],[3/2],z) is
 327  *
 328  * -2*(1/(sqrt(z)+1)^(3/4)-1/(1-sqrt(z))^(3/4))/(3*sqrt(z))
 329  *
 330  * Applying all of this gives the expected answer below.
 331  */
 332
 333 specint(t^(3/4)*bessel_j(1/2,t)*bessel_j(1/4,t)*%e^(-p*t),t);
 334 (2^(1/4)*%i*(1/((2*%i)/p+1)^(3/4)-1/(1-(2*%i)/p)^(3/4)))/(gamma(1/4)*p^(3/2))$
 335
 336 /*
 337  * Not sure this is right.  We can convert bessel_y to bessel_j,
 338  * multiply them together and use the results for products of bessel_j
 339  * functions.
 340  *
 341  * bessel_y(1/2,sqrt(t)) = -bessel_j(-1/2,sqrt(t))
 342  *
 343  * And maxima should be able to compute the transform of
 344  *
 345  * t^(5/2)*bessel_j(-1/2,sqrt(t))^2
 346  *
 347  * Or note that bessel_y(1/2,sqrt(t)) =
 348  * -sqrt(2/%pi)*cos(sqrt(t))/t^(1/4).  Then the integrand becomes
 349  *
 350  * 2/%pi*t^2*cos(sqrt(t))^2
 351  *
 352  * And maxima should know how to compute the transform of this.
 353  *
 354  * Unfortunately, the transforms of these two approaches don't agree.
 355  * Yuck!
 356  * After revision 1.65 of hypgeo.lisp it works as expected.
 357  */
 358 result:factor(ratsimp(specint(t^(5/2)*bessel_y(1/2,t^(1/2))^2*%e^(-p*t),t)));
 359 %e^-(1/p)*(16*p^3*%e^(1/p)-18*p^2*%e^(1/p)+4*p*%e^(1/p)
 360                                 +15*sqrt(%pi)*%i*erf(%i/sqrt(p))*p^(5/2)
 361                                 -20*sqrt(%pi)*%i*erf(%i/sqrt(p))*p^(3/2)
 362                                 +4*sqrt(%pi)*%i*erf(%i/sqrt(p))*sqrt(p))
 363        /(4*%pi*p^6);
 364
 365 /* This is equal to the Laplace transform of 2/%pi*t^2*cos(sqt(t))^2 */
 366 expand(result-specint(t^(5/2)*bessel_y(1/2,t^(1/2))^2*%e^(-p*t),t)),
 367    besselexpand:true;
 368  0;
 369
 370 /*
 371  * See formula (42), p. 187:
 372  *
 373  * t^(lam-1)*bessel_j(2*u,2*sqrt(a)*sqrt(t))*bessel_j(2*v,2*sqrt(a)*sqrt(t))
 374  *
 375  *    -> 2*gamma(lam+u+v)*a^(u+v)/gamma(2*u+1)/gamma(2*v+1)/p^(lam+u+v)
 376  *        *%f[3,3]([u+v+1/2,u+v+1,lam+u+v],[2*u+1,2*v+1,2*u+2*v+1],-4*a/p)
 377  *
 378  * with Re(lam + u + v) > 0.
 379  *
 380  * So, we have lam = 3/2, u=v=1/4, a = 1/4, we get
 381  *
 382  * 4/%pi/p^2*%f[3,3]([1,3/2,2],[3/2,3/2,2],-1/p)
 383  *  = 4/%pi/p^2*%f[1,1]([1],[3/2],-1/p)
 384  *
 385  * And %f[1,1]([1],[3/2],-1/p) is
 386  *
 387  *     -sqrt(%pi)*%i*erf(%i*sqrt(1/p))*%e^-(1/p)/(2*sqrt(1/p))
 388  *
 389  * So, the final result is:
 390  *
 391  * -2*%i*erf(%i/sqrt(p))*%e^-(1/p)/(sqrt(%pi)*p^(3/2))
 392  *
 393  * But we also have
 394  *
 395  * bessel_j(u,t)*bessel_j(v,t)
 396  *    = (z/2)^(u+v)/gamma(u+1)/gamma(v+1)
 397  *        * %f[2,3]([(u+v+1)/2,(u+v+2)/2],[u+1,v+1,u+v+1],-z^2)
 398  *
 399  * So bessel_j(1/2,sqrt(t))^2 is
 400  *
 401  *    2/%pi*%f[2,3]([1,3/2],[3/2,3/2,2],-t)*sqrt(t)
 402  *
 403  * So the integrand is
 404  *
 405  *    2/%pi*t*%f[2,3]([1,3/2],[3/2,3/2,2],-t)
 406  *     = 2/%pi*t*%f[1,2]([1],[3/2,2],-t)
 407  *
 408  * f19p220 then gives us the desired transform:
 409  *
 410  *    t*%f[1,2]([1],[3/2,2],-t)
 411  *      -> gamma(2)*p^(-2)*%f[2,2]([1,2],[3/2,2],-1/p)
 412  *
 413  *      = p^(-2)*%f[1,1]([1],[3/2],-1/p)
 414  *
 415  * So the final answer is
 416  *
 417  *    -%i*erf(%i/sqrt(p))*%e^-(1/p)/(sqrt(%pi)*p^(3/2))
 418  *
 419  * Hmm.  This differs from formula 42 above.  I think there's a bug in
 420  * formula 42, and it should be divided by 2.
 421  *
 422  * If we use the expression for the product of Bessel functions and
 423  * f19p220, we can easily derive the result of formula 42, except,
 424  * we're missing the factor of 2.  So, I think formula 42 is wrong.
 425  */
 426
 427 specint(t^(1/2)*bessel_j(1/2,t^(1/2))^2*%e^(-p*t),t);
 428 -%i*erf(%i/sqrt(p))*%e^-(1/p)/(sqrt(%pi)*p^(3/2)) $
 429
 430 /*
 431  * See formula (8), section 4.16 of Table of Integral Transforms:
 432  *
 433  * t^u*bessel_i(v,a*t) -> gamma(u+v+1)*s^(-u-1)*assoc_legendre_p(u,-v,p/s)
 434  *
 435  * where s = sqrt(p^2-a^2);
 436  */
 437 factor(ratsimp(specint(t^(1/2)*bessel_i(1,t)*%e^(-p*t),t)));
 438 3*sqrt(%pi)*assoc_legendre_p(1/2,-1,p/sqrt(p^2-1))/(4*(p^2-1)^(3/4))$
 439
 440 /*
 441  * hankel_1(2/3,sqrt(t)) = bessel_j(2/3,sqrt(t))+%i*bessel_y(2/3,sqrt(t))
 442  *
 443  * Formula (34) below gives:
 444  *
 445  * t^(u-1/2)*bessel_j(2*v,2*sqrt(a)*sqrt(t)) ->
 446  *   gamma(u+v+1/2)/sqrt(a)/gamma(2*v+1)*p^(-u)*exp(-a/p/2)*%m[u,v](a/p)
 447  *
 448  * Formula (50) gives
 449  *
 450  * t^(u-1/2)*bessel_y(2*v,2*sqrt(a)*sqrt(t)) ->
 451  *   1/sqrt(a)*p^(-u)*exp(-a/p/2)*
 452  *     (tan((u-v)*%pi)*gamma(u+v-1/2)/gamma(2*v+1) * %m[u,v](a/p)
 453  *       - sec((u-v)*%pi)*%w[u,v](a/p))
 454  *
 455  * But A&S 13.1.34 says
 456  *
 457  * %w[k,u](z) = gamma(-2*u)/gamma(1/2-u-k)*%m[k,u](z)
 458  *               + gamma(2*u)/gamma(1/2+u-k)*%m[k,-u](z)
 459  *
 460  */
 461
 462 ratsimp(specint(t*hankel_1(2/3,t^(1/2))*%e^(-p*t),t));
 463 /* Because of revision 1.110 of hyp.lisp Maxima knows in addition
 464  *    hgfred([7/3],[5/3],-1/(4*x)),
 465  * the result is in terms of the bessel_i function.
 466
 467 -4*%i*gamma(1/3)*%m[-3/2,1/3](-1/(4*p))*%e^-(1/(8*p))/(3*(-1)^(5/6)*sqrt(3)
 468  *gamma(2/3)*p^(3/2))+4*gamma(1/3)*%m[-3/2,1/3](-1/(4*p))*%e^-(1/(8*p))/(3*(
 469  -1)^(5/6)*gamma(2/3)*p^(3/2))-8*%i*gamma(2/3)*%m[-3/2,-1/3](-1/(4*p))*%e^-
 470  (1/(8*p))/(3*(-1)^(1/6)*sqrt(3)*gamma(1/3)*p^(3/2)) $ */
 471
 472 (((-1)^(1/6)*2^(2/3)*sqrt(3)*%i-3*(-1)^(1/6)*2^(2/3))
 473         *gamma(1/3)^2*gamma(5/6)*bessel_i(11/6,-1/(8*p))
 474         +10*(-1)^(5/6)*sqrt(3)*4^(2/3)*%i*gamma(1/6)*gamma(2/3)^2
 475            *bessel_i(7/6,-1/(8*p))
 476         +(45*(-1)^(1/6)*2^(2/3)-5*(-1)^(1/6)*2^(2/3)*3^(3/2)*%i)
 477          *gamma(1/3)^2*gamma(5/6)*bessel_i(5/6,-1/(8*p))
 478         +((9*(-1)^(1/6)*2^(5/3)-(-1)^(1/6)*2^(5/3)*3^(3/2)*%i)
 479          *bessel_i(-1/6,-1/(8*p))
 480          +(5*(-1)^(1/6)*2^(5/3)*sqrt(3)*%i-15*(-1)^(1/6)*2^(5/3))
 481           *bessel_i(-7/6,-1/(8*p)))
 482          *gamma(1/3)^2*gamma(5/6)
 483         +(((-1)^(5/6)*sqrt(3)*4^(2/3)*%i*bessel_i(-11/6,-1/(8*p))
 484          -5*(-1)^(5/6)*3^(3/2)*4^(2/3)*%i*bessel_i(-5/6,-1/(8*p)))
 485          *gamma(1/6)
 486          -2*(-1)^(5/6)*3^(3/2)*4^(2/3)*%i*gamma(1/6)*bessel_i(1/6,-1/(8*p)))
 487          *gamma(2/3)^2)
 488         *%e^-(1/(8*p))
 489         /(15*2^(13/3)*4^(1/3)*gamma(1/3)*gamma(2/3)*p^(7/2))/-1;
 490
 491 /*
 492  * hankel_2(3/4,t) = bessel_j(3/4,t)-%i*bessel_y(3/4,t)
 493  *
 494  * Sec 4.14, formula (9):
 495  *
 496  * t^u*bessel_j(v,a*t) -> gamma(u+v+1)*r^(-u-1)*assoc_legendre_p(u,-v,p/r)
 497  *
 498  * where r = sqrt(p^2+a^2)
 499  *
 500  * Sec 4.14, formula (48)
 501  *
 502  * t^u*bessel_y(v,a*t)
 503  *    -> r^(-u-1)*(gamma(u+v+1)*cot(v*%pi)*assoc_legendre_p(u,-v,p/r)
 504  *                  -gamma(u-v+1)*csc(v*%pi)*assoc_legendre_p(u,v,p/r))
 505  *
 506  * So,
 507  *
 508  * t^(1/2)*bessej_j(3/4,t)
 509  *    -> gamma(9/4)*r^(-3/2)*assoc_legendre_p(1/2,-3/4,p/r)
 510  *    =  5*gamma(1/4)/16*r^(-3/2)*assoc_legendre_p(1/2,-3/4,p/r)
 511  *
 512  * t^(1/2)*bessel_y(3/4,t)
 513  *    -> r^(-3/2)*(gamma(9/4)*cot(3/4*%pi)*assoc_legendre_p(1/2,-3/4,p/r)
 514  *                  -gamma(3/4)*csc(3/4*%pi)*assoc_legendre_p(1/2,3/4,p/r))
 515  *    =  r^(-3/2)*(-gamma(9/4)*assoc_legendre_p(1/2,-3/4,p/r)
 516  *                  -gamma(3/4)*sqrt(2)*assoc_legendre_p(1/2,3/4,p/r))
 517  */
 518
 519 ratsimp(specint(t^(1/2)*hankel_2(3/4,t)*%e^(-p*t),t));
 520 ''(ratsimp(5*gamma(1/4)/16/(p^2+1)^(3/4)*assoc_legendre_p(1/2,-3/4,p/sqrt(p^2+1))
 521 +sqrt(2)*%i*assoc_legendre_p(1/2,3/4,p/sqrt(p^2+1))*gamma(3/4)
 522        /(p^2+1)^(3/4)
 523        +5*%i*gamma(1/4)*assoc_legendre_p(1/2,-3/4,p/sqrt(p^2+1))
 524    /(16*(p^2+1)^(3/4))))$
 525
 526 /*
 527  * hankel_1(1/2,t) = bessel_j(1/2,t)+%i*bessel_y(1/2,t)
 528  *
 529  * So,
 530  *
 531  * t^(3/2)*bessel_j(1/2,t)
 532  *     -> gamma(3/2+1/2+1)*r^(-5/2)*assoc_legendre_p(3/2,-1/2,p/r)
 533  *     =  2*r^(-5/2)*assoc_legendre_p(3/2,-1/2,p/r)
 534  * t^(3/2)*bessel_y(1/2,t)
 535  *     -> r^(-5/2)*(gamma(3/2+1/2+1)*cot(%pi/2)*assoc_legendre_p(3/2,-1/2,p/r)
 536  *                  -gamma(3/2-1/2+1)*csc(%pi/2)*assoc_legendre_p(3/2,1/2,p/r))
 537  *     =  -r^(-5/2)*assoc_legendre_p(3/2,1/2,p/r))
 538  *
 539  * assoc_legendre_p(3/2,+/-1/2,z) can be expressed in terms of
 540  * hypergeometric functions (A&S 8.1.2).
 541  *
 542  * assoc_legendre_p(3/2,-1/2,z)
 543  *   = 1/gamma(3/2)*((z+1)/(z-1))^(-1/4)*F(-3/2,5/2;3/2;(1-z)/2)
 544  *   = sqrt(2)*(z-1)^(1/4)*z*(z+1)^(1/4)/sqrt(%pi)
 545  *
 546  * assoc_legendre_p(3/2,1/2,z)
 547  *   = 1/gamma(-1/2)*((z+1)/(z-1))^(1/4)*F(-3/2,5/2;-1/2;(1-z)/2)
 548  *   = sqrt(2)*z*(2*z^2-3)/(sqrt(%pi)*(z-1)^(1/4)*(z+1)^(5/4))
 549  *
 550  *
 551  * So the result should be
 552  *
 553  * t^(3/2)*bessel_j(1/2,t)
 554  *    -> 4*p/(sqrt(2)*sqrt(%pi)*(p^2+1)^2)
 555  *
 556  * t^(3/2)*bessel_y(1/2,t)
 557  *    -> -sqrt(2)*(p-1)*(p+1)/(sqrt(%pi)*(p^2+1)^2)
 558  */
 559
 560 ratsimp(specint(t^(3/2)*hankel_1(1/2,t)*%e^(-p*t),t));
 561 -((sqrt(%pi)*(sqrt(2)*%i*p^2-2*sqrt(2)*p-sqrt(2)*%i))/(%pi*p^4+2*%pi*p^2+%pi))$
 562
 563 /*
 564  * Formula 2, p 105:
 565  *
 566  * t^(u-1)*bessel_y(v,a*t)
 567  *    -> -2/%pi*gamma(u+v)*(a^2+p^2)^(-u/2)
 568  *         *assoc_legendre_q(u-1,-v,p/sqrt(a^2+p^2))
 569  *
 570  * for a > 0, Re u > |Re v|
 571  *
 572  * We have u = 5/2, v = 1, so the result is
 573  *
 574  *    -4/%pi/(p^2+a^2)*assoc_legendre_q(3/2,-1,p/sqrt(p^2+a^2))
 575  *
 576  * The expected result is not correct.
 577  * With gamma(1+5/2) = 15*sqrt(%pi)/8 we get:
 578  *
 579  * 15/(4*sqrt(%pi))*(p^2+a^2)^(-5/4)*assoc_legendre_q(3/3,1,p/sqrt(p^2+a^2))
 580  *
 581  * That is the result of Maxima too. The example is correct.
 582  */
 583 factor(specint(t^(5/2-1)*bessel_y(1,a*t)*%e^(-p*t),t));
 584 -15*assoc_legendre_q(3/2,-1,p/sqrt(p^2+a^2))/(4*sqrt(%pi)*(p^2+a^2)^(5/4));
 585
 586 /*
 587  * A&S 13.1.32:
 588  *
 589  *   %m[k,u](t) = exp(-t/2)*t^(u+1/2)*M(1/2+u-k,1+2*u,t)
 590  *
 591  * So
 592  *
 593  *   %m[1/2,1](t) = exp(-t/2)*t^(3/2)*M(1,3,t)
 594  *
 595  * and the integrand is:
 596  *
 597  *   t^(3/2)*%m[1/2,1](t)*exp(-p*t)
 598  *      = t^3*M(1,3,t)*exp(-(p+1/2)*t)
 599  *      = t^3*M(1,3,t)*exp(-p'*t)
 600  *
 601  * f19p220 will give us the Laplace transform of t^3*M(1,3,t):
 602  *
 603  *   gamma(4)/p'^4*F(1,4;3;1/p')
 604  *
 605  * But p' = p+1/2, so the final result is
 606  *
 607  *   32*(6*p-1)/(2*p-1)^2/(2*p+1)^3
 608  *
 609  */
 610 ratsimp(specint( t^(3/2)*%m[1/2,1](t)*%e^(-p*t),t));
 611 ''(ratsimp(32*(6*p-1)/(2*p-1)^2/(2*p+1)^3));
 612
 613 (assume(p>a),true);
 614 true;
 615 /*
 616  * exp(a*t)*t^2*erf(sqrt(t))*exp(-p*t)
 617  *
 618  * A&S 7.1.21 gives erf(z) = 2/sqrt(%pi)*z*M(1/2,3/2,-z^2) so
 619  * erf(sqrt(t)) = 2/sqrt(%pi)*sqrt(t)*M(1/2,3/2,-t)
 620  *
 621  * Therefore, the integrand, with p' = p-a, is
 622  *
 623  *   2/sqrt(%pi)*t^(5/2)*M(1/2,3/2,-t)*exp(-p'*t)
 624  *
 625  * Applying f19p220, the Laplace transform is
 626  *
 627  *   2/sqrt(%pi)*gamma(7/2)/p'^(7/2)*F(1/2,7/2;3/2;-1/p')
 628  *
 629  * Maxima can compute F(1/2,7/2;3/2;-1/p') and substituting p'=p-a
 630  * gives us the desired answer.
 631  *
 632  * 15*(1/sqrt(1/(p-a)+1)-2/(3*(p-a)*(1/(p-a)+1)^(3/2))
 633  *                      +1/(5*(p-a)^2*(1/(p-a)+1)^(5/2)))
 634  *  /(4*(p-a)^(7/2))
 635  */
 636 specint(%e^(a*t)*t^2*erf(t^(1/2))*%e^(-p*t),t);
 637 15*(1/sqrt(1/(p-a)+1)-2/(3*(p-a)*(1/(p-a)+1)^(3/2))
 638                      +1/(5*(p-a)^2*(1/(p-a)+1)^(5/2)))
 639  /(4*(p-a)^(7/2)) $
 640
 641 /*
 642  * Laplace transforms from Tables of Integral Transforms
 643  */
 644
 645 /*
 646  * p 182, (1)
 647  *
 648  * bessel_j(v,a*t) -> r^(-1)*(a/R)^v
 649  *
 650  * where r = sqrt(p^2+a^2) and R = p + r
 651  */
 652 (assume(v>0),true);
 653 true$
 654
 655 radcan(specint(bessel_j(v,a*t)*exp(-p*t),t));
 656 a^v/(sqrt(p^2+a^2)*(sqrt(p^2+a^2)+p)^v)$
 657
 658 /*
 659  * (5)
 660  * bessel_j(v,a*t)/t -> v^(-1)*(a/R)^v
 661  *
 662  * (Maxima doesn't recognize that gamma(v)/gamma(v+1) is 1/v.)
 663  */
 664 radcan(specint(bessel_j(v,a*t)/t*exp(-p*t),t));
 665 a^v*gamma(v)/((sqrt(p^2+a^2)+p)^v*gamma(v+1))$
 666
 667 /*
 668  * (7)
 669  * t^v*bessel_j(v,a*t) -> 2^v/sqrt(%pi)*gamma(v+1/2)*a^v*r^(-2*v-1)
 670  *
 671  * Maxima doesn't recognize the relationship between gamma(2*v+1) and
 672  * gamma(v+1).
 673  */
 674 radcan(specint(t^v*bessel_j(v,a*t)*exp(-p*t),t));
 675 a^v*gamma(2*v+1)/((p^2+a^2)^((2*v+1)/2)*2^v*gamma(v+1))$
 676
 677 /*
 678  * (9)
 679  * t^u*bessel_j(v,a*t) -> gamma(u+v+1)*r^(-u-1)*P(u,-v,p/r)
 680  */
 681 (assume(v+u+1>0),true);
 682 true$
 683 (assume(a>0),true);
 684 true$
 685
 686 radcan(specint(t^u*bessel_j(v,a*t)*exp(-p*t),t));
 687 /* This is not the correct answer: see the formula above which is correct.
 688    We had a bug in the routine legf24 in hyp.lisp.
 689 assoc_legendre_p(-u-1,-v,p/sqrt(p^2+a^2))*gamma(v+u+1)
 690         /((p^2+a^2)^((u+1)/2))$
 691  */
 692
 693 assoc_legendre_p(u,-v,p/sqrt(p^2+a^2))*gamma(v+u+1)
 694    /((p^2+a^2)^((u+1)/2))$
 695
 696 /*
 697  * (25)
 698  * bessel_j(0,2*sqrt(a)*sqrt(t)) -> exp(-a/p)/p
 699  */
 700 specint(bessel_j(0,2*sqrt(a)*sqrt(t))*exp(-p*t),t);
 701 %e^-(a/p)/p$
 702
 703 /*
 704  * (27)
 705  * sqrt(t)*bessel_j(1,2*sqrt(a)*sqrt(t)) -> sqrt(a)*p^(-2)*exp(-a/p)
 706  */
 707 specint(sqrt(t)*bessel_j(1,2*sqrt(a)*sqrt(t))*exp(-p*t),t);
 708 sqrt(a)*%e^-(a/p)/p^2$
 709
 710 /*
 711  * (29)
 712  * t^(-1/2)*bessel_j(1,2*sqrt(a)*sqrt(t)) ->
 713  *    sqrt(%pi)/sqrt(p)*exp(-a/2/p)*bessel_i(v/2,a/2/p)
 714  */
 715 specint(t^(-1/2)*bessel_j(1,2*sqrt(a)*sqrt(t))*exp(-p*t),t);
 716 sqrt(%pi)*bessel_i(1/2,a/(2*p))*%e^-(a/(2*p))/sqrt(p)$
 717
 718 /*
 719  * (30)
 720  * t^(v/2)*bessel_j(v,2*sqrt(a)*sqrt(t)) ->
 721  *    a^(v/2)/p^(v+1)*exp(-a/p)
 722  */
 723 specint(t^(v/2)*bessel_j(v,2*sqrt(a)*sqrt(t))*exp(-p*t),t);
 724 a^(v/2)*p^(-v-1)*%e^-(a/p)$
 725
 726 /*
 727  * (31)
 728  * t^(-v/2)*bessel_j(v,2*sqrt(a)*sqrt(t)) ->
 729  *    exp(%i*v*%pi)*p^(v-1)/a^(v/2)/gamma(v)*exp(-a/p)*
 730  *     gamma_greek(v,a/p*exp(-%i*%pi)
 731  *
 732  * gamma_greek is the incomplete gamma function.
 733  */
 734 specint(t^(-v/2)*bessel_j(v,2*sqrt(a)*sqrt(t))*exp(-p*t),t);
 735 p^(v-1)*%e^-(a/p)*v*gamma_greek(v,-a/p)/(a^(v/2)*(-1)^v*gamma(v+1))$
 736
 737 /*
 738  * (32)
 739  * t^(v/2-1)*bessel_j(v,2*sqrt(a)*sqrt(t)) ->
 740  *    a^(-v/2)*gamma_greek(v,a/p)
 741  */
 742 specint(t^(v/2-1)*bessel_j(v,2*sqrt(a)*sqrt(t))*exp(-p*t),t);
 743 v*gamma(v)*gamma_greek(v,a/p)/(a^(v/2)*gamma(v+1))$
 744
 745 /*
 746  * (34)
 747  * t^(u-1/2)*bessel_j(2*v,2*sqrt(a)*sqrt(t)) ->
 748  *   gamma(u+v+1/2)/sqrt(a)/gamma(2*v+1)*p^(-u)*exp(-a/p/2)*%m[u,v](a/p)
 749  *
 750  * A&S 13.1.32 gives
 751  *
 752  *   %m[k,u](z) = exp(-z/2)*z^(u+1/2)*M(1/2+u-k,1+2*u,z)
 753  *
 754  * A&S 13.1.27 (Kummer Transformation):
 755  *
 756  *   M(a,b,z) = exp(z)*M(b-a,b,-z)
 757  *
 758  * So
 759  *
 760  *   %m[k,u](z) = exp(z/2)*z^(u+1/2)*M(1/2+u+k,1+2*u,-z)
 761  *
 762  * But %m[-k,u](-z) = exp(z/2)*(-z)^(u+1/2)*M(1/2+u+k,1+2*u,-z)
 763  *
 764  * Therefore
 765  *
 766  *   %m[k,u](z) = (-1)^(u+1/2)*%m[-k,u](-z)
 767  *
 768  * So the Laplace transform can also be written as
 769  *
 770  *   gamma(u+v+1/2)/sqrt(a)/gamma(2*v+1)*p^(-u)*exp(-a/p/2)
 771  *     *%m[-u,v](-a/p)*(-1)^(v+1/2)
 772  *
 773  * Which is the answer we produce.
 774  */
 775 prefer_whittaker:true;
 776 true$
 777 (assume(2*v+2*u+1>0),true);
 778 true$
 779
 780 specint(t^(u-1/2)*bessel_j(2*v,2*sqrt(a)*sqrt(t))*exp(-p*t),t);
 781 %m[-u,v](-a/p)*%e^-(a/(2*p))*(-1)^(-v-1/2)*gamma(v+u+1/2)
 782          /(sqrt(a)*p^u*gamma(2*v+1))$
 783
 784 /*
 785  * (35)
 786  * t^(u-1)*bessel_j(2*v,2*sqrt(a)*sqrt(t)) ->
 787  *    gamma(u+v)*a^v/gamma(2*v+1)/p^(u+v)*%f[1,1](u+v,2*v+1,-a/p)
 788  */
 789 prefer_whittaker:false;
 790 false$
 791 (assume(v+u>0),true);
 792 true$
 793
 794 specint(t^(u-1)*bessel_j(2*v,2*sqrt(a)*sqrt(t))*exp(-p*t),t);
 795 a^v*p^(-v-u)*gamma(v+u)*%f[1,1]([v+u],[2*v+1],-a/p)/gamma(2*v+1)$
 796
 797 /*
 798  * (45)
 799  * bessel_y(v,a*t) ->
 800  *    a^v*cot(v*%pi)/r*R^(-v)-a^(-v)*csc(v*%pi)/r*R^v
 801  * For |Re v| < 1.
 802  *
 803  */
 804 expand(factor(radcan(specint(exp(-p*t)*bessel_y(1/6,a*t),t))));
 805 sqrt(3)*a^(1/6)/(sqrt(p^2+a^2)*(sqrt(p^2+a^2)+p)^(1/6))
 806         -2*(sqrt(p^2+a^2)+p)^(1/6)/(a^(1/6)*sqrt(p^2+a^2))$
 807
 808 (assume(v1 > 0, v1 < 1), true);
 809 true$
 810 expand(factor(radcan(specint(exp(-p*t)*bessel_y(v1,a*t),t))));
 811 a^v1*cot(%pi*v1)/(sqrt(p^2+a^2)*(sqrt(p^2+a^2)+p)^v1)
 812        -(sqrt(p^2+a^2)+p)^v1/(a^v1*sqrt(p^2+a^2)*sin(%pi*v1)) $
 813
 814
 815 /*
 816  * (42)
 817  *
 818  * t^(lam-1)*bessel_j(2*u,2*sqrt(a)*sqrt(t))*bessel_j(2*v,2*sqrt(a)*sqrt(t)) ->
 819  *    gamma(lam+u+v)/gamma(2*u+1)/gamma(2*v+1)*a^(u+v)/p^(lam+u+v)
 820  *      *%f[3,3]([u+v+1/2,u+v+1,lam+u+v],[2*u+1,2*v+1,2*u+2*v+1],-4*a/p)
 821  *
 822  */
 823 (assume(u>0,v>0,lam>0),true);
 824 true$
 825 specint(t^(lam-1)*bessel_j(2*u,2*sqrt(a)*sqrt(t))*bessel_j(2*v,2*sqrt(a)*sqrt(t))*exp(-p*t),t);
 826 a^(v+u)*p^(-v-u-lam)*gamma(v+u+lam)
 827               *%f[3,3]([v+u+1/2,v+u+1,v+u+lam],[2*u+1,2*v+1,2*v+2*u+1],-4*a/p)
 828         /(gamma(2*u+1)*gamma(2*v+1))$
 829
 830
 831 /*
 832  * (44)
 833  *
 834  * bessel_y(0,a*t) -> -2/%pi/sqrt(p^2+a^2)*asinh(p/a)
 835  *
 836  * Maxima returns
 837  *
 838  * -2/%pi/sqrt(p^2+a^2)*legendre_q(0,p/sqrt(p^2+a^2))
 839  *
 840  * But legendre_q(0,p/r) = log((1+p/r)/(1-p/r))/2, where r = sqrt(p^2+a^2).
 841  * This simplifies to log((1+p/r)/a) = log(p/a+sqrt(1+(p/a)^2)) = asinh(p/a).
 842  *
 843  * So we have -2/%pi/sqrt(p^2+a^2)*asinh(p/a).
 844  *
 845  * With revision 1.64 of hypgeo.lisp we simplify the Legendre Q function.
 846  * The result is equivalent to the above formula.
 847  */
 848 specint(bessel_y(0,a*t)*exp(-p*t),t);
 849
 850 /*-2/%pi/sqrt(p^2+a^2)*legendre_q(0,p/sqrt(p^2+a^2)) $*/
 851
 852  -log((sqrt(p^2+a^2)+p)/(sqrt(p^2+a^2)-p))/(%pi*sqrt(p^2+a^2));
 853
 854 /*
 855  * (46)
 856  *
 857  * t*bessel_y(0,a*t)
 858  *     -> 2/%pi/(p^2+a^2)*(1-p/sqrt(p^2+a^2)*log((p+sqrt(p^2+a^2))/a))
 859  *
 860  * Maxima returns
 861  *
 862  *    -2/%pi/(p^2+a^2)*legendre_q(1,p/sqrt(p^2+a^2))
 863  *
 864  * But
 865  *     legendre_q(1,p/r) = p/r/2*log((r+p)/(r-p)) - 1
 866  *                       = p/r*log((p+r)/a) - 1
 867  *
 868  * So, the transform is
 869  *
 870  *     -2/%pi/(p^2+a^2)*(p/r*log((p+r)/a) - 1)
 871  *
 872  *       = 2/%pi/(p^2+a^2)*(1-p/sqrt(p^2+a^2)*log((p+sqrt(p^2+a^2))/a))
 873  *
 874  * The Legendre Q function simplifes accordingly.
 875  */
 876 factor(specint(t*bessel_y(0,a*t)*exp(-p*t),t));
 877 /*-2/%pi/(p^2+a^2)*legendre_q(1,p/sqrt(p^2+a^2)) $*/
 878
 879 (p*log((sqrt(p^2+a^2)+p)/(sqrt(p^2+a^2)-p))-2*sqrt(p^2+a^2))
 880     /(-%pi*(p^2+a^2)^(3/2));
 881
 882 /*
 883  * (47)
 884  *
 885  * t*bessel_y(1,a*t)
 886  *     -> -2/%pi/(p^2+a^2)*(p/a+a/sqrt(p^2+a^2)*log((p+sqrt(p^2+a^2))/a)
 887  *
 888  * Maxima returns
 889  *   -4/%pi/(p^2+a^2)*assoc_legendre_q(1,-1,p/r)
 890  *
 891  * But
 892  *
 893  *   assoc_legendre_q(1,-1,z)
 894  *      = sqrt(1-z^2)/2/(z^2-1)*((z^2-1)*log((1+z)/(1-z)) - 2*z)
 895  *
 896  * So
 897  *
 898  *   assoc_legendre_q(1,-1,p/r)
 899  *      = (a/r)/2*(-(r/a)^2)*(-(a/r)^2*log((R/a)^2)-2*p/r)
 900  *      = 1/2*(p/a+a/r*log(R/a))
 901  *
 902  * where R = p + r
 903  *
 904  * Finally, the transform is
 905  *
 906  * -2/%pi/(p^2+a^2)*(p/a+a/r*log(R/a))
 907  *
 908  * as expected.
 909  *
 910  */
 911 factor(specint(t*bessel_y(1,a*t)*exp(-p*t),t));
 912 /*-4/%pi/(p^2+a^2)*assoc_legendre_q(1,-1,p/sqrt(p^2+a^2))$*/
 913
 914 (a^2*log((sqrt(p^2+a^2)+p)/(sqrt(p^2+a^2)-p))+2*p*sqrt(p^2+a^2))
 915         /(%pi*a*(p^2+a^2)^(3/2));
 916
 917
 918 /*
 919  * Some tests for step7
 920  */
 921
 922 /*
 923  * F(s,s+1/2;2*s+1;z) can be transformed to F(s,s+1/2;2*s+2;z) via
 924  * A&S 15.2.6.  And we know that F(s,s+1/2;2*s+1;z) =
 925  * 2^(2*s)/(1+sqrt(1-z))^(2*s).
 926  *
 927  * A&S 15.2.6 says
 928  * F(a,b;c+n;z) = poch(c,n)/poch(c-a,n)/poch(c-b,n)*(1-z)^(c+n-a-b)
 929  *                 *diff((1-z)^(a+b-c)*F(a,b;c;z),z,n)
 930  *
 931  * F(s,s+1/2;2*s+2;z)
 932  *    = poch(2*s+1,1)/poch(s+1,1)/poch(s+1/2,1)*(1-z)^(3/2)
 933  *       *diff((1-z)^(-1/2)*F(s,s+1/2;2*s+1;z),z,1)
 934  *
 935  *
 936  */
 937
 938 hgfred([s,s+1/2],[2*s+2],z)
 939   -(2*s+1)/(s+1)/(s+1/2)*(1-z)^(3/2)*diff((1-z)^(-1/2)
 940     *hgfred([s,s+1/2],[2*s+1],z),z);
 941 0$
 942
 943 /*
 944  * F(s,s+1/2;2*s+1;z) can be transformed to F(s+2,s+1/2;2*s+1;z) via
 945  * A&S 15.2.3:
 946  *    F(a+n,b;c;z) = z^(1-a)/poch(a,n)*diff(z^(a+n-1)*F(a,b;c;z),z,n)
 947  *
 948  * F(s+2,s+1/2;2*s+1,z)
 949  *    = z^(1-s)/s/(s+1)*diff(z^(s+1)*F(s,s+1/2;2*s+1;z),z,2)
 950  */
 951
 952 hgfred([s+2,s+1/2],[2*s+1],z)
 953   - z^(1-s)/s/(s+1)*diff(z^(s+1)*hgfred([s,s+1/2],[2*s+1],z),z,2);
 954 0$
 955
 956 /* Tests for Airy functions */
 957 closeto(e,tol):=block([numer:true,abse],abse:abs(e),if(abse<tol) then true else abse);
 958 closeto(e,tol):=block([numer:true,abse],abse:abs(e),if(abse<tol) then true else abse);
 959
 960 /* Derivatives of Airy functions */
 961 diff(airy_ai(x),x);
 962 airy_dai(x);
 963 diff(airy_dai(x),x);
 964 x*airy_ai(x);
 965 diff(airy_bi(x),x);
 966 airy_dbi(x);
 967 diff(airy_dbi(x),x);
 968 x*airy_bi(x);
 969
 970 /* Integrals of Airy functions */
 971 integrate(airy_ai(z),z);
 972 hypergeometric([1/3],[2/3,4/3],z^3/9)*z/(3^(2/3)*gamma(2/3))
 973    -3^(1/6)*gamma(2/3)*hypergeometric([2/3],[4/3,5/3],z^3/9)*z^2/(4*%pi);
 974 integrate(airy_dai(z),z);
 975 airy_ai(z);
 976 integrate(airy_bi(z),z);
 977 3^(2/3)*gamma(2/3)*hypergeometric([2/3],[4/3,5/3],z^3/9)*z^2/(4*%pi)
 978    +hypergeometric([1/3],[2/3,4/3],z^3/9)*z/(3^(1/6)*gamma(2/3));
 979 integrate(airy_dbi(z),z);
 980 airy_bi(z);
 981
 982 /* A&S 10.4.1  Airy functions satisfy the Airy differential equation */
 983 diff(airy_ai(x),x,2)-x*airy_ai(x);
 984 0;
 985 diff(airy_bi(x),x,2)-x*airy_bi(x);
 986 0;
 987
 988 /* A&S 10.4.4  Normalization of Airy Ai function */
 989 (c1:3^(-2/3)/gamma(2/3), closeto(airy_ai(0)-c1,1.0e-15));
 990 true;
 991 closeto(airy_bi(0)/sqrt(3)-c1,1.0e-15);
 992 true;
 993
 994 /* A&S 10.4.5  Normalization of Airy Bi function */
 995 (c2:3^(-1/3)/gamma(1/3),closeto(-airy_dai(0)-c2,1.0e-15));
 996 true;
 997 closeto(airy_dbi(0)/sqrt(3)-c2,1.0e-14);
 998 true;
 999
1000 /* Exact values A&S 10.4.4 and 10.4.5 */
1001 airy_ai(0);
1002 1/(3^(2/3)*gamma(2/3));
1003 airy_dai(0);
1004 -(1/(3^(1/3)*gamma(1/3)));
1005 airy_bi(0);
1006 1/(3^(1/6)*gamma(2/3));
1007 airy_dbi(0);
1008 3^(1/6)/gamma(1/3);
1009
1010 /* A&S 10.4.10 - Wronskian */
1011 AS_10_4_10(z):=expand(airy_ai(z)*airy_dbi(z)-airy_bi(z)*airy_dai(z)-1/%pi);
1012 AS_10_4_10(z):=expand(airy_ai(z)*airy_dbi(z)-airy_bi(z)*airy_dai(z)-1/%pi);
1013 closeto(AS_10_4_10(1),1.0e-15);
1014 true;
1015 closeto(AS_10_4_10(1+%i),1.0e-15);
1016 true;
1017 closeto(AS_10_4_10(%i),1.0e-15);
1018 true;
1019 closeto(AS_10_4_10(-1+%i),2.0e-15);
1020 true;
1021 closeto(AS_10_4_10(-1),1.0e-15);
1022 true;
1023 closeto(AS_10_4_10(-1-%i),2.0e-15);
1024 true;
1025 closeto(AS_10_4_10(-%i),1.0e-15);
1026 true;
1027 closeto(AS_10_4_10(1-%i),1.0e-15);
1028 true;
1029
1030 /* A&S 10.4.14 - only for z>0 ? */
1031 AS_10_4_14(z):=block([y:(2/3)*(z)^(3/2)],airy_ai(z)-(sqrt(z/3)*bessel_k(1/3,y)/%pi));
1032 AS_10_4_14(z):=block([y:(2/3)*(z)^(3/2)],airy_ai(z)-(sqrt(z/3)*bessel_k(1/3,y)/%pi));
1033 closeto(AS_10_4_14(1),1.0e-15);
1034 true;
1035 closeto(AS_10_4_14(2),1.0e-15);
1036 true;
1037 closeto(AS_10_4_14(5),1.0e-15);
1038 true;
1039 closeto(AS_10_4_14(10),1.0e-15);
1040 true;
1041
1042 /* A&S 10.4.15 - only for z<0 ? */
1043 AS_10_4_15(z):=block([y:(2/3)*(-z)^(3/2)],airy_ai(z)-(1/3)*sqrt(-z)*(bessel_j(1/3,y)+bessel_j(-1/3,y)));
1044 AS_10_4_15(z):=block([y:(2/3)*(-z)^(3/2)],airy_ai(z)-(1/3)*sqrt(-z)*(bessel_j(1/3,y)+bessel_j(-1/3,y)));
1045 closeto(AS_10_4_15(-1),1.0e-15);
1046 true;
1047 closeto(AS_10_4_15(-2),1.0e-15);
1048 true;
1049 closeto(AS_10_4_15(-5),1.0e-14);
1050 true;
1051 closeto(AS_10_4_15(-10),1.0e-15);
1052 true;
1053
1054 /* A&S 10.4.16 - only for z>0 ? */
1055 AS_10_4_16(z):=block([y:(2/3)*(z)^(3/2)],airy_dai(z)+(z/sqrt(3))*bessel_k(2/3,y)/%pi);
1056 AS_10_4_16(z):=block([y:(2/3)*(z)^(3/2)],airy_dai(z)+(z/sqrt(3))*bessel_k(2/3,y)/%pi);
1057 closeto(AS_10_4_16(1),1.0e-15);
1058 true;
1059 closeto(AS_10_4_16(2),1.0e-15);
1060 true;
1061 closeto(AS_10_4_16(5),1.0e-15);
1062 true;
1063 closeto(AS_10_4_16(10),1.0e-15);
1064 true;
1065
1066 /* A&S 10.4.17 - only for z<0 ?, Appears to be a sign error in A&S */
1067 AS_10_4_17(z):=block([y:(2/3)*(-z)^(3/2)],airy_dai(z)-(z/3)*(bessel_j(-2/3,y)-bessel_j(2/3,y)));
1068 AS_10_4_17(z):=block([y:(2/3)*(-z)^(3/2)],airy_dai(z)-(z/3)*(bessel_j(-2/3,y)-bessel_j(2/3,y)));
1069 closeto(AS_10_4_17(-1),1.0e-15);
1070 true;
1071 closeto(AS_10_4_17(-2),1.0e-15);
1072 true;
1073 closeto(AS_10_4_17(-5),1.0e-14);
1074 true;
1075 closeto(AS_10_4_17(-10),1.0e-15);
1076 true;
1077
1078 /* Test that complex float arguments are evaluated */
1079 airy_ai(%i);
1080 airy_ai(%i);
1081 floatnump(realpart(airy_ai(1.0*%i)));
1082 true;
1083
1084 kill(c1,c2,AS_10_4_10,AS_10_4_14,AS_10_4_15,AS_10_4_16,AS_10_4_17);
1085 done;
1086
1087 /* End of Airy function tests */
1088
1089 /* Numerical tests of gamma function. */
1090
1091 /* A&S Table 1.1, to 15 DP */
1092 closeto(gamma(1/2)-1.772453850905516,2e-15);
1093 true;
1094 closeto(gamma(1/3)-2.678938534707748,3.0e-15);
1095 true;
1096 closeto(gamma(7/4)-0.919062526848883,1e-15);
1097 true;
1098
1099 /* Complex values.  Checked against A&S Table 6.7 to 12 DP */
1100 closeto(gamma(1+%i)-(0.49801566811836-0.15494982830181*%i),1e-14);
1101 true;
1102 closeto(gamma(1+5*%i)-(-0.00169966449436-0.00135851941753*%i),1e-14);
1103 true;
1104 closeto(gamma(2+3*%i)-(-0.08239527266561+0.09177428743526*%i),1e-14);
1105 true;
1106
1107 /* Test numerical evaluation of Bessel functions
1108  * When the order is 0, and the arg is a float, we should produce a number.
1109  */
1110 closeto(bessel_j(0,1.0) - .7651976865579666, 1e-14);
1111 true;
1112
1113 closeto(bessel_y(0,1.0) - .08825696421567691, 1e-14);
1114 true;
1115
1116 closeto(bessel_i(0,1.0) - 1.266065877752009, 1e-14);
1117 true;
1118
1119 closeto(bessel_k(0,1.0) - .4210244382407085, 1e-14);
1120 true;
1121
1122 /*
1123  * Tests for failed cases to see if we're returning the noun forms
1124  */
1125 /* fail-on-f24p146test */
1126 specint(t^(-3/2)*exp(-t^2/8/a)*exp(-p*t),t);
1127 'specint(t^(-3/2)*exp(-t^2/8/a)*exp(-p*t),t);
1128
1129 /* fail-on-f35p147test
1130    Because we modified the construction of a noun form, we get a sligthly
1131    different noun form as result. DK. */
1132 specint((2*t)^(-3/2)*exp(-2*sqrt(a)*sqrt(t))*exp(-p*t),t);
1133 /*'specint(exp(-p*t -2*sqrt(a)*sqrt(t))/(8*t^(3/2)),t);*/
1134 'specint(%e^(-p*t-2*sqrt(a)*sqrt(t))/(2*sqrt(2)*t^(3/2)),t);
1135
1136 /* fail-on-f29p146test */
1137 specint(t^(v-1)*exp(1/8/t)*exp(-p*t),t);
1138 'specint(t^(v-1)*exp(1/8/t)*exp(-p*t),t);
1139
1140 /* fail-in-arbpow (aka f1p137test) */
1141 specint(t^(-1)*exp(-p*t),t);
1142 'specint(exp(-p*t)/t,t);
1143
1144 /* fail-in-f2p105vcond */
1145 (assume(p>3),0);
1146 0;
1147
1148 specint(8*t^1*exp(3*t)*bessel_y(3,t)*exp(-p*t),t);
1149 'specint(8*bessel_y(3,t)*t*%e^(3*t-p*t),t)$
1150
1151 /* fail-in-f50cond */
1152 specint(8*t^1*exp(3*t)*bessel_y(8,7*sqrt(t))*exp(-p*t),t);
1153 'specint(8*bessel_y(8,7*sqrt(t))*t*%e^(3*t-p*t),t);
1154
1155 /* fail-in-dionarghyp-y */
1156 specint(8*t^1*exp(3*t)*bessel_y(8,7*t^(3/2))*exp(-p*t),t);
1157 'specint(8*bessel_y(8,7*t^(3/2))*t*%e^(3*t-p*t),t);
1158
1159 /* The additionally phase factor in the calculation vanish, because of
1160    the modificaton of the transformation to Bessel J in the code. DK.
1161 */
1162 (assume(t>0,v2>1), radcan(specint(bessel_i(-v2,t)*exp(-p*t),t)));
1163 /*'specint((%e^((%i*%pi*v2-2*p*t)/2)*bessel_i(-v2,t))/(-1)^(v2/2),t);*/
1164 '(specint(bessel_i(-v2,t)*exp(-p*t),t));
1165
1166 /* Verify fix for [ 1877522 ] erf(1.0),nouns wrong; causes plot2d(erf) to fail
1167  */
1168 erf(1.0), nouns;
1169 0.8427007929497148;
1170
1171 erf(1.0), erf;
1172 0.8427007929497148;
1173
1174 erf(1.0);
1175 0.8427007929497148;
1176
1177 /* [ 789059 ] simpexpt problem: sign called on imag arg */
1178 (-(-1)^(1/6))^(1/2);
1179 sqrt(-(-1)^(1/6));
1180
1181 /* Further tests of bessel functions */
1182 /* The following numerical values are evaluated with the evaluation tool
1183    of the website www.functions.wolfram.com with a precision of 16 digits
1184    1998-2014 Wolfram Research, Inc. */
1185
1186 (test_bessel(actual, ref, digits) := closeto(realpart(actual)-realpart(ref), 10^(-digits)) and closeto(imagpart(actual)-imagpart(ref), 10^(-digits)), 0);
1187 0;
1188
1189 /* Numerical values for the bessel function J with negative order */
1190
1191 test_bessel(bessel_j(-1,-2.0),0.5767248077568734, 15);
1192 true;
1193
1194 test_bessel(bessel_j(-1,2.0), -0.5767248077568734, 15);
1195 true;
1196
1197 test_bessel(bessel_j(-1,-1.5), 0.5579365079100996, 15);
1198 true;
1199
1200 test_bessel(bessel_j(-1,1.5), -0.5579365079100996, 15);
1201 true;
1202
1203 test_bessel(bessel_j(-1.5, -2.0), -0.3956232813587035*%i, 15);
1204 true;
1205
1206 test_bessel(bessel_j(-1.5, 2.0), -0.3956232813587035, 15);
1207 true;
1208
1209 test_bessel
1210   (bessel_j(-1.8, -1.5),- 0.2033279902093184 - 0.1477264320209275 * %i,15);
1211 true;
1212
1213 test_bessel(bessel_j(-1.8, 1.5), -0.251327217627129314,15);
1214 true;
1215
1216 test_bessel(bessel_j(-2,-1.5), 0.2320876721442147, 15);
1217 true;
1218
1219 test_bessel(bessel_j(-2,1.5), 0.2320876721442147, 15);
1220 true;
1221
1222 test_bessel(bessel_j(-2.5,-1.5), -1.315037204805194 * %i,15);
1223 true;
1224
1225 test_bessel(bessel_j(-2.5,1.5), 1.315037204805194,15);
1226 true;
1227
1228 test_bessel
1229   (bessel_j(-2.3,-1.5), 0.5949438455752484 - 0.8188699527453657 * %i,14);
1230 true;
1231
1232 test_bessel(bessel_j(-2.3,1.5), 1.012178926325313, 14);
1233 true;
1234
1235 /* Numerical values for the bessel function J with positive order */
1236
1237 test_bessel(bessel_j(1.5,1.0), 0.2402978391234270, 15);
1238 true;
1239
1240 test_bessel(bessel_j(1.5,-1.0), -0.2402978391234270 * %i, 15);
1241 true;
1242
1243 test_bessel(bessel_j(1.8,1.0), 0.1564953153109239, 14);
1244 true;
1245
1246 test_bessel
1247   (bessel_j(1.8,-1.0), 0.1266073696266034 - 0.0919856383926216 * %i, 15);
1248 true;
1249
1250 test_bessel(bessel_j(2.0,1.0), 0.1149034849319005,15);
1251 true;
1252
1253 test_bessel(bessel_j(2.0,-1.0),0.1149034849319005,15);
1254 true;
1255
1256 test_bessel(bessel_j(2.5,1.0),0.04949681022847794,15);
1257 true;
1258
1259 test_bessel(bessel_j(2.5,-1.0),0.04949681022847794 * %i,15);
1260 true;
1261
1262 /* Numerical values for the bessel function J with complex arg
1263    and positive or negative order*/
1264
1265 test_bessel
1266   (bessel_j(0,1.0+%i),0.9376084768060293 - 0.4965299476091221 * %i,15);
1267 true;
1268
1269 test_bessel
1270   (bessel_j(1,1.0+%i),0.6141603349229036 + 0.3650280288270878 * %i,15);
1271 true;
1272
1273 test_bessel
1274   (bessel_j(-1,1.0+%i),-0.6141603349229036 - 0.3650280288270878 * %i,14);
1275 true;
1276
1277 test_bessel
1278   (bessel_j(2,1.0+%i),0.0415798869439621 + 0.2473976415133063 * %i,15);
1279 true;
1280
1281 test_bessel
1282   (bessel_j(-2,1.0+%i),0.0415798869439621 + 0.2473976415133063 * %i,15);
1283 true;
1284
1285 test_bessel
1286   (bessel_j(2.3,1.0+%i),-0.0141615213034667 + 0.1677798241687935 * %i,15);
1287 true;
1288
1289 test_bessel
1290   (bessel_j(-2.3,1.0+%i),0.1920598664138632 - 0.5158676904105332 * %i,14);
1291 true;
1292
1293 /* Numerical values for the bessel function J with complex order */
1294
1295 test_bessel
1296   (bessel_j(%i,1.0),1.641024179495082 - 0.437075010213683*%i,15);
1297 true;
1298
1299 test_bessel
1300   (bessel_j(%i,1.5),1.401883276281807 + 0.473362399311655*%i,15);
1301 true;
1302
1303 test_bessel
1304   (bessel_j(1.0+%i,-1.0),-0.01142279482478010 + 0.02390070064911069*%i,15);
1305 true;
1306
1307 test_bessel
1308   (bessel_j(1.5*%i,-2.0),0.01925195427338360 + 0.01442616961986814*%i,15);
1309 true;
1310
1311 /******************************************************************/
1312 /* Numerical values for the bessel function Y with negative order */
1313
1314 test_bessel
1315   (bessel_y(-1,-2.0),-0.1070324315409375 + 1.1534496155137468 * %i,14);
1316 true;
1317
1318 test_bessel(bessel_y(-1,2.0),0.1070324315409375,15);
1319 true;
1320
1321 test_bessel
1322   (bessel_y(-1,-1.5),-0.4123086269739113 + 1.1158730158201993 * %i,15);
1323 true;
1324
1325 test_bessel(bessel_y(-1,1.5),0.4123086269739113,15);
1326 true;
1327
1328 test_bessel(bessel_y(-1.5, -2.0),0.4912937786871623 * %i,15);
1329 true;
1330
1331 test_bessel(bessel_y(-1.5, 2.0),-0.4912937786871623,15);
1332 true;
1333
1334 test_bessel
1335   (bessel_y(-1.8, -1.5),-0.6777414340388011 + 0.0857519944018923 * %i,14);
1336 true;
1337
1338 test_bessel(bessel_y(-1.8, 1.5),-0.8377344836401481,14);
1339 true;
1340
1341 test_bessel
1342   (bessel_y(-2,-1.5),-0.9321937597629739 + 0.4641753442884295 * %i,15);
1343 true;
1344
1345 test_bessel(bessel_y(-2,1.5),-0.9321937597629739,14);
1346 true;
1347
1348 test_bessel(bessel_y(-2.5,-1.5),0.1244463597983876 * %i,14);
1349 true;
1350
1351 test_bessel(bessel_y(-2.5,1.5),0.1244463597983876,15);
1352 true;
1353
1354 test_bessel
1355   (bessel_y(-2.3,-1.5),-0.3148570865836879 + 0.7565240896444820 * %i,14);
1356 true;
1357
1358 test_bessel(bessel_y(-2.3,1.5),-0.5356668704355646,14);
1359 true;
1360
1361 /* Numerical values for the bessel function Y with positive order */
1362
1363 test_bessel(bessel_y(1.5,1.0),-1.102495575160179,14);
1364 true;
1365
1366 test_bessel(bessel_y(1.5,-1.0),-1.102495575160179 * %i,14);
1367 true;
1368
1369 test_bessel(bessel_y(1.8,1.0),-1.382351995367631,14);
1370 true;
1371
1372 test_bessel
1373   (bessel_y(1.8,-1.0),-1.1183462564605324 - 0.5593113771009602 * %i,14);
1374 true;
1375
1376 test_bessel(bessel_y(2.0,1.0),-1.650682606816254,14);
1377 true;
1378
1379 test_bessel
1380   (bessel_y(2.0,-1.0),-1.650682606816254 + 0.229806969863801 * %i,14);
1381 true;
1382
1383 test_bessel(bessel_y(2.5,1.0),-2.876387857462161,14);
1384 true;
1385
1386 test_bessel(bessel_y(2.5,-1.0),2.876387857462161 * %i,14);
1387 true;
1388
1389
1390 /* Numerical values for the bessel function Y with complex arg
1391    and positive or negative order*/
1392
1393 test_bessel
1394   (bessel_y(0,1.0+%i),0.4454744889360325 + 0.7101585820037345 * %i,15);
1395 true;
1396
1397 test_bessel
1398   (bessel_y(1,1.0+%i),-0.6576945355913452 + 0.6298010039928844 * %i,15);
1399 true;
1400
1401 test_bessel
1402 (bessel_y(-1,1.0+%i),0.6576945355913452 - 0.6298010039928844 * %i,15);
1403 true;
1404
1405 test_bessel
1406   (bessel_y(2,1.0+%i),-0.4733680205344934 + 0.5773369575804951 * %i,14);
1407 true;
1408
1409 test_bessel
1410   (bessel_y(-2,1.0+%i),-0.4733680205344934 + 0.5773369575804951 * %i,14);
1411 true;
1412
1413 test_bessel
1414   (bessel_y(2.3,1.0+%i),-0.2476879981252862 + 0.7595467103431256 * %i,15);
1415 true;
1416
1417 test_bessel
1418   (bessel_y(-2.3,1.0+%i),-0.1570442638685963 + 0.5821870838327466 * %i,14);
1419 true;
1420
1421 /* Numerical values for the bessel function Y with complex order */
1422
1423 test_bessel
1424   (bessel_y(%i,1.0),-0.476556612479964 - 1.505069159110387*%i,14);
1425 true;
1426
1427 test_bessel
1428   (bessel_y(%i,1.5),0.5161218926267688 - 1.2857405211747503*%i,14);
1429 true;
1430
1431 test_bessel
1432   (bessel_y(1.0+%i,-1.0),7.708594946281541 + 1.233384674244926*%i,14);
1433 true;
1434
1435 test_bessel
1436   (bessel_y(1.5*%i,-2.0),3.226466016458932 + 4.267260420563194*%i,13);
1437 true;
1438
1439 /******************************************************************/
1440 /* Numerical values for the bessel function I with negative order */
1441
1442 test_bessel(bessel_i(-1,-2.0),-1.590636854637329,15);
1443 true;
1444
1445 test_bessel(bessel_i(-1,2.0),1.590636854637329,15);
1446 true;
1447
1448 test_bessel(bessel_i(-1,-1.5),-0.9816664285779076,15);
1449 true;
1450
1451 test_bessel(bessel_i(-1,1.5),0.9816664285779076,15);
1452 true;
1453
1454 test_bessel(bessel_i(-1.5, -2.0),0.9849410530002364 * %i,14);
1455 true;
1456
1457 test_bessel(bessel_i(-1.5, 2.0),0.9849410530002364,14);
1458 true;
1459
1460 test_bessel
1461   (bessel_i(-1.8, -1.5),0.2026903980307014 + 0.1472631941876387 * %i,15);
1462 true;
1463
1464 test_bessel(bessel_i(-1.8, 1.5),0.2505391103524365,15);
1465 true;
1466
1467 test_bessel(bessel_i(-2,-1.5),0.3378346183356807,15);
1468 true;
1469
1470 test_bessel(bessel_i(-2,1.5),0.3378346183356807,15);
1471 true;
1472
1473 test_bessel(bessel_i(-2.5,-1.5),-0.8015666610717216*%i,14);
1474 true;
1475
1476 test_bessel(bessel_i(-2.5,1.5),0.8015666610717216,14);
1477 true;
1478
1479 test_bessel
1480   (bessel_i(-2.3,-1.5),0.3733945265830230 - 0.5139334755917659 * %i,15);
1481 true;
1482
1483 test_bessel(bessel_i(-2.3,1.5),0.6352567117441516,15);
1484 true;
1485
1486 /* Numerical values for the bessel function I with positive order */
1487
1488 test_bessel(bessel_i(1.5,1.0),0.2935253263474798,15);
1489 true;
1490
1491 test_bessel(bessel_i(1.5,-1.0),-0.2935253263474798*%i,13);
1492 true;
1493
1494 test_bessel(bessel_i(1.8,1.0),0.1871011888310777,15);
1495 true;
1496
1497 test_bessel
1498   (bessel_i(1.8,-1.0),0.1513680414320980 - 0.1099753194812967 * %i,14);
1499 true;
1500
1501 test_bessel(bessel_i(2.0,1.0),0.1357476697670383,15);
1502 true;
1503
1504 test_bessel(bessel_i(2.0,-1.0),0.1357476697670383,15);
1505 true;
1506
1507 test_bessel(bessel_i(2.5,1.0),0.05709890920304825,15);
1508 true;
1509
1510 test_bessel(bessel_i(2.5,-1.0),0.05709890920304825 * %i,15);
1511 true;
1512
1513
1514 /* Numerical values for the bessel function I with complex arg
1515    and positive or negative order*/
1516
1517 test_bessel
1518   (bessel_i(0,1.0+%i),0.9376084768060293 + 0.4965299476091221 * %i,15);
1519 true;
1520
1521 test_bessel
1522   (bessel_i(1,1.0+%i),0.3650280288270878 + 0.6141603349229036 * %i,15);
1523 true;
1524
1525 test_bessel
1526   (bessel_i(-1,1.0+%i),0.3650280288270878 + 0.6141603349229036 * %i,15);
1527 true;
1528
1529 test_bessel
1530   (bessel_i(2,1.0+%i),-0.0415798869439621 + 0.2473976415133063 * %i,15);
1531 true;
1532
1533 test_bessel
1534   (bessel_i(-2,1.0+%i),-0.0415798869439621 + 0.2473976415133063 * %i,15);
1535 true;
1536
1537 test_bessel
1538   (bessel_i(2.3,1.0+%i),-0.0635524383467825 + 0.1559221140952053 * %i,15);
1539 true;
1540
1541 test_bessel
1542   (bessel_i(-2.3,1.0+%i),-0.4053256245784623 - 0.3724481230406298 * %i,14);
1543 true;
1544
1545 /* Numerical values for the bessel function I with complex order */
1546
1547 test_bessel
1548   (bessel_i(%i,1.0),1.900799675819425 - 1.063960013554441*%i,14);
1549 true;
1550
1551 test_bessel
1552   (bessel_i(%i,1.5),2.495473638417463 - 0.601347501920535*%i,14);
1553 true;
1554
1555 test_bessel
1556   (bessel_i(1.0+%i,-1.0),-0.01096313515009349 + 0.03043920114776303*%i,15);
1557 true;
1558
1559 test_bessel
1560   (bessel_i(1.5*%i,-2.0),0.04238259669782487 - 0.01125055344512197*%i,15);
1561 true;
1562
1563 /******************************************************************/
1564 /* Numerical values for the bessel function K with negative order */
1565
1566 test_bessel
1567   (bessel_k(-1,-2.0),-0.139865881816522-4.997133057057809*%i,14);
1568 true;
1569
1570 test_bessel(bessel_k(-1,2.0),0.139865881816522,14);
1571 true;
1572
1573 test_bessel
1574   (bessel_k(-1,-1.5),-0.277387800456844 - 3.083996040296084*%i,14);
1575 true;
1576
1577 test_bessel(bessel_k(-1,1.5),0.2773878004568438,15);
1578 true;
1579
1580 test_bessel(bessel_k(-1.5, -2.0),-3.2741902342766302*%i,13);
1581 true;
1582
1583 test_bessel(bessel_k(-1.5, 2.0),0.1799066579520922,15);
1584 true;
1585
1586 test_bessel
1587   (bessel_k(-1.8, -1.5),0.3929372194683435 - 1.0725774293000646*%i,15);
1588 true;
1589
1590 test_bessel(bessel_k(-1.8, 1.5),0.4856971141526263,15);
1591 true;
1592
1593 test_bessel
1594   (bessel_k(-2,-1.5),0.5836559632566508 - 1.0613387550916862*%i,14);
1595 true;
1596
1597 test_bessel(bessel_k(-2,1.5),0.5836559632566508,15);
1598 true;
1599
1600 test_bessel
1601   (bessel_k(-2.3,-1.5),0.465598659425186 - 1.354876241730594*%i,14);
1602 true;
1603
1604 test_bessel(bessel_k(-2.3,1.5),0.7921237520153218,15);
1605 true;
1606
1607 /* Numerical values for the bessel function K with positive order */
1608
1609 test_bessel(bessel_k(1.5,1.0),0.9221370088957891,14);
1610 true;
1611
1612 /* kein Wert von function-site !? Ist das eine Nullstelle??? */
1613 test_bessel(bessel_k(1.5,-1.0),0,13);
1614 true;
1615
1616 test_bessel(bessel_k(1.8,1.0),1.275527037541854,15);
1617 true;
1618
1619 test_bessel
1620   (bessel_k(1.8,-1.0),1.0319230501560911 + 0.1619402612577788*%i,14);
1621 true;
1622
1623 test_bessel(bessel_k(2.0,1.0),1.624838898635177,14);
1624 true;
1625
1626 test_bessel(bessel_k(2.0,-1.0),1.624838898635177 - 0.426463882082061*%i,14);
1627 true;
1628
1629 test_bessel(bessel_k(2.5,1.0),3.227479531135262,14);
1630 true;
1631
1632 test_bessel(bessel_k(2.5,-1.0),-3.406861044815549*%i,14);
1633 true;
1634
1635 /* Numerical values for the bessel function K with complex arg
1636    and positive or negative order*/
1637
1638 test_bessel
1639   (bessel_k(0,1.0+%i),0.0801977269465178 - 0.3572774592853303*%i,15);
1640 true;
1641
1642 test_bessel
1643   (bessel_k(1,1.0+%i),0.0245683055237403 - 0.4597194738011894 * %i,15);
1644 true;
1645
1646 test_bessel
1647   (bessel_k(-1,1.0+%i),0.0245683055237403 - 0.4597194738011894 * %i,15);
1648 true;
1649
1650 test_bessel
1651   (bessel_k(2,1.0+%i),-0.3549534413309312 - 0.8415652386102600 * %i,13);
1652 true;
1653
1654 test_bessel
1655   (bessel_k(-2,1.0+%i),-0.3549534413309312 - 0.8415652386102600 * %i,13);
1656 true;
1657
1658 test_bessel
1659   (bessel_k(2.3,1.0+%i),-0.6635905911278042 - 1.0258894849569300 * %i,15);
1660 true;
1661
1662 test_bessel
1663   (bessel_k(-2.3,1.0+%i),-0.6635905911278042 - 1.0258894849569300 * %i,13);
1664 true;
1665
1666 /* Numerical values for the bessel function K with complex order */
1667
1668 test_bessel
1669   (bessel_k(%i,1.0),0.2894280370259921,14);
1670 true;
1671
1672 test_bessel
1673   (bessel_k(%i,1.5),0.1635839926633096,14);
1674 true;
1675
1676 test_bessel
1677   (bessel_k(1.0+%i,-1.0),-9.744252766030894 - 7.494570149760043*%i,14);
1678 true;
1679
1680 test_bessel
1681   (bessel_k(1.5*%i,-2.0),3.93512366711118 - 14.82183468316553*%i,13);
1682 true;
1683
1684 kill(closeto, test_bessel);
1685 done;
1686
1687 /* Numerical tests of the Bessel functions using the Wronskians
1688
1689    The Wronskians combines different types of Bessel functions and
1690    Bessel functions with negative and positve order.
1691    The results are very simple. Therefore the numerical calculation of the
1692    Wronskians is a good test for the different parts of the algorithmen.
1693
1694    Based on code by Dieter Kaiser.
1695 */
1696
1697 /* Test the Wronskian.  wf is the Wronskian function.  w_true is simplified Wronskian.
1698  * eps is the max absolute error allowed, and the Wronskian is tested
1699  * for values of the arg between -zlimit and zlimit.
1700  */
1701 (test_wronskian(wf, w_true, eps, zlimit) :=
1702 block([badpoints : [],
1703        ratprint : false,
1704        abserr : 0,
1705        maxerr : -1],
1706   for order:-1 thru 1 step 1/10 do
1707   (
1708     for z: -zlimit thru zlimit step 1 do
1709     (
1710       if notequal(z,0.0) then
1711       (
1712         result : float(rectform(wf(float(order),z))),
1713         answer : float(rectform(w_true(float(order),z))),
1714         abserr : abs(result - answer),
1715         maxerr : max(maxerr, abserr),
1716         if abserr > eps then
1717         (
1718           badpoints : cons([[order, z], result, answer, abserr], badpoints)
1719         )
1720       )
1721     )
1722   ),
1723   /*
1724    * For debugging, if there are any bad points, return the maximum error
1725    * found as the first element.
1726    */
1727   if badpoints # [] then
1728     cons(maxerr, badpoints)
1729   else
1730     badpoints
1731 ), 0);
1732 0;
1733
1734 /*******************************************************************************
1735
1736    Wronskian w_jj
1737
1738    A&S 9.1.15 : J[n+1](z)*J[-n](z)+J[n](z)*J[-(n+1)](z) = -2*sin(n*pi)/(pi*z)
1739
1740 *******************************************************************************/
1741
1742 w_jj(n,z) := bessel_j(n+1,z)*bessel_j(-n,z) + bessel_j(n,z)*bessel_j(-n-1,z);
1743 w_jj(n,z) := bessel_j(n+1,z) *bessel_j(-n,z) + bessel_j(n,z)*bessel_j(-n-1,z);
1744
1745 /* Calculation of w_jj for real argument */
1746
1747 test_wronskian('w_jj, lambda([n,z], -2.0*sin(n*%pi)/(z*%pi)), 1e-14, 10);
1748 [];
1749
1750
1751 /* Calculation of w_jj for complex argument */
1752
1753 test_wronskian(lambda([n,z], expand(w_jj(n,%i*z))),
1754                lambda([n,z],-2.0*sin(n*%pi)/(%i*z*%pi)),
1755                1e-8, 10);
1756 [];
1757
1758 /*******************************************************************************
1759
1760    Wronskian w_jy
1761
1762    A&S 9.1.16: J[n+1](z)*Y[n](z)-J[n](z)*Y[n+1,z] = 2/(pi*z)
1763
1764 *******************************************************************************/
1765
1766 w_jy(n,z) := bessel_j(n+1,z)*bessel_y(n,z) - bessel_j(n,z)*bessel_y(n+1,z);
1767 w_jy(n,z) := bessel_j(n+1,z)*bessel_y(n,z) - bessel_j(n,z)*bessel_y(n+1,z);
1768
1769 /* Calculation of w_yj for real argument */
1770
1771 test_wronskian(w_jy, lambda([n,z], 2.0/(z*%pi)), 1e-14, 10);
1772 [];
1773
1774
1775 /* Calculation of w_jy for complex argument */
1776
1777 test_wronskian(lambda([n,z], w_jy(n,z*%i)), lambda([n,z],2.0/(z*%i*%pi)), 1e-8, 10);
1778 [];
1779
1780 /*******************************************************************************
1781
1782    Wronskian w_ii
1783
1784    A&S 9.6.14: I[n](z)*I[-(n+1)](z)-I[n+1](z)*I[-n](z) = -2*sin(n*pi)/(pi*z)
1785
1786 *******************************************************************************/
1787
1788 w_ii(n,z) := bessel_i(n,z)*bessel_i(-n-1,z) - bessel_i(n+1,z)*bessel_i(-n,z);
1789 w_ii(n,z) := bessel_i(n,z)*bessel_i(-n-1,z) - bessel_i(n+1,z)*bessel_i(-n,z);
1790
1791 /* Calculation of w_ii for real argument */
1792
1793 test_wronskian(w_ii, lambda([n,z],-2.0*sin(n*%pi)/(z*%pi)), 1e-10, 5);
1794 [];
1795
1796 /* Calculation of w_ii for complex argument */
1797
1798 test_wronskian(lambda([n,z], w_ii(n,z*%i)),
1799                lambda([n,z], -2.0*sin(n*%pi)/(z*%i*%pi)),
1800                1e-10, 5);
1801 [];
1802
1803 /*******************************************************************************
1804
1805    Test Wronskian w_ik
1806
1807    A&S 9.6.15: I[n](z)*K[n+1](z)+I[n+1](z)*K[n](z) = 1/z
1808
1809 *******************************************************************************/
1810
1811 w_ik(n,z) := bessel_i(n,z)*bessel_k(n+1,z) + bessel_i(n+1,z)*bessel_k(n,z);
1812 w_ik(n,z) := bessel_i(n,z)*bessel_k(n+1,z) + bessel_i(n+1,z)*bessel_k(n,z);
1813
1814 /* Calculation of w_ik for real argument */
1815
1816 test_wronskian(w_ik, lambda([n,z], 1/z), 1e-10, 5);
1817 [];
1818
1819 /* Calculation of w_ik for complex argument */
1820
1821 test_wronskian(lambda([n,z], w_ik(n,z*%i)), lambda([n,z], 1/(z*%i)), 1e-12, 5);
1822 [];
1823
1824 /*******************************************************************************
1825
1826    Test Wronskian w_h1h2
1827
1828    A&S 9.1.17: H1[v+1](z)*H2[v](z)-H1[v](z)*H2[v+1](z) = -4*%i/(%pi*z)
1829
1830 *******************************************************************************/
1831
1832 w_h1h2(v,z) := hankel_1(v+1,z)*hankel_2(v,z) - hankel_1(v,z)*hankel_2(v+1,z);
1833 w_h1h2(v,z) := hankel_1(v+1,z)*hankel_2(v,z) - hankel_1(v,z)*hankel_2(v+1,z);
1834
1835 /* Calculation of w_h1h2 for real argument */
1836
1837 test_wronskian(w_h1h2, lambda([v,z], -4*%i/(%pi*z)), 1e-14, 5);
1838 [];
1839
1840 /* Calculation of w_h1h2 for complex argument */
1841
1842 test_wronskian(lambda([v,z], w_h1h2(v,z*%i)),
1843                lambda([v,z], -4/(%pi*z)),
1844                1e-13, 5);
1845 [];
1846
1847 /* Calculation of w_h1h2 for complex order and argument */
1848
1849 test_wronskian(lambda([v,z], w_h1h2(v*%i,z*%i)),
1850                lambda([v,z], -4/(%pi*z)),
1851                1e-10, 5);
1852 [];
1853
1854 /******************************************************************************
1855
1856     Integrals of Bessel functions
1857
1858 *******************************************************************************/
1859
1860 integrate(bessel_j(0,x),x);
1861 x*(bessel_j(0,x)*(2-%pi*struve_h(1,x))+%pi*bessel_j(1,x)*struve_h(0,x))/2;
1862
1863 integrate(bessel_j(1,x),x);
1864 -bessel_j(0,x);
1865
1866 integrate(bessel_j(2,x),x);
1867 hypergeometric([3/2],[5/2,3],-x^2/4)*x^3/24;
1868
1869 integrate(bessel_j(1/2,x),x);
1870  2^(3/2)*hypergeometric([3/4],[3/2,7/4],-x^2/4)*x^(3/2)/(3*sqrt(%pi));
1871
1872 /* http://functions.wolfram.com/Bessel-TypeFunctions/BesselY/21/01/01/ */
1873 integrate(bessel_y(0,x),x);
1874 %pi*x*(bessel_y(1,x)*struve_h(0,x)+bessel_y(0,x)*struve_h(-1,x))/2;
1875
1876 integrate(bessel_y(1,x),x);
1877 -bessel_y(0,x);
1878
1879 integrate(bessel_y(2,x),x);
1880 %pi*x*(bessel_y(1,x)*struve_h(0,x)+bessel_y(0,x)*struve_h(-1,x))/2-2*bessel_y(1,x);
1881
1882 integrate(bessel_y(3,x),x);
1883 -2*bessel_y(2,x)-bessel_y(0,x);
1884
1885 integrate(bessel_y(10,x),x);
1886 %pi*x*(bessel_y(1,x)*struve_h(0,x)+bessel_y(0,x)*struve_h(-1,x))/2
1887   -2*'sum(bessel_y(2*i+1,x),i,0,4);
1888
1889 integrate(bessel_y(11,x),x);
1890 -2*'sum(bessel_y(2*i,x),i,1,5)-bessel_y(0,x);
1891
1892 integrate(bessel_i(0,x),x);
1893 x*(bessel_i(0,x)*(%pi*struve_l(1,x)+2)-%pi*bessel_i(1,x)*struve_l(0,x))/2;
1894
1895 integrate(bessel_i(1,x),x);
1896 bessel_i(0,x);
1897
1898 integrate(bessel_i(2,x),x);
1899 hypergeometric([3/2],[5/2,3],x^2/4)*x^3/24;
1900
1901 integrate(bessel_i(1/2,x),x);
1902 2^(3/2)*hypergeometric([3/4],[3/2,7/4],x^2/4)*x^(3/2)/(3*sqrt(%pi));
1903
1904 /* http://functions.wolfram.com/Bessel-TypeFunctions/BesselK/21/01/01/ */
1905 integrate(bessel_k(0,x),x);
1906 %pi*x*(bessel_k(1,x)*struve_l(0,x)+bessel_k(0,x)*struve_l(-1,x))/2;
1907
1908 integrate(bessel_k(1,x),x);
1909 -bessel_k(0,x);
1910
1911 integrate(bessel_k(7,x),x);
1912 2*(-bessel_k(6,x)+bessel_k(4,x)-bessel_k(2,x))+bessel_k(0,x);
1913
1914 integrate(bessel_k(8,x),x);
1915 %pi*(struve_l(0,x)*bessel_k(1,x)+struve_l(-1,x)*bessel_k(0,x))*x/2
1916  +2*(-bessel_k(7,x)+bessel_k(5,x)-bessel_k(3,x)+bessel_k(1,x))$
1917
1918 /******************************************************************************
1919
1920     Check the handling of realpart and impagpart for special case
1921     of the order and arg of the Bessel functions.
1922
1923 *******************************************************************************/
1924
1925 /* Check for the Bessel J function */
1926
1927 (f(n,x):=[realpart(bessel_j(n,x)),imagpart(bessel_j(n,x))],done);
1928 done;
1929
1930 (declare(n,integer),assume(x>0),done);
1931 done;
1932
1933 f(n,x);
1934 [bessel_j(n,x),0];
1935 f(n,-x);
1936 [bessel_j(n,-x),0];
1937
1938 f(n+1/2,x);
1939 [bessel_j(n+1/2,x),0];
1940 f(n+1/2,-x);
1941 [0,-%i*bessel_j(n+1/2,-x)];
1942
1943 (declare(n_even,even),declare(n_odd,odd),%iargs:false,done);
1944 done;
1945
1946 f(n_even,%i);
1947 [bessel_j(n_even,%i),0];
1948 f(n_odd,%i);
1949 [0,-%i*bessel_j(n_odd,%i)];
1950
1951 f(n_even,x*%i);
1952 [bessel_j(n_even,x*%i),0];
1953 f(n_odd,x*%i);
1954 [0,-%i*bessel_j(n_odd,x*%i)];
1955
1956 f(n_even,(x+1)^2*%i);
1957 [bessel_j(n_even,(x+1)^2*%i),0];
1958 f(n_odd,(x+1)^2*%i);
1959 [0,-%i*bessel_j(n_odd,(x+1)^2*%i)];
1960
1961 (declare(j,imaginary),done);
1962 done;
1963
1964 f(n_even,(x+1)^2*j);
1965 [bessel_j(n_even,(x+1)^2*j),0];
1966 f(n_odd,(x+1)^2*j);
1967 [0,-%i*bessel_j(n_odd,(x+1)^2*j)];
1968
1969 /* Check the handling of realpart and imagpart for the Bessel I function */
1970
1971 (f(n,x):=[realpart(bessel_i(n,x)),imagpart(bessel_i(n,x))],done);
1972 done;
1973
1974 f(n,x);
1975 [bessel_i(n,x),0];
1976 f(n,-x);
1977 [bessel_i(n,-x),0];
1978
1979 f(n+1/2,x);
1980 [bessel_i(n+1/2,x),0];
1981 f(n+1/2,-x);
1982 [0,-%i*bessel_i(n+1/2,-x)];
1983
1984 f(n_even,%i);
1985 [bessel_i(n_even,%i),0];
1986 f(n_odd,%i);
1987 [0,-%i*bessel_i(n_odd,%i)];
1988
1989 f(n_even,x*%i);
1990 [bessel_i(n_even,x*%i),0];
1991 f(n_odd,x*%i);
1992 [0,-%i*bessel_i(n_odd,x*%i)];
1993
1994 f(n_even,(x+1)^2*%i);
1995 [bessel_i(n_even,(x+1)^2*%i),0];
1996 f(n_odd,(x+1)^2*%i);
1997 [0,-%i*bessel_i(n_odd,(x+1)^2*%i)];
1998
1999 f(n_even,(x+1)^2*j);
2000 [bessel_i(n_even,(x+1)^2*j),0];
2001 f(n_odd,(x+1)^2*j);
2002 [0,-%i*bessel_i(n_odd,(x+1)^2*j)];
2003
2004 /* Check the handling of realpart and imagpart for the Bessel K function */
2005
2006 (f(n,x):=[realpart(bessel_k(n,x)),imagpart(bessel_k(n,x))],done);
2007 done;
2008
2009 f(n,x);
2010 [bessel_k(n,x),0];
2011 f(n,-x);
2012 [realpart(bessel_k(n,-x)),imagpart(bessel_k(n,-x))];
2013
2014 f(n+1/2,x);
2015 [bessel_k(n+1/2,x),0];
2016 f(n+1/2,-x);
2017 [0,-%i*bessel_k(n+1/2,-x)];
2018
2019 f(n_even,%i);
2020 [realpart(bessel_k(n_even,%i)),imagpart(bessel_k(n_even,%i))];
2021 f(n_odd,%i);
2022 [realpart(bessel_k(n_odd,%i)),imagpart(bessel_k(n_odd,%i))];
2023
2024
2025 /* Check the handling of realpart and imagpart for the Bessel Y function */
2026
2027 (f(n,x):=[realpart(bessel_y(n,x)),imagpart(bessel_y(n,x))],done);
2028 done;
2029
2030 f(n,x);
2031 [bessel_y(n,x),0];
2032 f(n,-x);
2033 [realpart(bessel_y(n,-x)),imagpart(bessel_y(n,-x))];
2034
2035 f(n+1/2,x);
2036 [bessel_y(n+1/2,x),0];
2037 f(n+1/2,-x);
2038 [0,-%i*bessel_y(n+1/2,-x)];
2039
2040 f(n_even,%i);
2041 [realpart(bessel_y(n_even,%i)),imagpart(bessel_y(n_even,%i))];
2042 f(n_odd,%i);
2043 [realpart(bessel_y(n_odd,%i)),imagpart(bessel_y(n_odd,%i))];
2044
2045 /* %iargs is false, so transformation disabled */
2046 bessel_j(v, %i*x);
2047 bessel_j(v, %i*x);
2048
2049 bessel_i(v, %i*x);
2050 bessel_i(v, %i*x);
2051
2052 /* Set %iargs to true to enable transformation */
2053 (%iargs:true, bessel_j(v, %i*x));
2054 (%i*x)^v/x^v*bessel_i(v,x);
2055
2056 bessel_i(v, %i*x);
2057 (%i*x)^v/(x^v)*bessel_j(v, x);
2058
2059 imagpart(bessel_j(2,%i*3.0));
2060 0;
2061 realpart(bessel_j(3,%i*3.0));
2062 0;
2063
2064 imagpart(bessel_i(2,%i*3.0));
2065 0;
2066 realpart(bessel_i(3,%i*3.0));
2067 0;
2068
2069 /*******************************************************/
2070 /* Check the handling of conjugate on Bessel functions */
2071 /*******************************************************/
2072
2073 declare([w,z],complex,n,integer);
2074 done;
2075 assume(nonneg>=0,notequal(nonzero,0));
2076 [nonneg>=0,notequal(nonzero,0)];
2077
2078 conjugate(bessel_j(w,z));
2079 '(conjugate(bessel_j(w,z)));
2080
2081 conjugate(bessel_y(w,z));
2082 '(conjugate(bessel_y(w,z)));
2083
2084 conjugate(bessel_i(w,z));
2085 '(conjugate(bessel_i(w,z)));
2086
2087 conjugate(bessel_k(w,z));
2088 '(conjugate(bessel_k(w,z)));
2089
2090 conjugate(hankel_1(w,z));
2091 '(conjugate(hankel_1(w,z)));
2092
2093 conjugate(hankel_2(w,z));
2094 '(conjugate(hankel_2(w,z)));
2095
2096 /* Bessel functions with arguments off of the negative real axis
2097  * commute with conjugate
2098  */
2099 conjugate(bessel_j(z,x+%i*nonzero));
2100 '(bessel_j(conjugate(z),x-%i*nonzero));
2101
2102 conjugate(bessel_j(z,nonneg));
2103 '(bessel_j(conjugate(z),nonneg));
2104
2105 conjugate(bessel_y(z,x+%i*nonzero));
2106 '(bessel_y(conjugate(z),x-%i*nonzero));
2107
2108 conjugate(bessel_y(z,nonneg));
2109 '(bessel_y(conjugate(z),nonneg));
2110
2111 conjugate(bessel_i(z,x+%i*nonzero));
2112 '(bessel_i(conjugate(z),x-%i*nonzero));
2113
2114 conjugate(bessel_i(z,nonneg));
2115 '(bessel_i(conjugate(z),nonneg));
2116
2117 conjugate(bessel_k(z,x+%i*nonzero));
2118 '(bessel_k(conjugate(z),x-%i*nonzero));
2119
2120 conjugate(bessel_k(z,nonneg));
2121 '(bessel_k(conjugate(z),nonneg));
2122
2123 conjugate(hankel_1(z,x+%i*nonzero));
2124 '(hankel_2(conjugate(z),x-%i*nonzero));
2125
2126 conjugate(hankel_1(z,nonneg));
2127 '(hankel_2(conjugate(z),nonneg));
2128
2129 conjugate(hankel_2(z,x+%i*nonzero));
2130 '(hankel_1(conjugate(z),x-%i*nonzero));
2131
2132 conjugate(hankel_2(z,nonneg));
2133 '(hankel_1(conjugate(z),nonneg));
2134
2135 /* Bessel functions J and I of integral order commute with conjugate */
2136 conjugate(bessel_j(n,z));
2137 '(bessel_j(n,conjugate(z)));
2138
2139 conjugate(bessel_i(n,z));
2140 '(bessel_i(n,conjugate(z)));
2141
2142 remove([w,z],complex,n,integer);
2143 done;
2144 forget(nonneg>=0,notequal(nonzero,0));
2145 [nonneg>=0,notequal(nonzero,0)];
2146