1;;; -*- Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;; 2;;; 3;;; ** (c) Copyright 1976, 1983 Massachusetts Institute of Technology ** 4;;; 5;;; 6;;; These are the main routines for finding the Laplace Transform 7;;; of special functions --- written by Yannis Avgoustis 8;;; --- modified by Edward Lafferty 9;;; Latest mod by jpg 8/21/81 10;;; 11;;; This program uses the programs on ELL;HYP FASL. 12;;; 13;;; References: 14;;; 15;;; Definite integration using the generalized hypergeometric functions 16;;; Avgoustis, Ioannis Dimitrios 17;;; Thesis. 1977. M.S.--Massachusetts Institute of Technology. Dept. 18;;; of Electrical Engineering and Computer Science 19;;; http://dspace.mit.edu/handle/1721.1/16269 20;;; 21;;; Avgoustis, I. D., Symbolic Laplace Transforms of Special Functions, 22;;; Proceedings of the 1977 MACSYMA Users' Conference, pp 21-41 23 24(in-package :maxima) 25 26(macsyma-module hypgeo) 27 28(declare-top (special checkcoefsignlist $exponentialize $radexpand $logexpand 29 $expintrep)) 30 31(defmvar $prefer_d nil 32 "When NIL express a parabolic cylinder function as hypergeometric function.") 33 34;; The properties NOUN and VERB give correct linear display 35(defprop $specint %specint verb) 36(defprop %specint $specint noun) 37 38(defvar *hyp-return-noun-form-p* t 39 "Return noun form instead of internal Lisp symbols for integrals 40 that we don't support.") 41 42(defvar *hyp-return-noun-flag* nil 43 "Flag to signal that we have no result and to return a noun.") 44 45(defvar *debug-hypgeo* nil 46 "Print debug information if enabled.") 47 48;; The variables *var* and *par* are global to this file only. 49;; They are initialized in the routine defexec. The values are never changed. 50;; These globals are introduced to avoid passing the values of *par* and *var* 51;; through all functions of this file. 52 53(defvar *var* nil 54 "Variable of integration of Laplace transform.") 55(defvar *par* nil 56 "Parameter of Laplace transform.") 57 58;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 59 60;;; Helper function for this file 61 62(defun substl (p1 p2 p3) 63 (cond ((eq p1 p2) p3) 64 (t (maxima-substitute p1 p2 p3)))) 65 66;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 67 68;;; Test functions for pattern match, which use the globals var and *par* 69 70(defun parp (a) 71 (eq a *par*)) 72 73(defun freevar0 (m) 74 (cond ((equal m 0) nil) 75 (t (freevar m)))) 76 77;;; Test functions which do not depend on globals 78 79;; Test if EXP is 1 or %e. 80(defun expor1p (expr) 81 (or (equal expr 1) 82 (eq expr '$%e))) 83 84(defun has (expr x) 85 (not (free expr x))) 86 87(defun free-not-zero-p (expr x) 88 (and (not (equal expr 0)) (free expr x))) 89 90(defun free2 (expr x y) 91 (and (free expr x) (free expr y))) 92 93(defun has-not-alike1-p (expr x) 94 (and (not (alike1 expr x)) (has expr x))) 95 96;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 97 98;;; Some shortcuts for special functions. 99 100;; Lommel's little s[u,v](z) function. 101(defun littleslommel (m n z) 102 (list '(mqapply simp) (list '($%s simp array) m n) z)) 103 104;; Whittaker's M function 105(defun mwhit (a i1 i2) 106 (list '(mqapply simp) (list '($%m simp array) i1 i2) a)) 107 108;; Whittaker's W function 109(defun wwhit (a i1 i2) 110 (list '(mqapply simp) (list '($%w simp array) i1 i2) a)) 111 112;; Jacobi P 113(defun pjac (x n a b) 114 (list '(mqapply simp) (list '($%p simp array) n a b) x)) 115 116;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 117 118;;; Two general pattern for the routine lt-sf-log. 119 120;; Recognize c*u^v + a and a=0. 121(defun m2-arbpow1 (expr var) 122 (m2 expr 123 `((mplus) 124 ((coeffpt) 125 (c free ,var) ; more special to ensure that c is constant 126 ((mexpt) (u has ,var) (v free ,var))) 127 ((coeffpp) (a zerp))))) 128 129;; Recognize c*u^v*(a+b*u)^w+d and d=0. This is a generalization of arbpow1. 130(defun m2-arbpow2 (expr var) 131 (m2 expr 132 `((mplus) 133 ((mtimes) 134 ((coefftt) (c free ,var)) 135 ((mexpt) (u equal ,var) (v free ,var)) 136 ((mexpt) 137 ((mplus) (a free ,var) ((coefft) (b free ,var) (u equal ,var))) 138 (w free-not-zero-p ,var))) 139 ((coeffpp) (d zerp))))) 140 141;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 142 143;;; The pattern to match special functions in the routine lt-sf-log. 144 145;; Recognize asin(w) 146(defun m2-asin (expr var) 147 (m2 expr 148 `((mplus) 149 ((coeffpt) (u nonzerp) 150 ((%asin) (w has ,var))) 151 ((coeffpp) (a equal 0))))) 152 153;; Recognize atan(w) 154(defun m2-atan (expr var) 155 (m2 expr 156 `((mplus) 157 ((coeffpt) (u nonzerp) 158 ((%atan) (w has ,var))) 159 ((coeffpp) (a equal 0))))) 160 161;; Recognize bessel_j(v,w) 162(defun m2-onej (expr var) 163 (m2 expr 164 `((mplus) 165 ((coeffpt) 166 (u nonzerp) 167 ((%bessel_j) (v free ,var) (w has ,var))) 168 ((coeffpp) (a zerp))))) 169 170;; Recognize bessel_j(v1,w1)*bessel_j(v2,w2) 171(defun m2-twoj (expr var) 172 (m2 expr 173 `((mplus) 174 ((coeffpt) 175 (u nonzerp) 176 ((%bessel_j) (v1 free ,var) (w1 has ,var)) 177 ((%bessel_j) (v2 free ,var) (w2 has ,var))) 178 ((coeffpp) (a zerp))))) 179 180;; Recognize bessel_y(v1,w1)*bessel_y(v2,w2) 181(defun m2-twoy (expr var) 182 (m2 expr 183 `((mplus) 184 ((coeffpt) 185 (u nonzerp) 186 ((%bessel_y) (v1 free ,var) (w1 has ,var)) 187 ((%bessel_y) (v2 free ,var) (w2 has ,var))) 188 ((coeffpp) (a zerp))))) 189 190;; Recognize bessel_k(v1,w1)*bessel_k(v2,w2) 191(defun m2-twok (expr var) 192 (m2 expr 193 `((mplus) 194 ((coeffpt) 195 (u nonzerp) 196 ((%bessel_k) (v1 free ,var) (w1 has ,var)) 197 ((%bessel_k) (v2 free ,var) (w2 has ,var))) 198 ((coeffpp) (a zerp))))) 199 200;; Recognize bessel_k(v1,w1)*bessel_y(v2,w2) 201(defun m2-onekoney (expr var) 202 (m2 expr 203 `((mplus) 204 ((coeffpt) 205 (u nonzerp) 206 ((%bessel_k) (v1 free ,var) (w1 has ,var)) 207 ((%bessel_y) (v2 free ,var) (w2 has ,var))) 208 ((coeffpp) (a zerp))))) 209 210;; Recognize bessel_j(v,w)^2 211(defun m2-onej^2 (expr var) 212 (m2 expr 213 `((mplus) 214 ((coeffpt) 215 (u nonzerp) 216 ((mexpt) 217 ((%bessel_j) (v free ,var) (w has ,var)) 218 2)) 219 ((coeffpp) (a zerp))))) 220 221;; Recognize bessel_y(v,w)^2 222(defun m2-oney^2 (expr var) 223 (m2 expr 224 `((mplus) 225 ((coeffpt) 226 (u nonzerp) 227 ((mexpt) 228 ((%bessel_y) (v free ,var) (w has ,var)) 229 2)) 230 ((coeffpp) (a zerp))))) 231 232;; Recognize bessel_k(v,w)^2 233(defun m2-onek^2 (expr var) 234 (m2 expr 235 `((mplus) 236 ((coeffpt) 237 (u nonzerp) 238 ((mexpt) 239 ((%bessel_k) (v free ,var) (w has ,var)) 240 2)) 241 ((coeffpp) (a zerp))))) 242 243;; Recognize bessel_i(v,w) 244(defun m2-onei (expr var) 245 (m2 expr 246 `((mplus) 247 ((coeffpt) 248 (u nonzerp) 249 ((%bessel_i) (v free ,var) (w has ,var))) 250 ((coeffpp) (a zerp))))) 251 252;; Recognize bessel_i(v1,w1)*bessel_i(v2,w2) 253(defun m2-twoi (expr var) 254 (m2 expr 255 `((mplus) 256 ((coeffpt) 257 (u nonzerp) 258 ((%bessel_i) (v1 free ,var) (w1 has ,var)) 259 ((%bessel_i) (v2 free ,var) (w2 has ,var))) 260 ((coeffpp) (a zerp))))) 261 262;; Recognize hankel_1(v1,w1)*hankel_1(v2,w2), product of 2 Hankel 1 functions. 263(defun m2-two-hankel_1 (expr var) 264 (m2 expr 265 `((mplus) 266 ((coeffpt) 267 (u nonzerp) 268 ((%hankel_1) (v1 free ,var) (w1 has ,var)) 269 ((%hankel_1) (v2 free ,var) (w2 has ,var))) 270 ((coeffpp) (a zerp))))) 271 272;; Recognize hankel_2(v1,w1)*hankel_2(v2,w2), product of 2 Hankel 2 functions. 273(defun m2-two-hankel_2 (expr var) 274 (m2 expr 275 `((mplus) 276 ((coeffpt) 277 (u nonzerp) 278 ((%hankel_2) (v1 free ,var) (w1 has ,var)) 279 ((%hankel_2) (v2 free ,var) (w2 has ,var))) 280 ((coeffpp) (a zerp))))) 281 282;; Recognize hankel_1(v1,w1)*hankel_2(v2,w2), product of 2 Hankel functions. 283(defun m2-hankel_1*hankel_2 (expr var) 284 (m2 expr 285 `((mplus) 286 ((coeffpt) 287 (u nonzerp) 288 ((%hankel_1) (v1 free ,var) (w1 has ,var)) 289 ((%hankel_2) (v2 free ,var) (w2 has ,var))) 290 ((coeffpp) (a zerp))))) 291 292;; Recognize bessel_y(v1,w1)*bessel_j(v2,w2) 293(defun m2-oneyonej (expr var) 294 (m2 expr 295 `((mplus) 296 ((coeffpt) 297 (u nonzerp) 298 ((%bessel_y) (v1 free ,var) (w1 has ,var)) 299 ((%bessel_j) (v2 free ,var) (w2 has ,var))) 300 ((coeffpp) (a zerp))))) 301 302;; Recognize bessel_k(v1,w1)*bessel_j(v2,w2) 303(defun m2-onekonej (expr var) 304 (m2 expr 305 `((mplus) 306 ((coeffpt) 307 (u nonzerp) 308 ((%bessel_k) (v1 free ,var) (w1 has ,var)) 309 ((%bessel_j) (v2 free ,var) (w2 has ,var))) 310 ((coeffpp) (a zerp))))) 311 312;; Recognize bessel_y(v1,w1)*hankel_1(v2,w2) 313(defun m2-bessel_y*hankel_1 (expr var) 314 (m2 expr 315 `((mplus) 316 ((coeffpt) 317 (u nonzerp) 318 ((%bessel_y) (v1 free ,var) (w1 has ,var)) 319 ((%hankel_1) (v2 free ,var) (w2 has ,var))) 320 ((coeffpp) (a zerp))))) 321 322;; Recognize bessel_y(v1,w1)*hankel_2(v2,w2) 323(defun m2-bessel_y*hankel_2 (expr var) 324 (m2 expr 325 `((mplus) 326 ((coeffpt) 327 (u nonzerp) 328 ((%bessel_y) (v1 free ,var) (w1 has ,var)) 329 ((%hankel_2) (v2 free ,var) (w2 has ,var))) 330 ((coeffpp) (a zerp))))) 331 332;; Recognize bessel_k(v1,w1)*hankel_1(v2,w2) 333(defun m2-bessel_k*hankel_1 (expr var) 334 (m2 expr 335 `((mplus) 336 ((coeffpt) 337 (u nonzerp) 338 ((%bessel_k) (v1 free ,var) (w1 has ,var)) 339 ((%hankel_1) (v1 free ,var) (w2 has ,var))) 340 ((coeffpp) (a zerp))))) 341 342;; Recognize bessel_k(v1,w1)*hankel_2(v2,w2) 343(defun m2-bessel_k*hankel_2 (expr var) 344 (m2 expr 345 `((mplus) 346 ((coeffpt) 347 (u nonzerp) 348 ((%bessel_k) (v1 free ,var) (w1 has ,var)) 349 ((%hankel_2) (v2 free ,var) (w2 has ,var))) 350 ((coeffpp) (a zerp))))) 351 352;; Recognize bessel_i(v1,w1)*bessel_j(v2,w2) 353(defun m2-oneionej (expr var) 354 (m2 expr 355 `((mplus) 356 ((coeffpt) 357 (u nonzerp) 358 ((%bessel_i) (v1 free ,var) (w1 has ,var)) 359 ((%bessel_j) (v2 free ,var) (w2 has ,var))) 360 ((coeffpp) (a zerp))))) 361 362;; Recognize bessel_i(v1,w1)*hankel_1(v2,w2) 363(defun m2-bessel_i*hankel_1 (expr var) 364 (m2 expr 365 `((mplus) 366 ((coeffpt) 367 (u nonzerp) 368 ((%bessel_i) (v1 free ,var) (w1 has ,var)) 369 ((%hankel_1) (v2 free ,var) (w2 has ,var))) 370 ((coeffpp) (a zerp))))) 371 372;; Recognize bessel_i(v1,w1)*hankel_2(v2,w2) 373(defun m2-bessel_i*hankel_2 (expr var) 374 (m2 expr 375 `((mplus) 376 ((coeffpt) 377 (u nonzerp) 378 ((%bessel_i) (v1 free ,var) (w1 has ,var)) 379 ((%hankel_2) (v2 free ,var) (w2 has ,var))) 380 ((coeffpp) (a zerp))))) 381 382;; Recognize hankel_1(v1,w1)*bessel_j(v2,w2) 383(defun m2-hankel_1*bessel_j (expr var) 384 (m2 expr 385 `((mplus) 386 ((coeffpt) 387 (u nonzerp) 388 ((%hankel_1) (v1 free ,var) (w1 has ,var)) 389 ((%bessel_j) (v2 free ,var) (w2 has ,var))) 390 ((coeffpp) (a zerp))))) 391 392;; Recognize hankel_2(v1,w1)*bessel_j(v2,w2) 393(defun m2-hankel_2*bessel_j (expr var) 394 (m2 expr 395 `((mplus) 396 ((coeffpt) 397 (u nonzerp) 398 ((%hankel_2) (v1 free ,var) (w1 has ,var)) 399 ((%bessel_j) (v2 free ,var) (w2 has ,var))) 400 ((coeffpp) (a zerp))))) 401 402;; Recognize bessel_i(v1,w1)*bessel_y(v2,w2) 403(defun m2-oneioney (expr var) 404 (m2 expr 405 `((mplus) 406 ((coeffpt) 407 (u nonzerp) 408 ((%bessel_i) (v1 free ,var) (w1 has ,var)) 409 ((%bessel_y) (v2 free ,var) (w2 has ,var))) 410 ((coeffpp) (a zerp))))) 411 412;; Recognize bessel_i(v1,w1)*bessel_k(v2,w2) 413(defun m2-oneionek (expr var) 414 (m2 expr 415 `((mplus) 416 ((coeffpt) 417 (u nonzerp) 418 ((%bessel_i) (v1 free ,var) (w1 has ,var)) 419 ((%bessel_k) (v2 free ,var) (w2 has ,var))) 420 ((coeffpp) (a zerp))))) 421 422;; Recognize bessel_i(v,w)^2 423(defun m2-onei^2 (expr var) 424 (m2 expr 425 `((mplus) 426 ((coeffpt) 427 (u nonzerp) 428 ((mexpt) 429 ((%bessel_i) (v free ,var) (w has ,var)) 430 2)) 431 ((coeffpp) (a zerp))))) 432 433;; Recognize hankel_1(v,w)^2 434(defun m2-hankel_1^2 (expr var) 435 (m2 expr 436 `((mplus) 437 ((coeffpt) 438 (u nonzerp) 439 ((mexpt) 440 ((%hankel_1) (v free ,var) (w has ,var)) 441 2)) 442 ((coeffpp) (a zerp))))) 443 444;; Recognize hankel_2(v,w)^2 445(defun m2-hankel_2^2 (expr var) 446 (m2 expr 447 `((mplus) 448 ((coeffpt) 449 (u nonzerp) 450 ((mexpt) 451 ((%hankel_2) (v free ,var) (w has ,var)) 452 2)) 453 ((coeffpp) (a zerp))))) 454 455;; Recognize bessel_y(v,w) 456(defun m2-oney (expr var) 457 (m2 expr 458 `((mplus) 459 ((coeffpt) 460 (u nonzerp) 461 ((%bessel_y) (v free ,var) (w has ,var))) 462 ((coeffpp) (a zerp))))) 463 464;; Recognize bessel_k(v,w) 465(defun m2-onek (expr var) 466 (m2 expr 467 `((mplus) 468 ((coeffpt) 469 (u nonzerp) 470 ((%bessel_k) (v free ,var) (w has ,var))) 471 ((coeffpp) (a zerp))))) 472 473;; Recognize hankel_1(v,w) 474(defun m2-hankel_1 (expr var) 475 (m2 expr 476 `((mplus) 477 ((coeffpt) 478 (u nonzerp) 479 ((%hankel_1) (v free ,var) (w has ,var))) 480 ((coeffpp) (a zerp))))) 481 482;; Recognize hankel_2(v,w) 483(defun m2-hankel_2 (expr var) 484 (m2 expr 485 `((mplus) 486 ((coeffpt) 487 (u nonzerp) 488 ((%hankel_2) (v free ,var) (w has ,var))) 489 ((coeffpp) (a zerp))))) 490 491;; Recognize log(w) 492(defun m2-onelog (expr var) 493 (m2 expr 494 `((mplus) 495 ((coeffpt) 496 (u nonzerp) 497 ((%log) (w has ,var))) 498 ((coeffpp) (a zerp))))) 499 500;; Recognize erf(w) 501(defun m2-onerf (expr var) 502 (m2 expr 503 `((mplus) 504 ((coeffpt) 505 (u nonzerp) 506 ((%erf) (w has ,var))) 507 ((coeffpp) (a zerp))))) 508 509;; Recognize erfc(w) 510(defun m2-onerfc (expr var) 511 (m2 expr 512 `((mplus) 513 ((coeffpt) 514 (u nonzerp) 515 ((%erfc) (w has ,var))) 516 ((coeffpp) (a zerp))))) 517 518;; Recognize fresnel_s(w) 519(defun m2-onefresnel_s (expr var) 520 (m2 expr 521 `((mplus) 522 ((coeffpt) 523 (u nonzerp) 524 ((%fresnel_s) (w has ,var))) 525 ((coeffpp) (a zerp))))) 526 527;; Recognize fresnel_c(w) 528(defun m2-onefresnel_c (expr var) 529 (m2 expr 530 `((mplus) 531 ((coeffpt) 532 (u nonzerp) 533 ((%fresnel_c) (w has ,var))) 534 ((coeffpp) (a zerp))))) 535 536;; Recognize expintegral_e(v,w) 537(defun m2-oneexpintegral_e (expr var) 538 (m2 expr 539 `((mplus) 540 ((coeffpt) 541 (u nonzerp) 542 ((%expintegral_e) (v free ,var) (w has ,var))) 543 ((coeffpp) (a zerp))))) 544 545;; Recognize expintegral_ei(w) 546(defun m2-oneexpintegral_ei (expr var) 547 (m2 expr 548 `((mplus) 549 ((coeffpt) 550 (u nonzerp) 551 ((%expintegral_ei) (w has ,var))) 552 ((coeffpp) (a zerp))))) 553 554;; Recognize expintegral_e1(w) 555(defun m2-oneexpintegral_e1 (expr var) 556 (m2 expr 557 `((mplus) 558 ((coeffpt) 559 (u nonzerp) 560 ((%expintegral_e1) (w has ,var))) 561 ((coeffpp) (a zerp))))) 562 563;; Recognize expintegral_si(w) 564(defun m2-oneexpintegral_si (expr var) 565 (m2 expr 566 `((mplus) 567 ((coeffpt) 568 (u nonzerp) 569 ((%expintegral_si) (w has ,var))) 570 ((coeffpp) (a zerp))))) 571 572;; Recognize expintegral_shi(w) 573(defun m2-oneexpintegral_shi (expr var) 574 (m2 expr 575 `((mplus) 576 ((coeffpt) 577 (u nonzerp) 578 ((%expintegral_shi) (w has ,var))) 579 ((coeffpp) (a zerp))))) 580 581;; Recognize expintegral_ci(w) 582(defun m2-oneexpintegral_ci (expr var) 583 (m2 expr 584 `((mplus) 585 ((coeffpt) 586 (u nonzerp) 587 ((%expintegral_ci) (w has ,var))) 588 ((coeffpp) (a zerp))))) 589 590;; Recognize expintegral_chi(w) 591(defun m2-oneexpintegral_chi (expr var) 592 (m2 expr 593 `((mplus) 594 ((coeffpt) 595 (u nonzerp) 596 ((%expintegral_chi) (w has ,var))) 597 ((coeffpp) (a zerp))))) 598 599;; Recognize kelliptic(w), (new would be elliptic_kc) 600(defun m2-onekelliptic (expr var) 601 (m2 expr 602 `((mplus) 603 ((coeffpt) 604 (u nonzerp) 605 (($kelliptic) (w has ,var))) 606 ((coeffpp) (a zerp))))) 607 608;; Recognize elliptic_kc 609(defun m2-elliptic_kc (expr var) 610 (m2 expr 611 `((mplus) 612 ((coeffpt) 613 (u nonzerp) 614 ((%elliptic_kc) (w has ,var))) 615 ((coeffpp) (a zerp))))) 616 617;; Recognize %e(w), (new would be elliptic_ec) 618(defun m2-onee (expr var) 619 (m2 expr 620 `((mplus) 621 ((coeffpt) 622 (u nonzerp) 623 (($%e) (w has ,var))) 624 ((coeffpp) (a zerp))))) 625 626;; Recognize elliptic_ec 627(defun m2-elliptic_ec (expr var) 628 (m2 expr 629 `((mplus) 630 ((coeffpt) 631 (u nonzerp) 632 ((%elliptic_ec) (w has ,var))) 633 ((coeffpp) (a zerp))))) 634 635;; Recognize gamma_incomplete(w1, w2), Incomplete Gamma function 636(defun m2-onegammaincomplete (expr var) 637 (m2 expr 638 `((mplus) 639 ((coeffpt) 640 (u nonzerp) 641 ((%gamma_incomplete) (w1 free ,var) (w2 has ,var))) 642 ((coeffpp) (a zerp))))) 643 644;; Recognize gamma_incomplete_lower(w1,w2), gamma(a)-gamma_incomplete(w1,w2) 645(defun m2-onegamma-incomplete-lower (expr var) 646 (m2 expr 647 `((mplus) 648 ((coeffpt) 649 (u nonzerp) 650 (($gamma_incomplete_lower) (w1 free ,var) (w2 has ,var))) 651 ((coeffpp) (a zerp))))) 652 653;; Recognize Struve H function. 654(defun m2-struve_h (expr var) 655 (m2 expr 656 `((mplus) 657 ((coeffpt) 658 (u nonzerp) 659 ((%struve_h) (v free ,var) (w has ,var))) 660 ((coeffpp)(a zerp))))) 661 662;; Recognize Struve L function. 663(defun m2-struve_l (expr var) 664 (m2 expr 665 `((mplus) 666 ((coeffpt) 667 (u nonzerp) 668 ((%struve_l) (v free ,var) (w has ,var))) 669 ((coeffpp) (a zerp))))) 670 671;; Recognize Lommel s[v1,v2](w) function. 672(defun m2-ones (expr var) 673 (m2 expr 674 `((mplus) 675 ((coeffpt) 676 (u nonzerp) 677 ((mqapply) 678 (($%s array) (v1 free ,var) (v2 free ,var)) (w has ,var))) 679 ((coeffpp) (a zerp))))) 680 681;; Recognize S[v1,v2](w), Lommel function 682(defun m2-oneslommel (expr var) 683 (m2 expr 684 `((mplus) 685 ((coeffpt) 686 (u nonzerp) 687 ((mqapply) 688 (($slommel array) (v1 free ,var) (v2 free ,var)) (w has ,var))) 689 ((coeffpp) (a zerp))))) 690 691;; Recognize parabolic_cylinder_d function 692(defun m2-parabolic_cylinder_d (expr var) 693 (m2 expr 694 `((mplus) 695 ((coeffpt) 696 (u nonzerp) 697 (($parabolic_cylinder_d) (v free ,var) (w has ,var))) 698 ((coeffpp) (a zerp))))) 699 700;; Recognize kbatmann(v,w), Batemann function 701(defun m2-onekbateman (expr var) 702 (m2 expr 703 `((mplus) 704 ((coeffpt) 705 (u nonzerp) 706 ((mqapply) 707 (($kbateman array) (v free ,var)) (w has ,var))) 708 ((coeffpp) (a zerp))))) 709 710;; Recognize %l[v1,v2](w), Generalized Laguerre function 711(defun m2-onel (expr var) 712 (m2 expr 713 `((mplus) 714 ((coeffpt) 715 (u nonzerp) 716 ((mqapply) 717 (($%l array) (v1 free ,var) (v2 free ,var)) (w has ,var))) 718 ((coeffpp) (a zerp))))) 719 720;; Recognize gen_laguerre(v1,v2,w), Generalized Laguerre function 721(defun m2-one-gen-laguerre (expr var) 722 (m2 expr 723 `((mplus) 724 ((coeffpt) 725 (u nonzerp) 726 ((%gen_laguerre) (v1 free ,var) (v2 free ,var) (w has ,var))) 727 ((coeffpp) (a zerp))))) 728 729;; Recognize laguerre(v1,w), Laguerre function 730(defun m2-one-laguerre (expr var) 731 (m2 expr 732 `((mplus) 733 ((coeffpt) 734 (u nonzerp) 735 ((%laguerre) (v1 free ,var) (w has ,var))) 736 ((coeffpp) (a zerp))))) 737 738;; Recognize %c[v1,v2](w), Gegenbauer function 739(defun m2-onec (expr var) 740 (m2 expr 741 `((mplus) 742 ((coeffpt) 743 (u nonzerp) 744 ((mqapply) 745 (($%c array) (v1 free ,var) (v2 free ,var)) (w has ,var))) 746 ((coeffpp) (a zerp))))) 747 748;; Recognize %t[v1](w), Chebyshev function of the first kind 749(defun m2-onet (expr var) 750 (m2 expr 751 `((mplus) 752 ((coeffpt) 753 (u nonzerp) 754 ((mqapply) (($%t array) (v1 free ,var)) (w has ,var))) 755 ((coeffpp) (a zerp))))) 756 757;; Recognize %u[v1](w), Chebyshev function of the second kind 758(defun m2-oneu (expr var) 759 (m2 expr 760 `((mplus) 761 ((coeffpt) 762 (u nonzerp) 763 ((mqapply) (($%u array) (v1 free ,var)) (w has ,var))) 764 ((coeffpp) (a zerp))))) 765 766;; Recognize %p[v1,v2,v3](w), Jacobi function 767(defun m2-onepjac (expr var) 768 (m2 expr 769 `((mplus) 770 ((coeffpt) 771 (u nonzerp) 772 ((mqapply) 773 (($%p array) 774 (v1 free ,var) (v2 free ,var) (v3 free ,var)) (w has ,var))) 775 ((coeffpp) (a zerp))))) 776 777;; Recognize jacobi_p function 778(defun m2-jacobi_p (expr var) 779 (m2 expr 780 `((mplus) 781 ((coeffpt) 782 (u nonzerp) 783 (($jacobi_p) 784 (v1 free ,var) (v2 free ,var) (v3 free ,var) (w has ,var))) 785 ((coeffpp) (a zerp))))) 786 787;; Recognize %p[v1,v2](w), Associated Legendre P function 788(defun m2-hyp-onep (expr var) 789 (m2 expr 790 `((mplus) 791 ((coeffpt) 792 (u nonzerp) 793 ((mqapply) 794 (($%p array) (v1 free ,var) (v2 free ,var)) (w has ,var))) 795 ((coeffpp) (a zerp))))) 796 797;; Recognize assoc_legendre_p function 798(defun m2-assoc_legendre_p (expr var) 799 (m2 expr 800 `((mplus) 801 ((coeffpt) 802 (u nonzerp) 803 (($assoc_legendre_p) (v1 free ,var) (v2 free ,var) (w has ,var))) 804 ((coeffpp) (a zerp))))) 805 806;; Recognize %p[v1](w), Legendre P function 807(defun m2-onep0 (expr var) 808 (m2 expr 809 `((mplus) 810 ((coeffpt) 811 (u nonzerp) 812 ((mqapply)(($%p array) (v1 free ,var)) (w has ,var))) 813 ((coeffpp) (a zerp))))) 814 815;; Recognize %p[v1](w), Legendre P function 816(defun m2-legendre_p (expr var) 817 (m2 expr 818 `((mplus) 819 ((coeffpt) 820 (u nonzerp) 821 (($legendre_p) (v free ,var)) (w has ,var)) 822 ((coeffpp) (a zerp))))) 823 824;; Recognize hermite(v1,w), Hermite function 825(defun m2-one-hermite (expr var) 826 (m2 expr 827 `((mplus) 828 ((coeffpt) 829 (u nonzerp) 830 ((%hermite) (v1 free ,var) (w has ,var))) 831 ((coeffpp) (a zerp))))) 832 833;; Recognize %q[v1,v2](w), Associated Legendre function of the second kind 834(defun m2-oneq (expr var) 835 (m2 expr 836 `((mplus) 837 ((coeffpt) 838 (u nonzerp) 839 ((mqapply) 840 (($%q array) (v1 free ,var) (v2 free ,var)) (w has ,var))) 841 ((coeffpp) (a zerp))))) 842 843;; Recognize assoc_legendre_q function 844(defun m2-assoc_legendre_q (expr var) 845 (m2 expr 846 `((mplus) 847 ((coeffpt) 848 (u nonzerp) 849 (($assoc_legendre_q) (v1 free ,var) (v2 free ,var) (w has ,var))) 850 ((coeffpp) (a zerp))))) 851 852;; Recognize %w[v1,v2](w), Whittaker W function. 853(defun m2-onew (expr var) 854 (m2 expr 855 `((mplus) 856 ((coeffpt) 857 (u nonzerp) 858 ((mqapply) 859 (($%w array) (v1 free ,var) (v2 free ,var)) (w has ,var))) 860 ((coeffpp) (a zerp))))) 861 862;; Recognize whittaker_w function. 863(defun m2-whittaker_w (expr var) 864 (m2 expr 865 `((mplus) 866 ((coeffpt) 867 (u nonzerp) 868 (($whittaker_w) (v1 free ,var) (v2 free ,var) (w has ,var))) 869 ((coeffpp) (a zerp))))) 870 871;; Recognize %m[v1,v2](w), Whittaker M function 872(defun m2-onem (expr var) 873 (m2 expr 874 `((mplus) 875 ((coeffpt) 876 (u nonzerp) 877 ((mqapply) 878 (($%m array) (v1 free ,var) (v2 free ,var)) (w has ,var))) 879 ((coeffpp) (a zerp))))) 880 881;; Recognize whittaker_m function. 882(defun m2-whittaker_m (expr var) 883 (m2 expr 884 `((mplus) 885 ((coeffpt) 886 (u nonzerp) 887 (($whittaker_m) (v1 free ,var) (v2 free ,var) (w has ,var))) 888 ((coeffpp) (a zerp))))) 889 890;; Recognize %f[v1,v2](w1,w2,w3), Hypergeometric function 891(defun m2-onef (expr var) 892 (m2 expr 893 `((mplus) 894 ((coeffpt) 895 (u nonzerp) 896 ((mqapply) 897 (($%f array) (v1 free ,var) (v2 free ,var)) 898 (w1 free ,var) 899 (w2 free ,var) 900 (w3 has ,var))) 901 ((coeffpp) (a zerp))))) 902 903;; Recognize hypergeometric function 904(defun m2-hypergeometric (expr var) 905 (m2 expr 906 `((mplus) 907 ((coeffpt) 908 (u nonzerp) 909 (($hypergeometric) (w1 free ,var) (w2 free ,var) (w3 has ,var))) 910 ((coeffpp) (a zerp))))) 911 912;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 913 914;;; Pattern for the routine hypgeo-exec. 915;;; RECOGNIZES L.T.E. "U*%E^(A*X+E*F(X)-P*X+C)+D". 916 917(defun m2-ltep (expr var par) 918 (m2 expr 919 `((mplus) 920 ((coeffpt) 921 (u nonzerp) 922 ((mexpt) 923 $%e 924 ((mplus) 925 ((coeffpt) (a free2 ,var ,par) (x alike1 ,var)) 926 ((coeffpt) (e free2 ,var ,par) (f has ,var)) 927 ((mtimes) -1 (p alike1 ,par) (x alike1 ,var)) 928 ((coeffpp) (c free2 ,var ,par))))) 929 ((coeffpp) (d equal 0))))) 930 931;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 932 933;;; Pattern for the routine defexec. 934;;; This is trying to match EXP to u*%e^(a*x+e*f+c)+d 935;;; where a, c, and e are free of x, f is free of p, and d is 0. 936 937(defun m2-defltep (expr var) 938 (m2 expr 939 `((mplus) 940 ((coeffpt) 941 (u nonzerp) 942 ((mexpt) 943 $%e 944 ((mplus) 945 ((coeffpt) (a free ,var) (x alike1 ,var)) 946 ((coeffpt) (e free ,var) (f has-not-alike1-p ,var)) 947 ((coeffpp) (c free ,var))))) 948 ((coeffpp) (d equal 0))))) 949 950;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 951 952;;; $specint is the Maxima User function 953 954(defmfun $specint (expr var) 955 (prog ($radexpand checkcoefsignlist) 956 (setq $radexpand '$all) 957 (return (defintegrate expr var)))) 958 959(defun defintegrate (expr var) 960 ;; This used to have $exponentialize enabled for everything, but I 961 ;; don't think we should do that. If various routines want 962 ;; $exponentialize, let them set it themselves. So, for here, we 963 ;; want to expand the form with $exponentialize to convert trig and 964 ;; hyperbolic functions to exponential functions that we can handle. 965 (let ((form (let (($exponentialize t)) 966 ($factor (resimplify expr))))) ; At first factor the integrand. 967 968 ;; Because we call defintegrate recursively, we add code to end the 969 ;; recursion safely. 970 971 (when (atom form) 972 (cond ((and (numberp form) (zerop form)) (return-from defintegrate 0)) 973 (t (return-from defintegrate (list '(%specint simp) form var))))) 974 975 ;; We try to find a constant denominator. This is necessary to get results 976 ;; for integrands like u(t)/(a+b+c+...). 977 978 (let ((den ($denom form))) 979 (when (and (not (equal 1 den)) ($freeof var den)) 980 (return-from defintegrate 981 (div (defintegrate (mul den form) var) den)))) 982 983 ;; We search for a sum of Exponential functions which we can integrate. 984 ;; This code finds result for Trigonometric or Hyperbolic functions with 985 ;; a factor t^-1 or t^-2 e.g. t^-1*sin(a*t). 986 987 (let* ((l (m2-defltep form var)) 988 (s (mul -1 (cdras 'a l))) 989 (u ($expand (cdras 'u l))) 990 (l1)) 991 (cond 992 ((setq l1 (m2-sum-with-exp-case1 u var)) 993 ;; c * t^-1 * (%e^(-a*t) - %e^(-b*t)) + d 994 (let ((c (cdras 'c l1)) 995 (a (mul -1 (cdras 'a l1))) 996 (b (mul -1 (cdras 'b l1))) 997 (d (cdras 'd l1))) 998 (add (mul c (take '(%log) (div (add s b) (add s a)))) 999 (defintegrate (mul d (power '$%e (mul -1 s var))) var)))) 1000 1001 ((setq l1 (m2-sum-with-exp-case2 u var)) 1002 ;; c * t^(-3/2) * (%e^(-a*t) - %e^(-b*t)) + d 1003 (let ((c (cdras 'c l1)) 1004 (a (mul -1 (cdras 'a l1))) 1005 (b (mul -1 (cdras 'b l1))) 1006 (d (cdras 'd l1))) 1007 (add (mul 2 c 1008 (power '$%pi '((rat simp) 1 2)) 1009 (sub (power (add s b) '((rat simp) 1 2)) 1010 (power (add s a) '((rat simp) 1 2)))) 1011 (defintegrate (mul d (power '$%e (mul -1 s var))) var)))) 1012 1013 ((setq l1 (m2-sum-with-exp-case3 u var)) 1014 ;; c * t^-2 * (1 - 2 * %e^(-a*t) + %e^(2*a*t)) + d 1015 (let ((c (cdras 'c l1)) 1016 (a (div (cdras 'a l1) -2)) 1017 (d (cdras 'd l1))) 1018 (add (mul c 1019 (add (mul (add s a a) (take '(%log) (add s a a))) 1020 (mul s (take '(%log) s)) 1021 (mul -2 (add s a) (take '(%log) (add s a))))) 1022 (defintegrate (mul d (power '$%e (mul -1 s var))) var)))) 1023 1024 ((setq l1 (m2-sum-with-exp-case4 u var)) 1025 ;; c * t^-1 * (1 - 2 * %e^(-a*t) + %e^(2*a*t)) + d 1026 (let ((c (cdras 'c l1)) 1027 (a (div (cdras 'a l1) (mul 4 '$%i))) 1028 (d (cdras 'd l1))) 1029 (add (mul -1 c 1030 (take '(%log) 1031 (add 1 1032 (div (mul 4 a a) 1033 (mul (sub s (mul 2 '$%i a)) 1034 (sub s (mul 2 '$%i a))))))) 1035 (defintegrate (mul d (power '$%e (mul -1 s var))) var)))) 1036 1037 ((setq l1 (m2-sum-with-exp-case5 u var)) 1038 ;; c * t^-1 * (1 - %e^(2*a*t)) + d 1039 (let ((c (cdras 'c l1)) 1040 (a (cdras 'a l1)) 1041 (d (cdras 'd l1))) 1042 (add (mul c (take '(%log) (div (sub s a) s))) 1043 (defintegrate (mul d (power '$%e (mul -1 s var))) var)))) 1044 1045 (t 1046 ;; At this point we expand the integrand. 1047 (distrdefexecinit ($expand form) var)))))) 1048 1049;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1050 1051;;; Five pattern to find sums of Exponential functions which we can integrate. 1052 1053;; Case 1: c * u^-1 * (%e^(-a*u) - %e^(-b*u)) 1054(defun m2-sum-with-exp-case1 (expr var) 1055 (m2 expr 1056 `((mplus) 1057 ((coefft) 1058 (c free ,var) 1059 ((mexpt) (u alike1 ,var) -1) 1060 ((mexpt) $%e 1061 ((mplus) 1062 ((coeffpt) (a nonzerp) (u alike1 ,var)) 1063 ((coeffpp) (z1 zerp))))) 1064 ((coefft) 1065 (c2 equal-times-minus-one c) 1066 ((mexpt) (u alike1 ,var) -1) 1067 ((mexpt) $%e 1068 ((mplus) 1069 ((coeffpt) (b nonzerp) (u alike1 ,var)) 1070 ((coeffpp) (z2 zerp))))) 1071 ((coeffpp) (d true))))) 1072 1073;; Case 2: c * u^(-3/2) * (%e^(-a*u) - %e^(-b*u)) 1074(defun m2-sum-with-exp-case2 (expr var) 1075 (m2 expr 1076 `((mplus) 1077 ((coefft) 1078 (c free ,var) 1079 ((mexpt) (u alike1 ,var) ((rat) -3 2)) 1080 ((mexpt) $%e 1081 ((mplus) 1082 ((coeffpt) (a nonzerp) (u alike1 ,var)) 1083 ((coeffpp) (z1 zerp))))) 1084 ((coefft) 1085 (c2 equal-times-minus-one c) 1086 ((mexpt) (u alike1 ,var) ((rat) -3 2)) 1087 ((mexpt) $%e 1088 ((mplus) 1089 ((coeffpt) (b nonzerp) (u alike1 ,var)) 1090 ((coeffpp) (z2 zerp))))) 1091 ((coeffpp) (d true))))) 1092 1093;; Case 3: c * u^-2 * (1 - 2 * %e^(-a*u) + %e^(2*a*u)) 1094(defun m2-sum-with-exp-case3 (expr var) 1095 (m2 expr 1096 `((mplus) 1097 ((coefft) 1098 (c free ,var) 1099 ((mexpt) (u alike1 ,var) -2)) 1100 ((coefft) 1101 (c2 equal c) 1102 ((mexpt) (u alike1 ,var) -2) 1103 ((mexpt) $%e 1104 ((mplus) 1105 ((coeffpt) (a nonzerp) (u alike1 ,var)) 1106 ((coeffpp) (z1 zerp))))) 1107 ((coefft) 1108 (c3 equal-times-minus-two c) 1109 ((mexpt) (u alike1 ,var) -2) 1110 ((mexpt) $%e 1111 ((mplus) 1112 ((coeffpt) (b equal-div-two a) (u alike1 ,var)) 1113 ((coeffpp) (z2 zerp))))) 1114 ((coeffpp) (d true))))) 1115 1116;; Case 4: c * t^-1 * (1 - 2 * %e^(-a*t) + %e^(2*a*t)) 1117(defun m2-sum-with-exp-case4 (expr var) 1118 (m2 expr 1119 `((mplus) 1120 ((coefft) 1121 (c free ,var) 1122 ((mexpt) (u alike1 ,var) -1)) 1123 ((coefft) 1124 (c2 equal c) 1125 ((mexpt) (u alike1 ,var) -1) 1126 ((mexpt) $%e 1127 ((mplus) 1128 ((coeffpt) (a nonzerp) (u alike1 ,var)) 1129 ((coeffpp) (z1 zerp))))) 1130 ((coefft) 1131 (c3 equal-times-minus-two c) 1132 ((mexpt) (u alike1 ,var) -1) 1133 ((mexpt) $%e 1134 ((mplus) 1135 ((coeffpt) (b equal-div-two a) (u alike1 ,var)) 1136 ((coeffpp) (z2 zerp))))) 1137 ((coeffpp) (d true))))) 1138 1139;; Case 5: c* t^-1 * (1 - %e^(2*a*t)) 1140(defun m2-sum-with-exp-case5 (expr var) 1141 (m2 expr 1142 `((mplus) 1143 ((coefft) 1144 (c free ,var) 1145 ((mexpt) (u alike1 ,var) -1)) 1146 ((coefft) 1147 (c2 equal-times-minus-one c) 1148 ((mexpt) (u alike1 ,var) -1) 1149 ((mexpt) $%e 1150 ((mplus) 1151 ((coeffpt) (a nonzerp) (u alike1 ,var)) 1152 ((coeffpp) (z1 zerp))))) 1153 ((coeffpp) (d true))))) 1154 1155;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1156 1157;;; Test functions for the pattern m2-sum-with-exp-case<n> 1158 1159(defun equal-times-minus-one (a b) 1160 (equal a (mul -1 b))) 1161 1162(defun equal-times-minus-two (a b) 1163 (equal a (mul -2 b))) 1164 1165(defun equal-div-two (a b) 1166 (equal a (div b 2))) 1167 1168;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1169 1170;;; Called by defintegrate. 1171;;; Call for every term of a sum defexec and add up the results. 1172 1173;; Evaluate the transform of a sum as sum of transforms. 1174(defun distrdefexecinit (expr var) 1175 (cond ((equal (caar expr) 'mplus) 1176 (distrdefexec (cdr expr) var)) 1177 (t (defexec expr var)))) 1178 1179;; FUN is a list of addends. Compute the transform of each addend and 1180;; add them up. 1181(defun distrdefexec (expr var) 1182 (cond ((null expr) 0) 1183 (t (add (defexec (car expr) var) 1184 (distrdefexec (cdr expr) var))))) 1185 1186;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1187 1188;;; Called after transformation of a integrand to a new representation. 1189;;; Evalutate the tansform of a sum as sum of transforms. 1190 1191(defun sendexec (r a) 1192 (distrexecinit ($expand (mul (init r) a)))) 1193 1194;; Compute r*exp(-p*t), where t is the variable of integration and 1195;; p is the parameter of the Laplace transform. 1196(defun init (r) 1197 (mul r (power '$%e (mul -1 *var* *par*)))) 1198 1199(defun distrexecinit (expr) 1200 (cond ((and (consp expr) 1201 (consp (car expr)) 1202 (equal (caar expr) 'mplus)) 1203 (distrexec (cdr expr))) 1204 (t (hypgeo-exec expr)))) 1205 1206(defun distrexec (expr) 1207 (cond ((null expr) 0) 1208 (t (add (hypgeo-exec (car expr)) 1209 (distrexec (cdr expr)))))) 1210 1211;; It dispatches according to the kind of transform it matches. 1212(defun hypgeo-exec (expr) 1213 (prog (l u a c e f) 1214 (cond ((setq l (m2-ltep expr *var* *par*)) 1215 (setq u (cdras 'u l) 1216 a (cdras 'a l) 1217 c (cdras 'c l) 1218 e (cdras 'e l) 1219 f (cdras 'f l)) 1220 (return (ltscale u c a e f)))) 1221 (return (setq *hyp-return-noun-flag* 'other-trans-to-follow)))) 1222 1223;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1224 1225;;; Compute transform of EXP wrt the variable of integration VAR. 1226 1227(defun defexec (expr var) 1228 (let* ((*par* 'psey) ; Set parameter of Laplace transform 1229 (*var* var) ; Set variable of integration 1230 (*hyp-return-noun-flag* nil) ; Reset the flag 1231 (form expr) 1232 (l (m2-defltep expr var)) 1233 (s (cdras 'a l))) ; Get the parameter of the Laplace transform. 1234 1235 ;; If we have not found a parameter, we try to factor the integrand. 1236 1237 (when (and (numberp s) (equal s 0)) 1238 (setq l (m2-defltep ($factor form) *var*)) 1239 (setq s (cdras 'a l))) 1240 1241 (cond (l 1242 ;; EXP is an expression of the form u*%e^(s*t+e*f+c). So s 1243 ;; is basically the parameter of the Laplace transform. 1244 (let ((result (negtest l s))) 1245 ;; At this point we construct the noun form if one of the 1246 ;; called routines set the global flag. If the global flag 1247 ;; is not set, the noun form has been already constructed. 1248 (if (and *hyp-return-noun-form-p* *hyp-return-noun-flag*) 1249 (list '(%specint simp) expr *var*) 1250 result))) 1251 (t 1252 ;; If necessary we construct the noun form. 1253 (if *hyp-return-noun-form-p* 1254 (list '(%specint simp) expr *var*) 1255 'other-defint-to-follow-defexec))))) 1256 1257;; L is the integrand of the transform, after pattern matching. S is 1258;; the parameter (p) of the transform. 1259(defun negtest (l s) 1260 (prog (u e f c) 1261 (cond ((eq ($asksign ($realpart s)) '$neg) 1262 ;; The parameter of transform must have a negative 1263 ;; realpart. Break out the integrand into its various 1264 ;; components. 1265 (setq u (cdras 'u l) 1266 e (cdras 'e l) 1267 f (cdras 'f l) 1268 c (cdras 'c l)) 1269 (when (equal e 0) (setq f 1)) 1270 ;; To compute the transform, we replace A with PSEY for 1271 ;; simplicity. After the transform is computed, replace 1272 ;; PSEY with A. 1273 ;; 1274 ;; FIXME: Sometimes maxima will ask for the sign of PSEY. 1275 ;; But that doesn't occur in the original expression, so 1276 ;; it's very confusing. What should we do? 1277 1278 ;; We know psey must be positive. psey is a substitution 1279 ;; for the paratemter a and we have checked the sign. 1280 ;; So it is the best to add a rule for the sign of psey. 1281 1282 (mfuncall '$assume `((mgreaterp) ,*par* 0)) 1283 1284 (return 1285 (prog1 1286 (maxima-substitute 1287 (mul -1 s) 1288 *par* 1289 (ltscale u c 0 e f)) 1290 1291 ;; We forget the rule after finishing the calculation. 1292 (mfuncall '$forget `((mgreaterp) ,*par* 0)))))) 1293 1294 (return 1295 (setq *hyp-return-noun-flag* 'other-defint-to-follow-negtest)))) 1296 1297;; Compute the transform of 1298;; 1299;; U * %E^(-VAR * (*PAR* - PAR0) + E*F + C) 1300(defun ltscale (u c par0 e f) 1301 (mul (power '$%e c) 1302 (substl (sub *par* par0) *par* (lt-exec u e f)))) 1303 1304;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1305 1306;;; Compute the transform of u*%e^(-p*t+e*f) 1307 1308(defun lt-exec (u e f) 1309 (let (l a) 1310 (cond ((setq l (m2-sum u *var*)) 1311 ;; We have found a summation. 1312 (mul (cdras 'c l) 1313 (take '(%sum) 1314 (sendexec 1 (cdras 'u l)) 1315 (cdras 'i l) 1316 (cdras 'l l) 1317 (cdras 'h l)))) 1318 1319 ((setq l (m2-unit_step u *var*)) 1320 ;; We have found the Unit Step function. 1321 (setq u (cdras 'u l) 1322 a (cdras 'a l)) 1323 (mul (power '$%e (mul a *par*)) 1324 (sendexec (cond (($freeof *var* u) u) 1325 (t (maxima-substitute (sub *var* a) *var* u))) 1326 1))) 1327 1328 ((equal e 0) 1329 ;; The simple case of u*%e^(-p*t) 1330 (lt-sf-log u)) 1331 ((and (not (equal e 0)) 1332 (setq l (m2-c*t^v u *var*))) 1333 ;; We have u*%e^(-p*t+e*f). Try to see if U is of the form 1334 ;; c*t^v. If so, we can handle it here. 1335 (lt-exp l e f)) 1336 (t 1337 ;; The complicated case. Remove the e*f term and move it to u. 1338 (lt-sf-log (mul u (power '$%e (mul e f)))))))) 1339 1340;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1341 1342;;; Pattern for the routine lt-exec 1343 1344;; Recognize c*sum(u,index,low,high) 1345(defun m2-sum (expr var) 1346 (m2 expr 1347 `((mplus) 1348 ((coeffpt) 1349 (c free ,var) 1350 ((%sum) (u true) (i true) (l true) (h true))) 1351 ((coeffpp) (d zerp))))) 1352 1353;; Recognize u(t)*unit_step(x-a) 1354(defun m2-unit_step (expr var) 1355 (m2 expr 1356 `((mplus) 1357 ((coeffpt) 1358 (u nonzerp) 1359 (($unit_step) ((mplus) (x alike1 ,var) ((coeffpp) (a true))))) 1360 ((coeffpp) (d zerp))))) 1361 1362;; Recognize c*t^v. 1363;; This is a duplicate of m2-arbpow1. Look if we can use it. 1364(defun m2-c*t^v (expr var) 1365 (m2 expr 1366 `((mtimes) 1367 ((coefftt) (c free ,var)) 1368 ((mexpt) (u alike1 ,var) (v free ,var))))) 1369 1370;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1371;;; 1372;;; Algorithm 1: Laplace transform of c*t^v*exp(-s*t+e*f) 1373;;; 1374;;; L contains the pattern for c*t^v. 1375;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1376 1377(defun lt-exp (l e f) 1378 (let ((c (cdras 'c l)) 1379 (v (cdras 'v l))) 1380 (cond ((m2-t^2 f *var*) 1381 (setq e (inv (mul -8 e)) v (add v 1)) 1382 (f24p146test c v e)) 1383 ((m2-sqroott f *var*) 1384 ;; We don't do the transformation at this place. Because we take the 1385 ;; square of e we lost the sign and get wrong results. 1386 ;(setq e (mul* e e (inv 4)) v (add v 1)) 1387 (f35p147test c v e)) 1388 ((m2-t^-1 f *var*) 1389 (setq e (mul -4 e) v (add v 1)) 1390 (f29p146test c v e)) ; We have to call with the constant c. 1391 ((and (equal v 0) ; We have to test for v=0 and to call 1392 (m2-e^-t f *var*)) 1393 (f36p147 c e)) ; with the constant c. 1394 ((and (equal v 0) (m2-e^t f *var*)) 1395 (f37p147 c (mul -1 e))) 1396 (t 1397 (setq *hyp-return-noun-flag* 'other-lt-exponential-to-follow))))) 1398 1399;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1400 1401;;; Pattern for the routine lt-exp 1402 1403;; Recognize t^2 1404(defun m2-t^2 (expr var) 1405 (m2 expr `((mexpt) (u alike1 ,var) 2))) 1406 1407;; Recognize sqrt(t) 1408(defun m2-sqroott (expr var) 1409 (m2 expr `((mexpt) (u alike1 ,var) ((rat) 1 2)))) 1410 1411;; Recognize t^-1 1412(defun m2-t^-1 (expr var) 1413 (m2 expr `((mexpt) (u alike1 ,var) -1))) 1414 1415;; Recognize %e^-t 1416(defun m2-e^-t (expr var) 1417 (m2 expr `((mexpt) $%e ((mtimes) -1 (u alike1 ,var))))) 1418 1419;; Recognize %e^t 1420(defun m2-e^t (expr var) 1421 (m2 expr `((mexpt) $%e (u alike1 ,var)))) 1422 1423;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1424;;; 1425;;; Algorithm 1.1: Laplace transform of c*t^v*exp(-a*t^2) 1426;;; 1427;;; Table of Integral Transforms 1428;;; 1429;;; p. 146, formula 24: 1430;;; 1431;;; t^(v-1)*exp(-t^2/8/a) 1432;;; -> gamma(v)*2^v*a^(v/2)*exp(a*p^2)*D[-v](2*p*sqrt(a)) 1433;;; 1434;;; Re(a) > 0, Re(v) > 0 1435;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1436 1437(defun f24p146test (c v a) 1438 (cond ((and (eq ($asksign a) '$pos) 1439 (eq ($asksign v) '$pos)) 1440 ;; Both a and v must be positive 1441 (f24p146 c v a)) 1442 (t 1443 (setq *hyp-return-noun-flag* 'fail-on-f24p146test)))) 1444 1445(defun f24p146 (c v a) 1446 (mul c 1447 (take '(%gamma) v) 1448 (power 2 v) 1449 (power a (div v 2)) 1450 (power '$%e (mul a *par* *par*)) 1451 (dtford (mul 2 *par* (power a '((rat simp) 1 2))) 1452 (mul -1 v)))) 1453 1454;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1455;;; 1456;;; Algorithm 1.2: Laplace transform of c*t^v*exp(-a*sqrt(t)) 1457;;; 1458;;; Table of Integral Transforms 1459;;; 1460;;; p. 147, formula 35: 1461;;; 1462;;; (2*t)^(v-1)*exp(-2*sqrt(a)*sqrt(t)) 1463;;; -> gamma(2*v)*p^(-v)*exp(a/p/2)*D[-2*v](sqrt(2*a/p)) 1464;;; 1465;;; Re(v) > 0, Re(p) > 0 1466;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1467 1468;; Check if conditions for f35p147 hold 1469(defun f35p147test (c v a) 1470 (cond ((eq ($asksign (add v 1)) '$pos) 1471 ;; v must be positive 1472 (f35p147 c v a)) 1473 (t 1474 ;; Set a global flag. When *hyp-return-noun-form-p* is T the noun 1475 ;; form will be constructed in the routine DEFEXEC. 1476 (setq *hyp-return-noun-flag* 'fail-on-f35p147test)))) 1477 1478(defun f35p147 (c v a) 1479 ;; We have not done the calculation v->v+1 and a-> a^2/4 1480 ;; and subsitute here accordingly. 1481 (let ((v (add v 1))) 1482 (mul c 1483 (take '(%gamma) (add v v)) 1484 (power 2 (sub 1 v)) ; Is this supposed to be here? 1485 (power *par* (mul -1 v)) 1486 (power '$%e (mul a a '((rat simp) 1 8) (inv *par*))) 1487 ;; We need an additional factor -1 to get the expected results. 1488 ;; What is the mathematically reason? 1489 (dtford (mul -1 a (inv (power (mul 2 *par*) '((rat simp) 1 2)))) 1490 (mul -2 v))))) 1491 1492;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1493 1494;; Express a parabolic cylinder function as either a parabolic 1495;; cylinder function or as hypergeometric function. 1496;; 1497;; Tables of Integral Transforms, p. 386 gives this definition: 1498;; 1499;; D[v](z) = 2^(v/2+1/4)*z^(-1/2)*W[v/2+1/4,1/4](z^2/2) 1500 1501(defun dtford (z v) 1502 (let ((inv4 (inv 4))) 1503 (cond ((or $prefer_d 1504 (whittindtest (add (div v 2) inv4) inv4)) 1505 ;; At this time the Parabolic Cylinder D function is not implemented 1506 ;; as a simplifying function. We call nevertheless the simplifer 1507 ;; to simplify the arguments. When we implement the function 1508 ;; The symbol has to be changed to the noun form. 1509 (take '($parabolic_cylinder_d) v z)) 1510 (t (simpdtf z v))))) 1511 1512(defun whittindtest (i1 i2) 1513 (or (maxima-integerp (add i2 i2)) 1514 (neginp (sub (sub '((rat simp) 1 2) i2) i1)) 1515 (neginp (sub (add '((rat simp) 1 2) i2) i1)))) 1516 1517;; Return T if a is a non-positive integer. 1518;; (Do we really want maxima-integerp or hyp-integerp here?) 1519(defun neginp (a) 1520 (cond ((maxima-integerp a) 1521 (or (zerop1 a) 1522 (eq ($sign a) '$neg))))) 1523 1524;; Express parabolic cylinder function as a hypergeometric function. 1525;; 1526;; A&S 19.3.1 says 1527;; 1528;; U(a,x) = D[-a-1/2](x) 1529;; 1530;; and A&S 19.12.3 gives 1531;; 1532;; U(a,+/-x) = sqrt(%pi)*2^(-1/4-a/2)*exp(-x^2/4)/gamma(3/4+a/2) 1533;; *M(a/2+1/4,1/2,x^2/2) 1534;; -/+ sqrt(%pi)*2^(1/4-a/2)*x*exp(-x^2/4)/gamma(1/4+a/2) 1535;; *M(a/2+3/4,3/2,x^2/2) 1536;; 1537;; So 1538;; 1539;; D[v](z) = U(-v-1/2,z) 1540;; = sqrt(%pi)*2^(v/2+1/2)*x*exp(-x^2/4) 1541;; *M(1/2-v/2,3/2,x^2/2)/gamma(-v/2) 1542;; + sqrt(%pi)*2^(v/2)*exp(-x^2/4)/gamma(1/2-v/2) 1543;; *M(-v/2,1/2,x^2/2) 1544 1545(defun simpdtf (z v) 1546 (let ((inv2 '((rat simp) 1 2)) 1547 (pow (power '$%e (mul z z '((rat simp) -1 4))))) 1548 (add (mul (power 2 (div (sub v 1) 2)) 1549 z 1550 -2 (power '$%pi inv2) ; gamma(-1/2) 1551 (inv (take '(%gamma) (mul v -1 inv2))) 1552 pow 1553 (hgfsimp-exec (list (sub inv2 (div v 2))) 1554 (list '((rat simp) 3 2)) 1555 (mul z z inv2))) 1556 (mul (power 2 (div v 2)) 1557 (power '$%pi inv2) ; gamma(1/2) 1558 pow 1559 (inv (take '(%gamma) (sub inv2 (mul v inv2)))) 1560 (hgfsimp-exec (list (mul v -1 inv2)) 1561 (list inv2) 1562 (mul z z inv2)))))) 1563 1564;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1565;;; 1566;;; Algorithm 1.3: Laplace transform of t^v*exp(1/t) 1567;;; 1568;;; Table of Integral Transforms 1569;;; 1570;;; p. 146, formula 29: 1571;;; 1572;;; t^(v-1)*exp(-a/t/4) 1573;;; -> 2*(a/p/4)^(v/2)*bessel_k(v, sqrt(a)*sqrt(p)) 1574;;; 1575;;; Re(a) > 0 1576;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1577 1578;; Check if conditions for f29p146test hold 1579(defun f29p146test (c v a) 1580 (cond ((eq ($asksign a) '$pos) 1581 (f29p146 c v a)) 1582 (t 1583 (setq *hyp-return-noun-flag* 'fail-on-f29p146test)))) 1584 1585(defun f29p146 (c v a) 1586 (mul 2 c 1587 (power (mul a '((rat simp) 1 4) (inv *par*)) 1588 (div v 2)) 1589 (ktfork a v))) 1590 1591;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1592 1593;; bessel_k(v, sqrt(a)*sqrt(p)) in terms of bessel_k or in terms of 1594;; hypergeometric functions. 1595;; 1596;; Choose bessel_k if the order v is an integer. (Why?) 1597 1598(defun ktfork (a v) 1599 (let ((z (power (mul a *par*) '((rat simp) 1 2)))) 1600 (cond ((maxima-integerp v) 1601 (take '(%bessel_k) v z)) 1602 (t 1603 (simpktf z v))))) 1604 1605;; Express the Bessel K function in terms of hypergeometric functions. 1606;; 1607;; K[v](z) = %pi/2*(bessel_i(-v,z)-bessel(i,z))/sin(v*%pi) 1608;; 1609;; and 1610;; 1611;; bessel_i(v,z) = (z/2)^v/gamma(v+1)*0F1(;v+1;z^2/4) 1612 1613(defun simpktf (z v) 1614 (let ((dz2 (div z 2))) 1615 (mul '$%pi 1616 '((rat simp) 1 2) 1617 (inv (take '(%sin) (mul v '$%pi))) 1618 (sub (mul (power dz2 (mul -1 v)) 1619 (inv (take '(%gamma) (sub 1 v))) 1620 (hgfsimp-exec nil 1621 (list (sub 1 v)) 1622 (mul z z '((rat simp) 1 4)))) 1623 (mul (power dz2 v) 1624 (inv (take '(%gamma) (add v 1))) 1625 (hgfsimp-exec nil 1626 (list (add v 1)) 1627 (mul z z '((rat simp) 1 4)))))))) 1628 1629;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1630;;; 1631;;; Algorithm 1.4: Laplace transform of exp(exp(-t)) 1632;;; 1633;;; Table of Integral Transforms 1634;;; 1635;;; p. 147, formula 36: 1636;;; 1637;;; exp(-a*exp(-t)) 1638;;; -> a^(-p)*gamma(p,a) 1639;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1640 1641(defun f36p147 (c a) 1642 (let ((-a (mul -1 a))) 1643 (mul c 1644 (power -a (mul -1 *par*)) 1645 `(($gamma_incomplete_lower simp) ,*par* ,-a)))) 1646 1647;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1648;;; 1649;;; Algorithm 1.5: Laplace transform of exp(exp(t)) 1650;;; 1651;;; Table of Integral Transforms 1652;;; 1653;;; p. 147, formula 36: 1654;;; 1655;;; exp(-a*exp(t)) 1656;;; -> a^(-p)*gamma_incomplete(-p,a) 1657;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1658 1659(defun f37p147 (c a) 1660 (mul c 1661 (power a *par*) 1662 (take '(%gamma_incomplete) (mul -1 *par*) a))) 1663 1664;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1665;;; 1666;;; Algorithm 2: Laplace transform of u(t)*%e^(-p*t). 1667;;; 1668;;; u contains one or more special functions. Dispatches according to the 1669;;; special functions involved in the Laplace transformable expression. 1670;;; 1671;;; We have three general types of integrands: 1672;;; 1673;;; 1. Call a function to return immediately the Laplace transform, 1674;;; e.g. call lt-arbpow, lt-arbpow2, lt-log, whittest to return the 1675;;; Laplace transform. 1676;;; 2. Call lt-ltp directly or via an "expert function on Laplace transform", 1677;;; transform the special function to a representation in terms of one 1678;;; hypergeometric function and do the integration 1679;;; e.g. for a direct call of lt-ltp asin, atan or via lt2j for 1680;;; an integrand with two bessel function. 1681;;; 3. Call fractest, fractest1, ... which transform the involved special 1682;;; function to a new representation. Send the transformed expression with 1683;;; the routine sendexec to the integrator and try to integrate the new 1684;;; representation, e.g. gamma_incomplete is first transformed to a new 1685;;; representation. 1686;;; 1687;;; The ordering of the calls to match a pattern is important. 1688;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 1689 1690(defun lt-sf-log (u) 1691 (prog (l index1 index11 index2 index21 arg1 arg2 rest) 1692 1693 ;; Laplace transform of asin(w) 1694 (cond ((setq l (m2-asin u *var*)) 1695 (setq arg1 (cdras 'w l) 1696 rest (cdras 'u l)) 1697 (return (lt-ltp 'asin rest arg1 nil)))) 1698 1699 ;; Laplace transform of atan(w) 1700 (cond ((setq l (m2-atan u *var*)) 1701 (setq arg1 (cdras 'w l) 1702 rest (cdras 'u l)) 1703 (return (lt-ltp 'atan rest arg1 nil)))) 1704 1705 ;; Laplace transform of two Bessel J functions 1706 (cond ((setq l (m2-twoj u *var*)) 1707 (setq index1 (cdras 'v1 l) 1708 index2 (cdras 'v2 l) 1709 arg1 (cdras 'w1 l) 1710 arg2 (cdras 'w2 l) 1711 rest (cdras 'u l)) 1712 (return (lt2j rest arg1 arg2 index1 index2)))) 1713 1714 ;; Laplace transform of two hankel_1 functions 1715 (cond ((setq l (m2-two-hankel_1 u *var*)) 1716 (setq index1 (cdras 'v1 l) 1717 index2 (cdras 'v2 l) 1718 arg1 (cdras 'w1 l) 1719 arg2 (cdras 'w2 l) 1720 rest (cdras 'u l)) 1721 ;; Call the code for the symbol %h. 1722 (return (fractest rest arg1 arg2 index1 1 index2 1 '2htjory)))) 1723 1724 ;; Laplace transform of two hankel_2 functions 1725 (cond ((setq l (m2-two-hankel_2 u *var*)) 1726 (setq index1 (cdras 'v1 l) 1727 index2 (cdras 'v2 l) 1728 arg1 (cdras 'w1 l) 1729 arg2 (cdras 'w2 l) 1730 rest (cdras 'u l)) 1731 ;; Call the code for the symbol %h. 1732 (return (fractest rest arg1 arg2 index1 2 index2 2 '2htjory)))) 1733 1734 ;; Laplace transform of hankel_1 * hankel_2 1735 (cond ((setq l (m2-hankel_1*hankel_2 u *var*)) 1736 (setq index1 (cdras 'v1 l) 1737 index2 (cdras 'v2 l) 1738 arg1 (cdras 'w1 l) 1739 arg2 (cdras 'w2 l) 1740 rest (cdras 'u l)) 1741 ;; Call the code for the symbol %h. 1742 (return (fractest rest arg1 arg2 index1 1 index2 2 '2htjory)))) 1743 1744 ;; Laplace transform of two Bessel Y functions 1745 (cond ((setq l (m2-twoy u *var*)) 1746 (setq index1 (cdras 'v1 l) 1747 index2 (cdras 'v2 l) 1748 arg1 (cdras 'w1 l) 1749 arg2 (cdras 'w2 l) 1750 rest (cdras 'u l)) 1751 (return (fractest rest arg1 arg2 index1 nil index2 nil '2ytj)))) 1752 1753 ;; Laplace transform of two Bessel K functions 1754 (cond ((setq l (m2-twok u *var*)) 1755 (setq index1 (cdras 'v1 l) 1756 index2 (cdras 'v2 l) 1757 arg1 (cdras 'w1 l) 1758 arg2 (cdras 'w2 l) 1759 rest (cdras 'u l)) 1760 (return (fractest rest arg1 arg2 index1 nil index2 nil '2kti)))) 1761 1762 ;; Laplace transform of Bessel K and Bessel Y functions 1763 (cond ((setq l (m2-onekoney u *var*)) 1764 (setq index1 (cdras 'v1 l) 1765 index2 (cdras 'v2 l) 1766 arg1 (cdras 'w1 l) 1767 arg2 (cdras 'w2 l) 1768 rest (cdras 'u l)) 1769 (return (fractest rest arg1 arg2 index1 nil index2 nil 'ktiytj)))) 1770 1771 ;; Laplace transform of Bessel I and Bessel J functions 1772 (cond ((setq l (m2-oneionej u *var*)) 1773 (setq index1 (cdras 'v1 l) 1774 index2 (cdras 'v2 l) 1775 index21 (cdras 'v21 l) 1776 arg1 (mul '$%i (cdras 'w1 l)) 1777 arg2 (cdras 'w2 l) 1778 rest (mul (power '$%i (neg index1)) (cdras 'u l))) 1779 (return (lt2j rest arg1 arg2 index1 index2)))) 1780 1781 ;; Laplace transform of Bessel I and Hankel 1 functions 1782 (cond ((setq l (m2-bessel_i*hankel_1 u *var*)) 1783 (setq index1 (cdras 'v1 l) 1784 index2 (cdras 'v2 l) 1785 arg1 (mul '$%i (cdras 'w1 l)) 1786 arg2 (cdras 'w2 l) 1787 rest (mul (power '$%i (neg index1)) (cdras 'u l))) 1788 (return (fractest1 rest arg1 arg2 index1 index2 1 'besshtjory)))) 1789 1790 ;; Laplace transform of Bessel I and Hankel 2 functions 1791 (cond ((setq l (m2-bessel_i*hankel_2 u *var*)) 1792 (setq index1 (cdras 'v1 l) 1793 index2 (cdras 'v2 l) 1794 arg1 (mul '$%i (cdras 'w1 l)) 1795 arg2 (cdras 'w2 l) 1796 rest (mul (power '$%i (neg index1)) (cdras 'u l))) 1797 (return (fractest1 rest arg1 arg2 index1 index2 2 'besshtjory)))) 1798 1799 ;; Laplace transform of Bessel Y and Bessel J functions 1800 (cond ((setq l (m2-oneyonej u *var*)) 1801 (setq index1 (cdras 'v1 l) 1802 index2 (cdras 'v2 l) 1803 arg1 (cdras 'w1 l) 1804 arg2 (cdras 'w2 l) 1805 rest (cdras 'u l)) 1806 (return (fractest1 rest arg2 arg1 index2 index1 nil 'bessytj)))) 1807 1808 ;; Laplace transform of Bessel K and Bessel J functions 1809 (cond ((setq l (m2-onekonej u *var*)) 1810 (setq index1 (cdras 'v1 l) 1811 index2 (cdras 'v2 l) 1812 arg1 (cdras 'w1 l) 1813 arg2 (cdras 'w2 l) 1814 rest (cdras 'u l)) 1815 (return (fractest1 rest arg2 arg1 index2 index1 nil 'besskti)))) 1816 1817 ;; Laplace transform of Hankel 1 and Bessel J functions 1818 (cond ((setq l (m2-hankel_1*bessel_j u *var*)) 1819 (setq index1 (cdras 'v1 l) 1820 index2 (cdras 'v2 l) 1821 arg1 (cdras 'w1 l) 1822 arg2 (cdras 'w2 l) 1823 rest (cdras 'u l)) 1824 (return (fractest1 rest arg2 arg1 index2 index1 1 'besshtjory)))) 1825 1826 ;; Laplace transform of Hankel 2 and Bessel J functions 1827 (cond ((setq l (m2-hankel_2*bessel_j u *var*)) 1828 (setq index1 (cdras 'v1 l) 1829 index2 (cdras 'v2 l) 1830 arg1 (cdras 'w1 l) 1831 arg2 (cdras 'w2 l) 1832 rest (cdras 'u l)) 1833 (return (fractest1 rest arg2 arg1 index2 index1 2 'besshtjory)))) 1834 1835 ;; Laplace transform of Bessel Y and Hankel 1 functions 1836 (cond ((setq l (m2-bessel_y*hankel_1 u *var*)) 1837 (setq index1 (cdras 'v1 l) 1838 index2 (cdras 'v2 l) 1839 arg1 (cdras 'w1 l) 1840 arg2 (cdras 'w2 l) 1841 rest (cdras 'u l)) 1842 (return (fractest1 rest arg2 arg1 index2 index1 1 'htjoryytj)))) 1843 1844 ;; Laplace transform of Bessel Y and Hankel 2 functions 1845 (cond ((setq l (m2-bessel_y*hankel_2 u *var*)) 1846 (setq index1 (cdras 'v1 l) 1847 index2 (cdras 'v2 l) 1848 arg1 (cdras 'w1 l) 1849 arg2 (cdras 'w2 l) 1850 rest (cdras 'u l)) 1851 (return (fractest1 rest arg2 arg1 index2 index1 2 'htjoryytj)))) 1852 1853 ;; Laplace transform of Bessel K and Hankel 1 functions 1854 (cond ((setq l (m2-bessel_k*hankel_1 u *var*)) 1855 (setq index1 (cdras 'v1 l) 1856 index2 (cdras 'v2 l) 1857 arg1 (cdras 'w1 l) 1858 arg2 (cdras 'w2 l) 1859 rest (cdras 'u l)) 1860 (return (fractest1 rest arg2 arg1 index2 index1 1 'htjorykti)))) 1861 1862 ;; Laplace transform of Bessel K and Hankel 2 functions 1863 (cond ((setq l (m2-bessel_k*hankel_2 u *var*)) 1864 (setq index1 (cdras 'v1 l) 1865 index2 (cdras 'v2 l) 1866 arg1 (cdras 'w1 l) 1867 arg2 (cdras 'w2 l) 1868 rest (cdras 'u l)) 1869 (return (fractest1 rest arg2 arg1 index2 index1 2 'htjorykti)))) 1870 1871 ;; Laplace transform of Bessel I and Bessel Y functions 1872 (cond ((setq l (m2-oneioney u *var*)) 1873 (setq index1 (cdras 'v1 l) 1874 index2 (cdras 'v2 l) 1875 arg1 (mul '$%i (cdras 'w1 l)) 1876 arg2 (cdras 'w2 l) 1877 rest (mul (power '$%i (neg index1)) (cdras 'u l))) 1878 (return (fractest1 rest arg1 arg2 index1 index2 nil 'bessytj)))) 1879 1880 ;; Laplace transform of Bessel I and Bessel K functions 1881 (cond ((setq l (m2-oneionek u *var*)) 1882 (setq index1 (cdras 'v1 l) 1883 index2 (cdras 'v2 l) 1884 arg1 (mul '$%i (cdras 'w1 l)) 1885 arg2 (cdras 'w2 l) 1886 rest (mul (power '$%i (neg index1)) (cdras 'u l))) 1887 (return (fractest1 rest arg1 arg2 index1 index2 nil 'besskti)))) 1888 1889 ;; Laplace transform of Struve H function 1890 (cond ((setq l (m2-struve_h u *var*)) 1891 (setq index1 (cdras 'v l) 1892 arg1 (cdras 'w l) 1893 rest (cdras 'u l)) 1894 (return (lt1hstruve rest arg1 index1)))) 1895 1896 ;; Laplace transform of Struve L function 1897 (cond ((setq l (m2-struve_l u *var*)) 1898 (setq index1 (cdras 'v l) 1899 arg1 (cdras 'w l) 1900 rest (cdras 'u l)) 1901 (return (lt1lstruve rest arg1 index1)))) 1902 1903 ;; Laplace transform of little Lommel s function 1904 (cond ((setq l (m2-ones u *var*)) 1905 (setq index1 (cdras 'v1 l) 1906 index2 (cdras 'v2 l) 1907 arg1 (cdras 'w l) 1908 rest (cdras 'u l)) 1909 (return (lt1s rest arg1 index1 index2)))) 1910 1911 ;; Laplace transform of Lommel S function 1912 (cond ((setq l (m2-oneslommel u *var*)) 1913 (setq index1 (cdras 'v1 l) 1914 index2 (cdras 'v2 l) 1915 arg1 (cdras 'w l) 1916 rest (cdras 'u l)) 1917 (return (fractest2 rest arg1 index1 index2 'slommel)))) 1918 1919 ;; Laplace transform of Bessel Y function 1920 (cond ((setq l (m2-oney u *var*)) 1921 (setq index1 (cdras 'v l) 1922 arg1 (cdras 'w l) 1923 rest (cdras 'u l)) 1924 (return (lt1yref rest arg1 index1)))) 1925 1926 ;; Laplace transform of Bessel K function 1927 (cond ((setq l (m2-onek u *var*)) 1928 (setq index1 (cdras 'v l) 1929 arg1 (cdras 'w l) 1930 rest (cdras 'u l)) 1931 (cond ((zerop1 index1) 1932 ;; Special case for a zero index 1933 (return (lt-bessel_k0 rest arg1))) 1934 (t 1935 (return (fractest2 rest arg1 index1 nil 'kti)))))) 1936 1937 ;; Laplace transform of Parabolic Cylinder function 1938 (cond ((setq l (m2-parabolic_cylinder_d u *var*)) 1939 (setq index1 (cdras 'v l) 1940 arg1 (cdras 'w l) 1941 rest (cdras 'u l)) 1942 (return (fractest2 rest arg1 index1 nil 'd)))) 1943 1944 ;; Laplace transform of Incomplete Gamma function 1945 (cond ((setq l (m2-onegammaincomplete u *var*)) 1946 (setq arg1 (cdras 'w1 l) 1947 arg2 (cdras 'w2 l) 1948 rest (cdras 'u l)) 1949 (return (fractest2 rest arg1 arg2 nil 'gamma_incomplete)))) 1950 1951 ;; Laplace transform of Batemann function 1952 (cond ((setq l (m2-onekbateman u *var*)) 1953 (setq index1 (cdras 'v l) 1954 arg1 (cdras 'w l) 1955 rest (cdras 'u l)) 1956 (return (fractest2 rest arg1 index1 nil 'kbateman)))) 1957 1958 ;; Laplace transform of Bessel J function 1959 (cond ((setq l (m2-onej u *var*)) 1960 (setq index1 (cdras 'v l) 1961 arg1 (cdras 'w l) 1962 rest (cdras 'u l)) 1963 (return (lt1j rest arg1 index1)))) 1964 1965 ;; Laplace transform of lower incomplete Gamma function 1966 (cond ((setq l (m2-onegamma-incomplete-lower u *var*)) 1967 (setq arg1 (cdras 'w1 l) 1968 arg2 (cdras 'w2 l) 1969 rest (cdras 'u l)) 1970 (return (lt1gamma-incomplete-lower rest arg1 arg2)))) 1971 1972 ;; Laplace transform of Hankel 1 function 1973 (cond ((setq l (m2-hankel_1 u *var*)) 1974 (setq index1 (cdras 'v l) 1975 arg1 (cdras 'w l) 1976 rest (cdras 'u l)) 1977 (return (fractest2 rest arg1 index1 1 'htjory)))) 1978 1979 ;; Laplace transform of Hankel 2 function 1980 (cond ((setq l (m2-hankel_2 u *var*)) 1981 (setq index1 (cdras 'v l) 1982 arg1 (cdras 'w l) 1983 rest (cdras 'u l)) 1984 (return (fractest2 rest arg1 index1 2 'htjory)))) 1985 1986 ;; Laplace transform of Whittaker M function 1987 (cond ((setq l (m2-onem u *var*)) 1988 (setq index1 (cdras 'v1 l) 1989 index11 (cdras 'v2 l) 1990 arg1 (cdras 'w l) 1991 rest (cdras 'u l)) 1992 (return (lt1m rest arg1 index1 index11)))) 1993 1994 ;; Laplace transform of Whittaker M function 1995 (cond ((setq l (m2-whittaker_m u *var*)) 1996 (setq index1 (cdras 'v1 l) 1997 index2 (cdras 'v2 l) 1998 arg1 (cdras 'w l) 1999 rest (cdras 'u l)) 2000 (return (lt1m rest arg1 index1 index2)))) 2001 2002 ;; Laplace transform of the Generalized Laguerre function, %l[v1,v2](w) 2003 (cond ((setq l (m2-onel u *var*)) 2004 (setq index1 (cdras 'v1 l) 2005 index11 (cdras 'v2 l) 2006 arg1 (cdras 'w l) 2007 rest (cdras 'u l)) 2008 (return (integertest rest arg1 index1 index11 'l)))) 2009 2010 ;; Laplace transform for the Generalized Laguerre function 2011 ;; We call the routine for %l[v1,v2](w). 2012 (cond ((setq l (m2-one-gen-laguerre u *var*)) 2013 (setq index1 (cdras 'v1 l) 2014 index2 (cdras 'v2 l) 2015 arg1 (cdras 'w l) 2016 rest (cdras 'u l)) 2017 (return (integertest rest arg1 index1 index2 'l)))) 2018 2019 ;; Laplace transform for the Laguerre function 2020 ;; We call the routine for %l[v1,0](w). 2021 (cond ((setq l (m2-one-laguerre u *var*)) 2022 (setq index1 (cdras 'v1 l) 2023 arg1 (cdras 'w l) 2024 rest (cdras 'u l)) 2025 (return (integertest rest arg1 index1 0 'l)))) 2026 2027 ;; Laplace transform of Gegenbauer function 2028 (cond ((setq l (m2-onec u *var*)) 2029 (setq index1 (cdras 'v1 l) 2030 index11 (cdras 'v2 l) 2031 arg1 (cdras 'w l) 2032 rest (cdras 'u l)) 2033 (return (integertest rest arg1 index1 index11 'c)))) 2034 2035 ;; Laplace transform of Chebyshev function of the first kind 2036 (cond ((setq l (m2-onet u *var*)) 2037 (setq index1 (cdras 'v1 l) 2038 arg1 (cdras 'w l) 2039 rest (cdras 'u l)) 2040 (return (integertest rest arg1 index1 nil 't)))) 2041 2042 ;; Laplace transform of Chebyshev function of the second kind 2043 (cond ((setq l (m2-oneu u *var*)) 2044 (setq index1 (cdras 'v1 l) 2045 arg1 (cdras 'w l) 2046 rest (cdras 'u l)) 2047 (return (integertest rest arg1 index1 nil 'u)))) 2048 2049 ;; Laplace transform for the Hermite function, hermite(index1,arg1) 2050 (cond ((setq l (m2-one-hermite u *var*)) 2051 (setq index1 (cdras 'v1 l) 2052 arg1 (cdras 'w l) 2053 rest (cdras 'u l)) 2054 (return 2055 (cond ((and (maxima-integerp index1) 2056 (or (mevenp index1) 2057 (moddp index1))) 2058 ;; When index1 is an even or odd integer, we transform 2059 ;; directly to a hypergeometric function. For this case we 2060 ;; get a Laplace transform when the arg is the 2061 ;; square root of the variable. 2062 (sendexec rest (hermite-to-hypergeometric index1 arg1))) 2063 (t 2064 (integertest rest arg1 index1 nil 'he)))))) 2065 2066 ;; Laplace transform of %p[v1,v2](w), Associated Legendre P function 2067 (cond ((setq l (m2-hyp-onep u *var*)) 2068 (setq index1 (cdras 'v1 l) 2069 index11 (cdras 'v2 l) 2070 arg1 (cdras 'w l) 2071 rest (cdras 'u l)) 2072 (return (lt1p rest arg1 index1 index11)))) 2073 2074 ;; Laplace transform of Associated Legendre P function 2075 (cond ((setq l (m2-assoc_legendre_p u *var*)) 2076 (setq index1 (cdras 'v1 l) 2077 index2 (cdras 'v2 l) 2078 arg1 (cdras 'w l) 2079 rest (cdras 'u l)) 2080 (return (lt1p rest arg1 index1 index2)))) 2081 2082 ;; Laplace transform of %p[v1,v2,v3](w), Jacobi function 2083 (cond ((setq l (m2-onepjac u *var*)) 2084 (setq index1 (cdras 'v1 l) 2085 index2 (cdras 'v2 l) 2086 index21 (cdras 'v3 l) 2087 arg1 (cdras 'w l) 2088 rest (cdras 'u l)) 2089 (return (pjactest rest arg1 index1 index2 index21)))) 2090 2091 ;; Laplace transform of Jacobi P function 2092 (cond ((setq l (m2-jacobi_p u *var*)) 2093 (setq index1 (cdras 'v1 l) 2094 index2 (cdras 'v2 l) 2095 index21 (cdras 'v3 l) 2096 arg1 (cdras 'w l) 2097 rest (cdras 'u l)) 2098 (return (pjactest rest arg1 index1 index2 index21)))) 2099 2100 ;; Laplace transform of Associated Legendre function of the second kind 2101 (cond ((setq l (m2-oneq u *var*)) 2102 (setq index1 (cdras 'v1 l) 2103 index11 (cdras 'v2 l) 2104 arg1 (cdras 'w l) 2105 rest (cdras 'u l)) 2106 (return (lt1q rest arg1 index1 index11)))) 2107 2108 ;; Laplace transform of Associated Legendre function of the second kind 2109 (cond ((setq l (m2-assoc_legendre_q u *var*)) 2110 (setq index1 (cdras 'v1 l) 2111 index2 (cdras 'v2 l) 2112 arg1 (cdras 'w l) 2113 rest (cdras 'u l)) 2114 (return (lt1q rest arg1 index1 index2)))) 2115 2116 ;; Laplace transform of %p[v1](w), Legendre P function 2117 (cond ((setq l (m2-onep0 u *var*)) 2118 (setq index1 (cdras 'v1 l) 2119 index11 0 2120 arg1 (cdras 'w l) 2121 rest (cdras 'u l)) 2122 (return (lt1p rest arg1 index1 index11)))) 2123 2124 ;; Laplace transform of Legendre P function 2125 (cond ((setq l (m2-legendre_p u *var*)) 2126 (setq index1 (cdras 'v1 l) 2127 arg1 (cdras 'w l) 2128 rest (cdras 'u l)) 2129 (return (lt1p rest arg1 index1 0)))) 2130 2131 ;; Laplace transform of Whittaker W function 2132 (cond ((setq l (m2-onew u *var*)) 2133 (setq index1 (cdras 'v1 l) 2134 index11 (cdras 'v2 l) 2135 arg1 (cdras 'w l) 2136 rest (cdras 'u l)) 2137 (return (whittest rest arg1 index1 index11)))) 2138 2139 ;; Laplace transform of Whittaker W function 2140 (cond ((setq l (m2-whittaker_w u *var*)) 2141 (setq index1 (cdras 'v1 l) 2142 index2 (cdras 'v2 l) 2143 arg1 (cdras 'w l) 2144 rest (cdras 'u l)) 2145 (return (whittest rest arg1 index1 index2)))) 2146 2147 ;; Laplace transform of square of Bessel J function 2148 (cond ((setq l (m2-onej^2 u *var*)) 2149 (setq index1 (cdras 'v l) 2150 arg1 (cdras 'w l) 2151 rest (cdras 'u l)) 2152 (return (lt1j^2 rest arg1 index1)))) 2153 2154 ;; Laplace transform of square of Hankel 1 function 2155 (cond ((setq l (m2-hankel_1^2 u *var*)) 2156 (setq index1 (cdras 'v l) 2157 arg1 (cdras 'w l) 2158 rest (cdras 'u l)) 2159 (return (fractest rest arg1 arg1 index1 1 index1 1 '2htjory)))) 2160 2161 ;; Laplace transform of square of Hankel 2 function 2162 (cond ((setq l (m2-hankel_2^2 u *var*)) 2163 (setq index1 (cdras 'v l) 2164 arg1 (cdras 'w l) 2165 rest (cdras 'u l)) 2166 (return (fractest rest arg1 arg1 index1 2 index1 2 '2htjory)))) 2167 2168 ;; Laplace transform of square of Bessel Y function 2169 (cond ((setq l (m2-oney^2 u *var*)) 2170 (setq index1 (cdras 'v l) 2171 arg1 (cdras 'w l) 2172 rest (cdras 'u l)) 2173 (return (fractest rest arg1 arg1 index1 nil index1 nil '2ytj)))) 2174 2175 ;; Laplace transform of square of Bessel K function 2176 (cond ((setq l (m2-onek^2 u *var*)) 2177 (setq index1 (cdras 'v l) 2178 arg1 (cdras 'w l) 2179 rest (cdras 'u l)) 2180 (return (fractest rest arg1 arg1 index1 nil index1 nil '2kti)))) 2181 2182 ;; Laplace transform of two Bessel I functions 2183 (cond ((setq l (m2-twoi u *var*)) 2184 (setq index1 (cdras 'v1 l) 2185 index2 (cdras 'v2 l) 2186 arg1 (mul '$%i (cdras 'w1 l)) 2187 arg2 (mul '$%i (cdras 'w2 l)) 2188 rest (mul (power '$%i (neg index1)) 2189 (power '$%i (neg index1)) 2190 (cdras 'u l))) 2191 (return (lt2j rest arg1 arg2 index1 index2)))) 2192 2193 ;; Laplace transform of Bessel I. We use I[v](w)=%i^n*J[n](%i*w). 2194 (cond ((setq l (m2-onei u *var*)) 2195 (setq index1 (cdras 'v l) 2196 arg1 (mul '$%i (cdras 'w l)) 2197 rest (mul (power '$%i (neg index1)) (cdras 'u l))) 2198 (return (lt1j rest arg1 index1)))) 2199 2200 ;; Laplace transform of square of Bessel I function 2201 (cond ((setq l (m2-onei^2 u *var*)) 2202 (setq index1 (cdras 'v l) 2203 arg1 (mul '$%i (cdras 'w l)) 2204 rest (mul (power '$%i (neg index1)) 2205 (power '$%i (neg index1)) 2206 (cdras 'u l))) 2207 (return (lt1j^2 rest arg1 index1)))) 2208 2209 ;; Laplace transform of Erf function 2210 (cond ((setq l (m2-onerf u *var*)) 2211 (setq arg1 (cdras 'w l) 2212 rest (cdras 'u l)) 2213 (return (lt1erf rest arg1)))) 2214 2215 ;; Laplace transform of the logarithmic function. 2216 ;; We add an algorithm for the Laplace transform and call the routine 2217 ;; lt-log. The old code is still present, but isn't called. 2218 (cond ((setq l (m2-onelog u *var*)) 2219 (setq arg1 (cdras 'w l) 2220 rest (cdras 'u l)) 2221 (return (lt-log rest arg1)))) 2222 2223 ;; Laplace transform of Erfc function 2224 (cond ((setq l (m2-onerfc u *var*)) 2225 (setq arg1 (cdras 'w l) 2226 rest (cdras 'u l)) 2227 (return (fractest2 rest arg1 nil nil 'erfc)))) 2228 2229 ;; Laplace transform of expintegral_ei. 2230 ;; Maxima uses the build in transformation to the gamma_incomplete 2231 ;; function and simplifies the log functions of the transformation. We do 2232 ;; not use the dispatch mechanism of fractest2, but call sendexec directly 2233 ;; with the transformed function. 2234 (cond ((setq l (m2-oneexpintegral_ei u *var*)) 2235 (setq arg1 (cdras 'w l) 2236 rest (cdras 'u l)) 2237 (let (($expintrep '%gamma_incomplete) 2238 ($logexpand '$all)) 2239 (return (sratsimp (sendexec rest ($expintegral_ei arg1))))))) 2240 2241 ;; Laplace transform of expintegral_e1 2242 (cond ((setq l (m2-oneexpintegral_e1 u *var*)) 2243 (setq arg1 (cdras 'w l) 2244 rest (cdras 'u l)) 2245 (let (($expintrep '%gamma_incomplete) 2246 ($logexpand '$all)) 2247 (return (sratsimp (sendexec rest ($expintegral_e1 arg1))))))) 2248 2249 ;; Laplace transform of expintegral_e 2250 (cond ((setq l (m2-oneexpintegral_e u *var*)) 2251 (setq arg1 (cdras 'v l) 2252 arg2 (cdras 'w l) 2253 rest (cdras 'u l)) 2254 (let (($expintrep '%gamma_incomplete) 2255 ($logexpand '$all)) 2256 (return (sratsimp (sendexec rest ($expintegral_e arg1 arg2))))))) 2257 2258 ;; Laplace transform of expintegral_si 2259 (cond ((setq l (m2-oneexpintegral_si u *var*)) 2260 (setq arg1 (cdras 'w l) 2261 rest (cdras 'u l)) 2262 ;; We transform to the hypergeometric representation. 2263 (return 2264 (sendexec rest (expintegral_si-to-hypergeometric arg1))))) 2265 2266 ;; Laplace transform of expintegral_shi 2267 (cond ((setq l (m2-oneexpintegral_shi u *var*)) 2268 (setq arg1 (cdras 'w l) 2269 rest (cdras 'u l)) 2270 ;; We transform to the hypergeometric representation. 2271 (return 2272 (sendexec rest (expintegral_shi-to-hypergeometric arg1))))) 2273 2274 ;; Laplace transform of expintegral_ci 2275 (cond ((setq l (m2-oneexpintegral_ci u *var*)) 2276 (setq arg1 (cdras 'w l) 2277 rest (cdras 'u l)) 2278 ;; We transform to the hypergeometric representation. 2279 ;; Because we have Logarithmic terms in the transformation, 2280 ;; we switch on the flag $logexpand and do a ratsimp. 2281 (let (($logexpand '$super)) 2282 (return 2283 (sratsimp 2284 (sendexec rest (expintegral_ci-to-hypergeometric arg1))))))) 2285 2286 ;; Laplace transform of expintegral_chi 2287 (cond ((setq l (m2-oneexpintegral_chi u *var*)) 2288 (setq arg1 (cdras 'w l) 2289 rest (cdras 'u l)) 2290 ;; We transform to the hypergeometric representation. 2291 ;; Because we have Logarithmic terms in the transformation, 2292 ;; we switch on the flag $logexpand and do a ratsimp. 2293 (let (($logexpand '$super)) 2294 (return 2295 (sratsimp 2296 (sendexec rest (expintegral_chi-to-hypergeometric arg1))))))) 2297 2298 ;; Laplace transform of Complete elliptic integral of the first kind 2299 (cond ((setq l (m2-onekelliptic u *var*)) 2300 (setq arg1 (cdras 'w l) 2301 rest (cdras 'u l)) 2302 (return (lt1kelliptic rest arg1)))) 2303 2304 ;; Laplace transform of Complete elliptic integral of the first kind 2305 (cond ((setq l (m2-elliptic_kc u *var*)) 2306 (setq arg1 (cdras 'w l) 2307 rest (cdras 'u l)) 2308 (return (lt1kelliptic rest arg1)))) 2309 2310 ;; Laplace transform of Complete elliptic integral of the second kind 2311 (cond ((setq l (m2-onee u *var*)) 2312 (setq arg1 (cdras 'w l) 2313 rest (cdras 'u l)) 2314 (return (lt1e rest arg1)))) 2315 2316 ;; Laplace transform of Complete elliptic integral of the second kind 2317 (cond ((setq l (m2-elliptic_ec u *var*)) 2318 (setq arg1 (cdras 'w l) 2319 rest (cdras 'u l)) 2320 (return (lt1e rest arg1)))) 2321 2322 ;; Laplace transform of %f[v1,v2](w1,w2,w3), Hypergeometric function 2323 ;; We support the Laplace transform of the build in symbol %f. We do 2324 ;; not use the mechanism of defining an "Expert on Laplace transform", 2325 ;; the expert function does a call to lt-ltp. We do this call directly. 2326 (cond ((setq l (m2-onef u *var*)) 2327 (setq rest (cdras 'u l) 2328 arg1 (cdras 'w3 l) 2329 index1 (list (cdras 'w1 l) (cdras 'w2 l))) 2330 (return (lt-ltp 'f rest arg1 index1)))) 2331 2332 ;; Laplace transform of Hypergeometric function 2333 (cond ((setq l (m2-hypergeometric u *var*)) 2334 (setq rest (cdras 'u l) 2335 arg1 (cdras 'w3 l) 2336 index1 (list (cdras 'w1 l) (cdras 'w2 l))) 2337 (return (lt-ltp 'f rest arg1 index1)))) 2338 2339 ;; Laplace transform of c * t^v * (a+t)^w 2340 ;; It is possible to combine arbpow2 and arbpow. 2341 (cond ((setq l (m2-arbpow2 u *var*)) 2342 (setq rest (cdras 'c l) 2343 arg1 (cdras 'a l) 2344 arg2 (cdras 'b l) 2345 index1 (cdras 'v l) 2346 index2 (cdras 'w l)) 2347 (return (lt-arbpow2 rest arg1 arg2 index1 index2)))) 2348 2349 ;; Laplace transform of c * t^v 2350 (cond ((setq l (m2-arbpow1 u *var*)) 2351 (setq arg1 (cdras 'u l) 2352 arg2 (cdras 'c l) 2353 index1 (cdras 'v l)) 2354 (return (mul arg2 (lt-arbpow arg1 index1))))) 2355 2356 ;; We have specialized the pattern for arbpow1. Now a lot of integrals 2357 ;; will fail correctly and we have to return a noun form. 2358 (return (setq *hyp-return-noun-flag* 'other-j-cases-next)))) 2359 2360;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 2361;;; 2362;;; Algorithm 2.1: Laplace transform of c*t^v*%e(-p*t) 2363;;; 2364;;; Table of Integral Transforms 2365;;; 2366;;; p. 137, formula 1: 2367;;; 2368;;; t^u*exp(-p*t) 2369;;; -> gamma(u+1)*p^(-u-1) 2370;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 2371 2372(defun lt-arbpow (expr pow) 2373 (cond ((or (eq expr *var*) (equal pow 0)) 2374 (f1p137test pow)) 2375 (t 2376 (setq *hyp-return-noun-flag* 'lt-arbow-failed)))) 2377 2378;; Check if conditions for f1p137 hold 2379(defun f1p137test (pow) 2380 (cond ((eq ($asksign (add pow 1)) '$pos) 2381 (f1p137 pow)) 2382 (t 2383 (setq *hyp-return-noun-flag* 'fail-in-arbpow)))) 2384 2385(defun f1p137 (pow) 2386 (mul (take '(%gamma) (add pow 1)) 2387 (power *par* (sub (mul -1 pow) 1)))) 2388 2389;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 2390;;; 2391;;; Algorithm 2.2: Laplace transform of c*t^v*(1+t)^w 2392;;; 2393;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 2394 2395(defun lt-arbpow2 (c a b pow1 pow2) 2396 (if (eq ($asksign (add pow1 1)) '$pos) 2397 (cond 2398 ((equal pow1 0) 2399 ;; The Laplace transform is an Incomplete Gamma function. 2400 (mul c 2401 (power a (add pow2 1)) 2402 (inv b) 2403 (power (mul *par* a (inv b)) (mul -1 (add pow2 1))) 2404 (power '$%e (mul *par* a (inv b))) 2405 (take '(%gamma_incomplete) (add pow2 1) (mul *par* a (inv b))))) 2406 ((not (maxima-integerp (add pow1 pow2 2))) 2407 ;; The general result is a Hypergeometric U function U(a,b,z) which can 2408 ;; be represented by two Hypergeometic 1F1 functions for the special 2409 ;; case that the index b is not an integer value. 2410 (add (mul c 2411 (power a (add pow1 pow2 1)) 2412 (inv (power b (add pow1 1))) 2413 (take '(%gamma) (add pow1 pow2 1)) 2414 (power (mul *par* a (inv b)) (mul -1 (add pow1 pow2 1))) 2415 (hgfsimp-exec (list (mul -1 pow2)) 2416 (list (mul -1 (add pow1 pow2))) 2417 (mul *par* a (inv b)))) 2418 (mul c 2419 (power a (add pow1 pow2 1)) 2420 (inv (power b (add pow1 1))) 2421 (take '(%gamma) (add pow1 1)) 2422 (take '(%gamma) (mul -1 (add pow1 pow2 1))) 2423 (inv (take '(%gamma) (mul -1 pow2))) 2424 (hgfsimp-exec (list (add pow1 1)) 2425 (list (add pow1 pow2 2)) 2426 (mul *par* a (inv b)))))) 2427 (t 2428 ;; The most general case is a result with the Hypergeometric U function. 2429 (mul c 2430 (power a (add pow1 pow2 1)) 2431 (inv (power b (add pow1 1))) 2432 (take '(%gamma) (add pow1 1)) 2433 (list '(%hypergeometric_u simp) 2434 (add pow1 1) 2435 (add pow1 pow2 2) 2436 (mul *par* a (inv b)))))) 2437 (setq *hyp-return-noun-flag* 'lt-arbpow2-failed))) 2438 2439;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 2440;;; 2441;;; Algorithm 2.3: Laplace transform of the Logarithmic function 2442;;; 2443;;; c*t^(v-1)*log(a*t) 2444;;; -> c*gamma(v)*s^(-v)*(psi[0](v)-log(s/a)) 2445;;; 2446;;; This is the formula for an expression with log(t) scaled like 1/a*F(s/a). 2447;;; 2448;;; For the following cases we have to add further algorithm: 2449;;; log(1+a*x), log(x+a), log(x)^2. 2450;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 2451 2452(defun lt-log (rest arg) 2453 (let* ((l (m2-c*t^v rest *var*)) 2454 (c (cdras 'c l)) 2455 (v (add (cdras 'v l) 1))) ; because v -> v-1 2456 (cond 2457 ((and l (eq ($asksign v) '$pos)) 2458 (let* ((l1 (m2-a*t arg *var*)) 2459 (a (cdras 'a l1))) 2460 (cond (l1 2461 (mul c 2462 (take '(%gamma) v) 2463 (inv (power *par* v)) 2464 (sub (take '(mqapply) (list '($psi array) 0) v) 2465 (take '(%log) (div *par* a))))) 2466 (t 2467 (setq *hyp-return-noun-flag* 'lt-log-failed))))) 2468 (t 2469 (setq *hyp-return-noun-flag* 'lt-log-failed))))) 2470 2471;; Pattern for lt-log. 2472;; Extract the argument of a function: a*t+c for c=0. 2473(defun m2-a*t (expr var) 2474 (m2 expr 2475 `((mplus) 2476 ((mtimes) (u alike1 ,var) (a free ,var)) 2477 ((coeffpp) (c equal 0))))) 2478 2479;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 2480;;; Algorithm 2.4: Laplace transfom of the Whittaker function 2481;;; 2482;;; Test for Whittaker W function. Simplify this if possible, or 2483;;; convert to Whittaker M function. 2484;;; 2485;;; We have r * %w[i1,i2](a) 2486;;; 2487;;; Formula 16, p. 217 2488;;; 2489;;; t^(v-1)*%w[k,u](a*t) 2490;;; -> gamma(u+v+1/2)*gamma(v-u+1/2)*a^(u+1/2)/ 2491;;; (gamma(v-k+1)*(p+a/2)^(u+v+1/2) 2492;;; *2f1(u+v+1/2,u-k+1/2;v-k+1;(p-a/2)/(p+a/2)) 2493;;; 2494;;; For Re(v +/- mu) > -1/2 2495;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 2496 2497(defun whittest (r a i1 i2) 2498 (let (n m) 2499 (cond ((f16p217test r a i1 i2)) 2500 ((and 2501 (not 2502 (and (maxima-integerp (setq n (sub (sub '((rat simp) 1 2) i2) i1))) 2503 (member ($sign n) '($zero $neg $nz)))) 2504 (not 2505 (and (maxima-integerp (setq m (sub (add '((rat simp) 1 2) i2) i1))) 2506 (member ($sign m) '($zero $neg $nz))))) 2507 ;; 1/2-u-k and 1/2+u-k are not zero or a negative integer 2508 ;; Transform to Whittaker M and try again. 2509 (distrexecinit ($expand (mul (init r) (wtm a i1 i2))))) 2510 (t 2511 ;; Both conditions fails, return a noun form. 2512 (setq *hyp-return-noun-flag* 'whittest-failed))))) 2513 2514(defun f16p217test (r a i1 i2) 2515 ;; We have r*%w[i1,i2](a) 2516 (let ((l (m2-c*t^v r *var*))) 2517 ;; Make sure r is of the form c*t^v 2518 (when l 2519 (let* ((v (add (cdras 'v l) 1)) 2520 (c (cdras 'c l))) 2521 ;; Check that v + i2 + 1/2 > 0 and v - i2 + 1/2 > 0. 2522 (when (and (eq ($asksign (add (add v i2) '((rat simp) 1 2))) '$pos) 2523 (eq ($asksign (add (sub v i2) '((rat simp) 1 2))) '$pos)) 2524 ;; Ok, we satisfy the conditions. Now extract the arg. 2525 ;; The transformation is only valid for an argument a*t. We have 2526 ;; to special the pattern to make sure that we satisfy the condition. 2527 (let ((l (m2-a*t a *var*))) 2528 (when l 2529 (let ((a (cdras 'a l))) 2530 ;; We're ready now to compute the transform. 2531 (mul c 2532 (power a (add i2 '((rat simp) 1 2))) 2533 (take '(%gamma) (add (add v i2) '((rat simp) 1 2))) 2534 (take '(%gamma) (add (sub v i2) '((rat simp) 1 2))) 2535 (inv (mul (take '(%gamma) (add (sub v i1) 1)) 2536 (power (add *par* (div a 2)) 2537 (add (add i2 v) '((rat simp) 1 2))))) 2538 (hgfsimp-exec (list (add (add i2 v '((rat simp) 1 2))) 2539 (add (sub i2 i1) '((rat simp) 1 2))) 2540 (list (add (sub v i1) 1)) 2541 (div (sub *par* (div a 2)) 2542 (add *par* (div a 2))))))))))))) 2543 2544;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 2545;;; Algorithm 2.5: Laplace transfom of bessel_k(0,a*t) 2546;;; 2547;;; The general algorithm handles the Bessel K function for an order |v|<1. 2548;;; but does not include the special case v=0. Return the Laplace transform: 2549;;; 2550;;; bessel_k(0,a*t) --> acosh(s/a)/sqrt(s^2-a^2) 2551;;; 2552;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 2553 2554(defun lt-bessel_k0 (rest arg) 2555 (let* ((l (m2-c*t^v rest *var*)) 2556 (c (cdras 'c l)) 2557 (v (cdras 'v l)) 2558 (l (m2-a*t arg *var*)) 2559 (a (cdras 'a l))) 2560 (cond ((and l (zerop1 v)) 2561 (mul c 2562 (take '(%acosh) (div *par* a)) 2563 (inv (power (sub (mul *par* *par*) (mul a a)) 2564 '((rat simp) 1 2))))) 2565 (t 2566 (setq *hyp-return-noun-flag* 'lt-bessel_k-failed))))) 2567 2568;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 2569;;; 2570;;; DISPATCH FUNCTIONS TO CHANGE THE REPRESENTATION OF SPECIAL FUNCTIONS 2571;;; 2572;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 2573 2574;; Laplace transform of a product of Bessel functions. A1, A2 are 2575;; the args of the two functions. I1, I2 are the indices of each 2576;; function. I11, I21 are secondary indices of each function, if any. 2577;; FLG is a symbol indicating how we should handle the special 2578;; functions (and also indicates what the special functions are.) 2579;; 2580;; I11 and I21 are for the Hankel functions. 2581 2582(defun fractest (r a1 a2 i1 i11 i2 i21 flg) 2583 (cond ((or (and (listp i1) (equal (caar i1) 'rat) 2584 (listp i2) (equal (caar i2) 'rat)) 2585 (eq flg '2htjory)) 2586 ;; We have to Bessel or Hankel functions. Both indizes have to be 2587 ;; rational numbers or we have two Hankel functions. 2588 (sendexec r 2589 (cond ((eq flg '2ytj) 2590 (mul (ytj i1 a1) 2591 (ytj i2 a2))) 2592 ((eq flg '2htjory) 2593 (mul (htjory i1 i11 a1) 2594 (htjory i2 i21 a2))) 2595 ((eq flg 'ktiytj) 2596 (mul (kti i1 a1) 2597 (ytj i2 a2))) 2598 ((eq flg '2kti) 2599 (mul (kti i1 a1) 2600 (kti i2 a2)))))) 2601 (t 2602 (setq *hyp-return-noun-flag* 'product-of-y-with-nofract-indices)))) 2603 2604;; Laplace transform of a product of Bessel functions. A1, A2 are 2605;; the args of the two functions. I1, I2 are the indices of each 2606;; function. I is a secondary index to one function, if any. 2607;; FLG is a symbol indicating how we should handle the special 2608;; functions (and also indicates what the special functions are.) 2609;; 2610;; I is for the kind of Hankel function. 2611 2612(defun fractest1 (r a1 a2 i1 i2 i flg) 2613 (cond ((or (and (listp i2) 2614 (equal (caar i2) 'rat)) 2615 (eq flg 'besshtjory)) 2616 ;; We have two Bessel or Hankel functions. The second index has to 2617 ;; be a rational number or one of the functions is a Hankel function 2618 ;; and the second function is Bessel J or Bessel I 2619 (sendexec r 2620 (cond ((eq flg 'bessytj) 2621 (mul (take '(%bessel_j) i1 a1) 2622 (ytj i2 a2))) 2623 ((eq flg 'besshtjory) 2624 (mul (take '(%bessel_j) i1 a1) 2625 (htjory i2 i a2))) 2626 ((eq flg 'htjoryytj) 2627 (mul (htjory i1 i a1) 2628 (ytj i2 a2))) 2629 ((eq flg 'besskti) 2630 (mul (take '(%bessel_j) i1 a1) 2631 (kti i2 a2))) 2632 ((eq flg 'htjorykti) 2633 (mul (htjory i1 i a1) 2634 (kti i2 a2)))))) 2635 (t 2636 (setq *hyp-return-noun-flag* 'product-of-i-y-of-nofract-index)))) 2637 2638;; Laplace transform of a single special function. A is the arg of 2639;; the special function. I1, I11 are the indices of the function. FLG 2640;; is a symbol indicating how we should handle the special functions 2641;; (and also indicates what the special functions are.) 2642;; 2643;; I11 is the kind of Hankel function 2644 2645(defun fractest2 (r a1 i1 i11 flg) 2646 (cond ((or (and (listp i1) 2647 (equal (caar i1) 'rat)) 2648 (eq flg 'd) 2649 (eq flg 'kbateman) 2650 (eq flg 'gamma_incomplete) 2651 (eq flg 'htjory) 2652 (eq flg 'erfc) 2653 (eq flg 'slommel) 2654 (eq flg 'ytj)) 2655 ;; The index is a rational number or flg has the value of one of the 2656 ;; above special functions. 2657 (sendexec r 2658 (cond ((eq flg 'ytj) 2659 (ytj i1 a1)) 2660 ((eq flg 'htjory) 2661 (htjory i1 i11 a1)) 2662 ((eq flg 'd) 2663 (dtw i1 a1)) 2664 ((eq flg 'kbateman) 2665 (kbatemantw i1 a1)) 2666 ((eq flg 'gamma_incomplete) 2667 (gamma_incomplete-to-gamma-incomplete-lower a1 i1)) 2668 ((eq flg 'kti) 2669 (kti i1 a1)) 2670 ((eq flg 'erfc) 2671 (erfctd a1)) 2672 ((eq flg 'slommel) 2673 (slommeltjandy i1 i11 a1))))) 2674 (t 2675 (setq *hyp-return-noun-flag* 'y-of-nofract-index)))) 2676 2677;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 2678 2679(defun integertest (r arg i1 i2 flg) 2680 (cond ((maxima-integerp i1) 2681 (dispatchpoltrans r arg i1 i2 flg)) 2682 (t 2683 (setq *hyp-return-noun-flag* 'index-should-be-an-integer-in-polys)))) 2684 2685(defun dispatchpoltrans (r x i1 i2 flg) 2686 (sendexec r 2687 (cond ((eq flg 'l)(ltw x i1 i2)) 2688 ((eq flg 'he)(hetd x i1)) 2689 ((eq flg 'c)(ctpjac x i1 i2)) 2690 ((eq flg 't)(ttpjac x i1)) 2691 ((eq flg 'u)(utpjac x i1))))) 2692 2693;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 2694;; (-1)^n*n!*laguerre(n,a,x) = U(-n,a+1,x) 2695;; 2696;; W[k,u](z) = exp(-z/2)*z^(u+1/2)*U(1/2+u-k,1+2*u,z) 2697;; 2698;; So 2699;; 2700;; laguerre(n,a,x) = (-1)^n*U(-n,a+1,x)/n! 2701;; 2702;; U(-n,a+1,x) = exp(z/2)*z^(-a/2-1/2)*W[1/2+a/2+n,a/2](z) 2703;; 2704;; Finally, 2705;; 2706;; laguerre(n,a,x) = (-1)^n/n!*exp(z/2)*z^(-a/2-1/2)*M[1/2+a/2+n,a/2](z) 2707;; 2708;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 2709 2710(defun ltw (x n a) 2711 (let ((diva2 (div a 2))) 2712 (mul (take '(%gamma) (add n a 1)) 2713 (inv (take '(%gamma) (add a 1))) 2714 (inv (take '(%gamma) (add n 1))) 2715 (power x (sub '((rat simp) -1 2) diva2)) 2716 (power '$%e (div x 2)) 2717 (list '(mqapply simp) 2718 (list '($%m simp array) 2719 (add '((rat simp) 1 2) diva2 n) diva2) x)))) 2720 2721;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 2722;; Hermite He function as a parabolic cylinder function 2723;; 2724;; Tables of Integral Transforms 2725;; 2726;; p. 386 2727;; 2728;; D[n](z) = (-1)^n*exp(z^2/4)*diff(exp(-z^2/2),z,n); 2729;; 2730;; p. 369 2731;; 2732;; He[n](x) = (-1)^n*exp(x^2/2)*diff(exp(-x^2/2),x,n) 2733;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 2734 2735(defun hetd (x n) 2736 (mul (power '$%e (mul x x '((rat simp) 1 4))) 2737 ;; At this time the Parabolic Cylinder D function is not implemented 2738 ;; as a simplifying function. We call nevertheless the simplifer. 2739 (take '($parabolic_cylinder_d) n x))) 2740 2741;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 2742 2743;; Transform Gegenbauer function to Jacobi P function 2744;; ultraspherical(n,v,x) = gamma(2*v+n)*gamma(v+1/2)/gamma(2*v)/gamma(v+n+1/2) 2745;; *jacobi_p(n,v-1/2,v-1/2,x) 2746(defun ctpjac (x n v) 2747 (let ((inv2 '((rat simp) 1 2))) 2748 (mul (take '(%gamma) (add v v n)) 2749 (inv (take '(%gamma) (add v v))) 2750 (take '(%gamma) (add inv2 v)) 2751 (inv (take '(%gamma) (add v inv2 n))) 2752 (pjac x n (sub v inv2) (sub v inv2))))) 2753 2754;; Transform Chebyshev T function to Jacobi P function 2755;; chebyshev_t(n,x) = gamma(n+1)*sqrt(%pi)/gamma(n+1/2)*jacobi_p(n,-1/2,-1/2,x) 2756(defun ttpjac (x n) 2757 (let ((inv2 '((rat simp) 1 2))) 2758 (mul (take '(%gamma) n 1) 2759 (power '$%pi inv2) ; gamma(1/2) 2760 (inv (take '(%gamma) (add inv2 n))) 2761 (pjac x n (mul -1 inv2) (mul -1 inv2))))) 2762 2763;; Transform Chebyshev U function to Jacobi P function 2764;; chebyshev_u(n,x) = gamma(n+2)*sqrt(%pi)/2/gamma(3/2+n)*jacobi_p(n,1/2,1/2,x) 2765(defun utpjac (x n) 2766 (let ((inv2 '((rat simp) 1 2))) 2767 (mul (take '(%gamma) (add n 2)) 2768 inv2 2769 (power '$%pi inv2) ; gamma(1/2) 2770 (inv (take '(%gamma) (add inv2 n 1))) 2771 (pjac x n inv2 inv2)))) 2772 2773;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 2774 2775(defun pjactest (rest arg index1 index2 index3) 2776 (cond ((maxima-integerp index1) 2777 (lt-ltp 'onepjac 2778 rest 2779 arg 2780 (list index1 index2 index3))) 2781 (t 2782 (setq *hyp-return-noun-flag* 'ind-should-be-an-integer-in-polys)))) 2783 2784;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 2785;;; 2786;;; Laplace transform of a single Bessel Y function. 2787;;; 2788;;; REST is the multiplier, ARG1 is the arg, and INDEX1 is the order of 2789;;; the Bessel Y function. 2790;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 2791 2792(defun lt1yref (rest arg1 index1) 2793 ;; If the index is an integer, use LT1Y. Otherwise, convert Bessel 2794 ;; Y to Bessel J and compute the transform of that. We do this 2795 ;; because converting Y to J for an integer index doesn't work so 2796 ;; well without taking limits. 2797 (cond ((maxima-integerp index1) 2798 ;; Do not call lt1y but lty directly. 2799 ;; lt1y calls lt-ltp with the flag 'oney. lt-ltp checks this flag 2800 ;; and calls lty. So we can do it at this place and the algorithm is 2801 ;; more simple. 2802 (lty rest arg1 index1)) 2803 (t (fractest2 rest arg1 index1 nil 'ytj)))) 2804 2805;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 2806;;; 2807;;; TRANSFORMATIONS TO CHANGE THE REPRESENTATION OF SPECIAL FUNCTIONS 2808;;; 2809;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 2810 2811;; erfc in terms of D, parabolic cylinder function 2812;; 2813;; Tables of Integral Transforms 2814;; 2815;; p 387: 2816;; erfc(x) = (%pi*x)^(-1/2)*exp(-x^2/2)*W[-1/4,1/4](x^2) 2817;; 2818;; p 386: 2819;; D[v](z) = 2^(v/2+1/2)*z^(-1/2)*W[v/2+1/4,1/4](z^2/2) 2820;; 2821;; So 2822;; 2823;; erfc(x) = %pi^(-1/2)*2^(1/4)*exp(-x^2/2)*D[-1](x*sqrt(2)) 2824 2825(defun erfctd (x) 2826 (let ((inv2 '((rat simp) 1 2))) 2827 (mul (power 2 inv2) ; Should this be 2^(1/4)? 2828 (inv (power '$%pi inv2)) 2829 (power '$%e (mul -1 inv2 x x)) 2830 ;; At this time the Parabolic Cylinder D function is not implemented 2831 ;; as a simplifying function. We call nevertheless the simplifer 2832 ;; to simplify the arguments. When we implement the function 2833 ;; The symbol has to be changed to the noun form. 2834 (take '($parabolic_cylinder_d) -1 (mul (power 2 inv2) x))))) 2835 2836;; Lommel S function in terms of Bessel J and Bessel Y. 2837;; Luke gives 2838;; 2839;; S[u,v](z) = s[u,v](z) + {2^(u-1)*gamma((u-v+1)/2)*gamma((u+v+1)/2)} 2840;; * {sin[(u-v)*%pi/2]*bessel_j(v,z) 2841;; - cos[(u-v)*%pi/2]*bessel_y(v,z) 2842 2843(defun slommeltjandy (m n z) 2844 (let ((arg (mul '((rat simp) 1 2) '$%pi (sub m n)))) 2845 (add (littleslommel m n z) 2846 (mul (power 2 (sub m 1)) 2847 (take '(%gamma) (div (sub (add m 1) n) 2)) 2848 (take '(%gamma) (div (add m n 1) 2)) 2849 (sub (mul (take '(%sin) arg) 2850 (take '(%bessel_j) n z)) 2851 (mul (take '(%cos) arg) 2852 (take '(%bessel_y) n z))))))) 2853 2854;; Whittaker W function in terms of Whittaker M function 2855;; 2856;; A&S 13.1.34 2857;; 2858;; W[k,u](z) = gamma(-2*u)/gamma(1/2-u-k)*M[k,u](z) 2859;; + gamma(2*u)/gamma(1/2+u-k)*M[k,-u](z) 2860 2861(defun wtm (a i1 i2) 2862 (add (mul (take '(%gamma) (mul -2 i2)) 2863 (mwhit a i1 i2) 2864 (inv (take '(%gamma) (sub (sub '((rat simp) 1 2) i2) i1)))) 2865 (mul (take '(%gamma) (add i2 i2)) 2866 (mwhit a i1 (mul -1 i2)) 2867 (inv (take '(%gamma) (sub (add '((rat simp) 1 2) i2) i1)))))) 2868 2869;; Incomplete gamma function in terms of Whittaker W function 2870;; 2871;; Tables of Integral Transforms, p. 387 2872;; 2873;; gamma_incomplete(a,x) = x^((a-1)/2)*exp(-x/2)*W[(a-1)/2,a/2](x) 2874 2875(defun gammaincompletetw (a x) 2876 (mul (power x (div (sub a 1) 2)) 2877 (power '$%e (div x -2)) 2878 (wwhit x (div (sub a 1) 2)(div a 2)))) 2879 2880;;; Incomplete Gamma function in terms of lower incomplete Gamma function 2881;;; 2882;;; Only for a=0 we use the general representation as a Whittaker W function: 2883;;; 2884;;; gamma_incomplete(a,x) = x^((a-1)/2)*exp(-x/2)*W[(a-1)/2,a/2](x) 2885;;; 2886;;; In all other cases we transform to a lower incomplete Gamma function: 2887;;; 2888;;; gamma_incomplete(a,x) = gamma(a)- gamma_incomplete_lower(a,x) 2889;;; 2890;;; The lower incomplete Gamma function will be further transformed to a Hypergeometric 1F1 2891;;; representation. With this change we get more simple and correct results for 2892;;; the Laplace transform of the Incomplete Gamma function. 2893 2894(defun gamma_incomplete-to-gamma-incomplete-lower (a x) 2895 (if (or (eq ($sign a) '$zero) 2896 (and (integerp a) (< a 0))) 2897 ;; The representation as a Whittaker W function for a=0 or a negative 2898 ;; integer (The gamma function is not defined for this case.) 2899 (mul (power x (div (sub a 1) 2)) 2900 (power '$%e (div x -2)) 2901 (wwhit x (div (sub a 1) 2) (div a 2))) 2902 ;; In all other cases the representation as a lower incomplete Gamma function 2903 (sub (take '(%gamma) a) 2904 (list '($gamma_incomplete_lower simp) a x)))) 2905 2906;; Bessel Y in terms of Bessel J 2907;; 2908;; A&S 9.1.2: 2909;; 2910;; bessel_y(v,z) = bessel_j(v,z)*cot(v*%pi)-bessel_j(-v,z)/sin(v*%pi) 2911 2912(defun ytj (i a) 2913 (sub (mul (take '(%bessel_j) i a) 2914 (take '(%cot) (mul i '$%pi))) 2915 (mul (take '(%bessel_j) (mul -1 i) a) 2916 (inv (take '(%sin) (mul i '$%pi)))))) 2917 2918;; Parabolic cylinder function in terms of Whittaker W function. 2919;; 2920;; See Table of Integral Transforms, p.386: 2921;; 2922;; D[v](z) = 2^(v/2+1/4)*z^(-1/2)*W[v/2+1/4,1/4](z^2/2) 2923 2924(defun dtw (i a) 2925 (mul (power 2 (add (div i 2) '((rat simp) 1 4))) 2926 (power a '((rat simp) -1 2)) 2927 (wwhit (mul a a '((rat simp) 1 2)) 2928 (add (div i 2) '((rat simp) 1 4)) 2929 '((rat simp) 1 4)))) 2930 2931;; Bateman's function in terms of Whittaker W function 2932;; 2933;; See Table of Integral Transforms, p.386: 2934;; 2935;; k[2*v](z) = 1/gamma(v+1)*W[v,1/2](2*z) 2936 2937(defun kbatemantw (v a) 2938 (div (wwhit (add a a) (div v 2) '((rat simp) 1 2)) 2939 (take '(%gamma) (add (div v 2) 1)))) 2940 2941;; Bessel K in terms of Bessel I 2942;; 2943;; A&S 9.6.2 2944;; 2945;; bessel_k(v,z) = %pi/2*(bessel_i(-v,z)-bessel_i(v,z))/sin(v*%pi) 2946 2947(defun kti (i a) 2948 (mul '$%pi 2949 '((rat simp) 1 2) 2950 (inv (take '(%sin) (mul i '$%pi))) 2951 (sub (take '(%bessel_i) (mul -1 i) a) 2952 (take '(%bessel_i) i a)))) 2953 2954;; Express Hankel function in terms of Bessel J or Y function. 2955;; 2956;; A&S 9.1.3 2957;; 2958;; H[v,1](z) = %i*csc(v*%pi)*(exp(-v*%pi*%i)*bessel_j(v,z) - bessel_j(-v,z)) 2959;; 2960;; A&S 9.1.4: 2961;; H[v,2](z) = %i*csc(v*%pi)*(bessel_j(-v,z) - exp(-v*%pi*%i)*bessel_j(v,z)) 2962;; 2963;; Both formula are not valid for v an integer. 2964;; For this case use the definitions of the Hankel functions: 2965;; H[v,1](z) = bessel_j(v,z) + %i* bessel_y(v,z) 2966;; H[v,2](z) = bessel_j(v,z) - %i* bessel_y(v,z) 2967;; 2968;; All this can be implemented more simple. 2969;; We do not need the transformation to bessel_j for rational numbers, 2970;; because the correct transformation for bessel_y is already implemented. 2971;; It is enough to use the definitions for the Hankel functions. 2972 2973(defun htjory (v sort z) 2974 ;; V is the order, SORT is the kind of Hankel function (1 or 2), Z 2975 ;; is the arg. 2976 (cond ((and (consp v) 2977 (consp (car v)) 2978 (equal (caar v) 'rat)) 2979 ;; If the order is a rational number of some sort, 2980 ;; 2981 ;; (bessel_j(-v,z) - bessel_j(v,z)*exp(-v*%pi*%i))/(%i*sin(v*%pi*%i)) 2982 (div (numjory v sort z 'j) 2983 (mul '$%i (take '(%sin) (mul v '$%pi))))) 2984 ((equal sort 1) 2985 ;; Transform hankel_1(v,z) to bessel_j(v,z)+%i*bessel_y(v,z) 2986 (add (take '(%bessel_j) v z) 2987 (mul '$%i (take '(%bessel_y) v z)))) 2988 ((equal sort 2) 2989 ;; Transform hankel_2(v,z) to bessel_j(v,z)-%i*bessel_y(v,z) 2990 (sub (take '(%bessel_j) v z) 2991 (mul '$%i (take '(%bessel_y) v z)))) 2992 (t 2993 ;; We should never reach this point of code. 2994 ;; Problem: The user input for the symbol %h[v,sort](t) is not checked. 2995 ;; Therefore the user can generate this error as long as we do not cut 2996 ;; out the support for the Laplace transform of the symbol %h. 2997 (merror "htjory: Called with wrong sort of Hankel functions.")))) 2998 2999;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 3000 3001;;; Three helper functions only used by htjory 3002 3003;; Bessel J or Y, depending on if FLG is 'J or not. 3004(defun desjy (v z flg) 3005 (cond ((eq flg 'j) 3006 (take '(%bessel_j) v z)) 3007 (t 3008 (take '(%bessel_y) v z)))) 3009 3010(defun numjory (v sort z flg) 3011 (cond ((equal sort 1) 3012 ;; bessel(-v, z) - exp(-v*%pi*%i)*bessel(v, z) 3013 ;; 3014 ;; Where bessel is bessel_j if FLG is 'j. Otherwise, bessel 3015 ;; is bessel_y. 3016 ;; 3017 ;; bessel_y(-v, z) - exp(-v*%pi*%i)*bessel_y(v, z) 3018 (sub (desjy (mul -1 v) z flg) 3019 (mul (power '$%e (mul -1 v '$%pi '$%i)) 3020 (desjy v z flg)))) 3021 (t 3022 ;; exp(-v*%pi*%i)*bessel(v,z) - bessel(-v,z), where bessel is 3023 ;; bessel_j or bessel_y, depending on if FLG is 'j or not. 3024 (sub (mul (power '$%e (mul v '$%pi '$%i)) 3025 (desmjy v z flg)) 3026 (desmjy (mul -1 v) z flg))))) 3027 3028(defun desmjy (v z flg) 3029 (cond ((eq flg 'j) 3030 ;; bessel_j(v,z) 3031 (take '(%bessel_j) v z)) 3032 (t 3033 ;; -bessel_y(v,z) 3034 (mul -1 (take '(%bessel_y) v z))))) 3035 3036;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 3037;;; 3038;;; TRANSFORM TO HYPERGEOMETRIC FUNCTION WITHOUT USING THE ROUTINE REF 3039;;; 3040;;; This functions are called in the routine lt-sf-log to get the 3041;;; representation in terms of a hypergeometric function. 3042;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 3043 3044;; The algorithm of the implemented Hermite function %he does not work for 3045;; the known Laplace transforms. For an even or odd integer order, we 3046;; can represent the Hermite function by the Hypergeometric function 1F1. 3047;; With this representations we get the expected Laplace transforms. 3048(defun hermite-to-hypergeometric (order arg) 3049 (cond 3050 ((and (maxima-integerp order) 3051 (mevenp order)) 3052 ;; Transform to 1F1 for order an even integer 3053 (mul (power 2 order) 3054 (power '$%pi (div 1 2)) 3055 (inv (take '(%gamma) (div (sub 1 order) 2))) 3056 (list '(mqapply simp) 3057 (list '($%f array simp) 1 1) 3058 (list '(mlist) (div order -2)) 3059 (list '(mlist) '((rat simp) 1 2)) 3060 (mul arg arg)))) 3061 3062 ((and (maxima-integerp order) 3063 (moddp order)) 3064 ;; Transform to 1F1 for order an odd integer 3065 (mul -2 arg 3066 (power 2 order) 3067 (power '$%pi '((rat simp) 1 2)) 3068 (inv (take '(%gamma) (div order -2))) 3069 (list '(mqapply simp) 3070 (list '($%f simp array) 1 1) 3071 (list '(mlist simp) (div (sub 1 order) 2)) 3072 (list '(mlist simp) '((rat simp) 3 2)) 3073 (mul arg arg)))) 3074 (t 3075 ;; The general case, transform to 2F0 3076 ;; For this case we have no Laplace transform. 3077 (mul (power (mul 2 arg) order) 3078 (list '(mqapply simp) 3079 (list '($%f array simp) 2 0) 3080 (list '(mlist simp) (div order 2) (div (sub 1 order) 2)) 3081 (list '(mlist simp)) 3082 (div -1 (mul arg arg))))))) 3083 3084;;; Hypergeometric representation of the Exponential Integral Si 3085;;; Si(z) = z*1F2([1/2],[3/2,3/2],-z^2/4) 3086(defun expintegral_si-to-hypergeometric (arg) 3087 (mul arg 3088 (list '(mqapply simp) 3089 (list '($%f array simp) 1 2) 3090 (list '(mlist simp) '((rat simp) 1 2)) 3091 (list '(mlist simp) '((rat simp) 3 2) '((rat simp) 3 2)) 3092 (div (mul -1 arg arg) 4)))) 3093 3094;;; Hypergeometric representation of the Exponential Integral Shi 3095;;; Shi(z) = z*1F2([1/2],[3/2,3/2],z^2/4) 3096(defun expintegral_shi-to-hypergeometric (arg) 3097 (mul arg 3098 (list '(mqapply simp) 3099 (list '($%f simp array) 1 2) 3100 (list '(mlist simp) '((rat simp) 1 2)) 3101 (list '(mlist simp) '((rat simp) 3 2) '((rat simp) 3 2)) 3102 (div (mul arg arg) 4)))) 3103 3104;;; Hypergeometric representation of the Exponential Integral Ci 3105;;; Ci(z) = -z^2/4*2F3([1,1],[2,2,3/2],-z^2/4)+log(z)+%gamma 3106(defun expintegral_ci-to-hypergeometric (arg) 3107 (add (mul (div (mul -1 arg arg) 4) 3108 (list '(mqapply simp) 3109 (list '($%f simp array) 2 3) 3110 (list '(mlist simp) 1 1) 3111 (list '(mlist simp) 2 2 '((rat simp) 3 2)) 3112 (div (mul -1 arg arg) 4))) 3113 (take '(%log) arg) 3114 '$%gamma)) 3115 3116;;; Hypergeometric representation of the Exponential Integral Chi 3117;;; Chi(z) = z^2/4*2F3([1,1],[2,2,3/2],z^2/4)+log(z)+%gamma 3118(defun expintegral_chi-to-hypergeometric (arg) 3119 (add (mul (div (mul arg arg) 4) 3120 (list '(mqapply simp) 3121 (list '($%f array simp) 2 3) 3122 (list '(mlist simp) 1 1) 3123 (list '(mlist simp) 2 2 '((rat simp) 3 2)) 3124 (div (mul arg arg) 4))) 3125 (take '(%log) arg) 3126 '$%gamma)) 3127 3128;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 3129;;; 3130;;; EXPERTS ON LAPLACE TRANSFORMS 3131;;; 3132;;; LT<foo> functions are various experts on Laplace transforms of the 3133;;; function <foo>. The expression being transformed is 3134;;; r*<foo>(args). The first arg of each expert is r, The remaining 3135;;; args are the arg(s) and/or parameters. 3136;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 3137 3138;; Laplace transform of one Bessel J 3139(defun lt1j (rest arg index) 3140 (lt-ltp 'onej rest arg index)) 3141 3142;; Laplace transform of two Bessel J functions. 3143;; The argument of each must be the same, but the orders may be different. 3144(defun lt2j (rest arg1 arg2 index1 index2) 3145 (cond ((not (equal arg1 arg2)) 3146 (setq *hyp-return-noun-flag* 'product-of-bessel-with-different-args)) 3147 (t 3148 ;; Call lt-ltp two transform and integrate two Bessel J functions. 3149 (lt-ltp 'twoj rest arg1 (list 'list index1 index2))))) 3150 3151;; Laplace transform of a square of a Bessel J function 3152(defun lt1j^2 (rest arg index) 3153 (cond ((alike1 index '((rat) -1 2)) 3154 ;; Special case: Laplace transform of bessel_j(v,arg)^2, v = -1/2. 3155 ;; For this case the algorithm for the product of two Bessel functions 3156 ;; does not work. Call the integrator with the hypergeometric 3157 ;; representation: 2/%pi/arg*1/2*(1+0F1([], [1/2], -arg*arg)). 3158 (sendexec (mul (div 2 '$%pi) 3159 (inv arg) 3160 rest) 3161 (add '((rat simp) 1 2) 3162 (mul '((rat simp) 1 2) 3163 (list '(mqapply simp) 3164 (list '($%f simp array) 1 0) 3165 (list '(mlist simp) ) 3166 (list '(mlist simp) '((rat simp) 1 2)) 3167 (mul -1 (mul arg arg))))))) 3168 (t 3169 (lt-ltp 'twoj rest arg (list 'list index index))))) 3170 3171;; Laplace transform of Incomplete Gamma function 3172(defun lt1gamma-incomplete-lower (rest arg1 arg2) 3173 (lt-ltp 'gamma-incomplete-lower rest arg2 arg1)) 3174 3175;; Laplace transform of Whittaker M function 3176(defun lt1m (r a i1 i2) 3177 (lt-ltp 'onem r a (list i1 i2))) 3178 3179;; Laplace transform of Jacobi function 3180(defun lt1p (r a i1 i2) 3181 (lt-ltp 'hyp-onep r a (list i1 i2))) 3182 3183;; Laplace transform of Associated Legendre function of the second kind 3184(defun lt1q (r a i1 i2) 3185 (lt-ltp 'oneq r a (list i1 i2))) 3186 3187;; Laplace transform of Erf function 3188(defun lt1erf (rest arg) 3189 (lt-ltp 'onerf rest arg nil)) 3190 3191;; Laplace transform of Log function 3192(defun lt1log (rest arg) 3193 (lt-ltp 'onelog rest arg nil)) 3194 3195;; Laplace transform of Complete elliptic integral of the first kind 3196(defun lt1kelliptic (rest arg) 3197 (lt-ltp 'onekelliptic rest arg nil)) 3198 3199;; Laplace transform of Complete elliptic integral of the second kind 3200(defun lt1e (rest arg) 3201 (lt-ltp 'onee rest arg nil)) 3202 3203;; Laplace transform of Struve H function 3204(defun lt1hstruve (rest arg1 index1) 3205 (lt-ltp 'hs rest arg1 index1)) 3206 3207;; Laplace transform of Struve L function 3208(defun lt1lstruve (rest arg1 index1) 3209 (lt-ltp 'hl rest arg1 index1)) 3210 3211;; Laplace transform of Lommel s function 3212(defun lt1s (rest arg1 index1 index2) 3213 (lt-ltp 's rest arg1 (list index1 index2))) 3214 3215;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 3216;;; 3217;;; TRANSFORM TO HYPERGEOMETRIC FUNCTION AND DO THE INTEGRATION 3218;;; 3219;;; FLG = special function we're transforming 3220;;; REST = other stuff 3221;;; ARG = arg of special function 3222;;; INDEX = index of special function. 3223;;; 3224;;; So we're transforming REST*FLG(INDEX, ARG). 3225;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 3226 3227(defun lt-ltp (flg rest arg index) 3228 (let (l l1) 3229 (when (and (listp index) 3230 (eq (car index) 'list)) 3231 (setq index (list (cadr index) (caddr index)))) 3232 (cond ((setq l 3233 (m2-d*x^m*%e^a*x 3234 ($factor (mul rest 3235 (car (setq l1 (ref flg index arg))))) 3236 *var* *par*)) 3237 ;; Convert the special function to a hypgergeometric 3238 ;; function. L1 is the special function converted to the 3239 ;; hypergeometric function. d*x^m*%e^a*x looks for that 3240 ;; factor in the expanded form. 3241 (%$etest l l1)) 3242 (t 3243 ;; We currently don't know how to handle this yet. 3244 ;; We add the return of a noun form. 3245 (setq *hyp-return-noun-flag* 'other-ca-later))))) 3246 3247;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 3248;; Dispatch function to convert the given function to a hypergeometric 3249;; function. 3250;; 3251;; The first arg is a symbol naming the function; the last is the 3252;; argument to the function. The second arg is the index (or list of 3253;; indices) to the function. Not used if the function doesn't have 3254;; any indices 3255;; 3256;; The result is a list of 2 elements: The first element is a 3257;; multiplier; the second, the hypergeometric function itself. 3258;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 3259 3260(defun ref (flg index arg) 3261 (case flg 3262 (onej (j1tf index arg)) 3263 (twoj (j2tf (car index) (cadr index) arg)) 3264 (hs (hstf index arg)) 3265 (hl (lstf index arg)) 3266 (s (stf (car index) (cadr index) arg)) 3267 (onerf (erftf arg)) 3268 (onelog (logtf arg)) 3269 (onekelliptic (kelliptictf arg)) ; elliptic_kc 3270 (onee (etf arg)) ; elliptic_ec 3271 (onem (mtf (car index) (cadr index) arg)) 3272 (hyp-onep (ptf (car index) (cadr index) arg)) 3273 (oneq (qtf (car index) (cadr index) arg)) 3274 (gamma-incomplete-lower (gamma-incomplete-lower-tf index arg)) 3275 (onepjac (pjactf (car index) (cadr index) (caddr index) arg)) 3276 (asin (asintf arg)) 3277 (atan (atantf arg)) 3278 (f 3279 ;; Transform %f to internal representation FPQ 3280 (list 1 (ref-fpq (rest (car index)) (rest (cadr index)) arg))))) 3281 3282;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 3283;;; 3284;;; TRANSFORM FUNCTION IN TERMS OF HYPERGEOMETRIC FUNCTION 3285;;; 3286;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 3287 3288;; Create a hypergeometric form that we recognize. The function is 3289;; %f[n,m](p; q; arg). We represent this as a list of the form 3290;; (fpq (<length p> <length q>) <p> <q> <arg>) 3291 3292(defun ref-fpq (p q arg) 3293 (list 'fpq (list (length p) (length q)) 3294 p q arg)) 3295 3296;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 3297 3298;; Whittaker M function in terms of hypergeometric function 3299;; 3300;; A&S 13.1.32: 3301;; 3302;; M[k,u](z) = exp(-z/2)*z^(1/2+u)*M(1/2+u-k,1+2*u,z) 3303 3304(defun mtf (i1 i2 arg) 3305 (list (mul (power arg (add i2 '((rat simp) 1 2))) 3306 (power '$%e (div arg -2))) 3307 (ref-fpq (list (add '((rat simp) 1 2) i2 (mul -1 i1))) 3308 (list (add i2 i2 1)) 3309 arg))) 3310 3311;; Jacobi P in terms of hypergeometric function 3312;; 3313;; A&S 15.4.6: 3314;; 3315;; F(-n,a+1+b+n; a+1; x) = n!/poch(a+1,n)*jacobi_p(n,a,b,1-2*x) 3316;; 3317;; jacobi_p(n,a,b,x) = poch(a+1,n)/n!*F(-n,a+1+b+n; a+1; (1-x)/2) 3318;; = gamma(a+n+1)/gamma(a+1)/n!*F(-n,a+1+b+n; a+1; (1-x)/2) 3319;; 3320;; We have a problem: 3321;; We transform the argument x -> (1-x)/2. But an argument with a constant 3322;; part is not integrable with the implemented algorithm for the hypergeometric 3323;; function. It might be possible to get a result for an argument y=1-2*x 3324;; with x=t^-2. But for this case the routine lt-ltp fails. The routine 3325;; recognize a constant term in the argument, but does not take into account 3326;; that the constant term might vanish, when we transform to a hypergeometric 3327;; function. 3328;; Because of this the Laplace transform for the following functions does not 3329;; work too: Legendre P, Chebyshev T, Chebyshev U, and Gegenbauer. 3330 3331(defun pjactf (n a b x) 3332 (list (mul (take '(%gamma) (add n a 1)) 3333 (inv (take '(%gamma) (add a 1))) 3334 (inv (take '(%gamma) (add n 1)))) 3335 (ref-fpq (list (mul -1 n) (add n a b 1)) 3336 (list (add a 1)) 3337 (sub '((rat simp) 1 2) (div x 2))))) 3338 3339;; ArcSin in terms of hypergeometric function 3340;; 3341;; A&S 15.1.6: 3342;; 3343;; F(1/2,1/2; 3/2; z^2) = asin(z)/z 3344;; 3345;; asin(z) = z*F(1/2,1/2; 3/2; z^2) 3346 3347(defun asintf (arg) 3348 (let ((inv2 '((rat simp) 1 2))) 3349 (list arg 3350 (ref-fpq (list inv2 inv2) 3351 (list '((rat simp) 3 2)) 3352 (mul arg arg))))) 3353 3354;; ArcTan in terms of hypergeometric function 3355;; 3356;; A&S 15.1.5 3357;; 3358;; F(1/2,1; 3/2; -z^2) = atan(z)/z 3359;; 3360;; atan(z) = z*F(1/2,1; 3/2; -z^2) 3361 3362(defun atantf (arg) 3363 (list arg 3364 (ref-fpq (list '((rat simp) 1 2) 1) 3365 (list '((rat simp) 3 2)) 3366 (mul -1 arg arg)))) 3367 3368;; Associated Legendre function P in terms of hypergeometric function 3369;; 3370;; A&S 8.1.2 3371;; 3372;; assoc_legendre_p(v,u,z) = ((z+1)/(z-2))^(u/2)/gamma(1-u)*F(-v,v+1;1-u,(1-z)/2) 3373;; 3374;; FIXME: What about the branch cut? 8.1.2 is for z not on the real 3375;; line with -1 < z < 1. 3376 3377(defun ptf (n m z) 3378 (list (mul (inv (take '(%gamma) (sub 1 m))) 3379 (power (div (add z 1) 3380 (sub z 1)) 3381 (div m 2))) 3382 (ref-fpq (list (mul -1 n) (add n 1)) 3383 (list (sub 1 m)) 3384 (sub '((rat simp) 1 2) (div z 2))))) 3385 3386;; Associated Legendre function Q in terms of hypergeometric function 3387;; 3388;; A&S 8.1.3: 3389;; 3390;; assoc_legendre_q(v,u,z) 3391;; = exp(%i*u*%pi)*2^(-v-1)*sqrt(%pi) * 3392;; gamma(v+u+1)/gamma(v+3/2)*z^(-v-u-1)*(z^2-1)^(u/2) * 3393;; F(1+v/2+u/2, 1/2+v/2+u/2; v+3/2; 1/z^2) 3394;; 3395;; FIXME: What about the branch cut? 3396 3397(defun qtf (n m z) 3398 (list (mul (power '$%e (mul m '$%pi '$%i)) 3399 (power '$%pi '((rat simp) 1 2)) 3400 (take '(%gamma) (add m n 1)) 3401 (power 2 (sub -1 n)) 3402 (inv (take '(%gamma) (add n '((rat simp) 3 2)))) 3403 (power z (mul -1 (add m n 1))) 3404 (power (sub (mul z z) 1) (div m 2))) 3405 (ref-fpq (list (div (add m n 1) 2) (div (add m n 2) 2)) 3406 (list (add n '((rat simp) 3 2))) 3407 (power z -2)))) 3408 3409;; Incomplete lower Gamma in terms of hypergeometric function 3410;; 3411;; A&S 13.6.10: 3412;; 3413;; M(a,a+1,-x) = a*x^(-a)*gamma_incomplete_lower(a,x) 3414;; 3415;; gamma_incomplete_lower(a,x) = x^a/a*M(a,a+1,-x) 3416 3417(defun gamma-incomplete-lower-tf (a x) 3418 (list (mul (inv a) (power x a)) 3419 (ref-fpq (list a) 3420 (list (add a 1)) 3421 (mul -1 x)))) 3422 3423;; Complete elliptic K in terms of hypergeometric function 3424;; 3425;; A&S 17.3.9 3426;; 3427;; K(k) = %pi/2*2F1(1/2,1/2; 1; k^2) 3428 3429(defun kelliptictf (k) 3430 (let ((inv2 '((rat simp) 1 2))) 3431 (list (mul inv2 '$%pi) 3432 (ref-fpq (list inv2 inv2) 3433 (list 1) 3434 (mul k k))))) 3435 3436;; Complete elliptic E in terms of hypergeometric function 3437;; 3438;; A&S 17.3.10 3439;; 3440;; E(k) = %pi/2*2F1(-1/2,1/2;1;k^2) 3441 3442(defun etf (k) 3443 (let ((inv2 '((rat simp) 1 2))) 3444 (list (mul inv2 '$%pi) 3445 (ref-fpq (list (mul -1 inv2) inv2) 3446 (list 1) 3447 (mul k k))))) 3448 3449;; erf in terms of hypgeometric function. 3450;; 3451;; A&S 7.1.21 gives 3452;; 3453;; erf(z) = 2*z/sqrt(%pi)*M(1/2,3/2,-z^2) 3454;; = 2*z/sqrt(%pi)*exp(-z^2)*M(1,3/2,z^2) 3455 3456(defun erftf (arg) 3457 (list (mul 2 arg (power '$%pi '((rat simp) -1 2))) 3458 (ref-fpq (list '((rat simp) 1 2)) 3459 (list '((rat simp) 3 2)) 3460 (mul -1 arg arg)))) 3461 3462;; log in terms of hypergeometric function 3463;; 3464;; We know from A&S 15.1.3 that 3465;; 3466;; F(1,1;2;z) = -log(1-z)/z. 3467;; 3468;; So log(z) = (z-1)*F(1,1;2;1-z) 3469 3470(defun logtf (arg) 3471 (list (sub arg 1) 3472 (ref-fpq (list 1 1) 3473 (list 2) 3474 (sub 1 arg)))) 3475 3476;; Bessel J function expressed as a hypergeometric function. 3477;; 3478;; A&S 9.1.10: 3479;; inf 3480;; bessel_j(v,z) = (z/2)^v*sum (-z^2/4)^k/k!/gamma(v+k+1) 3481;; k=0 3482;; 3483;; = (z/2)^v/gamma(v+1)*sum 1/poch(v+1,k)*(-z^2/4)^k/k! 3484;; 3485;; = (z/2)^v/gamma(v+1) * 0F1(; v+1; -z^2/4) 3486 3487(defun j1tf (v z) 3488 (list (mul (inv (power 2 v)) 3489 (power z v) 3490 (inv (take '(%gamma) (add v 1)))) 3491 (ref-fpq nil 3492 (list (add v 1)) 3493 (mul '((rat simp) -1 4) (power z 2))))) 3494 3495;; Product of 2 Bessel J functions in terms of hypergeometric function 3496;; 3497;; See Y. L. Luke, formula 39, page 216: 3498;; 3499;; bessel_j(u,z)*bessel_j(v,z) 3500;; = (z/2)^(u+v)/gamma(u+1)/gamma(v+1) * 3501;; 2F3((u+v+1)/2, (u+v+2)/2; u+1, v+1, u+v+1; -z^2) 3502 3503(defun j2tf (n m arg) 3504 (list (mul (inv (take '(%gamma) (add n 1))) 3505 (inv (take '(%gamma) (add m 1))) 3506 (inv (power 2 (add n m))) 3507 (power arg (add n m))) 3508 (ref-fpq (list (add '((rat simp) 1 2) (div n 2) (div m 2)) 3509 (add 1 (div n 2) (div m 2))) 3510 (list (add 1 n) (add 1 m) (add 1 n m)) 3511 (mul -1 (power arg 2))))) 3512 3513;; Struve H function in terms of hypergeometric function. 3514;; 3515;; A&S 12.1.2 gives the following series for the Struve H function: 3516;; 3517;; inf 3518;; H[v](z) = (z/2)^(v+1)*sum (-1)^k*(z/2)^(2*k)/gamma(k+3/2)/gamma(k+v+3/2) 3519;; k=0 3520;; 3521;; We can write this in the form 3522;; 3523;; H[v](z) = 2/sqrt(%pi)*(z/2)^(v+1)/gamma(v+3/2) 3524;; 3525;; inf 3526;; * sum n!/poch(3/2,n)/poch(v+3/2,n)*(-z^2/4)^n/n! 3527;; n=0 3528;; 3529;; = 2/sqrt(%pi)*(z/2)^(v+1)/gamma(v+3/2) * 1F2(1;3/2,v+3/2;(-z^2/4)) 3530;; 3531;; See also A&S 12.1.21. 3532 3533(defun hstf (v z) 3534 (let ((d32 (div 3 2))) 3535 (list (mul (power (div z 2)(add v 1)) 3536 (inv (take '(%gamma) d32)) 3537 (inv (take '(%gamma) (add v d32)))) 3538 (ref-fpq (list 1) 3539 (list d32 (add v d32)) 3540 (mul '((rat simp) -1 4) z z))))) 3541 3542;; Struve L function in terms of hypergeometric function 3543;; 3544;; A&S 12.2.1: 3545;; 3546;; L[v](z) = -%i*exp(-v*%i*%pi/2)*H[v](%i*z) 3547;; 3548;; This function computes exactly this way. (But why is %i written as 3549;; exp(%i*%pi/2) instead of just %i) 3550;; 3551;; A&S 12.2.1 gives the series expansion as 3552;; 3553;; inf 3554;; L[v](z) = (z/2)^(v+1)*sum (z/2)^(2*k)/gamma(k+3/2)/gamma(k+v+3/2) 3555;; k=0 3556;; 3557;; It's quite easy to derive 3558;; 3559;; L[v](z) = 2/sqrt(%pi)*(z/2)^(v+1)/gamma(v+3/2) * 1F2(1;3/2,v+3/2;(z^2/4)) 3560 3561(defun lstf (v z) 3562 (let ((d32 (div 3 2))) 3563 (list (mul (power (div z 2) (add v 1)) 3564 (inv (take '(%gamma) d32)) 3565 (inv (take '(%gamma) (add v d32)))) 3566 (ref-fpq (list 1) 3567 (list d32 (add v d32)) 3568 (mul '((rat simp) 1 4) z z))))) 3569 3570;; Lommel s function in terms of hypergeometric function 3571;; 3572;; See Y. L. Luke, p 217, formula 1 3573;; 3574;; s(u,v,z) = z^(u+1)/(u-v+1)/(u+v+1)*1F2(1; (u-v+3)/2, (u+v+3)/2; -z^2/4) 3575 3576(defun stf (m n z) 3577 (list (mul (power z (add m 1)) 3578 (inv (sub (add m 1) n)) 3579 (inv (add m n 1))) 3580 (ref-fpq (list 1) 3581 (list (div (sub (add m 3) n) 2) 3582 (div (add m n 3) 2)) 3583 (mul '((rat simp) -1 4) z z)))) 3584 3585;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 3586;;; 3587;;; Algorithm 3: Laplace transform of a hypergeometric function 3588;;; 3589;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 3590 3591(defun %$etest (l l1) 3592 (prog(a q) 3593 (setq q (cdras 'q l)) 3594 (cond ((equal q 1)(setq a 0)(go loop))) 3595 (setq a (cdras 'a l)) 3596 loop 3597 (return (substl (sub *par* a) 3598 *par* 3599 (execf19 l (cadr l1)))))) 3600 3601(defun execf19 (l1 l2) 3602 (prog(ans) 3603 (setq ans (execargmatch (car (cddddr l2)))) 3604 (cond ((eq (car ans) 'dionimo) 3605 (return (dionarghyp l1 l2 (cadr ans))))) 3606 ;; We specialized the pattern for the argument. Now the test fails 3607 ;; correctly and we have to add the return of a noun form. 3608 (return 3609 (setq *hyp-return-noun-flag* 'next-for-other-args)))) 3610 3611;; Executive for recognizing the sort of argument to the 3612;; hypergeometric function. We look to see if the arg is of the form 3613;; a*x^m + c. Return a list of 'dionimo (what does that mean?) and 3614;; the match. 3615 3616(defun execargmatch (arg) 3617 (prog(l1) 3618 (cond ((setq l1 (m2-a*x^m+c ($factor arg) *var*)) 3619 (return (list 'dionimo l1)))) 3620 (cond ((setq l1 (m2-a*x^m+c ($expand arg) *var*)) 3621 (return (list 'dionimo l1)))) 3622 ;; The return value has to be a list. 3623 (return (list 'other-case-args-to-follow)))) 3624 3625;; We have hypergeometric function whose arg looks like a*x^m+c. L1 3626;; matches the d*x^m... part, L2 is the hypergeometric function and 3627;; arg is the match for a*x^m+c. 3628 3629(defun dionarghyp (l1 l2 arg) 3630 (prog(a m c) 3631 (setq a (cdras 'a arg) 3632 m (cdras 'm arg) 3633 c (cdras 'c arg)) 3634 (cond 3635 ((and (integerp m) ; The power of the argument has to be 3636 (plusp m) ; a positive integer. 3637 (equal 0 c)) ; No const term to the argument. 3638 (return (f19cond a m l1 l2))) 3639 (t 3640 (return 3641 (setq *hyp-return-noun-flag* 'prop4-and-other-cases-to-follow)))))) 3642 3643(defun f19cond (a m l1 l2) 3644 (prog(p q s d) 3645 (setq p (caadr l2) 3646 q (cadadr l2) 3647 s (cdras 'm l1) 3648 d (cdras 'd l1) 3649 l1 (caddr l2) 3650 l2 (cadddr l2)) 3651 ;; At this point, we have the function d*x^s*%f[p,q](l1, l2, (a*t)^m). 3652 ;; Check to see if Formula 19, p 220 applies. 3653 (cond ((and (not (eq ($asksign (sub (add p m -1) q)) '$pos)) 3654 (eq ($asksign (add s 1)) '$pos)) 3655 (return (mul d 3656 (f19p220-simp (add s 1) 3657 l1 3658 l2 3659 a 3660 m))))) 3661 (return (setq *hyp-return-noun-flag* 3662 'failed-on-f19cond-multiply-the-other-cases-with-d)))) 3663 3664;; Table of Laplace transforms, p 220, formula 19: 3665;; 3666;; If m + k <= n + 1, and Re(s) > 0, the Laplace transform of 3667;; 3668;; t^(s-1)*%f[m,n]([a1,...,am],[p1,...,pn],(c*t)^k) 3669;; is 3670;; 3671;; gamma(s)/p^s*%f[m+k,n]([a1,...,am,s/k,(s+1)/k,...,(s+k-1)/k],[p1,...,pm],(k*c/p)^k) 3672;; 3673;; with Re(p) > 0 if m + k <= n, Re(p+k*c*exp(2*%pi*%i*r/k)) > 0 for r 3674;; = 0, 1,...,k-1, if m + k = n + 1. 3675;; 3676;; The args below are s, [a's], [p's], c^k, k. 3677 3678(defun f19p220-simp (s l1 l2 cf k) 3679 (mul (take '(%gamma) s) 3680 (inv (power *par* s)) 3681 (hgfsimp-exec (append l1 (addarglist s k)) 3682 l2 3683 (mul cf 3684 (power k k) 3685 (power (inv *par*) k))))) 3686 3687;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 3688 3689;; Return a list of s/k, (s+1)/k, ..., (s+|k|-1)/k 3690(defun addarglist (s k) 3691 (let ((abs-k (abs k)) 3692 (res '())) 3693 (dotimes (n abs-k) 3694 (push (div (add s n) k) res)) 3695 (nreverse res))) 3696 3697;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 3698 3699;;; Pattern for the Laplace transform of the hypergeometric function 3700 3701;; Match d*x^m*%e^(a*x). If we match, Q is the e^(a*x) part, A is a, 3702;; M is M, and D is d. 3703(defun m2-d*x^m*%e^a*x (expr var par) 3704 (m2 expr 3705 `((mtimes) 3706 ((coefftt) (d free2 ,var ,par)) 3707 ((mexpt) (x alike1 ,var) (m free2 ,var ,par)) 3708 ((mexpt) 3709 (q expor1p) 3710 ((mtimes) 3711 ((coefftt) (a free2 ,var ,par)) 3712 (x alike1 ,var)))))) 3713 3714;; Match f(x)+c 3715(defun m2-f+c (expr var) 3716 (m2 expr 3717 `((mplus) 3718 ((coeffpt) (f has ,var)) 3719 ((coeffpp) (c free ,var))))) 3720 3721;; Match a*x^m+c. 3722;; The pattern was too general. We match also a*t^2+b*t. But that's not correct. 3723(defun m2-a*x^m+c (expr var) 3724 (m2 expr 3725 `((mplus) 3726 ((coefft) ; more special (not coeffpt) 3727 (a free ,var) 3728 ((mexpt) (x alike1 ,var) (m free-not-zero-p ,var))) 3729 ((coeffpp) (c free ,var))))) 3730 3731;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 3732;;; 3733;;; Algorithm 4: SPECIAL HANDLING OF Bessel Y for an integer order 3734;;; 3735;;; This is called for one Bessel Y function, when the order is an integer. 3736;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 3737 3738(defun lty (rest arg index) 3739 (prog(l) 3740 (cond ((setq l (m2-d*x^m*%e^a*x rest *var* *par*)) 3741 (return (execfy l arg index)))) 3742 (return (setq *hyp-return-noun-flag* 'fail-in-lty)))) 3743 3744(defun execfy (l arg index) 3745 (prog(ans) 3746 (setq ans (execargmatch arg)) 3747 (cond ((eq (car ans) 'dionimo) 3748 (return (dionarghyp-y l index (cadr ans))))) 3749 (return (setq *hyp-return-noun-flag* 'fail-in-execfy)))) 3750 3751(defun dionarghyp-y (l index arg) 3752 (prog (a m c) 3753 (setq a (cdras 'a arg) 3754 m (cdras 'm arg) 3755 c (cdras 'c arg)) 3756 (cond ((and (zerp c) (equal m 1)) 3757 (let ((ans (f2p105v2cond a l index))) 3758 (unless (symbolp ans) 3759 (return ans))))) 3760 (cond ((and (zerp c) (equal m '((rat simp) 1 2))) 3761 (let ((ans (f50cond a l index))) 3762 (unless (symbolp ans) 3763 (return ans))))) 3764 (return (setq *hyp-return-noun-flag* 'fail-in-dionarghyp-y)))) 3765 3766;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 3767;;; 3768;;; Algorithm 4.1: Laplace transform of t^n*bessel_y(v,a*t) 3769;;; v is an integer and n>=v 3770;;; 3771;;; Table of Integral Transforms 3772;;; 3773;;; Volume 2, p 105, formula 2 is a formula for the Y-transform of 3774;;; 3775;;; f(x) = x^(u-3/2)*exp(-a*x) 3776;;; 3777;;; where the Y-transform is defined by 3778;;; 3779;;; integrate(f(x)*bessel_y(v,x*y)*sqrt(x*y), x, 0, inf) 3780;;; 3781;;; which is 3782;;; 3783;;; -2/%pi*gamma(u+v)*sqrt(y)*(y^2+a^2)^(-u/2) 3784;;; *assoc_legendre_q(u-1,-v,a/sqrt(y^2+a^2)) 3785;;; 3786;;; with a > 0, Re u > |Re v|. 3787;;; 3788;;; In particular, with a slight change of notation, we have 3789;;; 3790;;; integrate(x^(u-1)*exp(-p*x)*bessel_y(v,a*x)*sqrt(a), x, 0, inf) 3791;:; 3792;;; which is the Laplace transform of x^(u-1/2)*bessel_y(v,x). 3793;;; 3794;;; Thus, the Laplace transform is 3795;;; 3796;;; -2/%pi*gamma(u+v)*sqrt(a)*(a^2+p^2)^(-u/2) 3797;;; *assoc_legendre_q(u-1,-v,p/sqrt(a^2+p^2)) 3798;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 3799 3800(defun f2p105v2cond (a l index) 3801 (prog (d m) 3802 (setq d (cdras 'd l) ; contains constant part of integrand 3803 m (cdras 'm l)) 3804 (setq m (add m 1.)) 3805 (cond ((eq ($asksign ($realpart (sub m index))) '$pos) 3806 (return (mul d (f2p105v2cond-simp m index a))))) 3807 (return (setq *hyp-return-noun-flag* 'fail-in-f2p105v2cond)))) 3808 3809(defun f2p105v2cond-simp (m v a) 3810 (mul -2 3811 (power '$%pi -1) 3812 (take '(%gamma) (add m v)) 3813 (power (add (mul a a) (mul *par* *par*)) 3814 (mul -1 '((rat simp) 1 2) m)) 3815 ;; Call Associated Legendre Q function, which simplifies accordingly. 3816 ;; We have to do a Maxima function call, because $assoc_legendre_q is 3817 ;; not in Maxima core and has to be autoloaded. 3818 (mfuncall '$assoc_legendre_q 3819 (sub m 1) 3820 (mul -1 v) 3821 (mul *par* 3822 (power (add (mul a a) (mul *par* *par*)) 3823 '((rat simp) -1 2)))))) 3824 3825;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 3826;;; 3827;;; Algorithm 4.2: Laplace transform of t^n*bessel_y(v, a*sqrt(t)) 3828;;; 3829;;; Table of Integral Transforms 3830;;; 3831;;; p. 188, formula 50: 3832;;; 3833;;; t^(u-1/2)*bessel_y(2*v,2*sqrt(a)*sqrt(t)) 3834;;; -> a^(-1/2)*p^(-u)*exp(-a/2/p) 3835;;; * [tan((u-v)*%pi)*gamma(u+v+1/2)/gamma(2*v+1)*M[u,v](a/p) 3836;;; -sec((u-v)*%pi)*W[u,v](a/p)] 3837;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 3838 3839(defun f50cond (a l v) 3840 (prog (d m) 3841 (setq d (cdras 'd l) 3842 m (add (cdras 'm l) '((rat simp) 1 2)) 3843 v (div v 2)) 3844 (cond 3845 ((and (eq ($asksign ($realpart (add m v '((rat simp) 1 2)))) '$pos) 3846 (eq ($asksign ($realpart (sub (add m '((rat simp) 1 2)) v))) '$pos) 3847 (not (maxima-integerp (mul (sub (add m m) (add v v 1)) 3848 '((rat simp) 1 2))))) 3849 (setq a (mul a a '((rat simp) 1 4))) 3850 (return (f50p188-simp d m v a)))) 3851 (return (setq *hyp-return-noun-flag* 'fail-in-f50cond)))) 3852 3853(defun f50p188-simp (d u v a) 3854 (mul d 3855 (power a '((rat simp) -1 2)) 3856 (power *par* (mul -1 u)) 3857 (power '$%e (div a (mul -2 *par*))) 3858 (sub (mul (take '(%tan) (mul '$%pi (sub u v))) 3859 (take '(%gamma) (add u v '((rat simp) 1 2))) 3860 (inv (take '(%gamma) (add v v 1))) 3861 (mwhit (div a *par*) u v)) 3862 (mul (take '(%sec) (mul '$%pi (sub u v))) 3863 (wwhit (div a *par*) u v))))) 3864 3865;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 3866