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 kill(c1,c2,AS_10_4_10,AS_10_4_14,AS_10_4_15,AS_10_4_16,AS_10_4_17);
1087 done;
1088
1089 /* End of Airy function tests */
1090
1091 /* Numerical tests of gamma function. */
1092
1093 /* A&S Table 1.1, to 15 DP */
1094 closeto(gamma(1/2)-1.772453850905516,2e-15);
1095 true;
1096 closeto(gamma(1/3)-2.678938534707748,3.0e-15);
1097 true;
1098 closeto(gamma(7/4)-0.919062526848883,1e-15);
1099 true;
1100
1101 /* Complex values.  Checked against A&S Table 6.7 to 12 DP */
1102 closeto(gamma(1+%i)-(0.49801566811836-0.15494982830181*%i),1e-14);
1103 true;
1104 closeto(gamma(1+5*%i)-(-0.00169966449436-0.00135851941753*%i),1e-14);
1105 true;
1106 closeto(gamma(2+3*%i)-(-0.08239527266561+0.09177428743526*%i),1e-14);
1107 true;
1108
1109 /* Test numerical evaluation of Bessel functions
1110  * When the order is 0, and the arg is a float, we should produce a number.
1111  */
1112 closeto(bessel_j(0,1.0) - .7651976865579666, 1e-14);
1113 true;
1114
1115 closeto(bessel_y(0,1.0) - .08825696421567691, 1e-14);
1116 true;
1117
1118 closeto(bessel_i(0,1.0) - 1.266065877752009, 1e-14);
1119 true;
1120
1121 closeto(bessel_k(0,1.0) - .4210244382407085, 1e-14);
1122 true;
1123
1124 /*
1125  * Tests for failed cases to see if we're returning the noun forms
1126  */
1127 /* fail-on-f24p146test */
1128 specint(t^(-3/2)*exp(-t^2/8/a)*exp(-p*t),t);
1129 'specint(t^(-3/2)*exp(-t^2/8/a)*exp(-p*t),t);
1130
1131 /* fail-on-f35p147test
1132    Because we modified the construction of a noun form, we get a sligthly
1133    different noun form as result. DK. */
1134 specint((2*t)^(-3/2)*exp(-2*sqrt(a)*sqrt(t))*exp(-p*t),t);
1135 /*'specint(exp(-p*t -2*sqrt(a)*sqrt(t))/(8*t^(3/2)),t);*/
1136 'specint(%e^(-p*t-2*sqrt(a)*sqrt(t))/(2*sqrt(2)*t^(3/2)),t);
1137
1138 /* fail-on-f29p146test */
1139 specint(t^(v-1)*exp(1/8/t)*exp(-p*t),t);
1140 'specint(t^(v-1)*exp(1/8/t)*exp(-p*t),t);
1141
1142 /* fail-in-arbpow (aka f1p137test) */
1143 specint(t^(-1)*exp(-p*t),t);
1144 'specint(exp(-p*t)/t,t);
1145
1146 /* fail-in-f2p105vcond */
1147 (assume(p>3),0);
1148 0;
1149
1150 specint(8*t^1*exp(3*t)*bessel_y(3,t)*exp(-p*t),t);
1151 'specint(8*bessel_y(3,t)*t*%e^(3*t-p*t),t)$
1152
1153 /* fail-in-f50cond */
1154 specint(8*t^1*exp(3*t)*bessel_y(8,7*sqrt(t))*exp(-p*t),t);
1155 'specint(8*bessel_y(8,7*sqrt(t))*t*%e^(3*t-p*t),t);
1156
1157 /* fail-in-dionarghyp-y */
1158 specint(8*t^1*exp(3*t)*bessel_y(8,7*t^(3/2))*exp(-p*t),t);
1159 'specint(8*bessel_y(8,7*t^(3/2))*t*%e^(3*t-p*t),t);
1160
1161 /* The additionally phase factor in the calculation vanish, because of
1162    the modificaton of the transformation to Bessel J in the code. DK.
1163 */
1164 (assume(t>0,v2>1), radcan(specint(bessel_i(-v2,t)*exp(-p*t),t)));
1165 /*'specint((%e^((%i*%pi*v2-2*p*t)/2)*bessel_i(-v2,t))/(-1)^(v2/2),t);*/
1166 '(specint(bessel_i(-v2,t)*exp(-p*t),t));
1167
1168 /* Verify fix for [ 1877522 ] erf(1.0),nouns wrong; causes plot2d(erf) to fail
1169  */
1170 erf(1.0), nouns;
1171 0.8427007929497148;
1172
1173 erf(1.0), erf;
1174 0.8427007929497148;
1175
1176 erf(1.0);
1177 0.8427007929497148;
1178
1179 /* [ 789059 ] simpexpt problem: sign called on imag arg */
1180 (-(-1)^(1/6))^(1/2);
1181 sqrt(-(-1)^(1/6));
1182
1183 /* Further tests of bessel functions */
1184 /* The following numerical values are evaluated with the evaluation tool
1185    of the website www.functions.wolfram.com with a precision of 16 digits
1186    1998-2014 Wolfram Research, Inc. */
1187
1188 (test_bessel(actual, ref, digits) := closeto(realpart(actual)-realpart(ref), 10^(-digits)) and closeto(imagpart(actual)-imagpart(ref), 10^(-digits)), 0);
1189 0;
1190
1191 /* Numerical values for the bessel function J with negative order */
1192
1193 test_bessel(bessel_j(-1,-2.0),0.5767248077568734, 15);
1194 true;
1195
1196 test_bessel(bessel_j(-1,2.0), -0.5767248077568734, 15);
1197 true;
1198
1199 test_bessel(bessel_j(-1,-1.5), 0.5579365079100996, 15);
1200 true;
1201
1202 test_bessel(bessel_j(-1,1.5), -0.5579365079100996, 15);
1203 true;
1204
1205 test_bessel(bessel_j(-1.5, -2.0), -0.3956232813587035*%i, 15);
1206 true;
1207
1208 test_bessel(bessel_j(-1.5, 2.0), -0.3956232813587035, 15);
1209 true;
1210
1211 test_bessel
1212   (bessel_j(-1.8, -1.5),- 0.2033279902093184 - 0.1477264320209275 * %i,15);
1213 true;
1214
1215 test_bessel(bessel_j(-1.8, 1.5), -0.251327217627129314,15);
1216 true;
1217
1218 test_bessel(bessel_j(-2,-1.5), 0.2320876721442147, 15);
1219 true;
1220
1221 test_bessel(bessel_j(-2,1.5), 0.2320876721442147, 15);
1222 true;
1223
1224 test_bessel(bessel_j(-2.5,-1.5), -1.315037204805194 * %i,15);
1225 true;
1226
1227 test_bessel(bessel_j(-2.5,1.5), 1.315037204805194,15);
1228 true;
1229
1230 test_bessel
1231   (bessel_j(-2.3,-1.5), 0.5949438455752484 - 0.8188699527453657 * %i,14);
1232 true;
1233
1234 test_bessel(bessel_j(-2.3,1.5), 1.012178926325313, 14);
1235 true;
1236
1237 /* Numerical values for the bessel function J with positive order */
1238
1239 test_bessel(bessel_j(1.5,1.0), 0.2402978391234270, 15);
1240 true;
1241
1242 test_bessel(bessel_j(1.5,-1.0), -0.2402978391234270 * %i, 15);
1243 true;
1244
1245 test_bessel(bessel_j(1.8,1.0), 0.1564953153109239, 14);
1246 true;
1247
1248 test_bessel
1249   (bessel_j(1.8,-1.0), 0.1266073696266034 - 0.0919856383926216 * %i, 15);
1250 true;
1251
1252 test_bessel(bessel_j(2.0,1.0), 0.1149034849319005,15);
1253 true;
1254
1255 test_bessel(bessel_j(2.0,-1.0),0.1149034849319005,15);
1256 true;
1257
1258 test_bessel(bessel_j(2.5,1.0),0.04949681022847794,15);
1259 true;
1260
1261 test_bessel(bessel_j(2.5,-1.0),0.04949681022847794 * %i,15);
1262 true;
1263
1264 /* Numerical values for the bessel function J with complex arg
1265    and positive or negative order*/
1266
1267 test_bessel
1268   (bessel_j(0,1.0+%i),0.9376084768060293 - 0.4965299476091221 * %i,15);
1269 true;
1270
1271 test_bessel
1272   (bessel_j(1,1.0+%i),0.6141603349229036 + 0.3650280288270878 * %i,15);
1273 true;
1274
1275 test_bessel
1276   (bessel_j(-1,1.0+%i),-0.6141603349229036 - 0.3650280288270878 * %i,14);
1277 true;
1278
1279 test_bessel
1280   (bessel_j(2,1.0+%i),0.0415798869439621 + 0.2473976415133063 * %i,15);
1281 true;
1282
1283 test_bessel
1284   (bessel_j(-2,1.0+%i),0.0415798869439621 + 0.2473976415133063 * %i,15);
1285 true;
1286
1287 test_bessel
1288   (bessel_j(2.3,1.0+%i),-0.0141615213034667 + 0.1677798241687935 * %i,15);
1289 true;
1290
1291 test_bessel
1292   (bessel_j(-2.3,1.0+%i),0.1920598664138632 - 0.5158676904105332 * %i,14);
1293 true;
1294
1295 /* Numerical values for the bessel function J with complex order */
1296
1297 test_bessel
1298   (bessel_j(%i,1.0),1.641024179495082 - 0.437075010213683*%i,15);
1299 true;
1300
1301 test_bessel
1302   (bessel_j(%i,1.5),1.401883276281807 + 0.473362399311655*%i,15);
1303 true;
1304
1305 test_bessel
1306   (bessel_j(1.0+%i,-1.0),-0.01142279482478010 + 0.02390070064911069*%i,15);
1307 true;
1308
1309 test_bessel
1310   (bessel_j(1.5*%i,-2.0),0.01925195427338360 + 0.01442616961986814*%i,15);
1311 true;
1312
1313 /******************************************************************/
1314 /* Numerical values for the bessel function Y with negative order */
1315
1316 test_bessel
1317   (bessel_y(-1,-2.0),-0.1070324315409375 + 1.1534496155137468 * %i,14);
1318 true;
1319
1320 test_bessel(bessel_y(-1,2.0),0.1070324315409375,15);
1321 true;
1322
1323 test_bessel
1324   (bessel_y(-1,-1.5),-0.4123086269739113 + 1.1158730158201993 * %i,15);
1325 true;
1326
1327 test_bessel(bessel_y(-1,1.5),0.4123086269739113,15);
1328 true;
1329
1330 test_bessel(bessel_y(-1.5, -2.0),0.4912937786871623 * %i,15);
1331 true;
1332
1333 test_bessel(bessel_y(-1.5, 2.0),-0.4912937786871623,15);
1334 true;
1335
1336 test_bessel
1337   (bessel_y(-1.8, -1.5),-0.6777414340388011 + 0.0857519944018923 * %i,14);
1338 true;
1339
1340 test_bessel(bessel_y(-1.8, 1.5),-0.8377344836401481,14);
1341 true;
1342
1343 test_bessel
1344   (bessel_y(-2,-1.5),-0.9321937597629739 + 0.4641753442884295 * %i,15);
1345 true;
1346
1347 test_bessel(bessel_y(-2,1.5),-0.9321937597629739,14);
1348 true;
1349
1350 test_bessel(bessel_y(-2.5,-1.5),0.1244463597983876 * %i,14);
1351 true;
1352
1353 test_bessel(bessel_y(-2.5,1.5),0.1244463597983876,15);
1354 true;
1355
1356 test_bessel
1357   (bessel_y(-2.3,-1.5),-0.3148570865836879 + 0.7565240896444820 * %i,14);
1358 true;
1359
1360 test_bessel(bessel_y(-2.3,1.5),-0.5356668704355646,14);
1361 true;
1362
1363 /* Numerical values for the bessel function Y with positive order */
1364
1365 test_bessel(bessel_y(1.5,1.0),-1.102495575160179,14);
1366 true;
1367
1368 test_bessel(bessel_y(1.5,-1.0),-1.102495575160179 * %i,14);
1369 true;
1370
1371 test_bessel(bessel_y(1.8,1.0),-1.382351995367631,14);
1372 true;
1373
1374 test_bessel
1375   (bessel_y(1.8,-1.0),-1.1183462564605324 - 0.5593113771009602 * %i,14);
1376 true;
1377
1378 test_bessel(bessel_y(2.0,1.0),-1.650682606816254,14);
1379 true;
1380
1381 test_bessel
1382   (bessel_y(2.0,-1.0),-1.650682606816254 + 0.229806969863801 * %i,14);
1383 true;
1384
1385 test_bessel(bessel_y(2.5,1.0),-2.876387857462161,14);
1386 true;
1387
1388 test_bessel(bessel_y(2.5,-1.0),2.876387857462161 * %i,14);
1389 true;
1390
1391
1392 /* Numerical values for the bessel function Y with complex arg
1393    and positive or negative order*/
1394
1395 test_bessel
1396   (bessel_y(0,1.0+%i),0.4454744889360325 + 0.7101585820037345 * %i,15);
1397 true;
1398
1399 test_bessel
1400   (bessel_y(1,1.0+%i),-0.6576945355913452 + 0.6298010039928844 * %i,15);
1401 true;
1402
1403 test_bessel
1404 (bessel_y(-1,1.0+%i),0.6576945355913452 - 0.6298010039928844 * %i,15);
1405 true;
1406
1407 test_bessel
1408   (bessel_y(2,1.0+%i),-0.4733680205344934 + 0.5773369575804951 * %i,14);
1409 true;
1410
1411 test_bessel
1412   (bessel_y(-2,1.0+%i),-0.4733680205344934 + 0.5773369575804951 * %i,14);
1413 true;
1414
1415 test_bessel
1416   (bessel_y(2.3,1.0+%i),-0.2476879981252862 + 0.7595467103431256 * %i,15);
1417 true;
1418
1419 test_bessel
1420   (bessel_y(-2.3,1.0+%i),-0.1570442638685963 + 0.5821870838327466 * %i,14);
1421 true;
1422
1423 /* Numerical values for the bessel function Y with complex order */
1424
1425 test_bessel
1426   (bessel_y(%i,1.0),-0.476556612479964 - 1.505069159110387*%i,14);
1427 true;
1428
1429 test_bessel
1430   (bessel_y(%i,1.5),0.5161218926267688 - 1.2857405211747503*%i,14);
1431 true;
1432
1433 test_bessel
1434   (bessel_y(1.0+%i,-1.0),7.708594946281541 + 1.233384674244926*%i,14);
1435 true;
1436
1437 test_bessel
1438   (bessel_y(1.5*%i,-2.0),3.226466016458932 + 4.267260420563194*%i,13);
1439 true;
1440
1441 /******************************************************************/
1442 /* Numerical values for the bessel function I with negative order */
1443
1444 test_bessel(bessel_i(-1,-2.0),-1.590636854637329,15);
1445 true;
1446
1447 test_bessel(bessel_i(-1,2.0),1.590636854637329,15);
1448 true;
1449
1450 test_bessel(bessel_i(-1,-1.5),-0.9816664285779076,15);
1451 true;
1452
1453 test_bessel(bessel_i(-1,1.5),0.9816664285779076,15);
1454 true;
1455
1456 test_bessel(bessel_i(-1.5, -2.0),0.9849410530002364 * %i,14);
1457 true;
1458
1459 test_bessel(bessel_i(-1.5, 2.0),0.9849410530002364,14);
1460 true;
1461
1462 test_bessel
1463   (bessel_i(-1.8, -1.5),0.2026903980307014 + 0.1472631941876387 * %i,15);
1464 true;
1465
1466 test_bessel(bessel_i(-1.8, 1.5),0.2505391103524365,15);
1467 true;
1468
1469 test_bessel(bessel_i(-2,-1.5),0.3378346183356807,15);
1470 true;
1471
1472 test_bessel(bessel_i(-2,1.5),0.3378346183356807,15);
1473 true;
1474
1475 test_bessel(bessel_i(-2.5,-1.5),-0.8015666610717216*%i,14);
1476 true;
1477
1478 test_bessel(bessel_i(-2.5,1.5),0.8015666610717216,14);
1479 true;
1480
1481 test_bessel
1482   (bessel_i(-2.3,-1.5),0.3733945265830230 - 0.5139334755917659 * %i,15);
1483 true;
1484
1485 test_bessel(bessel_i(-2.3,1.5),0.6352567117441516,15);
1486 true;
1487
1488 /* Numerical values for the bessel function I with positive order */
1489
1490 test_bessel(bessel_i(1.5,1.0),0.2935253263474798,15);
1491 true;
1492
1493 test_bessel(bessel_i(1.5,-1.0),-0.2935253263474798*%i,13);
1494 true;
1495
1496 test_bessel(bessel_i(1.8,1.0),0.1871011888310777,15);
1497 true;
1498
1499 test_bessel
1500   (bessel_i(1.8,-1.0),0.1513680414320980 - 0.1099753194812967 * %i,14);
1501 true;
1502
1503 test_bessel(bessel_i(2.0,1.0),0.1357476697670383,15);
1504 true;
1505
1506 test_bessel(bessel_i(2.0,-1.0),0.1357476697670383,15);
1507 true;
1508
1509 test_bessel(bessel_i(2.5,1.0),0.05709890920304825,15);
1510 true;
1511
1512 test_bessel(bessel_i(2.5,-1.0),0.05709890920304825 * %i,15);
1513 true;
1514
1515
1516 /* Numerical values for the bessel function I with complex arg
1517    and positive or negative order*/
1518
1519 test_bessel
1520   (bessel_i(0,1.0+%i),0.9376084768060293 + 0.4965299476091221 * %i,15);
1521 true;
1522
1523 test_bessel
1524   (bessel_i(1,1.0+%i),0.3650280288270878 + 0.6141603349229036 * %i,15);
1525 true;
1526
1527 test_bessel
1528   (bessel_i(-1,1.0+%i),0.3650280288270878 + 0.6141603349229036 * %i,15);
1529 true;
1530
1531 test_bessel
1532   (bessel_i(2,1.0+%i),-0.0415798869439621 + 0.2473976415133063 * %i,15);
1533 true;
1534
1535 test_bessel
1536   (bessel_i(-2,1.0+%i),-0.0415798869439621 + 0.2473976415133063 * %i,15);
1537 true;
1538
1539 test_bessel
1540   (bessel_i(2.3,1.0+%i),-0.0635524383467825 + 0.1559221140952053 * %i,15);
1541 true;
1542
1543 test_bessel
1544   (bessel_i(-2.3,1.0+%i),-0.4053256245784623 - 0.3724481230406298 * %i,14);
1545 true;
1546
1547 /* Numerical values for the bessel function I with complex order */
1548
1549 test_bessel
1550   (bessel_i(%i,1.0),1.900799675819425 - 1.063960013554441*%i,14);
1551 true;
1552
1553 test_bessel
1554   (bessel_i(%i,1.5),2.495473638417463 - 0.601347501920535*%i,14);
1555 true;
1556
1557 test_bessel
1558   (bessel_i(1.0+%i,-1.0),-0.01096313515009349 + 0.03043920114776303*%i,15);
1559 true;
1560
1561 test_bessel
1562   (bessel_i(1.5*%i,-2.0),0.04238259669782487 - 0.01125055344512197*%i,15);
1563 true;
1564
1565 /******************************************************************/
1566 /* Numerical values for the bessel function K with negative order */
1567
1568 test_bessel
1569   (bessel_k(-1,-2.0),-0.139865881816522-4.997133057057809*%i,14);
1570 true;
1571
1572 test_bessel(bessel_k(-1,2.0),0.139865881816522,14);
1573 true;
1574
1575 test_bessel
1576   (bessel_k(-1,-1.5),-0.277387800456844 - 3.083996040296084*%i,14);
1577 true;
1578
1579 test_bessel(bessel_k(-1,1.5),0.2773878004568438,15);
1580 true;
1581
1582 test_bessel(bessel_k(-1.5, -2.0),-3.2741902342766302*%i,13);
1583 true;
1584
1585 test_bessel(bessel_k(-1.5, 2.0),0.1799066579520922,15);
1586 true;
1587
1588 test_bessel
1589   (bessel_k(-1.8, -1.5),0.3929372194683435 - 1.0725774293000646*%i,15);
1590 true;
1591
1592 test_bessel(bessel_k(-1.8, 1.5),0.4856971141526263,15);
1593 true;
1594
1595 test_bessel
1596   (bessel_k(-2,-1.5),0.5836559632566508 - 1.0613387550916862*%i,14);
1597 true;
1598
1599 test_bessel(bessel_k(-2,1.5),0.5836559632566508,15);
1600 true;
1601
1602 test_bessel
1603   (bessel_k(-2.3,-1.5),0.465598659425186 - 1.354876241730594*%i,14);
1604 true;
1605
1606 test_bessel(bessel_k(-2.3,1.5),0.7921237520153218,15);
1607 true;
1608
1609 /* Numerical values for the bessel function K with positive order */
1610
1611 test_bessel(bessel_k(1.5,1.0),0.9221370088957891,14);
1612 true;
1613
1614 /* kein Wert von function-site !? Ist das eine Nullstelle??? */
1615 test_bessel(bessel_k(1.5,-1.0),0,13);
1616 true;
1617
1618 test_bessel(bessel_k(1.8,1.0),1.275527037541854,15);
1619 true;
1620
1621 test_bessel
1622   (bessel_k(1.8,-1.0),1.0319230501560911 + 0.1619402612577788*%i,14);
1623 true;
1624
1625 test_bessel(bessel_k(2.0,1.0),1.624838898635177,14);
1626 true;
1627
1628 test_bessel(bessel_k(2.0,-1.0),1.624838898635177 - 0.426463882082061*%i,14);
1629 true;
1630
1631 test_bessel(bessel_k(2.5,1.0),3.227479531135262,14);
1632 true;
1633
1634 test_bessel(bessel_k(2.5,-1.0),-3.406861044815549*%i,14);
1635 true;
1636
1637 /* Numerical values for the bessel function K with complex arg
1638    and positive or negative order*/
1639
1640 test_bessel
1641   (bessel_k(0,1.0+%i),0.0801977269465178 - 0.3572774592853303*%i,15);
1642 true;
1643
1644 test_bessel
1645   (bessel_k(1,1.0+%i),0.0245683055237403 - 0.4597194738011894 * %i,15);
1646 true;
1647
1648 test_bessel
1649   (bessel_k(-1,1.0+%i),0.0245683055237403 - 0.4597194738011894 * %i,15);
1650 true;
1651
1652 test_bessel
1653   (bessel_k(2,1.0+%i),-0.3549534413309312 - 0.8415652386102600 * %i,13);
1654 true;
1655
1656 test_bessel
1657   (bessel_k(-2,1.0+%i),-0.3549534413309312 - 0.8415652386102600 * %i,13);
1658 true;
1659
1660 test_bessel
1661   (bessel_k(2.3,1.0+%i),-0.6635905911278042 - 1.0258894849569300 * %i,15);
1662 true;
1663
1664 test_bessel
1665   (bessel_k(-2.3,1.0+%i),-0.6635905911278042 - 1.0258894849569300 * %i,13);
1666 true;
1667
1668 /* Numerical values for the bessel function K with complex order */
1669
1670 test_bessel
1671   (bessel_k(%i,1.0),0.2894280370259921,14);
1672 true;
1673
1674 test_bessel
1675   (bessel_k(%i,1.5),0.1635839926633096,14);
1676 true;
1677
1678 test_bessel
1679   (bessel_k(1.0+%i,-1.0),-9.744252766030894 - 7.494570149760043*%i,14);
1680 true;
1681
1682 test_bessel
1683   (bessel_k(1.5*%i,-2.0),3.93512366711118 - 14.82183468316553*%i,13);
1684 true;
1685
1686 kill(closeto, test_bessel);
1687 done;
1688
1689 /* Numerical tests of the Bessel functions using the Wronskians
1690
1691    The Wronskians combines different types of Bessel functions and
1692    Bessel functions with negative and positve order.
1693    The results are very simple. Therefore the numerical calculation of the
1694    Wronskians is a good test for the different parts of the algorithmen.
1695
1696    Based on code by Dieter Kaiser.
1697 */
1698
1699 /* Test the Wronskian.  wf is the Wronskian function.  w_true is simplified Wronskian.
1700  * eps is the max absolute error allowed, and the Wronskian is tested
1701  * for values of the arg between -zlimit and zlimit.
1702  */
1703 (test_wronskian(wf, w_true, eps, zlimit) :=
1704 block([badpoints : [],
1705        ratprint : false,
1706        abserr : 0,
1707        maxerr : -1],
1708   for order:-1 thru 1 step 1/10 do
1709   (
1710     for z: -zlimit thru zlimit step 1 do
1711     (
1712       if notequal(z,0.0) then
1713       (
1714         result : float(rectform(wf(float(order),z))),
1715         answer : float(rectform(w_true(float(order),z))),
1716         abserr : abs(result - answer),
1717         maxerr : max(maxerr, abserr),
1718         if abserr > eps then
1719         (
1720           badpoints : cons([[order, z], result, answer, abserr], badpoints)
1721         )
1722       )
1723     )
1724   ),
1725   /*
1726    * For debugging, if there are any bad points, return the maximum error
1727    * found as the first element.
1728    */
1729   if badpoints # [] then
1730     cons(maxerr, badpoints)
1731   else
1732     badpoints
1733 ), 0);
1734 0;
1735
1736 /*******************************************************************************
1737
1738    Wronskian w_jj
1739
1740    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)
1741
1742 *******************************************************************************/
1743
1744 w_jj(n,z) := bessel_j(n+1,z)*bessel_j(-n,z) + bessel_j(n,z)*bessel_j(-n-1,z);
1745 w_jj(n,z) := bessel_j(n+1,z) *bessel_j(-n,z) + bessel_j(n,z)*bessel_j(-n-1,z);
1746
1747 /* Calculation of w_jj for real argument */
1748
1749 test_wronskian('w_jj, lambda([n,z], -2.0*sin(n*%pi)/(z*%pi)), 1e-14, 10);
1750 [];
1751
1752
1753 /* Calculation of w_jj for complex argument */
1754
1755 test_wronskian(lambda([n,z], expand(w_jj(n,%i*z))),
1756                lambda([n,z],-2.0*sin(n*%pi)/(%i*z*%pi)),
1757                1e-8, 10);
1758 [];
1759
1760 /*******************************************************************************
1761
1762    Wronskian w_jy
1763
1764    A&S 9.1.16: J[n+1](z)*Y[n](z)-J[n](z)*Y[n+1,z] = 2/(pi*z)
1765
1766 *******************************************************************************/
1767
1768 w_jy(n,z) := bessel_j(n+1,z)*bessel_y(n,z) - bessel_j(n,z)*bessel_y(n+1,z);
1769 w_jy(n,z) := bessel_j(n+1,z)*bessel_y(n,z) - bessel_j(n,z)*bessel_y(n+1,z);
1770
1771 /* Calculation of w_yj for real argument */
1772
1773 test_wronskian(w_jy, lambda([n,z], 2.0/(z*%pi)), 1e-14, 10);
1774 [];
1775
1776
1777 /* Calculation of w_jy for complex argument */
1778
1779 test_wronskian(lambda([n,z], w_jy(n,z*%i)), lambda([n,z],2.0/(z*%i*%pi)), 1e-8, 10);
1780 [];
1781
1782 /*******************************************************************************
1783
1784    Wronskian w_ii
1785
1786    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)
1787
1788 *******************************************************************************/
1789
1790 w_ii(n,z) := bessel_i(n,z)*bessel_i(-n-1,z) - bessel_i(n+1,z)*bessel_i(-n,z);
1791 w_ii(n,z) := bessel_i(n,z)*bessel_i(-n-1,z) - bessel_i(n+1,z)*bessel_i(-n,z);
1792
1793 /* Calculation of w_ii for real argument */
1794
1795 test_wronskian(w_ii, lambda([n,z],-2.0*sin(n*%pi)/(z*%pi)), 1e-10, 5);
1796 [];
1797
1798 /* Calculation of w_ii for complex argument */
1799
1800 test_wronskian(lambda([n,z], w_ii(n,z*%i)),
1801                lambda([n,z], -2.0*sin(n*%pi)/(z*%i*%pi)),
1802                1e-10, 5);
1803 [];
1804
1805 /*******************************************************************************
1806
1807    Test Wronskian w_ik
1808
1809    A&S 9.6.15: I[n](z)*K[n+1](z)+I[n+1](z)*K[n](z) = 1/z
1810
1811 *******************************************************************************/
1812
1813 w_ik(n,z) := bessel_i(n,z)*bessel_k(n+1,z) + bessel_i(n+1,z)*bessel_k(n,z);
1814 w_ik(n,z) := bessel_i(n,z)*bessel_k(n+1,z) + bessel_i(n+1,z)*bessel_k(n,z);
1815
1816 /* Calculation of w_ik for real argument */
1817
1818 test_wronskian(w_ik, lambda([n,z], 1/z), 1e-10, 5);
1819 [];
1820
1821 /* Calculation of w_ik for complex argument */
1822
1823 test_wronskian(lambda([n,z], w_ik(n,z*%i)), lambda([n,z], 1/(z*%i)), 1e-12, 5);
1824 [];
1825
1826 /*******************************************************************************
1827
1828    Test Wronskian w_h1h2
1829
1830    A&S 9.1.17: H1[v+1](z)*H2[v](z)-H1[v](z)*H2[v+1](z) = -4*%i/(%pi*z)
1831
1832 *******************************************************************************/
1833
1834 w_h1h2(v,z) := hankel_1(v+1,z)*hankel_2(v,z) - hankel_1(v,z)*hankel_2(v+1,z);
1835 w_h1h2(v,z) := hankel_1(v+1,z)*hankel_2(v,z) - hankel_1(v,z)*hankel_2(v+1,z);
1836
1837 /* Calculation of w_h1h2 for real argument */
1838
1839 test_wronskian(w_h1h2, lambda([v,z], -4*%i/(%pi*z)), 1e-14, 5);
1840 [];
1841
1842 /* Calculation of w_h1h2 for complex argument */
1843
1844 test_wronskian(lambda([v,z], w_h1h2(v,z*%i)),
1845                lambda([v,z], -4/(%pi*z)),
1846                1e-13, 5);
1847 [];
1848
1849 /* Calculation of w_h1h2 for complex order and argument */
1850
1851 test_wronskian(lambda([v,z], w_h1h2(v*%i,z*%i)),
1852                lambda([v,z], -4/(%pi*z)),
1853                1e-10, 5);
1854 [];
1855
1856 /******************************************************************************
1857
1858     Integrals of Bessel functions
1859
1860 *******************************************************************************/
1861
1862 integrate(bessel_j(0,x),x);
1863 x*(bessel_j(0,x)*(2-%pi*struve_h(1,x))+%pi*bessel_j(1,x)*struve_h(0,x))/2;
1864
1865 integrate(bessel_j(1,x),x);
1866 -bessel_j(0,x);
1867
1868 integrate(bessel_j(2,x),x);
1869 hypergeometric([3/2],[5/2,3],-x^2/4)*x^3/24;
1870
1871 integrate(bessel_j(1/2,x),x);
1872  2^(3/2)*hypergeometric([3/4],[3/2,7/4],-x^2/4)*x^(3/2)/(3*sqrt(%pi));
1873
1874 /* http://functions.wolfram.com/Bessel-TypeFunctions/BesselY/21/01/01/ */
1875 integrate(bessel_y(0,x),x);
1876 %pi*x*(bessel_y(1,x)*struve_h(0,x)+bessel_y(0,x)*struve_h(-1,x))/2;
1877
1878 integrate(bessel_y(1,x),x);
1879 -bessel_y(0,x);
1880
1881 integrate(bessel_y(2,x),x);
1882 %pi*x*(bessel_y(1,x)*struve_h(0,x)+bessel_y(0,x)*struve_h(-1,x))/2-2*bessel_y(1,x);
1883
1884 integrate(bessel_y(3,x),x);
1885 -2*bessel_y(2,x)-bessel_y(0,x);
1886
1887 integrate(bessel_y(10,x),x);
1888 %pi*x*(bessel_y(1,x)*struve_h(0,x)+bessel_y(0,x)*struve_h(-1,x))/2
1889   -2*'sum(bessel_y(2*i+1,x),i,0,4);
1890
1891 integrate(bessel_y(11,x),x);
1892 -2*'sum(bessel_y(2*i,x),i,1,5)-bessel_y(0,x);
1893
1894 integrate(bessel_i(0,x),x);
1895 x*(bessel_i(0,x)*(%pi*struve_l(1,x)+2)-%pi*bessel_i(1,x)*struve_l(0,x))/2;
1896
1897 integrate(bessel_i(1,x),x);
1898 bessel_i(0,x);
1899
1900 integrate(bessel_i(2,x),x);
1901 hypergeometric([3/2],[5/2,3],x^2/4)*x^3/24;
1902
1903 integrate(bessel_i(1/2,x),x);
1904 2^(3/2)*hypergeometric([3/4],[3/2,7/4],x^2/4)*x^(3/2)/(3*sqrt(%pi));
1905
1906 /* http://functions.wolfram.com/Bessel-TypeFunctions/BesselK/21/01/01/ */
1907 integrate(bessel_k(0,x),x);
1908 %pi*x*(bessel_k(1,x)*struve_l(0,x)+bessel_k(0,x)*struve_l(-1,x))/2;
1909
1910 integrate(bessel_k(1,x),x);
1911 -bessel_k(0,x);
1912
1913 integrate(bessel_k(7,x),x);
1914 2*(-bessel_k(6,x)+bessel_k(4,x)-bessel_k(2,x))+bessel_k(0,x);
1915
1916 integrate(bessel_k(8,x),x);
1917 %pi*(struve_l(0,x)*bessel_k(1,x)+struve_l(-1,x)*bessel_k(0,x))*x/2
1918  +2*(-bessel_k(7,x)+bessel_k(5,x)-bessel_k(3,x)+bessel_k(1,x))$
1919
1920 /******************************************************************************
1921
1922     Check the handling of realpart and impagpart for special case
1923     of the order and arg of the Bessel functions.
1924
1925 *******************************************************************************/
1926
1927 /* Check for the Bessel J function */
1928
1929 (f(n,x):=[realpart(bessel_j(n,x)),imagpart(bessel_j(n,x))],done);
1930 done;
1931
1932 (declare(n,integer),assume(x>0),done);
1933 done;
1934
1935 f(n,x);
1936 [bessel_j(n,x),0];
1937 f(n,-x);
1938 [bessel_j(n,-x),0];
1939
1940 f(n+1/2,x);
1941 [bessel_j(n+1/2,x),0];
1942 f(n+1/2,-x);
1943 [0,-%i*bessel_j(n+1/2,-x)];
1944
1945 (declare(n_even,even),declare(n_odd,odd),%iargs:false,done);
1946 done;
1947
1948 f(n_even,%i);
1949 [bessel_j(n_even,%i),0];
1950 f(n_odd,%i);
1951 [0,-%i*bessel_j(n_odd,%i)];
1952
1953 f(n_even,x*%i);
1954 [bessel_j(n_even,x*%i),0];
1955 f(n_odd,x*%i);
1956 [0,-%i*bessel_j(n_odd,x*%i)];
1957
1958 f(n_even,(x+1)^2*%i);
1959 [bessel_j(n_even,(x+1)^2*%i),0];
1960 f(n_odd,(x+1)^2*%i);
1961 [0,-%i*bessel_j(n_odd,(x+1)^2*%i)];
1962
1963 (declare(j,imaginary),done);
1964 done;
1965
1966 f(n_even,(x+1)^2*j);
1967 [bessel_j(n_even,(x+1)^2*j),0];
1968 f(n_odd,(x+1)^2*j);
1969 [0,-%i*bessel_j(n_odd,(x+1)^2*j)];
1970
1971 /* Check the handling of realpart and imagpart for the Bessel I function */
1972
1973 (f(n,x):=[realpart(bessel_i(n,x)),imagpart(bessel_i(n,x))],done);
1974 done;
1975
1976 f(n,x);
1977 [bessel_i(n,x),0];
1978 f(n,-x);
1979 [bessel_i(n,-x),0];
1980
1981 f(n+1/2,x);
1982 [bessel_i(n+1/2,x),0];
1983 f(n+1/2,-x);
1984 [0,-%i*bessel_i(n+1/2,-x)];
1985
1986 f(n_even,%i);
1987 [bessel_i(n_even,%i),0];
1988 f(n_odd,%i);
1989 [0,-%i*bessel_i(n_odd,%i)];
1990
1991 f(n_even,x*%i);
1992 [bessel_i(n_even,x*%i),0];
1993 f(n_odd,x*%i);
1994 [0,-%i*bessel_i(n_odd,x*%i)];
1995
1996 f(n_even,(x+1)^2*%i);
1997 [bessel_i(n_even,(x+1)^2*%i),0];
1998 f(n_odd,(x+1)^2*%i);
1999 [0,-%i*bessel_i(n_odd,(x+1)^2*%i)];
2000
2001 f(n_even,(x+1)^2*j);
2002 [bessel_i(n_even,(x+1)^2*j),0];
2003 f(n_odd,(x+1)^2*j);
2004 [0,-%i*bessel_i(n_odd,(x+1)^2*j)];
2005
2006 /* Check the handling of realpart and imagpart for the Bessel K function */
2007
2008 (f(n,x):=[realpart(bessel_k(n,x)),imagpart(bessel_k(n,x))],done);
2009 done;
2010
2011 f(n,x);
2012 [bessel_k(n,x),0];
2013 f(n,-x);
2014 [realpart(bessel_k(n,-x)),imagpart(bessel_k(n,-x))];
2015
2016 f(n+1/2,x);
2017 [bessel_k(n+1/2,x),0];
2018 f(n+1/2,-x);
2019 [0,-%i*bessel_k(n+1/2,-x)];
2020
2021 f(n_even,%i);
2022 [realpart(bessel_k(n_even,%i)),imagpart(bessel_k(n_even,%i))];
2023 f(n_odd,%i);
2024 [realpart(bessel_k(n_odd,%i)),imagpart(bessel_k(n_odd,%i))];
2025
2026
2027 /* Check the handling of realpart and imagpart for the Bessel Y function */
2028
2029 (f(n,x):=[realpart(bessel_y(n,x)),imagpart(bessel_y(n,x))],done);
2030 done;
2031
2032 f(n,x);
2033 [bessel_y(n,x),0];
2034 f(n,-x);
2035 [realpart(bessel_y(n,-x)),imagpart(bessel_y(n,-x))];
2036
2037 f(n+1/2,x);
2038 [bessel_y(n+1/2,x),0];
2039 f(n+1/2,-x);
2040 [0,-%i*bessel_y(n+1/2,-x)];
2041
2042 f(n_even,%i);
2043 [realpart(bessel_y(n_even,%i)),imagpart(bessel_y(n_even,%i))];
2044 f(n_odd,%i);
2045 [realpart(bessel_y(n_odd,%i)),imagpart(bessel_y(n_odd,%i))];
2046
2047 /* %iargs is false, so transformation disabled */
2048 bessel_j(v, %i*x);
2049 bessel_j(v, %i*x);
2050
2051 bessel_i(v, %i*x);
2052 bessel_i(v, %i*x);
2053
2054 /* Set %iargs to true to enable transformation */
2055 (%iargs:true, bessel_j(v, %i*x));
2056 (%i*x)^v/x^v*bessel_i(v,x);
2057
2058 bessel_i(v, %i*x);
2059 (%i*x)^v/(x^v)*bessel_j(v, x);
2060
2061 imagpart(bessel_j(2,%i*3.0));
2062 0;
2063 realpart(bessel_j(3,%i*3.0));
2064 0;
2065
2066 imagpart(bessel_i(2,%i*3.0));
2067 0;
2068 realpart(bessel_i(3,%i*3.0));
2069 0;
2070
2071 /*******************************************************/
2072 /* Check the handling of conjugate on Bessel functions */
2073 /*******************************************************/
2074
2075 declare([w,z],complex,n,integer);
2076 done;
2077 assume(nonneg>=0,notequal(nonzero,0));
2078 [nonneg>=0,notequal(nonzero,0)];
2079
2080 conjugate(bessel_j(w,z));
2081 '(conjugate(bessel_j(w,z)));
2082
2083 conjugate(bessel_y(w,z));
2084 '(conjugate(bessel_y(w,z)));
2085
2086 conjugate(bessel_i(w,z));
2087 '(conjugate(bessel_i(w,z)));
2088
2089 conjugate(bessel_k(w,z));
2090 '(conjugate(bessel_k(w,z)));
2091
2092 conjugate(hankel_1(w,z));
2093 '(conjugate(hankel_1(w,z)));
2094
2095 conjugate(hankel_2(w,z));
2096 '(conjugate(hankel_2(w,z)));
2097
2098 /* Bessel functions with arguments off of the negative real axis
2099  * commute with conjugate
2100  */
2101 conjugate(bessel_j(z,x+%i*nonzero));
2102 '(bessel_j(conjugate(z),x-%i*nonzero));
2103
2104 conjugate(bessel_j(z,nonneg));
2105 '(bessel_j(conjugate(z),nonneg));
2106
2107 conjugate(bessel_y(z,x+%i*nonzero));
2108 '(bessel_y(conjugate(z),x-%i*nonzero));
2109
2110 conjugate(bessel_y(z,nonneg));
2111 '(bessel_y(conjugate(z),nonneg));
2112
2113 conjugate(bessel_i(z,x+%i*nonzero));
2114 '(bessel_i(conjugate(z),x-%i*nonzero));
2115
2116 conjugate(bessel_i(z,nonneg));
2117 '(bessel_i(conjugate(z),nonneg));
2118
2119 conjugate(bessel_k(z,x+%i*nonzero));
2120 '(bessel_k(conjugate(z),x-%i*nonzero));
2121
2122 conjugate(bessel_k(z,nonneg));
2123 '(bessel_k(conjugate(z),nonneg));
2124
2125 conjugate(hankel_1(z,x+%i*nonzero));
2126 '(hankel_2(conjugate(z),x-%i*nonzero));
2127
2128 conjugate(hankel_1(z,nonneg));
2129 '(hankel_2(conjugate(z),nonneg));
2130
2131 conjugate(hankel_2(z,x+%i*nonzero));
2132 '(hankel_1(conjugate(z),x-%i*nonzero));
2133
2134 conjugate(hankel_2(z,nonneg));
2135 '(hankel_1(conjugate(z),nonneg));
2136
2137 /* Bessel functions J and I of integral order commute with conjugate */
2138 conjugate(bessel_j(n,z));
2139 '(bessel_j(n,conjugate(z)));
2140
2141 conjugate(bessel_i(n,z));
2142 '(bessel_i(n,conjugate(z)));
2143
2144 remove([w,z],complex,n,integer);
2145 done;
2146 forget(nonneg>=0,notequal(nonzero,0));
2147 [nonneg>=0,notequal(nonzero,0)];
2148
2149 /* mailing list 2016-01-06: "coerce-float-fun and hashed arrays" */
2150
2151 (kill(a, s, t), a:make_array(hashed), a[s]:5, 0);
2152 0;
2153
2154 a[s];
2155 5;
2156
2157 a[s], nouns;
2158 5;
2159
2160 a[t];
2161 false;
2162
2163 a[t], nouns;
2164 false;
2165