tests/rtest_hypgeo.mac

   1 /******************************************************************************
   2
   3    Tests from A&S
   4
   5 ******************************************************************************/
   6
   7 /* Make assumptions about parameters which we will use for the tests */
   8
   9 (kill(all),declare(n,integer, k,real),assume(s>0,n>0, k>0),'done);
  10 done;
  11
  12 /******************************************************************************
  13    First part: A&S 29.3.1 to 29.3.33
  14
  15    Problem 2: A&S: 29.3.1 OK
  16 ******************************************************************************/
  17
  18 specint(%e^(-s*t),t);
  19 1/s;
  20
  21 /*****************************************************************************
  22    Problem 3: A&S: 29.3.2 OK
  23 ******************************************************************************/
  24
  25 specint(%e^(-s*t)*t,t);
  26 1/s^2;
  27
  28 /*****************************************************************************
  29    Problem 4: A&S: 29.3.3 OK
  30 ******************************************************************************/
  31
  32 specint(%e^(-s*t)*t^(n-1)/gamma(n),t);  /* use gamma(n) instead of (n-1)! */
  33 1/s^n;
  34
  35 /*****************************************************************************
  36    Problem 5: A&S: 29.3.4 OK
  37 ******************************************************************************/
  38
  39 specint(%e^(-s*t)/sqrt(%pi*t),t);
  40 1/sqrt(s);
  41
  42 /*****************************************************************************
  43    Problem 6: A&S: 29.3.5 OK
  44 ******************************************************************************/
  45
  46 specint(%e^(-s*t)*2*sqrt(t/%pi),t);
  47 1/s^(3/2);
  48
  49 /*****************************************************************************
  50    Problem 7-10: A&S: 29.3.6 OK
  51 ******************************************************************************/
  52
  53 /* Define the function to be integrated */
  54
  55 f(n,t):= 2^n*t^(n-1/2)/(product(2*k-1,k,1,n)*sqrt(%pi));
  56 f(n,t):= 2^n*t^(n-1/2)/(product(2*k-1,k,1,n)*sqrt(%pi));
  57
  58 /* Test for 3 values of n.
  59    Can Maxima find the general solution, too? */
  60
  61 specint(%e^(-s*t)*f(3,t),t);
  62 s^(-3-1/2);
  63
  64 specint(%e^(-s*t)*f(7,t),t);
  65 s^(-7-1/2);
  66
  67 specint(%e^(-s*t)*f(100,t),t);
  68 s^(-100-1/2);
  69
  70 /*****************************************************************************
  71    Problem 11: A&S: 29.3.7 OK
  72 ******************************************************************************/
  73
  74 specint(%e^(-s*t)*t^(k-1),t);
  75 gamma(k)/s^k;
  76
  77 /*****************************************************************************
  78    Problem 12: A&S: 29.3.8 OK
  79 ******************************************************************************/
  80
  81 (assume(s+a>0), specint(%e^(-s*t)*%e^(-a*t),t));
  82 1/(s+a);
  83
  84 /*****************************************************************************
  85    Problem 13: A&S: 29.3.9 OK
  86 ******************************************************************************/
  87
  88 specint(%e^(-s*t)*t*%e^(-a*t),t);
  89 1/(s+a)^2;
  90
  91 /*****************************************************************************
  92    Problem 14: A&S: 29.3.10 OK
  93 ******************************************************************************/
  94
  95 specint(%e^(-s*t)*t^(n-1)*%e^(-a*t)/gamma(n),t);
  96 1/(s+a)^n;
  97
  98 /*****************************************************************************
  99    Problem 15: A&S: 29.3.11 OK
 100 ******************************************************************************/
 101
 102 specint(%e^(-s*t)*t^(k-1)*%e^(-a*t),t);
 103 gamma(k)/(s+a)^k;
 104
 105 /*****************************************************************************
 106    Problem 16: A&S: 29.3.12 OK New
 107    No solution with the original code because of the denominator.
 108 ******************************************************************************/
 109
 110 (assume(s+b>0), factor(specint(%e^(-s*t)*(%e^(-a*t)-%e^(-b*t))/(b-a),t)));
 111 1/((s+a)*(s+b));
 112
 113 /*****************************************************************************
 114    Problem 17: A&S: 29.3.13 OK New
 115    No solution with the original code because of the denominator.
 116 ******************************************************************************/
 117
 118 factor(specint(%e^(-s*t)*(a*%e^(-a*t)-b*%e^(-b*t))/(a-b),t));
 119 s/((s+a)*(s+b));
 120
 121 /*****************************************************************************
 122    Problem 18: A&S: 29.3.14 OK New
 123    No solution with the original code because of the denominator.
 124 ******************************************************************************/
 125
 126 (assume(s+c>0),
 127   factor(specint(%e^(-s*t)*((b-c)*%e^(-a*t)+(c-a)*%e^(-b*t)+(a-b)*%e^(-c*t))/((b-a)*(b-c)*(c-a)),t)));
 128 1/((s+a)*(s+b)*(s+c));
 129
 130 /*****************************************************************************
 131    Problem 19: A&S: 29.3.15 OK
 132 ******************************************************************************/
 133
 134 ratsimp(specint(%e^(-s*t)*sin(a*t)/a,t));
 135 1/(s^2+a^2);
 136
 137 /*****************************************************************************
 138    Problem 20: A&S: 29.3.16 OK
 139 ******************************************************************************/
 140
 141 ratsimp(specint(%e^(-s*t)*cos(a*t),t));
 142 s/(s^2+a^2);
 143
 144 /*****************************************************************************
 145    Problem 21: A&S: 29.3.17 OK
 146 ******************************************************************************/
 147
 148 (assume(s-a>0),ratsimp(specint(%e^(-s*t)*sinh(a*t)/a,t)));
 149 1/(s^2-a^2);
 150
 151 /*****************************************************************************
 152    Problem 22: A&S: 29.3.18 OK
 153 ******************************************************************************/
 154
 155 ratsimp(specint(%e^(-s*t)*cosh(a*t),t));
 156 s/(s^2-a^2);
 157
 158 /*****************************************************************************
 159    Problem 23: A&S: 29.3.19 OK
 160 ******************************************************************************/
 161
 162 factor(ratsimp(specint(%e^(-s*t)*(1-cos(a*t))/a^2,t)));
 163 1/(s*(s^2+a^2));
 164
 165 /*****************************************************************************
 166    Problem 24: A&S: 29.3.20 OK
 167 ******************************************************************************/
 168
 169 factor(ratsimp(specint(%e^(-s*t)*(a*t-sin(a*t))/a^3,t)));
 170 1/(s^2*(s^2+a^2));
 171
 172 /*****************************************************************************
 173    Problem 25: A&S: 29.3.21 OK
 174 ******************************************************************************/
 175
 176 factor(ratsimp(specint(%e^(-s*t)*(sin(a*t)-a*t*cos(a*t))/(2*a^3),t)));
 177 1/(s^2+a^2)^2;
 178
 179 /*****************************************************************************
 180    Problem 26: A&S: 29.3.22 OK
 181 ******************************************************************************/
 182
 183 factor(ratsimp(specint(%e^(-s*t)*t*sin(a*t)/(2*a),t)));
 184 s/(s^2+a^2)^2;
 185
 186 /*****************************************************************************
 187    Problem 27: A&S: 29.3.23 OK
 188 ******************************************************************************/
 189
 190 factor(ratsimp(specint(%e^(-s*t)*(sin(a*t)+a*t*cos(a*t))/(2*a),t)));
 191 s^2/(s^2+a^2)^2;
 192
 193 /*****************************************************************************
 194    Problem 28: A&S: 29.3.24 OK
 195 ******************************************************************************/
 196
 197 factor(ratsimp(specint(%e^(-s*t)*t*cos(a*t),t)));
 198 (s-a)*(s+a)/(s^2+a^2)^2;
 199
 200 /*****************************************************************************
 201    Problem 29: A&S: 29.3.25 OK New
 202    After implementing a more general code to extract a constant denominator
 203    we get the correct result.
 204 ******************************************************************************/
 205
 206 factor(ratsimp(specint(%e^(-s*t)*(cos(a*t)-cos(b*t))/(b^2-a^2),t)));
 207 s/((s^2+a^2)*(s^2+b^2));
 208
 209 /*****************************************************************************
 210    Problem 30: A&S: 29.3.26 OK
 211 ******************************************************************************/
 212
 213 factor(ratsimp(specint(%e^(-s*t)*(%e^(-a*t)*sin(b*t))/b,t)));
 214 1/(s^2+2*a*s+b^2+a^2);
 215
 216 /*****************************************************************************
 217    Problem 31: A&S: 29.3.27 OK
 218 ******************************************************************************/
 219
 220 factor(ratsimp(specint(%e^(-s*t)*(%e^(-a*t)*cos(b*t)),t)));
 221 (s+a)/(s^2+2*a*s+b^2+a^2);
 222
 223 /*****************************************************************************
 224    Problem 32: A&S: 29.3.28 OK
 225 ******************************************************************************/
 226
 227 (assume(2*s-a>0),
 228 factor(ratsimp(specint(%e^(-s*t)*(%e^(-a*t)-%e^(a*t/2)*(cos(a*t*sqrt(3)/2)
 229                                            -sqrt(3)*sin(a*t*sqrt(3)/2))),t))));
 230 3*a^2/((s+a)*(s^2-a*s+a^2));  /* that's 3*a^2/(s^3+a^3) */
 231
 232 /*****************************************************************************
 233    Problem 33: A&S: 29.3.29 OK
 234 ******************************************************************************/
 235
 236 factor(ratsimp(specint(%e^(-s*t)*(sin(a*t)*cosh(a*t)-cos(a*t)*sinh(a*t)),t)));
 237 4*a^3/((s^2-2*a*s+2*a^2)*(s^2+2*a*s+2*a^2));
 238
 239 /*****************************************************************************
 240    Problem 34: A&S: 29.3.30 OK
 241 ******************************************************************************/
 242
 243 factor(ratsimp(specint(%e^(-s*t)*sin(a*t)*sinh(a*t)/(2*a^2),t)));
 244 s/((s^2-2*a*s+2*a^2)*(s^2+2*a*s+2*a^2));
 245
 246 /*****************************************************************************
 247    Problem 35: A&S: 29.3.31 OK
 248 ******************************************************************************/
 249
 250 factor(ratsimp(specint(%e^(-s*t)*(sinh(a*t)-sin(a*t))/(2*a^3),t)));
 251 1/((s-a)*(s+a)*(s^2+a^2));
 252
 253 /*****************************************************************************
 254    Problem 36: A&S: 29.3.32 OK
 255 ******************************************************************************/
 256
 257 factor(ratsimp(specint(%e^(-s*t)*(cosh(a*t)-cos(a*t))/(2*a^2),t)));
 258 s/((s-a)*(s+a)*(s^2+a^2));
 259
 260 /*****************************************************************************
 261    Problem 37: A&S: 29.3.33 OK
 262 ******************************************************************************/
 263
 264 factor(ratsimp(specint(%e^(-s*t)*((1+a^2*t^2)*sin(a*t)-a*t*cos(a*t)),t)));
 265 8*a^3*s^2/(s^2+a^2)^3;
 266
 267 /******************************************************************************
 268    Problem 38: A&S: 29.3.34
 269
 270    Laguerre Polynom
 271
 272 ******************************************************************************/
 273
 274 /* We test the integration of the Laguerre function for three different
 275    values for the order of the polynom.
 276 */
 277
 278 factor(specint(%e^(-s*t)*laguerre(1,t),t));
 279 (s-1)/s^2;
 280
 281 factor(specint(%e^(-s*t)*laguerre(2,t),t));
 282 (s-1)^2/s^3;
 283
 284 factor(specint(%e^(-s*t)*laguerre(10,t),t));
 285 (s-1)^10/s^(10+1);
 286
 287 /*****************************************************************************
 288    The code knows the general Laguerre function %l[n,a](x).
 289    There should be the relationship: laguerre(n,t) -> %l[n,0](t).
 290
 291    But it doesn't work. We got a lot of errors when we try to use %l[n,a](t)
 292    in a Laplace transformation.
 293
 294    First:  The function ltw seems not to be correct. There are differences
 295            beetween the code and the formula A&S 13.6.9.
 296            The code uses the Lisp-Function FACTORIAL. That's wrong for a
 297            symbol. We get a Lisp error. We have to use a general GAMMA
 298            expression.
 299    Second: In the transformation to a Hypergeomtric function there are
 300            similar problems. We get the expression gamma(1-n) which is not
 301            defined for all positive integers n. The code stops with the
 302            message gamma( .. ) is undefined.
 303            The complete transformation formula has to be verified.
 304    Last:   After some corrections to the routine ltw the result is OK.
 305 */
 306
 307 /* Problem 41: OK New
 308  * After correction of some bugs we get the correct result.
 309  * First a test with a special value for the order of the Laguerre function.
 310  */
 311
 312 factor(specint(%e^(-s*t)*'laguerre(10,t),t));
 313 (s-1)^10/s^11;
 314
 315 factor(specint(%e^(-s*t)*'laguerre(n,t),t));
 316 (s-1)^n/s^(n+1);
 317
 318 /* Problem 43: OK We get a result. New
 319    Because n is declared to be an integer, the implemented Laguerre function
 320    simplifies to a sum. The summand can be integrated.
 321    This summand is equivalent to the general solution above.
 322 */
 323
 324 (assume(i1>0), niceindices(specint(%e^(-s*t)*laguerre(n,t),t)));
 325 'sum(gamma(i+1)*pochhammer(-n,i)*s^(-i-1)/i!^2,i,0,n);
 326
 327 /*****************************************************************************
 328    Problem 44: A&S: 29.3.35 OK
 329 ******************************************************************************/
 330
 331 factor(specint(%e^(-s*t)*(1/(sqrt(%pi*t))*%e^(-a*t)*(1-2*a*t)),t));
 332 s/(s+a)^(3/2);
 333
 334 /*****************************************************************************
 335    Problem 45: A&S: 29.3.36 OK New
 336    sqrt(s+a)-sqrt(s+b);
 337 ******************************************************************************/
 338
 339 specint(%e^(-s*t)*(1/(2*sqrt(%pi*t^3)))*(%e^(-b*t)-%e^(-a*t)),t);
 340 -sqrt(s+b)+sqrt(s+a);
 341
 342 /******************************************************************************
 343  *
 344  * Error functions erf und erfc
 345  *
 346  ******************************************************************************/
 347
 348 assume(s-a^2>0,a>0);
 349 [s>a^2,a>0];
 350
 351 /*****************************************************************************
 352    Problem 47: A&S: 29.3.37 OK
 353 ******************************************************************************/
 354
 355 factor(ratsimp(specint(%e^(-s*t)*(1/sqrt(%pi*t)-a*%e^(a^2*t)
 356                                                 *erfc(a*sqrt(t))),t)));
 357 (sqrt(s)-a)/(s-a^2);
 358
 359 /*****************************************************************************
 360    Problem 48: A&S: 29.3.38 OK
 361 ******************************************************************************/
 362
 363 factor(ratsimp(specint(%e^(-s*t)*(1/sqrt(%pi*t)+a*%e^(a^2*t)
 364                                                 *erf(a*sqrt(t))),t)));
 365 sqrt(s)/(s-a^2);
 366
 367 /*****************************************************************************
 368    Problem 49: A&S: 29.3.39 OK
 369    In this Laplace transform the integral represents the erfi function.
 370    Maxima does the integration correct.
 371 ******************************************************************************/
 372
 373 (assume(t>0),
 374  factor(ratsimp(specint(%e^(-s*t)*(1/sqrt(%pi*t)-2*a/sqrt(%pi)*%e^(-a^2*t)
 375                             *integrate(%e^(lambda^2),lambda,0,a*sqrt(t))),t))));
 376 sqrt(s)/(s+a^2);
 377
 378 /*****************************************************************************
 379    Problem 50: A&S: 29.3.40 OK
 380 ******************************************************************************/
 381
 382 factor(ratsimp(specint(%e^(-s*t)*1/a*%e^(a^2*t)*erf(a*sqrt(t)),t)));
 383 1/(sqrt(s)*(s-a^2));
 384
 385 /*****************************************************************************
 386    Problem 51: A&S: 29.3.41 OK
 387 ******************************************************************************/
 388
 389 factor(ratsimp(specint(%e^(-s*t)*(2/(a*sqrt(%pi))*%e^(-a^2*t)
 390                              *integrate(%e^(lambda^2),lambda,0,a*sqrt(t))),t)));
 391 1/(sqrt(s)*(s+a^2));
 392
 393 /*****************************************************************************
 394    Problem 52: A&S: 29.3.42 OK
 395 ******************************************************************************/
 396
 397 (assume(s-b^2>0,b>0),
 398  factor(ratsimp(specint(%e^(-s*t)*(%e^(a^2*t)*(b-a*erf(a*sqrt(t)))
 399                                            -b*%e^(b^2*t)*erfc(b*sqrt(t))),t))));
 400 (b-a)*(b+a)*(sqrt(s)-b)/((s-a^2)*(s-b^2));
 401
 402 /*****************************************************************************
 403    Problem 53: A&S: 29.3.43 OK
 404 ******************************************************************************/
 405
 406 factor(ratsimp(specint(%e^(-s*t)*%e^(a^2*t)*erfc(a*sqrt(t)),t)));
 407 (sqrt(s)-a)/(sqrt(s)*(s-a^2));
 408
 409 /*****************************************************************************
 410    Problem 54: A&S: 29.3.44 OK
 411 ******************************************************************************/
 412
 413 factor(ratsimp(specint(%e^(-s*t)*(1/sqrt(b-a))*%e^(-a*t)
 414                                               *erf(sqrt(b-a)*sqrt(t)),t)));
 415 1/((s+a)*sqrt(s+b));
 416
 417 /*****************************************************************************
 418    Problem 55: A&S: 29.3.45 OK
 419 ******************************************************************************/
 420
 421 factor(ratsimp(specint(%e^(-s*t)*(%e^(a^2*t)*(b/a*erf(a*sqrt(t))-1)
 422                                               +%e^(b^2*t)*erfc(b*sqrt(t))),t)));
 423 ((b-a)*(b+a)*(s-b*sqrt(s)))/(s*(s-a^2)*(s-b^2));
 424
 425 /******************************************************************************
 426   Problem 56: A&S 29.3.46 and A&S 29.3.47
 427
 428   Hermite
 429 ******************************************************************************/
 430
 431 /* We test the integration of the Hermite Polynom for three different
 432    values for the order of the polynom.
 433 */
 434
 435 block (h1(n,x):= n!/((2*n)!*sqrt(%pi*t))*hermite(2*n,sqrt(t)),
 436    h1(n,x):= n!/((2*n)!*sqrt(%pi*t))*hermite(2*n,sqrt(t)),
 437    done
 438 );
 439 done;
 440
 441 /* Problem 57: OK */
 442
 443 factor(specint(%e^(-s*t)*h1(1,t),t));
 444 -((s-1)/(s^(3/2)));
 445
 446 /* Problem 58: OK */
 447
 448 factor(specint(%e^(-s*t)*h1(2,t),t));
 449 ((s-1)^2/(s^(5/2)));
 450
 451 /* Problem 59: OK */
 452
 453 factor(specint(%e^(-s*t)*h1(10,t),t));
 454 ((s-1)^10/(s^(10+1/2)));
 455
 456 /* AS 29.3.47 */
 457
 458 h2(n,x):= n!/((2*n+1)!*sqrt(%pi))*hermite(2*n+1,sqrt(t));
 459 h2(n,x):= n!/((2*n+1)!*sqrt(%pi))*hermite(2*n+1,sqrt(t));
 460
 461 /* Problem 61: OK */
 462
 463 factor(specint(%e^(-s*t)*h2(1,t),t));
 464 -((s-1)/(s^(5/2)));
 465
 466 /* Problem 62: OK */
 467
 468 factor(specint(%e^(-s*t)*h2(2,t),t));
 469 ((s-1)^2/(s^(7/2)));
 470
 471 /* Problem 63: OK */
 472
 473 factor(specint(%e^(-s*t)*h2(10,t),t));
 474 ((s-1)^10/(s^(10+3/2)));
 475
 476 /* Laplace transform of hermite for even and odd order:
 477  *
 478  * Problem 64: AS 29.3.46
 479  */
 480
 481 laplace(1/sqrt(t)*'hermite(2*n,sqrt(t)),t,s);
 482 %pi*2^(2*n)*(1-1/s)^n/(gamma((2*n-1)/-2)*sqrt(s));
 483
 484 /* Problem 65: AS 29.3.47 */
 485
 486 laplace('hermite(2*n+1,sqrt(t)),t,s);
 487 -%pi*2^(2*n+2)*(1-1/s)^n/(2*gamma((2*n+1)/-2)*s^(3/2));
 488
 489 /******************************************************************************
 490   A&S 29.3.48 to A&S 29.3.60
 491
 492   Bessel I und Bessel J
 493 ******************************************************************************/
 494
 495 /*****************************************************************************
 496    Problem 66: A&S: 29.3.48 OK
 497 ******************************************************************************/
 498
 499 factor(specint(%e^(-s*t)*a*%e^(-a*t)*(bessel_i(1,a*t)+bessel_i(0,a*t)),t));
 500 -1*((sqrt(s)*sqrt(s+2*a)-s-2*a)/(sqrt(s)*sqrt(s+2*a)));
 501
 502 /*****************************************************************************
 503    Problem 67: A&S: 29.3.49 OK
 504 ******************************************************************************/
 505
 506 factor(specint(%e^(-s*t)*%e^(-(a+b)/2*t)*bessel_i(0,(a-b)/2*t),t));
 507 1/sqrt(s^2+(b+a)*s+a*b);
 508
 509 /*****************************************************************************
 510    Problem 68: A&S: 29.3.50 OK New
 511    Maxima doesn't know the simplification for gamma(2*k), but the
 512    answer is correct.
 513    The original code gives additional phase factors.
 514 ******************************************************************************/
 515
 516 (assume(2*k-1>0),
 517  factor(radcan(specint(%e^(-s*t)*sqrt(%pi)*(t/(a-b))^(k-1/2)
 518                                *%e^(-(a+b)/2*t)*bessel_i(k-1/2,(a-b)/2*t),t))));
 519 2*sqrt(%pi)*gamma(2*k)/(2^(2*k)*gamma((2*k+1)/2)*(s+a)^k*(s+b)^k);
 520
 521 /*****************************************************************************
 522    Problem 69: A&S: 29.3.51 OK
 523 ******************************************************************************/
 524
 525 factor(specint(%e^(-s*t)*t*%e^(-(a+b)/2*t)*(bessel_i(0,(a-b)/2*t)
 526                                                     +bessel_i(1,(a-b)/2*t)),t));
 527 1/((s+b)*sqrt(s^2+(b+a)*s+a*b));
 528
 529 /*****************************************************************************
 530    Problem 70: A&S: 29.3.52 OK
 531 ******************************************************************************/
 532
 533 factor(specint(%e^(-s*t)*(1/t)*%e^(-a*t)*bessel_i(1,a*t),t));
 534 /*(sqrt(s+2*a)-sqrt(s))/(sqrt(s+2*a)+sqrt(s));*/
 535 -((sqrt(s)*sqrt(s+2*a)-s-a)/a);
 536
 537 /*****************************************************************************
 538    Problem 71: A&S: 29.3.53 OK New
 539    The additional phase factor vanish.
 540 ******************************************************************************/
 541
 542 /* The result isn't as simple as A&S but correct */
 543
 544 factor(radcan(specint(%e^(-s*t)*(k/t)*%e^(-(a+b)/2*t)*bessel_i(k,(a-b)/2*t),t)));
 545 /* (a-b)^k/(sqrt(s+a)+sqrt(s+b))^(2*k); */
 546 ((b-a)^k*k*(-1)^k*gamma(k))/(gamma(k+1)*(2*sqrt(s+a)*sqrt(s+b)+2*s+b+a)^k);
 547
 548 /*****************************************************************************
 549    Problem 72: A&S: 29.3.54 OK New
 550    The additional phase factor vanish.
 551 ******************************************************************************/
 552
 553 (assume (nu>0),
 554  factor(specint(%e^(-s*t)*(1/a^nu)*%e^(-a/2*t)*bessel_i(nu,a/2*t),t)));
 555 /*((sqrt(s+a)+sqrt(s))^(-2*nu))/(sqrt(s)*sqrt(s+a));*/
 556 1/(sqrt(s)*sqrt(s+a)*(2*sqrt(s)*sqrt(s+a)+2*s+a)^nu);
 557
 558 /*****************************************************************************
 559    Problem 73: A&S: 29.3.55 OK
 560 ******************************************************************************/
 561
 562 factor(specint(%e^(-s*t)*bessel_j(0,a*t),t));
 563 1/sqrt(s^2+a^2);
 564
 565 /*****************************************************************************
 566    Problem 74: A&S: 29.3.56 OK
 567 ******************************************************************************/
 568
 569 factor(specint(%e^(-s*t)*a^nu*bessel_j(nu,a*t),t));
 570 a^(2*nu)/(sqrt(s^2+a^2)*(sqrt(s^2+a^2)+s)^nu);
 571
 572 /*****************************************************************************
 573    Problem 75: A&S: 29.3.57 OK
 574    Maxima doesn't know the simplification for gamma(2*k)
 575 ******************************************************************************/
 576
 577 factor(specint(%e^(-s*t)*sqrt(%pi)/gamma(k)*(t/(2*a))^(k-1/2)*bessel_j(k-1/2,a*t),t));
 578 2*sqrt(%pi)*gamma(2*k)/(2^(2*k)*gamma(k)*gamma((2*k+1)/2)*(s^2+a^2)^k);
 579
 580 /*****************************************************************************
 581    Problem 76: A&S: 29.3.58 OK
 582    The result of Maxima is correct and can be further simplified to
 583    (sqrt(s^2+a^2)-s)^k
 584 ******************************************************************************/
 585
 586 factor(specint(%e^(-s*t)*(k*a^k/t)*bessel_j(k,a*t),t));
 587 a^(2*k)*k*gamma(k)/(gamma(k+1)*(sqrt(s^2+a^2)+s)^k);
 588
 589 /*****************************************************************************
 590    Problem 77: A&S: 29.3.59 OK New
 591    The additional phase factors vanish.
 592    The result can be further simplified to (s-sqrt(s^2-a^2))^nu/sqrt(s^2-a^2)
 593 ******************************************************************************/
 594
 595 factor(specint(%e^(-s*t)*a^nu*bessel_i(nu,a*t),t));
 596 a^(2*nu)/(sqrt(s^2-a^2)*(sqrt(s^2-a^2)+s)^nu);
 597
 598 /*****************************************************************************
 599    Problem 78: A&S: 29.3.60 OK New
 600    Use gamma(2*z)=2^(2*z-1)/sqrt(%pi)*gamma(z)*gamma(z+1/2) to show that
 601    the result of Maxima is equivalent to the result of A&S: 1/(s^2-a^2)^k;
 602    The additional phase factors vanish.
 603 ******************************************************************************/
 604
 605 radcan(specint(%e^(-s*t)*sqrt(%pi)/gamma(k)*(t/(2*a))^(k-1/2)
 606                                                        *bessel_i(k-1/2,a*t),t));
 607 2*sqrt(%pi)*gamma(2*k)/(2^(2*k)*gamma(k)*gamma((2*k+1)/2)*(s-a)^k*(s+a)^k);
 608
 609 /*****************************************************************************
 610    Problem 79: A&S: 29.3.61 OK New
 611    An algorithm for the Unit Step function has been implemented.
 612 ******************************************************************************/
 613
 614 specint(%e^(-s*t)*unit_step(t-k),t);
 615 1/s*%e^(-s*k);
 616
 617 /*****************************************************************************
 618    Problem 80: A&S: 29.3.62 OK New
 619 ******************************************************************************/
 620
 621 radcan(specint(%e^(-s*t)*(t-k)*unit_step(t-k),t));
 622 1/s^2*%e^(-s*k);
 623
 624 /*****************************************************************************
 625    Problem 81: A&S: 29.3.63 OK New
 626 ******************************************************************************/
 627
 628 ratsimp(specint(%e^(-s*t)*(t-k)^(nu-1)/gamma(nu)*unit_step(t-k),t));
 629 1/s^nu*%e^(-k*s);
 630
 631 /*****************************************************************************
 632    Problem 82: A&S: 29.3.64 OK New
 633 ******************************************************************************/
 634
 635 ratsimp(specint(%e^(-s*t)*(unit_step(t)-unit_step(t-k)),t));
 636 (%e^(-s*k)*(%e^(k*s)-1))/s;
 637
 638 /*****************************************************************************
 639    Problem 83: A&S: 29.3.65 OK New
 640 ******************************************************************************/
 641
 642 ratsimp(specint(%e^(-s*t)*sum(unit_step(t-n*k),n,0,inf),t)),simpsum;
 643 %e^(k*s)/(s*%e^(k*s)-s);
 644
 645 /*****************************************************************************
 646    Problem 84: A&S: 29.3.66 (OK NEW)
 647    $SPECINT gets a result with a sum. This sum simplifies with simpsum to
 648    the correct result: 1/(s*%e^(k*s)-a*s)
 649
 650    But Maxima ask in the simplifier of the sum for the sign of %e^(k*?psey)-a.
 651    We have a problem to give Maxima this as a rule. It doesn't work.
 652    But when we do it in two steps and give a rule for %e^(k*s)-a it works and
 653    we get the expected result.
 654 ******************************************************************************/
 655
 656 block([result],
 657  (assume(%e^(k*s)-a>0),
 658  result:ratsimp(specint(%e^(-s*t)*sum(a^(n-1)*unit_step(t-n*k),n,1,inf),t)),
 659  ev(result,simpsum)));
 660 %e^(-k*s)/(s*(1-a*%e^(-k*s)));
 661
 662 /*****************************************************************************
 663    Problem 85: A&S: 29.3.67 OK New
 664    The result of Maxima is equivalent to 1/s*tanh(k*s)
 665 ******************************************************************************/
 666
 667 factor(ratsimp(specint(%e^(-s*t)*(unit_step(t)+2*sum((-1)^n*unit_step(t-2*n*k),n,1,inf)),t))),simpsum;
 668 (%e^(k*s)-1)*(%e^(k*s)+1)/(s*(%e^(2*k*s)+1));
 669
 670 /*****************************************************************************
 671    Problem 86: A&S: 29.3.68 OK New
 672 ******************************************************************************/
 673
 674 ratsimp(specint(%e^(-s*t)*sum((-1)^n*unit_step(t-n*k),n,0,inf),t)),simpsum;
 675 %e^(k*s)/(s*%e^(k*s)+s);
 676
 677 /*****************************************************************************
 678    Problem 87: A&S: 29.3.69 OK New
 679 ******************************************************************************/
 680
 681 factor(ratsimp(specint(%e^(-s*t)*(t*unit_step(t)+2*sum((-1)^n*(t-2*n*k)
 682                                      *unit_step(t-2*n*k),n,1,inf)),t))),simpsum;
 683 (%e^(k*s)-1)*(%e^(k*s)+1)/(s^2*(%e^(2*k*s)+1));
 684
 685 /*****************************************************************************
 686    Problem 88: A&S: 29.3.70 OK New
 687 ******************************************************************************/
 688
 689 factor(ratsimp(specint(%e^(-s*t)
 690                             *2*sum(unit_step(t-(2*n+1)*k),n,0,inf),t))),simpsum;
 691 2*%e^(k*s)/(s*(%e^(k*s)-1)*(%e^(k*s)+1));
 692
 693 /*****************************************************************************
 694    Problem 89: A&S: 29.3.71 OK New
 695 ******************************************************************************/
 696
 697 factor(ratsimp(specint(%e^(-s*t)*2*sum((-1)^n
 698                                   *unit_step(t-(2*n+1)*k),n,0,inf),t))),simpsum;
 699 2*%e^(k*s)/(s*(%e^(2*k*s)+1));
 700
 701 /*****************************************************************************
 702    Problem 90: A&S: 29.3.72 OK New
 703 ******************************************************************************/
 704
 705 factor(ratsimp(specint(%e^(-s*t)*(unit_step(t)
 706                                +2*sum(unit_step(t-2*n*k),n,1,inf)),t))),simpsum;
 707 (%e^(2*k*s)+1)/(s*(%e^(k*s)-1)*(%e^(k*s)+1));
 708
 709 /*****************************************************************************
 710    Problem 91: A&S: 29.3.73 (OK Noun form New)
 711    We have no algorithm for the Abs function.
 712    The result is k/(s^2+k^2)*coth(%pi*s/(2*k))
 713    Because we exponentialize the integrand, we don't get
 714    a noun form with the original Sin function.
 715 ******************************************************************************/
 716
 717 specint(%e^(-s*t)*abs(sin(k*t)),t);
 718 'specint(%e^-(s*t)*abs(%e^(2*%i*k*t)-1),t)/2$
 719
 720 /*****************************************************************************
 721    Problem 92: A&S: 29.3.74 OK New
 722 ******************************************************************************/
 723
 724 factor(ratsimp(specint(%e^(-s*t)*sum((-1)^n*sin(t)
 725                                       *unit_step(t-n*%pi),n,0,inf),t))),simpsum;
 726 %e^(%pi*s)/((s^2+1)*(%e^(%pi*s)-1));
 727
 728 /*****************************************************************************
 729    Problem 93: A&S: 29.3.75 OK
 730 ******************************************************************************/
 731
 732 factor(specint(%e^(-s*t)*bessel_j(0,2*sqrt(k*t)),t));
 733 %e^(-k/s)/s;
 734
 735 /*****************************************************************************
 736    Problem 94: A&S: 29.3.76 OK New
 737 ******************************************************************************/
 738
 739 ratsimp(specint(%e^(-s*t)*1/sqrt(%pi*t)*cos(2*sqrt(k*t)),t));
 740 %e^(-k/s)/sqrt(s);
 741
 742 /*****************************************************************************
 743    Problem 95: A&S: 29.3.77 OK New
 744 ******************************************************************************/
 745
 746 ratsimp(specint(%e^(-s*t)*1/sqrt(%pi*t)*cosh(2*sqrt(k*t)),t));
 747 %e^(k/s)/sqrt(s);
 748
 749 /*****************************************************************************
 750    Problem 96: A&S: 29.3.78 OK New
 751 ******************************************************************************/
 752
 753 ratsimp(specint(%e^(-s*t)*1/sqrt(%pi*k)*sin(2*sqrt(k*t)),t));
 754 s^(-3/2)*%e^(-k/s);
 755
 756 /*****************************************************************************
 757    Problem 97: A&S: 29.3.79 OK New
 758 ******************************************************************************/
 759
 760 ratsimp(specint(%e^(-s*t)*1/sqrt(%pi*k)*sinh(2*sqrt(k*t)),t));
 761 s^(-3/2)*%e^(k/s);
 762
 763 /*****************************************************************************
 764    Problem 98: A&S: 29.3.80 OK
 765 ******************************************************************************/
 766
 767 (assume(nu-1>0),
 768  ratsimp(specint(%e^(-s*t)*(t/k)^((nu-1)/2)*bessel_j(nu-1,2*sqrt(k*t)),t)));
 769 (1/s^nu)*%e^(-k/s);
 770
 771 /*****************************************************************************
 772    Problem 99: A&S: 29.3.81 OK New
 773 ******************************************************************************/
 774
 775 (assume(nu-1>0),
 776  ratsimp(specint(%e^(-s*t)*(t/k)^((nu-1)/2)*bessel_i(nu-1,2*sqrt(k*t)),t)));
 777 (1/s^nu)*%e^(k/s);
 778
 779 /******************************************************************************
 780    Problem 100: A&S: 29.3.82 OK New
 781 ******************************************************************************/
 782
 783 (block [besselexpand:true],
 784   ratsimp(exponentialize(
 785     specint(%e^(-s*t)*k/(2*sqrt(%pi*t^3))*%e^(-k^2/(4*t)),t))));
 786 %e^(-k*sqrt(s));
 787
 788 /******************************************************************************
 789    Problem 101: A&S: 29.3.83 (OK noun form New)
 790    Maxima doesn't find a result. The noun form is correct.
 791    Tabulated result is: 1/s*%e^(-k*sqrt(s))
 792 ******************************************************************************/
 793
 794 specint(%e^(-s*t)*erfc(k/(2*sqrt(t))),t);
 795 'specint(%e^(-s*t)*erfc(k/(2*sqrt(t))),t);
 796
 797 /*****************************************************************************
 798    Problem 102: A&S: 29.3.84 OK New
 799 ******************************************************************************/
 800
 801 ratsimp(exponentialize(specint(%e^(-s*t)/sqrt(%pi*t)*%e^(-k^2/(4*t)),t)));
 802 1/sqrt(s)*%e^(-k*sqrt(s));
 803
 804 /*****************************************************************************
 805    Problem 103: A&S: 29.3.85 (OK noun form New)
 806    Tabulated result: 1/s^2*%e^(-k*sqrt(s))
 807 ******************************************************************************/
 808
 809 specint(%e^(-s*t)*2*sqrt(t)*%i*erfc(k/(2*sqrt(t))),t);
 810 'specint(%e^(-s*t)*2*sqrt(t)*%i*erfc(k/(2*sqrt(t))),t);
 811
 812 /*****************************************************************************
 813    Problem 104: A&S: 29.3.86 (OK noun form New)
 814    Tabulated result: 1/(1+n/2)*%e^(-k*sqrt(s))
 815 ******************************************************************************/
 816
 817 specint(%e^(-s*t)*(4*t)^(n/2)*%i*erfc(k/(2*sqrt(t))),t);
 818 'specint(%e^(-s*t)*(4*t)^(n/2)*%i*erfc(k/(2*sqrt(t))),t);
 819
 820 /*****************************************************************************
 821    Problem 105: A&S: 29.3.87 (OK noun form New)
 822    We have no algorithm for an argument 1/sqrt(t).
 823    Tabulated result: s^((n-1)/2)*%e^(-k*sqrt(s));
 824 ******************************************************************************/
 825
 826 specint(%e^(-s*t)*%e^(-k^2/(4*t))/(2^n*sqrt(%pi*t^(n+1)))
 827                                                   *hermite(n,(k/(2*sqrt(t)))),t);
 828 'specint(hermite(n,k/(2*sqrt(t)))*t^(-n/2-1/2)*%e^(-s*t-k^2/(4*t)),t)
 829  /(sqrt(%pi)*2^n);
 830
 831 /*****************************************************************************
 832    Problem 106: A&S: 29.3.88 (Ok noun form New)
 833    Tabulated result: %e^(-k*sqrt(s))/(a*sqrt(s));
 834    We have a sum with 3 terms. Maxima gives a partially result for 2 terms.
 835 ******************************************************************************/
 836
 837 factor(ratsimp(specint(%e^(-s*t)*(1/sqrt(%pi*t)*exp(-k^2/(4*t))
 838                      -a*%e^(a*k)*%e^(a^2*t)*erfc(a*sqrt(t)+k/(2*sqrt(t)))),t)));
 839 (sqrt(s)*'specint(-a*%e^(-s*t+a^2*t+a*k)*erfc((2*a*t+k)/(2*sqrt(t))),t)
 840                                       -sinh(k*sqrt(s))+cosh(k*sqrt(s)))/sqrt(s);
 841
 842 /*****************************************************************************
 843    Problem 107: A&S: 29.3.89 (OK noun form New)
 844    Tabulated result: a*e^(-k*sqrt(s))/(s*(a+sqrt(s)));
 845 ******************************************************************************/
 846
 847 specint(%e^(-s*t)*(-%e^(a*k)*%e^(a^2*t)*erfc(a*sqrt(t)+k/(2*sqrt(t)))
 848                                                        +erfc(k/(2*sqrt(t)))),t);
 849 'specint(-erfc(a*sqrt(t)+k/(2*sqrt(t)))*%e^(-s*t+a^2*t+a*k),t)
 850                                      +'specint(erfc(k/(2*sqrt(t)))*%e^-(s*t),t);
 851
 852 /*****************************************************************************
 853    Problem 108: A&S: 29.3.90 (OK noun form New)
 854    Tabulated result: %e^(-k*sqrt(s))/(sqrt(s)*(a+sqrt(s)));
 855 ******************************************************************************/
 856
 857 specint(%e^(-s*t)*%e^(a*k)*%e^(a^2*t)*erfc(a*sqrt(t)+k/(2*sqrt(t))),t);
 858 'specint(%e^(-s*t)*%e^(a*k)*%e^(a^2*t)*erfc(a*sqrt(t)+k/(2*sqrt(t))),t);
 859
 860 /*****************************************************************************
 861    Problem 109: A&S: 29.3.91 (Ok noun form New)
 862    Tabulated result: %e^(-k*sqrt(s*(s+a)))/sqrt(s*(s+a))
 863 ******************************************************************************/
 864
 865 specint(%e^(-s*t)*%e^(a*t/2)*bessel_i(0,1/2*a*sqrt(t^2-k^2))*unit_step(t-k),t);
 866 'specint(%e^(-s*t)*%e^(a*t/2)*bessel_i(0,1/2*a*sqrt(t^2-k^2))*unit_step(t-k),t);
 867
 868 /*****************************************************************************
 869    Problem 110: A&S: 29.3.92 (OK noun form New)
 870    Tabulated result: %e^(-k*sqrt(s^2+a^2))/sqrt(s^2+a^2)
 871 ******************************************************************************/
 872
 873 specint(%e^(-s*t)*bessel_j(0,a*sqrt(t^2-k^2))*unit_step(t-k),t);
 874 'specint(%e^(-s*t)*bessel_j(0,a*sqrt(t^2-k^2))*unit_step(t-k),t);
 875
 876 /*****************************************************************************
 877    Problem 111: A&S: 29.3.93 (OK noun form New)
 878    Tabulated result: %e^(-k*sqrt(s^2-a^2))/sqrt(s^2-a^2)
 879 ******************************************************************************/
 880
 881 specint(%e^(-s*t)*bessel_i(0,a*sqrt(t^2-k^2))*unit_step(t-k),t);
 882 'specint(%e^(-s*t)*bessel_i(0,a*sqrt(t^2-k^2))*unit_step(t-k),t);
 883
 884 /*****************************************************************************
 885    Problem 112: A&S: 29.3.94 (OK noun form New)
 886    Tabulated result: %e^(-k*(sqrt(s^2+a^2))-s)/sqrt(s^2+a^2)
 887 ******************************************************************************/
 888
 889 specint(%e^(-s*t)*bessel_j(0,a*sqrt(t^2+2*k*t)),t);
 890 'specint(%e^(-s*t)*bessel_j(0,a*sqrt(t)*sqrt(t+2*k)),t);
 891
 892 /*****************************************************************************
 893    Problem 113: A&S: 29.3.95 (OK noun form New)
 894    Tabulated result: %e^(-k*s)-%e^(-k*sqrt(s^2+a^2))
 895 ******************************************************************************/
 896
 897 specint(%e^(-s*t)*a*k/sqrt(t^2-k^2)*bessel_j(1,a*sqrt(t^2-k^2))*unit_step(t-k),t);
 898 'specint(%e^(-s*t)*a*k/sqrt(t^2-k^2)*bessel_j(1,a*sqrt(t^2-k^2))*unit_step(t-k),t);
 899
 900 /*****************************************************************************
 901    Problem 114: A&S: 29.3.96 (OK noun form New)
 902    Tabulated result: -%e^(-k*s)+%e^(-k*sqrt(s^2-a^2))
 903 ******************************************************************************/
 904
 905 specint(%e^(-s*t)*a*k/sqrt(t^2-k^2)*bessel_i(1,a*sqrt(t^2-k^2))*unit_step(t-k),t);
 906 'specint(%e^(-s*t)*a*k/sqrt(t^2-k^2)*bessel_i(1,a*sqrt(t^2-k^2))*unit_step(t-k),t);
 907
 908 /*****************************************************************************
 909    Problem 115: A&S: 29.3.97 (OK noun form New)
 910    Tabulated result:
 911       a^nu*%e^(-k*sqrt(s^2+a^2))/(sqrt(s^2+a^2)*a^2*(sqrt(s^2+a^2)+s)^nu)
 912 ******************************************************************************/
 913
 914 specint(%e^(-s*t)*((t-k)/(t+k))^(nu/2)*bessel_j(nu,a*sqrt(t^2-k^2))*unit_step(t-k),t);
 915 'specint(%e^(-s*t)*((t-k)/(t+k))^(nu/2)*bessel_j(nu,a*sqrt(t^2-k^2))*unit_step(t-k),t);
 916
 917 /*****************************************************************************
 918    Problem 116: A&S: 29.3.98 OK New
 919 ******************************************************************************/
 920
 921 ratsimp(specint(%e^(-s*t)*(-%gamma-log(t)),t));
 922 1/s*log(s);
 923
 924 /*****************************************************************************
 925    Problem 117: A&S: 29.3.99 OK New
 926 ******************************************************************************/
 927
 928 ratsimp(specint(%e^(-s*t)*t^(k-1)/gamma(k)*(psi[0](k)-log(t)),t));
 929 1/s^k*log(s);
 930
 931 /*****************************************************************************
 932    Problem 118: A&S: 29.3.100 OK NEW
 933    We support the expintegral_e1 and get a new result.
 934    Tabulated result: log(s)/(s-a);
 935 ******************************************************************************/
 936
 937 ratsimp(specint(%e^(-s*t)*%e^(a*t)*(log(a)+expintegral_e1(a*t)),t));
 938 log(s)/(s-a);
 939
 940 /*****************************************************************************
 941    Problem 119: A&S: 29.3.101 OK NEW
 942
 943    Tabulated result: log(s)/(s^2+1)
 944 ******************************************************************************/
 945
 946 radcan(specint(%e^(-s*t)*(cos(t)*expintegral_si(t)-sin(t)*expintegral_ci(t)),t)),
 947   logarc:true,logexpand:super;
 948 log(s)/(s^2+1);
 949
 950 /*****************************************************************************
 951    Problem 120: A&S: 29.3.102 OK NEW
 952    Tabulated result: s*log(s)/(s^2+1)
 953 ******************************************************************************/
 954
 955 radcan(specint(%e^(-s*t)*(-sin(t)*expintegral_si(t)-cos(t)*expintegral_ci(t)),t)),
 956    logarc:true,logexpand:super;
 957 s*log(s)/(s^2+1);
 958
 959 /*****************************************************************************
 960    Problem 121: A&S: 29.3.103 OK New
 961    We support the expintegral_e1 and get a new result.
 962    Tabulated result: 1/s*log(1+k*s)
 963 ******************************************************************************/
 964
 965 specint(%e^(-s*t)*expintegral_e1(t/k),t);
 966 log(1+k*s)/s;
 967
 968 /*****************************************************************************
 969    Problem 122: A&S: 29.3.104 OK New
 970 ******************************************************************************/
 971
 972 specint(%e^(-s*t)*1/t*(%e^(-b*t)-%e^(-a*t)),t);
 973 log((s+a)/(s+b));
 974
 975 /*****************************************************************************
 976    Problem 123: A&S: 29.3.105 OK NEW
 977    Tabulated result: 1/s*log(1+k^2*s^2)
 978 ******************************************************************************/
 979
 980 radcan(specint(%e^(-s*t)*(-2)*expintegral_ci(t/k),t));
 981 1/s*log(1+k^2*s^2);
 982
 983 /*****************************************************************************
 984    Problem 124: A&S: 29.3.106 OK NEW
 985    Tabulated result: 1/s*log(s^2+a^2)
 986 ******************************************************************************/
 987
 988 ratsimp(specint(%e^(-s*t)*(2*log(a)-2*expintegral_ci(a*t)),t));
 989 1/s*log(s^2+a^2);
 990
 991 /*****************************************************************************
 992    Problem 125: A&S: 29.3.107 OK NEW
 993    We get a correct noun form with a partially evaluated integral.
 994    Tabulated result: 1/s^2*log(s^2+a^2)
 995 ******************************************************************************/
 996
 997 ratsimp(specint(%e^(-s*t)*(2*t*log(a)+2/a*sin(a*t)-2*t*expintegral_ci(a*t)),t));
 998 1/s^2*log(s^2+a^2);
 999
1000 /*****************************************************************************
1001    Problem 126: A&S: 29.3.108 OK New
1002    Tabulated result: log((s^2+a^2)/s^2);
1003 ******************************************************************************/
1004
1005 specint(%e^(-s*t)*2/t*(1-cos(a*t)),t);
1006 log(a^2/s^2+1);
1007
1008 /*****************************************************************************
1009    Problem 127: A&S: 29.3.109 OK New
1010    Tabulated result: log((s^2-a^2)/s^2);
1011 ******************************************************************************/
1012
1013 specint(%e^(-s*t)*2*(1-cosh(a*t))/t,t);
1014 log(1-a^2/s^2);
1015
1016 /*****************************************************************************
1017    Problem 128: A&S: 29.3.110 OK New
1018    arctan(k/s);
1019 ******************************************************************************/
1020
1021 radcan(specint(%e^(-s*t)*sin(k*t)/t,t));
1022 -((%i*log(s+%i*k)-%i*log(s-%i*k)) /2);
1023
1024 /*****************************************************************************
1025    Problem 129: A&S: 29.3.111 (OK Noun form New)
1026    1/s*arctan(k/s);
1027 ******************************************************************************/
1028
1029 specint(%e^(-s*t)*expintegral_si(k*t),t);
1030 1/s*atan(k/s);
1031
1032 /*****************************************************************************
1033    Problem 130: A&S: 29.3.112 OK
1034 ******************************************************************************/
1035
1036 factor(ratsimp(specint(%e^(-s*t)*1/(k*sqrt(%pi))*exp(-t^2/(4*k^2)),t)));
1037 -%e^(k^2*s^2)*(erf(k*s)-1);
1038
1039 /*****************************************************************************
1040    Problem 131: A&S: 29.3.113 (Ok noun form New)
1041    Tabulated result: 1/s*%e^(k^2*s^2)*erfc(k*s);
1042 ******************************************************************************/
1043
1044 specint(%e^(-s*t)*erf(t/(2*k)),t);
1045 'specint(%e^(-s*t)*erf(t/(2*k)),t);
1046
1047 /*****************************************************************************
1048    Problem 132: A&S: 29.3.114 OK New
1049    Tabulated result is %e^(k*s)*erfc(sqrt(k*s));
1050 ******************************************************************************/
1051
1052 radcan(specint(%e^(-s*t)*sqrt(k)/(%pi*sqrt(t)*(t+k)),t));
1053 %e^(k*s)*(1-erf(sqrt(k*s)));
1054
1055 /*****************************************************************************
1056    Problem 133: A&S: 29.3.115 OK New
1057    gamma_incomplete(1/2,a*s)=sqrt(%pi)*erfc(sqrt(a*s))
1058    1/sqrt(s)*erfc(sqrt(k*s));
1059 ******************************************************************************/
1060
1061 specint(%e^(-s*t)/sqrt(%pi*t)*unit_step(t-k),t),gamma_expand:true;
1062 1/sqrt(s)*erfc(sqrt(k)*sqrt(s));
1063
1064 /*****************************************************************************
1065    Problem 134: A&S: 29.3.116 OK New
1066    1/sqrt(s)*%e^(k*s)*erfc(sqrt(k*s));
1067 ******************************************************************************/
1068
1069 specint(%e^(-s*t)*1/sqrt(%pi*(t+k)),t),gamma_expand:true;
1070 erfc(sqrt(k)*sqrt(s))*%e^(k*s)/sqrt(s);
1071
1072 /*****************************************************************************
1073    Problem 135: A&S: 29.3.117 (OK noun form)
1074    We have t^-1*sin(sqrt(t)) and F35P147 fails because of factor t^-1.
1075    Tabulated result: erf(k/sqrt(s));
1076 ******************************************************************************/
1077
1078 radcan(specint(%e^(-s*t)*sin(2*k*sqrt(t))/(%pi*t),t));
1079 'specint(-%i*%e^(2*%i*k*sqrt(t)-s*t)/(2*%pi*t),t)
1080    +'specint(%i*%e^(-s*t-2*%i*k*sqrt(t))/(2*%pi*t),t);
1081
1082 /*****************************************************************************
1083    Problem 136: A&S: 29.3.118 OK
1084 ******************************************************************************/
1085
1086 ratsimp(specint(%e^(-s*t)*1/sqrt(%pi*t)*%e^(-2*k*sqrt(t)),t));
1087 -((erf(k/sqrt(s))-1)*%e^(k^2/s)/sqrt(s));
1088
1089 /*****************************************************************************
1090    Problem 137: A&S: 29.3.119 (OK noun form New)
1091    bessel_k(0,k*s);
1092    We have no algorithm for 1/sqrt(t^2-k^2)
1093 ******************************************************************************/
1094
1095 specint(%e^(-s*t)*1/sqrt(t^2-k^2)*unit_step(t-k),t);
1096 'specint(%e^(-s*t)*1/sqrt(t^2-k^2)*unit_step(t-k),t);
1097
1098 /*****************************************************************************
1099    Problem 138: A&S: 29.3.120 OK New
1100 ******************************************************************************/
1101
1102 specint(%e^(-s*t)*1/(2*t)*exp(-k^2/(4*t)),t);
1103 bessel_k(0,k*sqrt(s));
1104
1105 /*****************************************************************************
1106    Problem 139: A&S: 29.3.121 OK New
1107    We know U(3/2,3,z)= %e^(k*s)*bessel_k(1,z/2)/(sqrt(%pi)*z) and we get
1108    the tabulated result: 1/s*%e^(k*s)*bessel_k(1,k*s);
1109 ******************************************************************************/
1110
1111 ratsimp(specint(%e^(-s*t)/k*sqrt(t*(t+2*k)),t));
1112 2*k*sqrt(%pi)*hypergeometric_u(3/2,3,(2*k*s));
1113
1114 /*****************************************************************************
1115    Problem 140: A&S: 29.3.122 OK New
1116    The original code had an extra factor.
1117 ******************************************************************************/
1118
1119 specint(%e^(-s*t)/k*%e^(-k^2/(4*t)),t);
1120 1/sqrt(s)*bessel_k(1,k*sqrt(s));
1121
1122 /*****************************************************************************
1123    Problem 141: A&S: 29.3.123 (OK noun formn New)
1124    1/sqrt(s)*%e^(k/s)*bessel_k(0,k/s);
1125 ******************************************************************************/
1126
1127 specint(%e^(-s*t)*2/sqrt(%pi*t)*bessel_k(0,2*sqrt(2*k*t)),t);
1128 'specint(2*bessel_k(0,2*sqrt(2)*sqrt(k)*sqrt(t))*%e^(-s*t)/(sqrt(%pi)*sqrt(t)),t);
1129
1130 /*****************************************************************************
1131    Problem 142: A&S: 29.3.124 OK New
1132    We get a result in parts. I think the algorithm to calculate t^v*(t+a)^w
1133    can be further generalized to get also a result for this integral.
1134    The problem is the pattern which doesn't match sqrt(2*k-t).
1135    The pattern is further generalizied. With
1136      U(1/2,1,z)=%e^(z/2)*K0(z/2)
1137      U(1/2,1-z)=%e^(-z/2)*(K0(z/2)+%i*%pi*I0(z/2))
1138    we get the expected result:
1139    %pi*%e^(-s*k)*bessel_i(0,k*s);
1140 ******************************************************************************/
1141
1142 specint(%e^(-s*t)*(unit_step(t)-unit_step(t-2*k))/sqrt(t*(2*k-t)),t);
1143 sqrt(%pi)*%i*'hypergeometric_u(1/2,1,2*k*s)*%e^-(2*k*s)-sqrt(%pi)*%i*'hypergeometric_u(1/2,1,-2*k*s);
1144
1145 /*****************************************************************************
1146    Problem 143: A&S: 29.3.125 Error
1147    Can we verrify the result: %e^(-s*k)*bessel_i(1,k*s)
1148    Here we have a problem: We get a result with four Hypergeometric U functions.
1149    But the result simplifies to something like
1150     (a * bessel_i(0,k*s)  + b * bessel_i(1,k*s))*%e^(-s*k)
1151    We have no bessel_i(0,k*s) in the expected result. What is wrong???
1152    We declare this test to be an error.
1153 ******************************************************************************/
1154
1155 radcan(specint(%e^(-s*t)*(unit_step(t)-unit_step(t-2*k))*(k-t)/(%pi*k*sqrt(t*(2*k-t))),t));
1156 done;
1157
1158 /*****************************************************************************
1159    Problem 144: A&S: 29.3.126 OK
1160    %e^(a*s)*E[1](a*s);
1161 ******************************************************************************/
1162
1163 specint(%e^(-s*t)/(t+a),t);
1164 gamma_incomplete(0,a*s)*%e^(a*s);
1165
1166 /*****************************************************************************
1167    Problem 145: A&S: 29.3.127 OK
1168    (1/a-s*%e^(a*s)*E[1](a*s));
1169 ******************************************************************************/
1170
1171 radcan(specint(%e^(-s*t)/(t+a)^2,t)),gamma_expand:true;
1172 (expintegral_ei(-a*s)*a*s*%e^(a*s)+1)/a;
1173
1174 /*****************************************************************************
1175    Problem 146: A&S: 29.3.128 OK
1176    (a^(1-n)*%e^(a*s)*E[n](a*s));
1177 ******************************************************************************/
1178
1179 specint(%e^(-s*t)/(t+a)^n,t);
1180 gamma_incomplete(1-n,a*s)*s^(n-1)*%e^(a*s);
1181
1182 /*****************************************************************************
1183    Problem 147: A&S: 29.3.129 (OK noun form)
1184    (%pi/2-Si(s))*cos(s)+Ci(s)*sin(s);
1185    For this case the noun form doesn't look nice.
1186 ******************************************************************************/
1187
1188 specint(%e^(-s*t)/(t^2+1),t);
1189 'specint(1/(t^2*%e^(s*t)+%e^(s*t)),t);
1190
1191 /*****************************************************************************
1192    End of tests from A&S
1193 ******************************************************************************/
1194 /*******************************************************************************
1195
1196   Do the tests with the tabulated Laplace transforms from EqWorld
1197
1198 *******************************************************************************/
1199
1200 (kill(all), assume(s>0, a>0, b>0, n>0, nu>-1), declare(n,integer, nu,real));
1201 done;
1202
1203 /*******************************************************************************
1204
1205    Expressions with Power-Law Functions
1206
1207 *******************************************************************************/
1208
1209 /* Problem 149: No. 1 OK */
1210
1211 specint(%e^(-s*t),t);
1212 1/s;
1213
1214 /* Problem 150: No. 2 OK New
1215    This is an example for the step function u(t).
1216    We have implemented an algorithm with unit_step.
1217 */
1218
1219 specint(%e^(-s*t)*(unit_step(t-a)-unit_step(t-b)),t);
1220 %e^-(a*s)/s-%e^-(b*s)/s;
1221
1222 /* Problem 151: No. 3 OK */
1223
1224 specint(%e^(-s*t)*t,t);
1225 1/s^2;
1226
1227 /* Problem 152: No. 4 OK New
1228    The result of Maxima is equivalent to the
1229    tabulated result: -%e^(a*s)*ei(-a*p)
1230 */
1231
1232 specint(%e^(-s*t)*(t+a)^-1,t);
1233 gamma_incomplete(0,a*s)*%e^(a*s);
1234
1235 /* Problem 153: No. 5 OK */
1236
1237 specint(%e^(-s*t)*t^n,t);
1238 gamma(n+1)/s^(1+n);
1239
1240 /* Problem 154: No. 6 OK */
1241
1242 specint(%e^(-s*t)*t^(n-1/2),t);
1243 gamma(n+1/2)/s^(1/2+n);
1244
1245 /* Problem 155: No. 7 OK New
1246    With gamma_incomplete(1/2,a*s)=sqrt(%pi)*erfc(sqrt(a*s)) the result of
1247    Maxima is equivalent to the
1248    tabulated result: sqrt(%pi/s)*%e^(a*s)*erfc(sqrt(a*s));
1249 */
1250
1251 specint(%e^(-s*t)/sqrt(t+a),t);
1252 gamma_incomplete(1/2,a*s)*%e^(a*s)/sqrt(s);
1253
1254 /* Problem 156: No. 8 OK New
1255    tabulated result: sqrt(%pi/s)-%pi*sqrt(a)*%e^(a*s)*erfc(sqrt(a*s))
1256 */
1257
1258 radcan(specint(%e^(-s*t)*sqrt(t)/(t+a),t));
1259 (sqrt(a)*(%pi*erf(sqrt(a)*sqrt(s))-%pi)*sqrt(s)*%e^(a*s)+sqrt(%pi))/sqrt(s);
1260
1261 /* Problem 157: No. OK New
1262    With gamma_incomplete(-1/2,z)=2*e^-z/sqrt(z)-2*sqrt(%pi)*erfc(sqrt(z))
1263    the result of Maxima is equivalent to the
1264    Tabulated result: 2*a^(-1/2)-2*(%pi*s)^(1/2)*%e^(a*s)*erfc(sqrt(a*s))
1265 */
1266
1267 specint(%e^(-s*t)*(t+a)^(-3/2),t);
1268 gamma_incomplete(-1/2,a*s)*sqrt(s)*%e^(a*s);
1269
1270 /* Problem 158: No. 10 OK New
1271    Tabulated result: (%pi/s)^(1/2)-%pi*a^(1/2)*e^(a*s)*erfc(sqrt(a*s));
1272 */
1273
1274 radcan(specint(%e^(-s*t)*t^(1/2)*(t+a)^(-1),t));
1275 (sqrt(a)*(%pi*erf(sqrt(a)*sqrt(s))-%pi)*sqrt(s)*%e^(a*s)+sqrt(%pi))/sqrt(s);
1276
1277 /* Problem 159: No. 11 OK New
1278    Tabulated result: %pi*a^(1/2)*e^(a*s)*erfc(sqrt(a*s));
1279 */
1280
1281 radcan(specint(%e^(-s*t)*t^(-1/2)*(t+a)^(-1),t));
1282 -((%pi*erf(sqrt(a)*sqrt(s))-%pi)*%e^(a*s)/sqrt(a));
1283
1284 /* Problem 160: No. 12 OK */
1285
1286 specint(%e^(-s*t)*t^nu,t);
1287 gamma(nu+1)/s^(nu+1);
1288
1289 /* Problem 161: No. 13 OK New
1290    Tabulated result: s^(-nu-1)*%e^(a*s)*gamma_incomplete(nu+1,a*s)
1291 */
1292
1293 specint(%e^(-s*t)*(t+a)^nu,t);
1294 gamma_incomplete(nu+1,a*s)*s^(-nu-1)*%e^(a*s);
1295
1296 /* Problem 162: No. 14 Ok New
1297    With
1298      gamma_incomplete(-v,a*s) = gamma(-v,a*s) - gamma_incomplete_lower(-v,a*s) and
1299      gamma(v+1) = v * gamma(v)
1300    the result of Maxima can be shown to be equivalent ot the
1301    tabulated result: a^nu*gamma(nu+1)*e^(a*s)*gamma_incomplete(-nu,a*s);
1302 */
1303
1304 radcan(specint(%e^(-s*t)*t^nu*(t+a)^(-1),t));
1305 (a^nu*gamma(-nu)*gamma(nu+1)-a^nu*nu*gamma_incomplete_lower(-nu,a*s)*gamma(nu))*%e^(a*s);
1306
1307 /*******************************************************************************
1308
1309    Expressions with Exponential Functions
1310
1311 *******************************************************************************/
1312
1313 /* Problem 163: No. 1 OK */
1314
1315 specint(%e^(-s*t)*%e^(-a*t),t);
1316 (s+a)^(-1);
1317
1318 /* Problem 164: No. 2 OK */
1319
1320 specint(%e^(-s*t)*t*%e^(-a*t),t);
1321 (s+a)^(-2);
1322
1323 /* Problem 165: No. 3 OK */
1324
1325 (assume(nu>0), specint(%e^(-s*t)*t^(nu-1)*%e^(-a*t),t));
1326 gamma(nu)*(s+a)^(-nu);
1327
1328 /* Problem 166: No. 4 OK New
1329    After extension of the algorithm we find the tabulated result.
1330 */
1331
1332 specint(%e^(-s*t)*(%e^(-a*t)-%e^(-b*t))/t,t);
1333 -log((s+a)/(s+b));
1334
1335 /* Problem 167: No. 5 OK New
1336    Tabulated result: (s+2*a)*log(s+2*a)+s*log(s)-2*(s+a)*log(s+a);
1337 */
1338
1339 specint(%e^(-s*t)*(1-%e^(-a*t))^2/t^2,t);
1340 (s+2*a)*log(s+2*a)+s*log(s)-2*(s+a)*log(s+a);
1341
1342 /* Problem 168: No. 6 OK */
1343
1344 factor(radcan(specint(%e^(-s*t)*%e^(-a*t^2),t)));
1345 - (sqrt(%pi)*%e^(s^2/(4*a))*(erf(s/(2*sqrt(a)))-1))/(2*sqrt(a));
1346
1347 /* Problem 169: No. 7 OK */
1348
1349 factor(expand(radcan(specint(%e^(-s*t)*t*%e^(-a*t^2),t))));
1350 (sqrt(%pi)*s*%e^(s^2/(4*a))*erf(s/(2*sqrt(a)))-sqrt(%pi)*s*%e^(s^2/(4*a))+2*sqrt(a))/(4*a^(3/2));
1351
1352 /* Problem 170: No. 8 OK */
1353
1354 specint(%e^(-s*t)*%e^(-a/t),t);
1355 2*sqrt(a/s)*bessel_k(1,2*sqrt(a*s));
1356
1357 /* Problem 171: No. 9 OK */
1358
1359 block([besselexpand:true,exponentialize:true],radcan(specint(%e^(-s*t)*sqrt(t)*%e^(-a/t),t)));
1360 ((2*sqrt(%pi)*sqrt(a)*sqrt(s)+sqrt(%pi))*%e^(-2*sqrt(a)*sqrt(s)))/(2*s^(3/2));
1361
1362 /* Problem 172: No. 10 OK */
1363
1364 block([besselexpand:true,exponentialize:true],radcan(specint(%e^(-s*t)/sqrt(t)*%e^(-a/t),t)));
1365 sqrt(%pi/s)*%e^(-2*sqrt(a*s));
1366
1367 /* Problem 173: No. 11 OK */
1368
1369 block([besselexpand:true,exponentialize:true],radcan(specint(%e^(-s*t)/(t*sqrt(t))*%e^(-a/t),t)));
1370 sqrt(%pi/a)*%e^(-2*sqrt(a*s));
1371
1372 /* Problem 174: No. 12 OK
1373    We calculate for an integer to prevent the conversion to bessel_i.
1374 */
1375
1376 factor(specint(%e^(-s*t)*t^(n-1)*%e^(-a/t),t));
1377 2*a^(n/2)*bessel_k(n,2*sqrt(a)*sqrt(s))/s^(n/2);
1378
1379 /* Problem 175: No. 13 OK
1380    Answer of EqWorld is
1381      s^(-1)-(%pi*a)^(1/2)*s^(-3/2)*e^(a/s)*erfc(sqrt(a/s))
1382 */
1383
1384 radcan(specint(%e^(-s*t)*%e^(-2*sqrt(a*t)),t));
1385 (sqrt(%pi)*sqrt(a)*(erf(sqrt(a)/sqrt(s))-1)*%e^(a/s)+sqrt(s))/s^(3/2);
1386
1387 /* Problem 176: No. 14 OK
1388    Result of EqWorld is (%pi/s)^(1/2)*%e^(a/s)*erfc(sqrt(a/s));
1389    Maxima doesn't convert erf(z)-1 -> erfc(z)
1390 */
1391
1392 radcan(specint(%e^(-s*t)/sqrt(t)*%e^(-2*sqrt(a*t)),t));
1393 -((sqrt(%pi)*(erf(sqrt(a)/sqrt(s))-1)*%e^(a/s))/sqrt(s));
1394
1395 /*******************************************************************************
1396
1397    Expressions with Trigonometric Functions
1398
1399 *******************************************************************************/
1400
1401 /* Problem 177: No. 1 OK */
1402
1403 radcan(specint(%e^(-s*t)*sin(a*t),t));
1404 a/(s^2+a^2);
1405
1406 /* Problem 178: No. 2 (OK noun form New)
1407    We have no algorithm for the abs function
1408    Tabulated result: a/(s^2+a^2)*coth(%pi*s/(2*a));
1409 */
1410
1411 specint(%e^(-s*t)*(abs(sin(a*t))),t);
1412 'specint(%e^-(s*t)*abs(%e^(2*%i*a*t)-1),t)/2$
1413
1414 /* Problem 179: No. 3 (OK noun form New)
1415    The result of EqWorld includes a product:
1416       a^(2*n)*(2*n)!/(s*product(s^2+(2*k*a)^2,k,1,n));
1417 */
1418
1419 specint(%e^(-s*t)*(sin(a*t))^(2*n),t);
1420 'specint((-1)^n*(%e^(2*%i*a*t)-1)^(2*n)*%e^(-s*t-2*%i*a*n*t),t)/2^(2*n);
1421
1422 /* Problem 180: No. 4 (OK noun form New)
1423    The result of EqWorld includes a product:
1424       a^(2*n+1)*(2*n+1)!/((s^2+a^2)*product(s^2+(2*k+1)^2*a^2,k,1,n));
1425 */
1426
1427 factor(specint(%e^(-s*t)*(sin(a*t))^(2*n+1),t));
1428
1429 /*
1430 specint(-%i*2^(-2*n-1)*(-1)^n*(%e^(2*%i*a*t)-1)^(2*n)*%e^(-s*t-2*%i*a*n*t+%i*a*t),t)+specint(%i*2^(-2*n-1)*(-1)^n*(%e^(2*%i*a*t)-1)^(2*n)*%e^(-s*t-2*%i*a*n*t-%i*a*t),t);
1431 */
1432
1433 'specint(-%i*2^(-2*n-1)*(-1)^n*(%e^(2*%i*a*t)-1)^(2*n)
1434             *%e^(-s*t-2*%i*a*n*t+%i*a*t),t)
1435  +'specint(%i*2^(-2*n-1)*(-1)^n*(%e^(2*%i*a*t)-1)^(2*n)
1436              *%e^(-s*t-2*%i*a*n*t-%i*a*t),t);
1437
1438 /* Problem 181: No. 5 OK New */
1439
1440 specint(%e^(-s*t)*t^n*sin(a*t),t);
1441 (%i*gamma(n+1)*(s+%i*a)^(-n-1)-%i*gamma(n+1)*(s-%i*a)^(-n-1))/2;
1442
1443 /* Problem 182: No. 6 OK New
1444    Tabulated result: atan(a/s);
1445    This is equivalent to the result of Maxima after extension of the algorithm.
1446 */
1447
1448 specint(%e^(-s*t)*sin(a*t)/t,t);
1449 %i*log((s-%i*a)/(s+%i*a))/2;
1450
1451 /* Problem 183: No. 7 OK New
1452    The result of Maxima is equivalent to the
1453    tabulated result: 1/4*log(1+4*a^2*s^(-2));
1454 */
1455
1456 radcan(specint(%e^(-s*t)/t*(sin(a*t))^2,t));
1457 (log(s^2+4*a^2)-2*log(s))/4;
1458
1459 /* Problem 184: No. 8 OK New
1460    The result of Maxima can be shown to be equivalent to the
1461    tabulated result: a*atan(2*a/s)-1/4*s*log(1+4*a^2*p^(-2));
1462 */
1463
1464 expand(radcan(specint(%e^(-s*t)/t^2*(sin(a*t))^2,t)));
1465 -s*log(s+2*%i*a)/4-%i*a*log(s+2*%i*a)/2-s*log(s-2*%i*a)/4+%i*a*log(s-2*%i*a)/2+s*log(s)/2;
1466
1467 /* Problem 185: No. 9 OK New
1468    We correct a bug in f35p147.
1469 */
1470
1471 radcan(specint(%e^(-s*t)*sin(2*sqrt(a*t)),t));
1472 sqrt(%pi)*sqrt(a)*%e^-(a/s)/s^(3/2);
1473
1474 /* Problem 186: No. 10 (OK noun form New)
1475    For this Problem the test f35p147test fail.
1476    The exponent v is negativ. We have no algorithm for this case.
1477    Tabulated answer is:
1478    %pi*erf(sqrt(a/s));
1479 */
1480
1481 radcan(specint(%e^(-s*t)*sin(2*sqrt(a*t))/t,t));
1482 specint(-%i*%e^(2*%i*sqrt(a)*sqrt(t)-s*t)/(2*t),t)+specint(%i*%e^(-s*t-2*%i*sqrt(a)*sqrt(t))/(2*t),t);
1483
1484 /* Problem 187: No. 11 OK */
1485
1486 radcan(specint(%e^(-s*t)*cos(a*t),t));
1487 s/(s^2+a^2);
1488
1489 /* Problem 188: No. 12 OK */
1490
1491 radcan(specint(%e^(-s*t)*(cos(a*t))^2,t));
1492 (s^2+2*a^2)/(s^3+4*a^2*s);
1493
1494 /* Problem 189: No. 13 (OK Result is not verified.)
1495    Maxima get a result, but we can't verify it. See Problem 33.
1496    Tabulated result:
1497    n!*s^(n+1)/(s^2+a^2)^(n+1)*sum((-1)^k*C[2*k+1,n+1]((a/s)^(2*k)),k,0,n);
1498 */
1499
1500 radcan(specint(%e^(-s*t)*t^n*cos(a*t),t));
1501 ((s+%i*a)^n*(gamma(n+1)*s+%i*a*gamma(n+1))+(s-%i*a)^n*(gamma(n+1)*s-%i*a*gamma(n+1)))/((s-%i*a)^n*(s+%i*a)^n*(2*s^2+2*a^2));
1502
1503 /* Problem 190: No. 14 OK New*/
1504
1505 specint(%e^(-s*t)*(1-cos(a*t))/t,t);
1506 1/2*log(1+a^2*s^(-2));
1507
1508 /* Problem 191: No. 15 OK New
1509    The result of Maxima is equivalent to the tabulated result:
1510       1/2*log((s^2+b^2)/(s^2+a^2));
1511 */
1512
1513 radcan(specint(%e^(-s*t)/t*(cos(a*t)-cos(b*t)),t));
1514 (log(s+%i*b)+log(s-%i*b)-log(s+%i*a)-log(s-%i*a))/2;
1515
1516 /* Problem 192: No. 16 OK New*/
1517
1518 radcan(specint(%e^(-s*t)*sqrt(t)*cos(2*sqrt(a*t)),t));
1519 1/2*%pi^(1/2)*s^(-5/2)*(s-2*a)*%e^(-a/s);
1520
1521 /* Problem 193: No. 17 OK New
1522     After correction of a bug Maxima gets the correct answer.
1523 */
1524
1525 radcan(specint(%e^(-s*t)/sqrt(t)*cos(2*sqrt(a*t)),t));
1526 sqrt(%pi/s)*%e^(-a/s);
1527
1528 /* Problem 194: No. 18 OK */
1529
1530 factor(radcan(specint(%e^(-s*t)*sin(a*t)*sin(b*t),t)));
1531 2*a*b*s/((s^2+(b^2-2*a*b+a^2))*(s^2+(b^2+2*a*b+a^2)));
1532
1533 /* Problem 195: No. 19 OK */
1534
1535 factor(radcan(specint(%e^(-s*t)*cos(a*t)*sin(b*t),t)));
1536 (b*(s^2+b^2-a^2))/((s^2+(b^2-2*a*b+a^2))*(s^2+(b^2+2*a*b+a^2)));
1537
1538 /* Problem 196: No. 20 OK */
1539
1540 factor(radcan(specint(%e^(-s*t)*cos(a*t)*cos(b*t),t)));
1541 (s*(s^2+a^2+b^2))/((s^2+b^2-2*a*b+a^2)*(s^2+b^2+2*a*b+a^2));
1542
1543 /* Problem 197: No. 21 OK */
1544
1545 (assume(s-b>0), factor(radcan(specint(%e^(-s*t)*%e^(b*t)*sin(a*t),t))));
1546 a/(s^2-2*b*s+b^2+a^2);
1547
1548 /* Problem 198: No. 22 OK */
1549
1550 factor(radcan(specint(%e^(-s*t)*%e^(b*t)*cos(a*t),t)));
1551 (s-b)/(s^2-2*b*s+b^2+a^2);
1552
1553 /*******************************************************************************
1554
1555    Expressions with Hyperbolic Functions
1556
1557 *******************************************************************************/
1558
1559 /* Problem 199: No. 01 OK */
1560
1561 (assume(s-a>0), radcan(specint(%e^(-s*t)*sinh(a*t),t)));
1562 a/(s^2-a^2);
1563
1564 /* Problem 200: No. 02 OK */
1565
1566 (assume(s-2*a>0),radcan(specint(%e^(-s*t)*(sinh(a*t))^2,t)));
1567 2*a^2/(s^3-4*a^2*s);
1568
1569 /* Problem 201: No. 03 OK New
1570    Correct after correction of the code.
1571 */
1572
1573 radcan(specint(%e^(-s*t)/t*sinh(a*t),t));
1574 (log(s+a)-log(s-a))/2;
1575
1576 /* Problem 202: No. 04 OK */
1577
1578 factor(radcan(specint(%e^(-s*t)*t^(nu-1)*sinh(a*t),t)));
1579 gamma(nu)*((s+a)^nu-(s-a)^nu)/(2*(s-a)^nu*(s+a)^nu);
1580
1581 /* Problem 203: No. 05 OK New
1582    Maxima gives the wrong result 0.
1583    Correct after bug fix.
1584 */
1585
1586 radcan(specint(%e^(-s*t)*sinh(2*sqrt(a*t)),t));
1587 sqrt(%pi)*sqrt(a)*%e^(a/s)/s^(3/2);
1588
1589 /* Problem 204: No. 06 (OK Sign? New)
1590    Tabulated result:
1591      %pi^(1/2)*s^(-5/2)*(1/2*s+a)*%e^(a/s)*erf(sqrt(a/s))-a^(1/2)*s^(-2);
1592    After bug fix nearly the tabulated answer.
1593    The last term of the result has a different sign. What's wrong???
1594    All other examples are correct with the corrected code!
1595    No equivalent example by A&S. But sinh(2*sqrt(a*t)) and a lot of
1596    other examples can be shown to be correct. I think EqWorld is wrong.
1597 */
1598
1599 factor(radcan(block([gamma_expand:true], specint(%e^(-s*t)*sqrt(t)*sinh(2*sqrt(a*t)),t))));
1600 (sqrt(%pi)*erf(sqrt(a)/sqrt(s))*s*%e^(a/s)+2*sqrt(%pi)*a*erf(sqrt(a)/sqrt(s))*%e^(a/s)+2*sqrt(a)*sqrt(s))/(2*s^(5/2));
1601
1602 /* Problem 205: No. 07 OK New
1603    Correct after bug fix
1604 */
1605
1606 factor(radcan(specint(%e^(-s*t)/sqrt(t)*sinh(2*sqrt(a*t)),t)));
1607 %pi^(1/2)*s^(-1/2)*%e^(a/s)*erf(sqrt(a/s));
1608
1609 /* Problem 206: No. 08 (OK extra factor New)
1610    The result looks perfect, but in the tabulated answer a factor of 4
1611    in the exponent is missing. What's wrong???
1612    Tabulated answer:
1613      1/2*%pi^(1/2)*s^(-1/2)*(%e^(a/s)-1);
1614    I think EqWorld is wrong.
1615 */
1616
1617 radcan(specint(%e^(-s*t)*(sinh(2*sqrt(a*t)))^2/sqrt(t),t));
1618 sqrt(%pi)*(%e^(4*a/s)-1)/(2*sqrt(s)); /* Extra factor 4 ??? */
1619
1620 /* Problem 207: No. 09 OK */
1621
1622 factor(radcan(specint(%e^(-s*t)*cosh(a*t)^2,t)));
1623 (s^2-2*a^2)/(s*(s-2*a)*(s+2*a));
1624
1625 /* Problem 208: No. 10 OK */
1626
1627 factor(radcan(specint(%e^(-s*t)*cosh(a*t)^2,t)));
1628 (s^2-2*a^2)/(s*(s-2*a)*(s+2*a));
1629
1630 /* Problem 209: No. 11 OK */
1631
1632 factor(radcan(specint(%e^(-s*t)*t^(nu-1)*cosh(a*t),t)));
1633 gamma(nu)*((s+a)^nu+(s-a)^nu)/(2*(s-a)^nu*(s+a)^nu);
1634
1635 /* Problem 210: No. 12 OK New*/
1636
1637 radcan(specint(%e^(-s*t)*cosh(2*sqrt(a*t)),t));
1638 (sqrt(%pi)*sqrt(a)*erf(sqrt(a)/sqrt(s))*%e^(a/s)+sqrt(s))/s^(3/2);
1639
1640 /* Problem 211: No. 13 OK New
1641    Maxima gives the result:
1642      -(sqrt(%pi)*((erf(sqrt(a)/sqrt(s))-1)*s+2*a*erf(sqrt(a)/sqrt(s))-2*a)
1643          *%e^(a/s)
1644              +2*sqrt(a)*sqrt(s))/(2*s^(5/2))
1645    We have again the erf function in the result.
1646    The result is now correct.
1647 */
1648
1649 radcan(specint(%e^(-s*t)*sqrt(t)*cosh(2*sqrt(a*t)),t));
1650 sqrt(%pi)*(s+2*a)*%e^(a/s)/(2*s^(5/2));
1651
1652 /* Problem 212: No. 14 OK New
1653    Maxima gives the result:
1654      -sqrt(%pi)*(erf(sqrt(a)/sqrt(s))-1)*%e^(a/s)/sqrt(s)
1655    Again the erf function.
1656    Now correct.
1657 */
1658
1659 radcan(specint(%e^(-s*t)/sqrt(t)*cosh(2*sqrt(a*t)),t));
1660 sqrt(%pi)*%e^(a/s)/sqrt(s);
1661
1662 /* Problem 213: No. 15 (OK but additional factor ? New)
1663    Maxima gives the result:
1664      (sqrt(s+4*a)*(2*sqrt(%pi)*sqrt(s-4*a)+sqrt(%pi)*sqrt(s))
1665         + sqrt(%pi)*sqrt(s)*sqrt(s-4*a))
1666            /(4*sqrt(s)*sqrt(s-4*a)*sqrt(s+4*a))$
1667    The result doesn't contain the requiered exponential function.
1668    Now correct. But there is an extra factor 4 in the exponent.
1669    I think EqWorld is wrong.
1670 */
1671
1672 radcan(specint(%e^(-s*t)*cosh(2*sqrt(a*t))^2/sqrt(t),t));
1673 1/2*%pi^(1/2)*s^(-1/2)*(%e^(4*a/s)+1);
1674
1675 /*******************************************************************************
1676
1677    Expressions with Logarithmic Functions
1678
1679 *******************************************************************************/
1680
1681 /* Problem 214: No. 01 OK New*/
1682
1683 radcan(specint(%e^(-s*t)*log(t),t));
1684 -(1/s)*(log(s)+%gamma);
1685
1686 /* Problem 215: No. 02 (Ok noun form)
1687    Tabulated result: -(1/s)*%e^(s/a)*ei(-s/a);
1688 */
1689
1690 radcan(specint(%e^(-s*t)*log(1+a*t),t));
1691 'specint(%e^-(s*t)*log(a*t+1),t);
1692
1693 /* Problem 216: No. 03 (OK noun form New)
1694    Tabulated result: 1/s*(log(a)-%e^(a*s)*ei(-a*s));
1695 */
1696
1697 radcan(specint(%e^(-s*t)*log(t+a),t));
1698 specint(%e^(-s*t)*log(t+a),t);
1699
1700 /* Problem 217: No. 04 OK New
1701    Tabulated result is equivalent to Maxima result.
1702       n!/s^(n+1)*(sum(1/k,1,n)-log(s)-%gamma);
1703 */
1704
1705 radcan(specint(%e^(-s*t)*t^n*log(t),t));
1706 -s^(-n-1)*(gamma(n+1)*log(s)-psi[0](n+1)*gamma(n+1));
1707
1708 /* Problem 218: No. 05 OK New*/
1709
1710 radcan(specint(%e^(-s*t)*log(t)/sqrt(t),t));
1711 -(sqrt(%pi)*(log(s)+2*log(2)+%gamma)/sqrt(s));
1712
1713 /* Problem 219: No. 06 OK New*/
1714
1715 radcan(specint(%e^(-s*t)*t^(n-1/2)*log(t),t));
1716 -((gamma((2*n+1)/2)*log(s)-psi[0]((2*n+1)/2)*gamma((2*n+1)/2))/s^((2*n+1)/2));
1717
1718 /* Problem 220: No. 07 OK New*/
1719
1720 radcan(specint(%e^(-s*t)*t^(nu-1)*log(t),t));
1721 -((gamma(nu)*log(s)-psi[0](nu)*gamma(nu))/s^nu);
1722
1723 /* Problem 221: No. 08 (OK noun form New)
1724    1/s*((log(x)+%gamma)^2+1/6*%pi^2);
1725 */
1726
1727 radcan(specint(%e^(-s*t)*log(t)^2,t));
1728 specint(%e^(-s*t)*log(t)^2,t);
1729
1730 /* Problem 222: No. 09 OK New*/
1731
1732 radcan(specint(%e^(-s*t)*%e^(-a*t)*log(t),t));
1733 -((log(s+a)+%gamma)/(s+a));
1734
1735 /*******************************************************************************
1736
1737    Expressions with Error Functions
1738
1739 *******************************************************************************/
1740
1741 /* Problem 223: No. 01 (OK noun form New)
1742    Tabulated result:
1743       1/s*%e^((1/(2*a))^2*s^2)*erfc((1/(2*a))*s);
1744    No result. Is the noun form correct. Why do we get the function abs in the
1745    answer.
1746 */
1747
1748 radcan(specint(%e^(-s*t)*erf(a*t),t));
1749 specint(%e^(-s*t)*erf(a*t),t);
1750
1751 /* Problem 224: No. 02 OK */
1752
1753 radcan(specint(%e^(-s*t)*erf(sqrt(a*t)),t));
1754 sqrt(a)/(s*sqrt(s+a));
1755
1756 /* Problem 225: No. 03 OK */
1757
1758 radcan(specint(%e^(-s*t)*%e^(a*t)*erf(sqrt(a*t)),t));
1759 sqrt(a)/(sqrt(s)*(s-a));
1760
1761 /* Problem 226: No. 04 NEW Noun form
1762    Maxima gives:
1763       (sqrt(a)*sqrt(s)*sqrt(4-a*s)+4*asin(sqrt(a)*sqrt(s)/2))/(4*s)
1764    But there is no asin function in the expected answer.
1765    I can't verify the answer of Maxima.
1766    We declare this to an Error until verification of the problem.
1767
1768    This is a bug. The argument of the hypergeometric function for this example
1769    is 1/t. For an argument with a negavtive power the alogrithm for the
1770    Laplace transform is not correct.
1771 */
1772
1773 specint(%e^(-s*t)*erf(sqrt(a)/sqrt(t)/2),t);
1774 /*1/s*(1-%e^(-sqrt(a*s)));*/
1775 'specint(%e^(-s*t)*erf(sqrt(a)/sqrt(t)/2),t);
1776
1777 /* Problem 227: No. 05 OK */
1778
1779 radcan(specint(%e^(-s*t)*erfc(sqrt(a*t)),t));
1780 (sqrt(s+a)-sqrt(a))/(s*sqrt(s+a));
1781
1782 /* Problem 228: No. 06 OK */
1783
1784 radcan(specint(%e^(-s*t)*%e^(a*t)*erfc(sqrt(a*t)),t));
1785 1/(s+sqrt(s*a));
1786
1787 /* Problem 229: No. 07 (OK noun form New)
1788    The same as Problem 226.
1789    For this case erfc is transformed to an expression with Whittaker M
1790    functions. The Laplace transform fails and we get "2 'other-ca-later".
1791    Extension of the code to return a noun form.
1792    EqWorld: 1/s*%e^(-sqrt(a*s));
1793 */
1794
1795 radcan(specint(%e^(-s*t)*%e^(a*t)*erfc(sqrt(a/t)/2),t));
1796 'specint(erfc(sqrt(a)/(2*sqrt(t)))*%e^((a-s)*t),t);
1797
1798 /*******************************************************************************
1799
1800    Expressions with Bessel and Modified Bessel Functions
1801
1802 *******************************************************************************/
1803
1804 /* Problem 230: No. 01 OK */
1805
1806 radcan(specint(%e^(-s*t)*bessel_j(0,a*t),t));
1807 1/sqrt(s^2+a^2);
1808
1809 /* Problem 231: No. 02 OK */
1810
1811 radcan(specint(%e^(-s*t)*bessel_j(nu,a*t),t));
1812 a^nu/(sqrt(s^2+a^2)*(sqrt(s^2+a^2)+s)^nu);
1813
1814 /* Problem 232: No. 03 (OK)
1815    Thats the same result as Problem 233, because Maxima don't distinguish
1816    real and integer order of the Bessel function.
1817    But we have a problem with a factor 1/2, see Problem 86.
1818    I think EqWorld is wrong.
1819 */
1820
1821 radcan(specint(%e^(-s*t)*t^n*bessel_j(n,a*t),t));
1822 a^n*gamma(2*n+1)/(2^n*gamma(n+1)*(s^2+a^2)^((2*n+1)/2));
1823
1824 /* Problem 233: No. 04 (OK)
1825    We can show that the Maxima result is equivalent to
1826        (1/2) * 2^nu/sqrt(%pi)*gamma(nu+1/2)*a^nu*(s^2+a^2)^(-nu-1/2)
1827    That differs by a factor of 1/2 from the tabulated result of EqWorld:
1828        2^nu*%pi^((-1)/2)*gamma(1/2+nu)*a^nu*(a^2+s^2)^(-1/2-nu);
1829    I think Eqworld is wrong. See also Problem: 241
1830 */
1831
1832 radcan(specint(%e^(-s*t)*t^nu*bessel_j(nu,a*t),t));
1833 a^nu*gamma(2*nu+1)/(2^nu*gamma(nu+1)*(s^2+a^2)^((2*nu+1)/2));
1834
1835 /* Problem 234: No. 05 OK
1836    For this case we can verify the Maxima result. It is equivalent to
1837    the tabulated solution:
1838       2^nu*%pi^(-1/2)*gamma(nu+1/2)*a^nu*(s^2+a^2)^(-nu-1/2)
1839 */
1840
1841 radcan(specint(%e^(-s*t)*t^(nu+1)*bessel_j(nu,a*t),t));
1842 a^nu*gamma(2*nu+2)*s/((s^2+a^2)^((2*nu+1)/2)*(2^nu*gamma(nu+1)*s^2+a^2*2^nu*gamma(nu+1)));
1843
1844 /* Problem 235: No. 06 OK */
1845
1846 radcan(specint(%e^(-s*t)*bessel_j(0,2*sqrt(a*t)),t));
1847 1/s*%e^(-a/s);
1848
1849 /* Problem 236: No. 07 OK */
1850
1851 radcan(specint(%e^(-s*t)*sqrt(t)*bessel_j(1,2*sqrt(a*t)),t));
1852 sqrt(a)/s^2*%e^(-a/s);
1853
1854 /* Problem 237: No. 08 OK */
1855
1856 radcan(specint(%e^(-s*t)*t^(nu/2)*bessel_j(nu,2*sqrt(a*t)),t));
1857 a^(nu/2)*s^(-nu-1)*%e^(-a/s);
1858
1859 /* Problem 238: No. 09 (OK noun form New)
1860    Maxima gives the result:
1861       1/sqrt(s^2+a^2*b+a^2)
1862    But the answer of EqWorld contains an exponential function.
1863    What's wrong?
1864    We have a bug in the routine dionarghyp and following. After transformation
1865    to a hypergeometric function we have an argument of the form (t^2+b*t).
1866    Maxima doesn't recognize the two terms of the argument correctly.
1867    After specialization of the pattern and correcting the return
1868    value of the routine execargmatch and execf19 we get a correct noun.
1869    1/sqrt(s^2+a^2)*%e^(b*s-b*sqrt(s^2+a^2));
1870 */
1871
1872 radcan(specint(%e^(-s*t)*bessel_j(0,a*sqrt(t^2+b*t)),t));
1873 'specint(bessel_j(0,a*sqrt(t)*sqrt(t+b))*%e^(-s*t),t);
1874
1875 /* Problem 239: No. 10 OK */
1876
1877 radcan(specint(%e^(-s*t)*bessel_i(0,a*t),t));
1878 1/(sqrt(s-a)*sqrt(s+a));
1879
1880 /* Problem 240: No. 11 OK New*/
1881
1882 radcan(specint(%e^(-s*t)*bessel_i(nu,a*t),t));
1883 a^nu/((sqrt(s-a)*sqrt(s+a))*(sqrt(s-a)*sqrt(s+a)+s)^nu);
1884
1885 /* Problem 241: No. 12 OK New
1886    But the result differs by a factor 1/2. See problem: 233.
1887    What's wrong.
1888    Tabulated result: 2^nu*%pi^(-1/2)*gamma(nu+1/2)*a^nu*(s^2-a^2)^(-nu-1/2);
1889    The tabuluated result of EqWorld is wrong. The first factor should be
1890    2^(nu-1). We have verified the result with the help of A&S 29.3.60.
1891    I think EqWorld is wrong.
1892 */
1893
1894 radcan(specint(%e^(-s*t)*t^nu*bessel_i(nu,a*t),t));
1895 a^nu*gamma(2*nu+1)/(2^nu*gamma(nu+1)*(s-a)^((2*nu+1)/2)*(s+a)^((2*nu+1)/2));
1896
1897 /* Problem 242: No. 13 OK New
1898    The result is equivalent to the tabulated solution:
1899       2^(nu+1)*%pi^(1/2)*gamma(nu+3/2)*a^nu*s*(s^2-a^2)^(-nu-3/2);
1900 */
1901
1902 radcan(specint(%e^(-s*t)*t^(nu+1)*bessel_i(nu,a*t),t));
1903 a^nu*gamma(2*nu+2)*s/((s-a)^((2*nu+1)/2)*(s+a)^((2*nu+1)/2)*(2^nu*gamma(nu+1)*s^2-a^2*2^nu*gamma(nu+1)));
1904
1905 /* Problem 243: No. 14 OK */
1906
1907 radcan(specint(%e^(-s*t)*bessel_i(0,2*sqrt(a*t)),t));
1908 1/s*%e^(a/s);
1909
1910 /* Problem 244: No. 15 OK */
1911
1912 block([besselexpand:true],expand(exponentialize(radcan(specint(%e^(-s*t)*bessel_i(1,2*sqrt(a*t))/sqrt(t),t)))));
1913 %e^(a/s)/sqrt(a)-1/sqrt(a);
1914
1915 /* Problem 245: No. 16 OK New*/
1916
1917 radcan(specint(%e^(-s*t)*t^(nu/2)*bessel_i(nu,2*sqrt(a*t)),t));
1918 a^(nu/2)*s^(-nu-1)*%e^(a/s);
1919
1920 /* Problem 246: No. 17 OK
1921    The tabulated result -2/%pi*asinh(s/a)/sqrt(s^2+a^2) is equivalent to
1922    the result of Maxima.
1923 */
1924
1925 specint(%e^(-s*t)*bessel_y(0,a*t),t);
1926  -log((sqrt(s^2+a^2)+s)/(sqrt(s^2+a^2)-s))/(%pi*sqrt(s^2+a^2));
1927
1928 /* Problem 247: No. 18 OK New (10/2009)
1929  * A special algorithm for bessel_k(0,a*t) has been added.
1930  */
1931
1932 specint(%e^(-s*t)*bessel_k(0,a*t),t);
1933 acosh(s/a)/sqrt(s^2-a^2);
1934
1935 /*******************************************************************************
1936  *
1937  * Bug ID: 2858045 - Laplace transform of asin or atan gives wrong noun form
1938  * Check the correct noun for the asin and atan function.
1939  */
1940
1941 specint(%e^(-s*t)*asin(t),t);
1942 'specint(%e^(-s*t)*asin(t),t);
1943
1944 specint(%e^(-s*t)*atan(t),t);
1945 'specint(%e^(-s*t)*atan(t),t);
1946
1947 /* The results for an argument 1/sqrt(t) have to be checked.
1948  *
1949  * The following results are wrong. The algorithm to integrate the
1950  * hypergeometric function only works for an argument t^n where n is positive
1951  * integer. We have asin(z) = z * F(1/2,1/2; 3/2; z^2), therefore the argument
1952  * of the hypergeometric function is 1/t.
1953  */
1954
1955 specint(%e^(-s*t)*asin(1/sqrt(t)),t);
1956 /*sqrt(%pi)*%f[3,1]([1/2,1/2,-1/2],[3/2],-s)/sqrt(s);*/
1957 'specint(%e^(-s*t)*asin(1/sqrt(t)),t);
1958
1959 specint(%e^(-s*t)*atan(1/sqrt(t)),t);
1960 /*sqrt(%pi)*%f[3,1]([1/2,1,-1/2],[3/2],s)/sqrt(s);*/
1961 'specint(%e^(-s*t)*atan(1/sqrt(t)),t);
1962 /*******************************************************************************
1963  *
1964  * Bug ID: 2864588 - specint(exp(-s*t)*bessel_i(1/2,t)^2,t) not correct
1965  *
1966  * The Laplace transform of the square of the bessel_i(1/2,t) function is equal
1967  * to the Laplace transform of the expanded form 2*sinh(t)^2/(%pi*t).
1968  */
1969
1970 specint(exp(-s*t)*bessel_i(1/2,t)^2,t),besselexpand:false;
1971  -log(1-4/s^2)/(2*%pi);
1972
1973 specint(exp(-s*t)*bessel_i(1/2,t)^2,t),besselexpand:true;
1974  -log(1-4/s^2)/(2*%pi);
1975
1976 /*******************************************************************************
1977  *
1978  * Bug ID: 2866802 - specint(exp(-s*t)*t^(5/2)*bessel_j(-1/2,sqrt(t))^2,t) wrong
1979  *
1980  * We add two tests to show that the square of bessel_j(v,t) works for
1981  * v=-1/2 (the originally bug) and v=1/2 as expected.
1982  */
1983
1984 res:factor(ratsimp(specint(exp(-s*t)*t*bessel_j(-1/2,t)^2,t))),
1985     besselexpand:false;
1986 2*(s^2+2)/(%pi*s*(s^2+4));
1987
1988 expand(res-ratsimp(specint(exp(-s*t)*t*bessel_j(-1/2,t)^2,t))),
1989    besselexpand:true;
1990 0;
1991
1992 res:factor(ratsimp(specint(exp(-s*t)*t*bessel_j(1/2,t)^2,t))),
1993     besselexpand:false;
1994 4/(%pi*s*(s^2+4));
1995
1996 expand(res-ratsimp(specint(exp(-s*t)*t*bessel_j(1/2,t)^2,t))),
1997    besselexpand:true;
1998 0;
1999
2000 /*******************************************************************************
2001  *
2002  * BUG ID: 2867499 -specint(exp(-s*t)*CONST*t^2*bessel_y(1,t),t) wrong
2003  *
2004  * Show that a constant term is not lost in the Laplace transform
2005  */
2006
2007 factor(ratsimp(specint(exp(-s*t)*CONST*t^2*bessel_y(1,t),t)));
2008 (3*s*log((sqrt(s^2+1)+s)/(sqrt(s^2+1)-s))
2009         +2*s^2*sqrt(s^2+1)-4*sqrt(s^2+1))
2010         *CONST
2011         /(%pi*(s^2+1)^(5/2));
2012
2013 /*******************************************************************************
2014  *
2015  * Bug ID: 2867434 - specint(exp(-s*t)*t^2*%h[1,1](t),t) does not work
2016  *
2017  * Show that the Laplace transform of Hankel functions work for an integer order
2018  */
2019
2020 res1:factor(ratsimp(specint(exp(-s*t)*t^2*hankel_1(1,t),t)));
2021 (3*%i*s*sqrt(s^2+1)*log((sqrt(s^2+1)+s)/(sqrt(s^2+1)-s))
2022         +3*%pi*s*sqrt(s^2+1)+2*%i*s^4-2*%i*s^2-4*%i)
2023         /(%pi*(s^2+1)^3);
2024
2025 res2:factor(ratsimp(specint(exp(-s*t)*t^2*hankel_2(1,t),t)));
2026 (3*%i*s*sqrt(s^2+1)*log((sqrt(s^2+1)+s)/(sqrt(s^2+1)-s))
2027  -3*%pi*s*sqrt(s^2+1)+2*%i*s^4-2*%i*s^2-4*%i)
2028  /(-%pi*(s^2+1)^3)$
2029
2030
2031 /* Proof that the results are equivalent to the Laplace transforms of
2032    bessel_j and bessel_y.
2033    */
2034
2035 ratsimp(specint(exp(-s*t)*t^2*bessel_j(1,t),t) - 1/2*(res1+res2));
2036 0;
2037 ratsimp(specint(exp(-s*t)*t^2*bessel_y(1,t),t) - 1/(2*%i)*(res1-res2));
2038 0;
2039
2040 /*******************************************************************************
2041  *
2042  * We know the following Laplace transform of the Parabolic Cylinder D function:
2043  * (functions.wolfram.com)
2044  *
2045  *   t^(a-1)*D[v](sqrt(t))
2046  *      ->   2^(-a+v/2+1)*sqrt(%pi)*gamma(2*a)/gamma(a+(1-v)/2)
2047  *              *2F1([a,a+1/2],[a+(1-v)/2],1/2-2*t)
2048  *
2049  * Check this formula for some values.
2050  */
2051 (assume(4*?psey-1>0, 4*s-1>0),done);
2052 done;
2053
2054 (lt(a,v,z):=2^(-a+v/2+1)*sqrt(%pi)*gamma(2*a)/gamma(a+(1-v)/2)
2055                *hgfred([a,a+1/2],[a+(1-v)/2],1/2-2*z),done);
2056 done;
2057
2058 res:factor(ratsimp(specint(exp(-s*t)*parabolic_cylinder_d(0,sqrt(t)),t)));
2059 4/(4*s+1);
2060
2061 res-factor(ratsimp(lt(1,0,s)));
2062 0;
2063
2064 res:factor(ratsimp(specint(exp(-s*t)*parabolic_cylinder_d(1,sqrt(t)),t)));
2065 4*sqrt(%pi)/(4*s+1)^(3/2);
2066
2067 res-factor(ratsimp(lt(1,1,s)));
2068 0;
2069
2070 res:factor(ratsimp(specint(exp(-s*t)*parabolic_cylinder_d(2,sqrt(t)),t)));
2071 -4*(4*s-3)/(4*s+1)^2;
2072
2073 res-factor(ratsimp(lt(1,2,s)));
2074 0;
2075
2076 res:factor(ratsimp(specint(exp(-s*t)*t*parabolic_cylinder_d(2,sqrt(t)),t)));
2077 -16*(4*s-7)/(4*s+1)^3;
2078
2079 res-factor(ratsimp(lt(2,2,s)));
2080 0;
2081
2082 res:factor(ratsimp(specint(exp(-s*t)*t^2*parabolic_cylinder_d(2,sqrt(t)),t)));
2083 -128*(4*s-11)/(4*s+1)^4;
2084
2085 ratsimp(res-factor(ratsimp(lt(3,2,s))));
2086 0;
2087
2088 res:factor(ratsimp(specint(exp(-s*t)*parabolic_cylinder_d(1/2,sqrt(t)),t)));
2089 2^(3/2)*sqrt(%pi)*assoc_legendre_p(3/4,-1/4,sqrt(2)/sqrt(4*s+1))
2090         /((4*s-1)^(1/8)*(4*s+1)^(7/8));
2091
2092 ratsimp(res-factor(ratsimp(lt(1,1/2,s))));
2093 0;
2094
2095 /*******************************************************************************
2096  *
2097  * BUG ID: 2867727 - specint: wrong result for Parabolic Cylinder D function
2098  *
2099  * The Laplace transform for the Parabolic Cylinder D function should only
2100  * work for an argument a*sqrt(t). Check this.
2101  *
2102  * But the Parabolic Cylinder D function has simple expansions for an
2103  * integer order and argument a*t:
2104  *    parabolic_cylinder_d(0,z):=exp(-z^2/4),
2105  *    parabolic_cylinder_d(1,z):=z*exp(-z^2/4)
2106  * We show the results for later improvement of the code.
2107  */
2108
2109 factor(ratsimp(specint(exp(-s*t)*exp(-t^2/4),t)));
2110 -sqrt(%pi)*%e^s^2*(erf(s)-1);
2111
2112 factor(ratsimp(specint(exp(-s*t)*parabolic_cylinder_d(0,t),t)));
2113 'specint(exp(-s*t)*parabolic_cylinder_d(0,t),t);
2114
2115 factor(ratsimp(specint(exp(-s*t)*t*exp(-t^2/4),t)));
2116 2*(sqrt(%pi)*s*%e^s^2*erf(s)-sqrt(%pi)*s*%e^s^2+1);
2117
2118 factor(ratsimp(specint(exp(-s*t)*parabolic_cylinder_d(1,t),t)));
2119 'specint(exp(-s*t)*parabolic_cylinder_d(1,t),t);
2120
2121 /*******************************************************************************
2122  *
2123  * Test Laplace transform of whittaker_m
2124  *
2125  *  Check it against two simple special cases of whittaker_m:
2126  *
2127  *    whittaker_m(-1/2-u,u,z) = exp(z/2) * z^(u+1/2)
2128  *    whittaker_m(1/2 +u,u,z) = exp(-z/2)* z^(u+1/2)
2129  */
2130
2131 (assume(s>0,2*u+3>0,2*s-1>0),done);
2132 done;
2133
2134 res1:factor(ratsimp(specint(exp(-s*t)*whittaker_m(-1/2-u,u,t),t)));
2135 (2*s-1)^(-u-3/2)*2^(u+3/2)*gamma((2*u+3)/2);
2136
2137 res2:factor(ratsimp(specint(exp(-s*t)*exp(t/2)*t^(u+1/2),t)));
2138 (2*s-1)^(-u-3/2)*2^(u+3/2)*gamma((2*u+3)/2);
2139
2140 ratsimp(res1-res2);
2141 0;
2142
2143 res1:factor(ratsimp(specint(exp(-s*t)*whittaker_m(1/2+u,u,t),t)));
2144 (2*s+1)^(-u-3/2)*2^(u+3/2)*gamma((2*u+3)/2);
2145
2146 res2:factor(ratsimp(specint(exp(-s*t)*exp(-t/2)*t^(u+1/2),t)));
2147 (2*s+1)^(-u-3/2)*2^(u+3/2)*gamma((2*u+3)/2);
2148
2149 ratsimp(res1-res2);
2150 0;
2151
2152 /* The general solution for an argument t, with 2*a+2*u+1>0 */
2153
2154 (assume(2*a+2*u+1>0),done);
2155 done;
2156
2157 factor(ratsimp(specint(exp(-s*t)*t^(a-1)*whittaker_m(v,u,t),t)));
2158 (2*s+1)^(-u-a-1/2)*2^(u+a+1/2)*gamma((2*u+2*a+1)/2)
2159                          *%f[2,1]([(2*v-2*u-1)/-2,(2*u+2*a+1)/2],[2*u+1],
2160                                   2/(2*s+1));
2161
2162 /* The general solution for an argument -t, with 2*a+2*u+1>0 */
2163
2164 factor(ratsimp(specint(exp(-s*t)*t^(a-1)*whittaker_m(v,u,-t),t)));
2165 %i*(2*s-1)^(-u-a-1/2)*(-1)^u*2^(u+a+1/2)*gamma((2*u+2*a+1)/2)
2166          *%f[2,1]([(2*v-2*u-1)/-2,(2*u+2*a+1)/2],[2*u+1],-2/(2*s-1));
2167
2168 /******************************************************************************/
2169
2170 kill(all);
2171 done;