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