1jn(3,4); 2jn(3,4); /*0.1320341839226408 ;*/ 3j0(1); 4j0(1); /*0.7651976865579665 ;*/ 5bessel_j(0,1.0); 60.7651976865579665 ; 7 8/* BESSEL is gone. This makes sure it's gone. */ 9bessel(2,3); 10bessel(2,3)$ 11 12bessel(2,3.0); 13bessel(2, 3.0); 14 15bessel_j(3,2.0); 160.1289432494744021; 17 18bessel_j(-5,x); 19-bessel_j(5,x); 20bessel_j(-n,x); 21bessel_j(-n,x); 22 23diff(bessel_j(0,x),x); 24-bessel_j(1,x); 25diff(bessel_j(1,x),x); 26(bessel_j(0,x)-bessel_j(2,x))/2; 27diff(bessel_y(0,x),x); 28-bessel_y(1,x); 29diff(bessel_y(1,x),x); 30(bessel_y(0,x)-bessel_y(2,x))/2; 31diff(bessel_i(0,x),x); 32bessel_i(1,x); 33diff(bessel_i(1,x),x); 34(bessel_i(2,x)+bessel_i(0,x))/2; 35diff(bessel_k(0,x),x); 36-bessel_k(1,x); 37diff(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 */ 41diff(bessel_y(v,z),v); 42cot(%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 44besselexpand:true; 45true$ 46 47bessel_j(1/2,x); 48sqrt(2/%pi)*sin(x)/sqrt(x); 49bessel_j(-1/2,x); 50sqrt(2/%pi)*cos(x)/sqrt(x); 51 52bessel_j(3/2,x); 53(sqrt(2)*sqrt(x)*(sin(x)/x^2-cos(x)/x))/sqrt(%pi); 54 55bessel_j(5/2,x); 56(sqrt(2)*sqrt(x)*((3/x^3-1/x)*sin(x)-(3*cos(x))/x^2))/sqrt(%pi)$ 57 58ratsimp(bessel_j(-3/2,x) - (2*(-1/2)/x*bessel_j(-1/2,x)-bessel_j(1/2,x))); 590$ 60 61ratsimp(bessel_j(-5/2,x) - (2*(-3/2)/x*bessel_j(-3/2,x)-bessel_j(-1/2,x))); 620$ 63 64ratsimp(bessel_y(1/2,x) + sqrt(2/%pi)*cos(x)/sqrt(x)); 650$ 66 67ratsimp(bessel_y(-1/2,x) - sqrt(2/%pi)*sin(x)/sqrt(x)); 680$ 69 70ratsimp(bessel_y(3/2,x) - (2*(1/2)/x*bessel_y(1/2,x)-bessel_y(-1/2,x))); 710$ 72 73ratsimp(bessel_y(5/2,x) - (2*(3/2)/x*bessel_y(3/2,x)-bessel_y(1/2,x))); 740$ 75 76ratsimp(bessel_y(-3/2,x) - (2*(-1/2)/x*bessel_y(-1/2,x)-bessel_y(1/2,x))); 770$ 78 79ratsimp(bessel_y(-5/2,x) - (2*(-3/2)/x*bessel_y(-3/2,x)-bessel_y(-1/2,x))); 800$ 81 82bessel_i(1/2,x); 83sqrt(2*x/%pi)*sinh(x)/x; 84 85bessel_i(-1/2,x); 86sqrt(2*x/%pi)*cosh(x)/x; 87 88ratsimp(bessel_i(3/2,x) - (-2*(1/2)/x*bessel_i(1/2,x)+bessel_i(-1/2,x))); 890$ 90 91ratsimp(bessel_i(5/2,x) -(-2*(3/2)/x*bessel_i(3/2,x)+bessel_i(1/2,x))); 920$ 93 94ratsimp(bessel_i(-3/2,x) -(2*(-1/2)/x*bessel_i(-1/2,x)+bessel_i(1/2,x))); 950$ 96 97ratsimp(bessel_i(-5/2,x) - (2*(-3/2)/x*bessel_i(-3/2,x)+bessel_i(-1/2,x))); 980$ 99 100ratsimp(bessel_k(1/2,x) - sqrt(%pi/2)*%e^(-x)/sqrt(x)); 1010$ 102 103ratsimp(bessel_k(-1/2,x)- sqrt(%pi/2)*%e^(-x)/sqrt(x)); 1040$ 105 106ratsimp(bessel_k(3/2,x) - (2*(1/2)/x*bessel_k(1/2,x)+bessel_k(-1/2,x))); 1070$ 108 109ratsimp(bessel_k(5/2,x) -(2*(3/2)/x*bessel_k(3/2,x)+bessel_k(1/2,x))); 1100$ 111 112ratsimp(bessel_k(-3/2,x) -(-2*(-1/2)/x*bessel_k(-1/2,x)+bessel_k(1/2,x))); 1130$ 114 115ratsimp(bessel_k(-5/2,x) -(-2*(-3/2)/x*bessel_k(-3/2,x)+bessel_k(-1/2,x))); 1160$ 117 118(assume(p>0),true); 119true$ 120(assume(4*p+a>0),true); 121true$ 122besselexpand:false; 123false$ 124 125specint(t^(1/2)*%e^(-a*t/4)*%e^(-p*t),t); 126sqrt(%pi)/(2*(p+a/4)^(3/2)); 127 128prefer_whittaker:true; 129true$ 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); 172true$ 173 174specint(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 1843*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 */ 196ratsimp(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/* 216radcan(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 232specint(t^(1/2)*bessel_j(1,2*a^(1/2)*t^(1/2))*%e^(-p*t),t); 233sqrt(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 */ 247factor(ratsimp(specint(t^2*bessel_j(1,a*t)*%e^(-p*t),t))); 2483*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(%%/%)); 2601$ 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 */ 293ratsimp(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 */ 304specint(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 $ */ 306specint(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 337specint(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 */ 362result: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 */ 370expand(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 431specint(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 */ 441factor(ratsimp(specint(t^(1/2)*bessel_i(1,t)*%e^(-p*t),t))); 4423*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 466ratsimp(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 523ratsimp(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 564ratsimp(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 */ 587factor(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 */ 614ratsimp(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); 618true; 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 */ 640specint(%e^(a*t)*t^2*erf(t^(1/2))*%e^(-p*t),t); 64115*(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); 657true$ 658 659radcan(specint(bessel_j(v,a*t)*exp(-p*t),t)); 660a^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 */ 668radcan(specint(bessel_j(v,a*t)/t*exp(-p*t),t)); 669a^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 */ 678radcan(specint(t^v*bessel_j(v,a*t)*exp(-p*t),t)); 679a^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); 686true$ 687(assume(a>0),true); 688true$ 689 690radcan(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. 693assoc_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 697assoc_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 */ 704specint(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 */ 711specint(sqrt(t)*bessel_j(1,2*sqrt(a)*sqrt(t))*exp(-p*t),t); 712sqrt(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 */ 719specint(t^(-1/2)*bessel_j(1,2*sqrt(a)*sqrt(t))*exp(-p*t),t); 720sqrt(%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 */ 727specint(t^(v/2)*bessel_j(v,2*sqrt(a)*sqrt(t))*exp(-p*t),t); 728a^(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 */ 736specint(t^(-v/2)*bessel_j(v,2*sqrt(a)*sqrt(t))*exp(-p*t),t); 737p^(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 */ 744specint(t^(v/2-1)*bessel_j(v,2*sqrt(a)*sqrt(t))*exp(-p*t),t); 745v*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 */ 777prefer_whittaker:true; 778true$ 779(assume(2*v+2*u+1>0),true); 780true$ 781 782specint(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 */ 791prefer_whittaker:false; 792false$ 793(assume(v+u>0),true); 794true$ 795 796specint(t^(u-1)*bessel_j(2*v,2*sqrt(a)*sqrt(t))*exp(-p*t),t); 797a^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 */ 806expand(factor(radcan(specint(exp(-p*t)*bessel_y(1/6,a*t),t)))); 807sqrt(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); 811true$ 812expand(factor(radcan(specint(exp(-p*t)*bessel_y(v1,a*t),t)))); 813a^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); 826true$ 827specint(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); 828a^(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 */ 850specint(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 */ 878factor(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 */ 913factor(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 940hgfred([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); 9430$ 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 954hgfred([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); 9560$ 957 958/* Tests for Airy functions */ 959closeto(e,tol):=block([numer:true,abse],abse:abs(e),if(abse<tol) then true else abse); 960closeto(e,tol):=block([numer:true,abse],abse:abs(e),if(abse<tol) then true else abse); 961 962/* Derivatives of Airy functions */ 963diff(airy_ai(x),x); 964airy_dai(x); 965diff(airy_dai(x),x); 966x*airy_ai(x); 967diff(airy_bi(x),x); 968airy_dbi(x); 969diff(airy_dbi(x),x); 970x*airy_bi(x); 971 972/* Integrals of Airy functions */ 973integrate(airy_ai(z),z); 974hypergeometric([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); 976integrate(airy_dai(z),z); 977airy_ai(z); 978integrate(airy_bi(z),z); 9793^(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)); 981integrate(airy_dbi(z),z); 982airy_bi(z); 983 984/* A&S 10.4.1 Airy functions satisfy the Airy differential equation */ 985diff(airy_ai(x),x,2)-x*airy_ai(x); 9860; 987diff(airy_bi(x),x,2)-x*airy_bi(x); 9880; 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)); 992true; 993closeto(airy_bi(0)/sqrt(3)-c1,1.0e-15); 994true; 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)); 998true; 999closeto(airy_dbi(0)/sqrt(3)-c2,1.0e-14); 1000true; 1001 1002/* Exact values A&S 10.4.4 and 10.4.5 */ 1003airy_ai(0); 10041/(3^(2/3)*gamma(2/3)); 1005airy_dai(0); 1006-(1/(3^(1/3)*gamma(1/3))); 1007airy_bi(0); 10081/(3^(1/6)*gamma(2/3)); 1009airy_dbi(0); 10103^(1/6)/gamma(1/3); 1011 1012/* A&S 10.4.10 - Wronskian */ 1013AS_10_4_10(z):=expand(airy_ai(z)*airy_dbi(z)-airy_bi(z)*airy_dai(z)-1/%pi); 1014AS_10_4_10(z):=expand(airy_ai(z)*airy_dbi(z)-airy_bi(z)*airy_dai(z)-1/%pi); 1015closeto(AS_10_4_10(1),1.0e-15); 1016true; 1017closeto(AS_10_4_10(1+%i),1.0e-15); 1018true; 1019closeto(AS_10_4_10(%i),1.0e-15); 1020true; 1021closeto(AS_10_4_10(-1+%i),2.0e-15); 1022true; 1023closeto(AS_10_4_10(-1),1.0e-15); 1024true; 1025closeto(AS_10_4_10(-1-%i),2.0e-15); 1026true; 1027closeto(AS_10_4_10(-%i),1.0e-15); 1028true; 1029closeto(AS_10_4_10(1-%i),1.0e-15); 1030true; 1031 1032/* A&S 10.4.14 - only for z>0 ? */ 1033AS_10_4_14(z):=block([y:(2/3)*(z)^(3/2)],airy_ai(z)-(sqrt(z/3)*bessel_k(1/3,y)/%pi)); 1034AS_10_4_14(z):=block([y:(2/3)*(z)^(3/2)],airy_ai(z)-(sqrt(z/3)*bessel_k(1/3,y)/%pi)); 1035closeto(AS_10_4_14(1),1.0e-15); 1036true; 1037closeto(AS_10_4_14(2),1.0e-15); 1038true; 1039closeto(AS_10_4_14(5),1.0e-15); 1040true; 1041closeto(AS_10_4_14(10),1.0e-15); 1042true; 1043 1044/* A&S 10.4.15 - only for z<0 ? */ 1045AS_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))); 1046AS_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))); 1047closeto(AS_10_4_15(-1),1.0e-15); 1048true; 1049closeto(AS_10_4_15(-2),1.0e-15); 1050true; 1051closeto(AS_10_4_15(-5),1.0e-14); 1052true; 1053closeto(AS_10_4_15(-10),1.0e-15); 1054true; 1055 1056/* A&S 10.4.16 - only for z>0 ? */ 1057AS_10_4_16(z):=block([y:(2/3)*(z)^(3/2)],airy_dai(z)+(z/sqrt(3))*bessel_k(2/3,y)/%pi); 1058AS_10_4_16(z):=block([y:(2/3)*(z)^(3/2)],airy_dai(z)+(z/sqrt(3))*bessel_k(2/3,y)/%pi); 1059closeto(AS_10_4_16(1),1.0e-15); 1060true; 1061closeto(AS_10_4_16(2),1.0e-15); 1062true; 1063closeto(AS_10_4_16(5),1.0e-15); 1064true; 1065closeto(AS_10_4_16(10),1.0e-15); 1066true; 1067 1068/* A&S 10.4.17 - only for z<0 ?, Appears to be a sign error in A&S */ 1069AS_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))); 1070AS_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))); 1071closeto(AS_10_4_17(-1),1.0e-15); 1072true; 1073closeto(AS_10_4_17(-2),1.0e-15); 1074true; 1075closeto(AS_10_4_17(-5),1.0e-14); 1076true; 1077closeto(AS_10_4_17(-10),1.0e-15); 1078true; 1079 1080/* Test that complex float arguments are evaluated */ 1081airy_ai(%i); 1082airy_ai(%i); 1083floatnump(realpart(airy_ai(1.0*%i))); 1084true; 1085 1086kill(c1,c2,AS_10_4_10,AS_10_4_14,AS_10_4_15,AS_10_4_16,AS_10_4_17); 1087done; 1088 1089/* End of Airy function tests */ 1090 1091/* Numerical tests of gamma function. */ 1092 1093/* A&S Table 1.1, to 15 DP */ 1094closeto(gamma(1/2)-1.772453850905516,2e-15); 1095true; 1096closeto(gamma(1/3)-2.678938534707748,3.0e-15); 1097true; 1098closeto(gamma(7/4)-0.919062526848883,1e-15); 1099true; 1100 1101/* Complex values. Checked against A&S Table 6.7 to 12 DP */ 1102closeto(gamma(1+%i)-(0.49801566811836-0.15494982830181*%i),1e-14); 1103true; 1104closeto(gamma(1+5*%i)-(-0.00169966449436-0.00135851941753*%i),1e-14); 1105true; 1106closeto(gamma(2+3*%i)-(-0.08239527266561+0.09177428743526*%i),1e-14); 1107true; 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 */ 1112closeto(bessel_j(0,1.0) - .7651976865579666, 1e-14); 1113true; 1114 1115closeto(bessel_y(0,1.0) - .08825696421567691, 1e-14); 1116true; 1117 1118closeto(bessel_i(0,1.0) - 1.266065877752009, 1e-14); 1119true; 1120 1121closeto(bessel_k(0,1.0) - .4210244382407085, 1e-14); 1122true; 1123 1124/* 1125 * Tests for failed cases to see if we're returning the noun forms 1126 */ 1127/* fail-on-f24p146test */ 1128specint(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. */ 1134specint((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 */ 1139specint(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) */ 1143specint(t^(-1)*exp(-p*t),t); 1144'specint(exp(-p*t)/t,t); 1145 1146/* fail-in-f2p105vcond */ 1147(assume(p>3),0); 11480; 1149 1150specint(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 */ 1154specint(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 */ 1158specint(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 */ 1170erf(1.0), nouns; 11710.8427007929497148; 1172 1173erf(1.0), erf; 11740.8427007929497148; 1175 1176erf(1.0); 11770.8427007929497148; 1178 1179/* [ 789059 ] simpexpt problem: sign called on imag arg */ 1180(-(-1)^(1/6))^(1/2); 1181sqrt(-(-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); 11890; 1190 1191/* Numerical values for the bessel function J with negative order */ 1192 1193test_bessel(bessel_j(-1,-2.0),0.5767248077568734, 15); 1194true; 1195 1196test_bessel(bessel_j(-1,2.0), -0.5767248077568734, 15); 1197true; 1198 1199test_bessel(bessel_j(-1,-1.5), 0.5579365079100996, 15); 1200true; 1201 1202test_bessel(bessel_j(-1,1.5), -0.5579365079100996, 15); 1203true; 1204 1205test_bessel(bessel_j(-1.5, -2.0), -0.3956232813587035*%i, 15); 1206true; 1207 1208test_bessel(bessel_j(-1.5, 2.0), -0.3956232813587035, 15); 1209true; 1210 1211test_bessel 1212 (bessel_j(-1.8, -1.5),- 0.2033279902093184 - 0.1477264320209275 * %i,15); 1213true; 1214 1215test_bessel(bessel_j(-1.8, 1.5), -0.251327217627129314,15); 1216true; 1217 1218test_bessel(bessel_j(-2,-1.5), 0.2320876721442147, 15); 1219true; 1220 1221test_bessel(bessel_j(-2,1.5), 0.2320876721442147, 15); 1222true; 1223 1224test_bessel(bessel_j(-2.5,-1.5), -1.315037204805194 * %i,15); 1225true; 1226 1227test_bessel(bessel_j(-2.5,1.5), 1.315037204805194,15); 1228true; 1229 1230test_bessel 1231 (bessel_j(-2.3,-1.5), 0.5949438455752484 - 0.8188699527453657 * %i,14); 1232true; 1233 1234test_bessel(bessel_j(-2.3,1.5), 1.012178926325313, 14); 1235true; 1236 1237/* Numerical values for the bessel function J with positive order */ 1238 1239test_bessel(bessel_j(1.5,1.0), 0.2402978391234270, 15); 1240true; 1241 1242test_bessel(bessel_j(1.5,-1.0), -0.2402978391234270 * %i, 15); 1243true; 1244 1245test_bessel(bessel_j(1.8,1.0), 0.1564953153109239, 14); 1246true; 1247 1248test_bessel 1249 (bessel_j(1.8,-1.0), 0.1266073696266034 - 0.0919856383926216 * %i, 15); 1250true; 1251 1252test_bessel(bessel_j(2.0,1.0), 0.1149034849319005,15); 1253true; 1254 1255test_bessel(bessel_j(2.0,-1.0),0.1149034849319005,15); 1256true; 1257 1258test_bessel(bessel_j(2.5,1.0),0.04949681022847794,15); 1259true; 1260 1261test_bessel(bessel_j(2.5,-1.0),0.04949681022847794 * %i,15); 1262true; 1263 1264/* Numerical values for the bessel function J with complex arg 1265 and positive or negative order*/ 1266 1267test_bessel 1268 (bessel_j(0,1.0+%i),0.9376084768060293 - 0.4965299476091221 * %i,15); 1269true; 1270 1271test_bessel 1272 (bessel_j(1,1.0+%i),0.6141603349229036 + 0.3650280288270878 * %i,15); 1273true; 1274 1275test_bessel 1276 (bessel_j(-1,1.0+%i),-0.6141603349229036 - 0.3650280288270878 * %i,14); 1277true; 1278 1279test_bessel 1280 (bessel_j(2,1.0+%i),0.0415798869439621 + 0.2473976415133063 * %i,15); 1281true; 1282 1283test_bessel 1284 (bessel_j(-2,1.0+%i),0.0415798869439621 + 0.2473976415133063 * %i,15); 1285true; 1286 1287test_bessel 1288 (bessel_j(2.3,1.0+%i),-0.0141615213034667 + 0.1677798241687935 * %i,15); 1289true; 1290 1291test_bessel 1292 (bessel_j(-2.3,1.0+%i),0.1920598664138632 - 0.5158676904105332 * %i,14); 1293true; 1294 1295/* Numerical values for the bessel function J with complex order */ 1296 1297test_bessel 1298 (bessel_j(%i,1.0),1.641024179495082 - 0.437075010213683*%i,15); 1299true; 1300 1301test_bessel 1302 (bessel_j(%i,1.5),1.401883276281807 + 0.473362399311655*%i,15); 1303true; 1304 1305test_bessel 1306 (bessel_j(1.0+%i,-1.0),-0.01142279482478010 + 0.02390070064911069*%i,15); 1307true; 1308 1309test_bessel 1310 (bessel_j(1.5*%i,-2.0),0.01925195427338360 + 0.01442616961986814*%i,15); 1311true; 1312 1313/******************************************************************/ 1314/* Numerical values for the bessel function Y with negative order */ 1315 1316test_bessel 1317 (bessel_y(-1,-2.0),-0.1070324315409375 + 1.1534496155137468 * %i,14); 1318true; 1319 1320test_bessel(bessel_y(-1,2.0),0.1070324315409375,15); 1321true; 1322 1323test_bessel 1324 (bessel_y(-1,-1.5),-0.4123086269739113 + 1.1158730158201993 * %i,15); 1325true; 1326 1327test_bessel(bessel_y(-1,1.5),0.4123086269739113,15); 1328true; 1329 1330test_bessel(bessel_y(-1.5, -2.0),0.4912937786871623 * %i,15); 1331true; 1332 1333test_bessel(bessel_y(-1.5, 2.0),-0.4912937786871623,15); 1334true; 1335 1336test_bessel 1337 (bessel_y(-1.8, -1.5),-0.6777414340388011 + 0.0857519944018923 * %i,14); 1338true; 1339 1340test_bessel(bessel_y(-1.8, 1.5),-0.8377344836401481,14); 1341true; 1342 1343test_bessel 1344 (bessel_y(-2,-1.5),-0.9321937597629739 + 0.4641753442884295 * %i,15); 1345true; 1346 1347test_bessel(bessel_y(-2,1.5),-0.9321937597629739,14); 1348true; 1349 1350test_bessel(bessel_y(-2.5,-1.5),0.1244463597983876 * %i,14); 1351true; 1352 1353test_bessel(bessel_y(-2.5,1.5),0.1244463597983876,15); 1354true; 1355 1356test_bessel 1357 (bessel_y(-2.3,-1.5),-0.3148570865836879 + 0.7565240896444820 * %i,14); 1358true; 1359 1360test_bessel(bessel_y(-2.3,1.5),-0.5356668704355646,14); 1361true; 1362 1363/* Numerical values for the bessel function Y with positive order */ 1364 1365test_bessel(bessel_y(1.5,1.0),-1.102495575160179,14); 1366true; 1367 1368test_bessel(bessel_y(1.5,-1.0),-1.102495575160179 * %i,14); 1369true; 1370 1371test_bessel(bessel_y(1.8,1.0),-1.382351995367631,14); 1372true; 1373 1374test_bessel 1375 (bessel_y(1.8,-1.0),-1.1183462564605324 - 0.5593113771009602 * %i,14); 1376true; 1377 1378test_bessel(bessel_y(2.0,1.0),-1.650682606816254,14); 1379true; 1380 1381test_bessel 1382 (bessel_y(2.0,-1.0),-1.650682606816254 + 0.229806969863801 * %i,14); 1383true; 1384 1385test_bessel(bessel_y(2.5,1.0),-2.876387857462161,14); 1386true; 1387 1388test_bessel(bessel_y(2.5,-1.0),2.876387857462161 * %i,14); 1389true; 1390 1391 1392/* Numerical values for the bessel function Y with complex arg 1393 and positive or negative order*/ 1394 1395test_bessel 1396 (bessel_y(0,1.0+%i),0.4454744889360325 + 0.7101585820037345 * %i,15); 1397true; 1398 1399test_bessel 1400 (bessel_y(1,1.0+%i),-0.6576945355913452 + 0.6298010039928844 * %i,15); 1401true; 1402 1403test_bessel 1404(bessel_y(-1,1.0+%i),0.6576945355913452 - 0.6298010039928844 * %i,15); 1405true; 1406 1407test_bessel 1408 (bessel_y(2,1.0+%i),-0.4733680205344934 + 0.5773369575804951 * %i,14); 1409true; 1410 1411test_bessel 1412 (bessel_y(-2,1.0+%i),-0.4733680205344934 + 0.5773369575804951 * %i,14); 1413true; 1414 1415test_bessel 1416 (bessel_y(2.3,1.0+%i),-0.2476879981252862 + 0.7595467103431256 * %i,15); 1417true; 1418 1419test_bessel 1420 (bessel_y(-2.3,1.0+%i),-0.1570442638685963 + 0.5821870838327466 * %i,14); 1421true; 1422 1423/* Numerical values for the bessel function Y with complex order */ 1424 1425test_bessel 1426 (bessel_y(%i,1.0),-0.476556612479964 - 1.505069159110387*%i,14); 1427true; 1428 1429test_bessel 1430 (bessel_y(%i,1.5),0.5161218926267688 - 1.2857405211747503*%i,14); 1431true; 1432 1433test_bessel 1434 (bessel_y(1.0+%i,-1.0),7.708594946281541 + 1.233384674244926*%i,14); 1435true; 1436 1437test_bessel 1438 (bessel_y(1.5*%i,-2.0),3.226466016458932 + 4.267260420563194*%i,13); 1439true; 1440 1441/******************************************************************/ 1442/* Numerical values for the bessel function I with negative order */ 1443 1444test_bessel(bessel_i(-1,-2.0),-1.590636854637329,15); 1445true; 1446 1447test_bessel(bessel_i(-1,2.0),1.590636854637329,15); 1448true; 1449 1450test_bessel(bessel_i(-1,-1.5),-0.9816664285779076,15); 1451true; 1452 1453test_bessel(bessel_i(-1,1.5),0.9816664285779076,15); 1454true; 1455 1456test_bessel(bessel_i(-1.5, -2.0),0.9849410530002364 * %i,14); 1457true; 1458 1459test_bessel(bessel_i(-1.5, 2.0),0.9849410530002364,14); 1460true; 1461 1462test_bessel 1463 (bessel_i(-1.8, -1.5),0.2026903980307014 + 0.1472631941876387 * %i,15); 1464true; 1465 1466test_bessel(bessel_i(-1.8, 1.5),0.2505391103524365,15); 1467true; 1468 1469test_bessel(bessel_i(-2,-1.5),0.3378346183356807,15); 1470true; 1471 1472test_bessel(bessel_i(-2,1.5),0.3378346183356807,15); 1473true; 1474 1475test_bessel(bessel_i(-2.5,-1.5),-0.8015666610717216*%i,14); 1476true; 1477 1478test_bessel(bessel_i(-2.5,1.5),0.8015666610717216,14); 1479true; 1480 1481test_bessel 1482 (bessel_i(-2.3,-1.5),0.3733945265830230 - 0.5139334755917659 * %i,15); 1483true; 1484 1485test_bessel(bessel_i(-2.3,1.5),0.6352567117441516,15); 1486true; 1487 1488/* Numerical values for the bessel function I with positive order */ 1489 1490test_bessel(bessel_i(1.5,1.0),0.2935253263474798,15); 1491true; 1492 1493test_bessel(bessel_i(1.5,-1.0),-0.2935253263474798*%i,13); 1494true; 1495 1496test_bessel(bessel_i(1.8,1.0),0.1871011888310777,15); 1497true; 1498 1499test_bessel 1500 (bessel_i(1.8,-1.0),0.1513680414320980 - 0.1099753194812967 * %i,14); 1501true; 1502 1503test_bessel(bessel_i(2.0,1.0),0.1357476697670383,15); 1504true; 1505 1506test_bessel(bessel_i(2.0,-1.0),0.1357476697670383,15); 1507true; 1508 1509test_bessel(bessel_i(2.5,1.0),0.05709890920304825,15); 1510true; 1511 1512test_bessel(bessel_i(2.5,-1.0),0.05709890920304825 * %i,15); 1513true; 1514 1515 1516/* Numerical values for the bessel function I with complex arg 1517 and positive or negative order*/ 1518 1519test_bessel 1520 (bessel_i(0,1.0+%i),0.9376084768060293 + 0.4965299476091221 * %i,15); 1521true; 1522 1523test_bessel 1524 (bessel_i(1,1.0+%i),0.3650280288270878 + 0.6141603349229036 * %i,15); 1525true; 1526 1527test_bessel 1528 (bessel_i(-1,1.0+%i),0.3650280288270878 + 0.6141603349229036 * %i,15); 1529true; 1530 1531test_bessel 1532 (bessel_i(2,1.0+%i),-0.0415798869439621 + 0.2473976415133063 * %i,15); 1533true; 1534 1535test_bessel 1536 (bessel_i(-2,1.0+%i),-0.0415798869439621 + 0.2473976415133063 * %i,15); 1537true; 1538 1539test_bessel 1540 (bessel_i(2.3,1.0+%i),-0.0635524383467825 + 0.1559221140952053 * %i,15); 1541true; 1542 1543test_bessel 1544 (bessel_i(-2.3,1.0+%i),-0.4053256245784623 - 0.3724481230406298 * %i,14); 1545true; 1546 1547/* Numerical values for the bessel function I with complex order */ 1548 1549test_bessel 1550 (bessel_i(%i,1.0),1.900799675819425 - 1.063960013554441*%i,14); 1551true; 1552 1553test_bessel 1554 (bessel_i(%i,1.5),2.495473638417463 - 0.601347501920535*%i,14); 1555true; 1556 1557test_bessel 1558 (bessel_i(1.0+%i,-1.0),-0.01096313515009349 + 0.03043920114776303*%i,15); 1559true; 1560 1561test_bessel 1562 (bessel_i(1.5*%i,-2.0),0.04238259669782487 - 0.01125055344512197*%i,15); 1563true; 1564 1565/******************************************************************/ 1566/* Numerical values for the bessel function K with negative order */ 1567 1568test_bessel 1569 (bessel_k(-1,-2.0),-0.139865881816522-4.997133057057809*%i,14); 1570true; 1571 1572test_bessel(bessel_k(-1,2.0),0.139865881816522,14); 1573true; 1574 1575test_bessel 1576 (bessel_k(-1,-1.5),-0.277387800456844 - 3.083996040296084*%i,14); 1577true; 1578 1579test_bessel(bessel_k(-1,1.5),0.2773878004568438,15); 1580true; 1581 1582test_bessel(bessel_k(-1.5, -2.0),-3.2741902342766302*%i,13); 1583true; 1584 1585test_bessel(bessel_k(-1.5, 2.0),0.1799066579520922,15); 1586true; 1587 1588test_bessel 1589 (bessel_k(-1.8, -1.5),0.3929372194683435 - 1.0725774293000646*%i,15); 1590true; 1591 1592test_bessel(bessel_k(-1.8, 1.5),0.4856971141526263,15); 1593true; 1594 1595test_bessel 1596 (bessel_k(-2,-1.5),0.5836559632566508 - 1.0613387550916862*%i,14); 1597true; 1598 1599test_bessel(bessel_k(-2,1.5),0.5836559632566508,15); 1600true; 1601 1602test_bessel 1603 (bessel_k(-2.3,-1.5),0.465598659425186 - 1.354876241730594*%i,14); 1604true; 1605 1606test_bessel(bessel_k(-2.3,1.5),0.7921237520153218,15); 1607true; 1608 1609/* Numerical values for the bessel function K with positive order */ 1610 1611test_bessel(bessel_k(1.5,1.0),0.9221370088957891,14); 1612true; 1613 1614/* kein Wert von function-site !? Ist das eine Nullstelle??? */ 1615test_bessel(bessel_k(1.5,-1.0),0,13); 1616true; 1617 1618test_bessel(bessel_k(1.8,1.0),1.275527037541854,15); 1619true; 1620 1621test_bessel 1622 (bessel_k(1.8,-1.0),1.0319230501560911 + 0.1619402612577788*%i,14); 1623true; 1624 1625test_bessel(bessel_k(2.0,1.0),1.624838898635177,14); 1626true; 1627 1628test_bessel(bessel_k(2.0,-1.0),1.624838898635177 - 0.426463882082061*%i,14); 1629true; 1630 1631test_bessel(bessel_k(2.5,1.0),3.227479531135262,14); 1632true; 1633 1634test_bessel(bessel_k(2.5,-1.0),-3.406861044815549*%i,14); 1635true; 1636 1637/* Numerical values for the bessel function K with complex arg 1638 and positive or negative order*/ 1639 1640test_bessel 1641 (bessel_k(0,1.0+%i),0.0801977269465178 - 0.3572774592853303*%i,15); 1642true; 1643 1644test_bessel 1645 (bessel_k(1,1.0+%i),0.0245683055237403 - 0.4597194738011894 * %i,15); 1646true; 1647 1648test_bessel 1649 (bessel_k(-1,1.0+%i),0.0245683055237403 - 0.4597194738011894 * %i,15); 1650true; 1651 1652test_bessel 1653 (bessel_k(2,1.0+%i),-0.3549534413309312 - 0.8415652386102600 * %i,13); 1654true; 1655 1656test_bessel 1657 (bessel_k(-2,1.0+%i),-0.3549534413309312 - 0.8415652386102600 * %i,13); 1658true; 1659 1660test_bessel 1661 (bessel_k(2.3,1.0+%i),-0.6635905911278042 - 1.0258894849569300 * %i,15); 1662true; 1663 1664test_bessel 1665 (bessel_k(-2.3,1.0+%i),-0.6635905911278042 - 1.0258894849569300 * %i,13); 1666true; 1667 1668/* Numerical values for the bessel function K with complex order */ 1669 1670test_bessel 1671 (bessel_k(%i,1.0),0.2894280370259921,14); 1672true; 1673 1674test_bessel 1675 (bessel_k(%i,1.5),0.1635839926633096,14); 1676true; 1677 1678test_bessel 1679 (bessel_k(1.0+%i,-1.0),-9.744252766030894 - 7.494570149760043*%i,14); 1680true; 1681 1682test_bessel 1683 (bessel_k(1.5*%i,-2.0),3.93512366711118 - 14.82183468316553*%i,13); 1684true; 1685 1686kill(closeto, test_bessel); 1687done; 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) := 1704block([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); 17340; 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 1744w_jj(n,z) := bessel_j(n+1,z)*bessel_j(-n,z) + bessel_j(n,z)*bessel_j(-n-1,z); 1745w_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 1749test_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 1755test_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 1768w_jy(n,z) := bessel_j(n+1,z)*bessel_y(n,z) - bessel_j(n,z)*bessel_y(n+1,z); 1769w_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 1773test_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 1779test_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 1790w_ii(n,z) := bessel_i(n,z)*bessel_i(-n-1,z) - bessel_i(n+1,z)*bessel_i(-n,z); 1791w_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 1795test_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 1800test_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 1813w_ik(n,z) := bessel_i(n,z)*bessel_k(n+1,z) + bessel_i(n+1,z)*bessel_k(n,z); 1814w_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 1818test_wronskian(w_ik, lambda([n,z], 1/z), 1e-10, 5); 1819[]; 1820 1821/* Calculation of w_ik for complex argument */ 1822 1823test_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 1834w_h1h2(v,z) := hankel_1(v+1,z)*hankel_2(v,z) - hankel_1(v,z)*hankel_2(v+1,z); 1835w_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 1839test_wronskian(w_h1h2, lambda([v,z], -4*%i/(%pi*z)), 1e-14, 5); 1840[]; 1841 1842/* Calculation of w_h1h2 for complex argument */ 1843 1844test_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 1851test_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 1862integrate(bessel_j(0,x),x); 1863x*(bessel_j(0,x)*(2-%pi*struve_h(1,x))+%pi*bessel_j(1,x)*struve_h(0,x))/2; 1864 1865integrate(bessel_j(1,x),x); 1866-bessel_j(0,x); 1867 1868integrate(bessel_j(2,x),x); 1869hypergeometric([3/2],[5/2,3],-x^2/4)*x^3/24; 1870 1871integrate(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/ */ 1875integrate(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 1878integrate(bessel_y(1,x),x); 1879-bessel_y(0,x); 1880 1881integrate(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 1884integrate(bessel_y(3,x),x); 1885-2*bessel_y(2,x)-bessel_y(0,x); 1886 1887integrate(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 1891integrate(bessel_y(11,x),x); 1892-2*'sum(bessel_y(2*i,x),i,1,5)-bessel_y(0,x); 1893 1894integrate(bessel_i(0,x),x); 1895x*(bessel_i(0,x)*(%pi*struve_l(1,x)+2)-%pi*bessel_i(1,x)*struve_l(0,x))/2; 1896 1897integrate(bessel_i(1,x),x); 1898bessel_i(0,x); 1899 1900integrate(bessel_i(2,x),x); 1901hypergeometric([3/2],[5/2,3],x^2/4)*x^3/24; 1902 1903integrate(bessel_i(1/2,x),x); 19042^(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/ */ 1907integrate(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 1910integrate(bessel_k(1,x),x); 1911-bessel_k(0,x); 1912 1913integrate(bessel_k(7,x),x); 19142*(-bessel_k(6,x)+bessel_k(4,x)-bessel_k(2,x))+bessel_k(0,x); 1915 1916integrate(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); 1930done; 1931 1932(declare(n,integer),assume(x>0),done); 1933done; 1934 1935f(n,x); 1936[bessel_j(n,x),0]; 1937f(n,-x); 1938[bessel_j(n,-x),0]; 1939 1940f(n+1/2,x); 1941[bessel_j(n+1/2,x),0]; 1942f(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); 1946done; 1947 1948f(n_even,%i); 1949[bessel_j(n_even,%i),0]; 1950f(n_odd,%i); 1951[0,-%i*bessel_j(n_odd,%i)]; 1952 1953f(n_even,x*%i); 1954[bessel_j(n_even,x*%i),0]; 1955f(n_odd,x*%i); 1956[0,-%i*bessel_j(n_odd,x*%i)]; 1957 1958f(n_even,(x+1)^2*%i); 1959[bessel_j(n_even,(x+1)^2*%i),0]; 1960f(n_odd,(x+1)^2*%i); 1961[0,-%i*bessel_j(n_odd,(x+1)^2*%i)]; 1962 1963(declare(j,imaginary),done); 1964done; 1965 1966f(n_even,(x+1)^2*j); 1967[bessel_j(n_even,(x+1)^2*j),0]; 1968f(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); 1974done; 1975 1976f(n,x); 1977[bessel_i(n,x),0]; 1978f(n,-x); 1979[bessel_i(n,-x),0]; 1980 1981f(n+1/2,x); 1982[bessel_i(n+1/2,x),0]; 1983f(n+1/2,-x); 1984[0,-%i*bessel_i(n+1/2,-x)]; 1985 1986f(n_even,%i); 1987[bessel_i(n_even,%i),0]; 1988f(n_odd,%i); 1989[0,-%i*bessel_i(n_odd,%i)]; 1990 1991f(n_even,x*%i); 1992[bessel_i(n_even,x*%i),0]; 1993f(n_odd,x*%i); 1994[0,-%i*bessel_i(n_odd,x*%i)]; 1995 1996f(n_even,(x+1)^2*%i); 1997[bessel_i(n_even,(x+1)^2*%i),0]; 1998f(n_odd,(x+1)^2*%i); 1999[0,-%i*bessel_i(n_odd,(x+1)^2*%i)]; 2000 2001f(n_even,(x+1)^2*j); 2002[bessel_i(n_even,(x+1)^2*j),0]; 2003f(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); 2009done; 2010 2011f(n,x); 2012[bessel_k(n,x),0]; 2013f(n,-x); 2014[realpart(bessel_k(n,-x)),imagpart(bessel_k(n,-x))]; 2015 2016f(n+1/2,x); 2017[bessel_k(n+1/2,x),0]; 2018f(n+1/2,-x); 2019[0,-%i*bessel_k(n+1/2,-x)]; 2020 2021f(n_even,%i); 2022[realpart(bessel_k(n_even,%i)),imagpart(bessel_k(n_even,%i))]; 2023f(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); 2030done; 2031 2032f(n,x); 2033[bessel_y(n,x),0]; 2034f(n,-x); 2035[realpart(bessel_y(n,-x)),imagpart(bessel_y(n,-x))]; 2036 2037f(n+1/2,x); 2038[bessel_y(n+1/2,x),0]; 2039f(n+1/2,-x); 2040[0,-%i*bessel_y(n+1/2,-x)]; 2041 2042f(n_even,%i); 2043[realpart(bessel_y(n_even,%i)),imagpart(bessel_y(n_even,%i))]; 2044f(n_odd,%i); 2045[realpart(bessel_y(n_odd,%i)),imagpart(bessel_y(n_odd,%i))]; 2046 2047/* %iargs is false, so transformation disabled */ 2048bessel_j(v, %i*x); 2049bessel_j(v, %i*x); 2050 2051bessel_i(v, %i*x); 2052bessel_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 2058bessel_i(v, %i*x); 2059(%i*x)^v/(x^v)*bessel_j(v, x); 2060 2061imagpart(bessel_j(2,%i*3.0)); 20620; 2063realpart(bessel_j(3,%i*3.0)); 20640; 2065 2066imagpart(bessel_i(2,%i*3.0)); 20670; 2068realpart(bessel_i(3,%i*3.0)); 20690; 2070 2071/*******************************************************/ 2072/* Check the handling of conjugate on Bessel functions */ 2073/*******************************************************/ 2074 2075declare([w,z],complex,n,integer); 2076done; 2077assume(nonneg>=0,notequal(nonzero,0)); 2078[nonneg>=0,notequal(nonzero,0)]; 2079 2080conjugate(bessel_j(w,z)); 2081'(conjugate(bessel_j(w,z))); 2082 2083conjugate(bessel_y(w,z)); 2084'(conjugate(bessel_y(w,z))); 2085 2086conjugate(bessel_i(w,z)); 2087'(conjugate(bessel_i(w,z))); 2088 2089conjugate(bessel_k(w,z)); 2090'(conjugate(bessel_k(w,z))); 2091 2092conjugate(hankel_1(w,z)); 2093'(conjugate(hankel_1(w,z))); 2094 2095conjugate(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 */ 2101conjugate(bessel_j(z,x+%i*nonzero)); 2102'(bessel_j(conjugate(z),x-%i*nonzero)); 2103 2104conjugate(bessel_j(z,nonneg)); 2105'(bessel_j(conjugate(z),nonneg)); 2106 2107conjugate(bessel_y(z,x+%i*nonzero)); 2108'(bessel_y(conjugate(z),x-%i*nonzero)); 2109 2110conjugate(bessel_y(z,nonneg)); 2111'(bessel_y(conjugate(z),nonneg)); 2112 2113conjugate(bessel_i(z,x+%i*nonzero)); 2114'(bessel_i(conjugate(z),x-%i*nonzero)); 2115 2116conjugate(bessel_i(z,nonneg)); 2117'(bessel_i(conjugate(z),nonneg)); 2118 2119conjugate(bessel_k(z,x+%i*nonzero)); 2120'(bessel_k(conjugate(z),x-%i*nonzero)); 2121 2122conjugate(bessel_k(z,nonneg)); 2123'(bessel_k(conjugate(z),nonneg)); 2124 2125conjugate(hankel_1(z,x+%i*nonzero)); 2126'(hankel_2(conjugate(z),x-%i*nonzero)); 2127 2128conjugate(hankel_1(z,nonneg)); 2129'(hankel_2(conjugate(z),nonneg)); 2130 2131conjugate(hankel_2(z,x+%i*nonzero)); 2132'(hankel_1(conjugate(z),x-%i*nonzero)); 2133 2134conjugate(hankel_2(z,nonneg)); 2135'(hankel_1(conjugate(z),nonneg)); 2136 2137/* Bessel functions J and I of integral order commute with conjugate */ 2138conjugate(bessel_j(n,z)); 2139'(bessel_j(n,conjugate(z))); 2140 2141conjugate(bessel_i(n,z)); 2142'(bessel_i(n,conjugate(z))); 2143 2144remove([w,z],complex,n,integer); 2145done; 2146forget(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); 21520; 2153 2154a[s]; 21555; 2156 2157a[s], nouns; 21585; 2159 2160a[t]; 2161false; 2162 2163a[t], nouns; 2164false; 2165 2166