1@menu 2* Introduction to orthogonal polynomials:: 3* Functions and Variables for orthogonal polynomials:: 4@end menu 5 6@node Introduction to orthogonal polynomials, Functions and Variables for orthogonal polynomials, orthopoly-pkg, orthopoly-pkg 7@section Introduction to orthogonal polynomials 8 9@code{orthopoly} is a package for symbolic and numerical evaluation of 10several kinds of orthogonal polynomials, including Chebyshev, 11Laguerre, Hermite, Jacobi, Legendre, and ultraspherical (Gegenbauer) 12polynomials. Additionally, @code{orthopoly} includes support for the spherical Bessel, 13spherical Hankel, and spherical harmonic functions. 14 15For the most part, @code{orthopoly} follows the conventions of Abramowitz and Stegun 16@i{Handbook of Mathematical Functions}, Chapter 22 (10th printing, December 1972); 17additionally, we use Gradshteyn and Ryzhik, 18@i{Table of Integrals, Series, and Products} (1980 corrected and 19enlarged edition), and Eugen Merzbacher @i{Quantum Mechanics} (2nd edition, 1970). 20 21@c INSTALLATION INSTRUCTIONS NO LONGER RELEVANT 22@c BUT MAYBE SOME OF THESE FILES SHOULD BE MENTIONED IN ANOTHER CONTEXT 23@c This will create a directory @code{orthopoly_x} (again x is the release 24@c identifier) that contains the source file @code{orthopoly.lisp}, user 25@c documentation in html and texi formats, a sample maxima initialization file 26@c @code{orthopoly-init.lisp}, a README file, a testing routine 27@c @code{test_orthopoly.mac}, and two demonstration files. 28 29@c Start Maxima and compile orthopoly. To do this, use the command 30@c 31@c (c1) compile_file("orthopoly.lisp"); 32 33Barton Willis of the University of Nebraska at Kearney (UNK) wrote 34the @code{orthopoly} package and its documentation. The package 35is released under the GNU General Public License (GPL). 36 37@opencatbox 38@category{Orthogonal polynomials} 39@category{Share packages} 40@category{Package orthopoly} 41@closecatbox 42 43@subsection Getting Started with orthopoly 44 45@code{load ("orthopoly")} loads the @code{orthopoly} package. 46 47To find the third-order Legendre polynomial, 48 49@c ===beg=== 50@c legendre_p (3, x); 51@c ===end=== 52@example 53(%i1) legendre_p (3, x); 54 3 2 55 5 (1 - x) 15 (1 - x) 56(%o1) - ---------- + ----------- - 6 (1 - x) + 1 57 2 2 58@end example 59 60To express this as a sum of powers of @var{x}, apply @var{ratsimp} or @var{rat} 61to the result. 62 63@c CONTINUING PREVIOUS EXAMPLE HERE 64@c ===beg=== 65@c [ratsimp (%), rat (%)]; 66@c ===end=== 67@example 68(%i2) [ratsimp (%), rat (%)]; 69 3 3 70 5 x - 3 x 5 x - 3 x 71(%o2)/R/ [----------, ----------] 72 2 2 73@end example 74 75Alternatively, make the second argument to @code{legendre_p} (its ``main'' variable) 76a canonical rational expression (CRE). 77 78@c ===beg=== 79@c legendre_p (3, rat (x)); 80@c ===end=== 81@example 82(%i1) legendre_p (3, rat (x)); 83 3 84 5 x - 3 x 85(%o1)/R/ ---------- 86 2 87@end example 88 89For floating point evaluation, @code{orthopoly} uses a running error analysis 90to estimate an upper bound for the error. For example, 91 92@c ===beg=== 93@c jacobi_p (150, 2, 3, 0.2); 94@c ===end=== 95@example 96(%i1) jacobi_p (150, 2, 3, 0.2); 97(%o1) interval(- 0.062017037936715, 1.533267919277521E-11) 98@end example 99 100Intervals have the form @code{interval (@var{c}, @var{r})}, where @var{c} is the 101center and @var{r} is the radius of the interval. Since Maxima 102does not support arithmetic on intervals, in some situations, such 103as graphics, you want to suppress the error and output only the 104center of the interval. To do this, set the option 105variable @code{orthopoly_returns_intervals} to @code{false}. 106 107@c ===beg=== 108@c orthopoly_returns_intervals : false; 109@c jacobi_p (150, 2, 3, 0.2); 110@c ===end=== 111@example 112(%i1) orthopoly_returns_intervals : false; 113(%o1) false 114(%i2) jacobi_p (150, 2, 3, 0.2); 115(%o2) - 0.062017037936715 116@end example 117 118Refer to the section @pxref{Floating point Evaluation} for more information. 119 120Most functions in @code{orthopoly} have a @code{gradef} property; thus 121 122@c ===beg=== 123@c diff (hermite (n, x), x); 124@c diff (gen_laguerre (n, a, x), x); 125@c ===end=== 126@example 127(%i1) diff (hermite (n, x), x); 128(%o1) 2 n H (x) 129 n - 1 130(%i2) diff (gen_laguerre (n, a, x), x); 131 (a) (a) 132 n L (x) - (n + a) L (x) unit_step(n) 133 n n - 1 134(%o2) ------------------------------------------ 135 x 136@end example 137 138The unit step function in the second example prevents an error that would 139otherwise arise by evaluating with @var{n} equal to 0. 140 141@c CONTINUING PREVIOUS EXAMPLE HERE 142@c ===beg=== 143@c ev (%, n = 0); 144@c ===end=== 145@example 146(%i3) ev (%, n = 0); 147(%o3) 0 148@end example 149 150The @code{gradef} property only applies to the ``main'' variable; derivatives with 151respect other arguments usually result in an error message; for example 152 153@c ===beg=== 154@c diff (hermite (n, x), x); 155@c diff (hermite (n, x), n); 156@c ===end=== 157@example 158(%i1) diff (hermite (n, x), x); 159(%o1) 2 n H (x) 160 n - 1 161(%i2) diff (hermite (n, x), n); 162 163Maxima doesn't know the derivative of hermite with respect the first 164argument 165 -- an error. Quitting. To debug this try debugmode(true); 166@end example 167 168Generally, functions in @code{orthopoly} map over lists and matrices. For 169the mapping to fully evaluate, the option variables 170@code{doallmxops} and @code{listarith} must both be @code{true} (the defaults). 171To illustrate the mapping over matrices, consider 172 173@c ===beg=== 174@c hermite (2, x); 175@c m : matrix ([0, x], [y, 0]); 176@c hermite (2, m); 177@c ===end=== 178@example 179(%i1) hermite (2, x); 180 2 181(%o1) - 2 (1 - 2 x ) 182(%i2) m : matrix ([0, x], [y, 0]); 183 [ 0 x ] 184(%o2) [ ] 185 [ y 0 ] 186(%i3) hermite (2, m); 187 [ 2 ] 188 [ - 2 - 2 (1 - 2 x ) ] 189(%o3) [ ] 190 [ 2 ] 191 [ - 2 (1 - 2 y ) - 2 ] 192@end example 193 194In the second example, the @code{i, j} element of the value 195is @code{hermite (2, m[i,j])}; this is not the same as computing 196@code{-2 + 4 m . m}, as seen in the next example. 197 198@c CONTINUING PREVIOUS EXAMPLE HERE 199@c ===beg=== 200@c -2 * matrix ([1, 0], [0, 1]) + 4 * m . m; 201@c ===end=== 202@example 203(%i4) -2 * matrix ([1, 0], [0, 1]) + 4 * m . m; 204 [ 4 x y - 2 0 ] 205(%o4) [ ] 206 [ 0 4 x y - 2 ] 207@end example 208 209If you evaluate a function at a point outside its domain, generally 210@code{orthopoly} returns the function unevaluated. For example, 211 212@c ===beg=== 213@c legendre_p (2/3, x); 214@c ===end=== 215@example 216(%i1) legendre_p (2/3, x); 217(%o1) P (x) 218 2/3 219@end example 220 221@code{orthopoly} supports translation into TeX; it also does two-dimensional 222output on a terminal. 223 224@c ===beg=== 225@c spherical_harmonic (l, m, theta, phi); 226@c tex (%); 227@c jacobi_p (n, a, a - b, x/2); 228@c tex (%); 229@c ===end=== 230@example 231(%i1) spherical_harmonic (l, m, theta, phi); 232 m 233(%o1) Y (theta, phi) 234 l 235(%i2) tex (%); 236$$Y_@{l@}^@{m@}\left(\vartheta,\varphi\right)$$ 237(%o2) false 238(%i3) jacobi_p (n, a, a - b, x/2); 239 (a, a - b) x 240(%o3) P (-) 241 n 2 242(%i4) tex (%); 243$$P_@{n@}^@{\left(a,a-b\right)@}\left(@{@{x@}\over@{2@}@}\right)$$ 244(%o4) false 245@end example 246 247@subsection Limitations 248 249When an expression involves several orthogonal polynomials with 250symbolic orders, it's possible that the expression actually 251vanishes, yet Maxima is unable to simplify it to zero. If you 252divide by such a quantity, you'll be in trouble. For example, 253the following expression vanishes for integers @var{n} greater than 1, yet Maxima 254is unable to simplify it to zero. 255 256@c ===beg=== 257@c (2*n - 1) * legendre_p (n - 1, x) * x - n * legendre_p (n, x) 258@c + (1 - n) * legendre_p (n - 2, x); 259@c ===end=== 260@example 261(%i1) (2*n - 1) * legendre_p (n - 1, x) * x - n * legendre_p (n, x) 262 + (1 - n) * legendre_p (n - 2, x); 263(%o1) (2 n - 1) P (x) x - n P (x) + (1 - n) P (x) 264 n - 1 n n - 2 265@end example 266 267For a specific @var{n}, we can reduce the expression to zero. 268 269@c CONTINUING PREVIOUS EXAMPLE HERE 270@c ===beg=== 271@c ev (% ,n = 10, ratsimp); 272@c ===end=== 273@example 274(%i2) ev (% ,n = 10, ratsimp); 275(%o2) 0 276@end example 277 278Generally, the polynomial form of an orthogonal polynomial is ill-suited 279for floating point evaluation. Here's an example. 280 281@c ACTUALLY NEEDS load(orthopoly); BEFORE ANYTHING ELSE 282@c ===beg=== 283@c p : jacobi_p (100, 2, 3, x)$ 284@c subst (0.2, x, p); 285@c jacobi_p (100, 2, 3, 0.2); 286@c float(jacobi_p (100, 2, 3, 2/10)); 287@c ===end=== 288@example 289(%i1) p : jacobi_p (100, 2, 3, x)$ 290 291(%i2) subst (0.2, x, p); 292(%o2) 3.4442767023833592E+35 293(%i3) jacobi_p (100, 2, 3, 0.2); 294(%o3) interval(0.18413609135169, 6.8990300925815987E-12) 295(%i4) float(jacobi_p (100, 2, 3, 2/10)); 296(%o4) 0.18413609135169 297@end example 298 299The true value is about 0.184; this calculation suffers from extreme 300subtractive cancellation error. Expanding the polynomial and then 301evaluating, gives a better result. 302@c CONTINUING PREVIOUS EXAMPLE HERE 303@c ===beg=== 304@c p : expand (p)$ 305@c subst (0.2, x, p); 306@c ===end=== 307@example 308(%i5) p : expand(p)$ 309(%i6) subst (0.2, x, p); 310(%o6) 0.18413609766122982 311@end example 312 313This isn't a general rule; expanding the polynomial does not always 314result in an expression that is better suited for numerical evaluation. 315By far, the best way to do numerical evaluation is to make one or more 316of the function arguments floating point numbers. By doing that, 317specialized floating point algorithms are used for evaluation. 318 319Maxima's @code{float} function is somewhat indiscriminate; if you apply 320@code{float} to an expression involving an orthogonal polynomial with a 321symbolic degree or order parameter, these parameters may be 322converted into floats; after that, the expression will not evaluate 323fully. Consider 324 325@c ===beg=== 326@c assoc_legendre_p (n, 1, x); 327@c float (%); 328@c ev (%, n=2, x=0.9); 329@c ===end=== 330@example 331(%i1) assoc_legendre_p (n, 1, x); 332 1 333(%o1) P (x) 334 n 335(%i2) float (%); 336 1.0 337(%o2) P (x) 338 n 339(%i3) ev (%, n=2, x=0.9); 340 1.0 341(%o3) P (0.9) 342 2 343@end example 344 345The expression in (%o3) will not evaluate to a float; @code{orthopoly} doesn't 346recognize floating point values where it requires an integer. Similarly, 347numerical evaluation of the @code{pochhammer} function for orders that 348exceed @code{pochhammer_max_index} can be troublesome; consider 349 350@c ===beg=== 351@c x : pochhammer (1, 10), pochhammer_max_index : 5; 352@c ===end=== 353@example 354(%i1) x : pochhammer (1, 10), pochhammer_max_index : 5; 355(%o1) (1) 356 10 357@end example 358 359Applying @code{float} doesn't evaluate @var{x} to a float 360 361@c CONTINUING PREVIOUS EXAMPLE HERE 362@c ===beg=== 363@c float (x); 364@c ===end=== 365@example 366(%i2) float (x); 367(%o2) (1.0) 368 10.0 369@end example 370 371To evaluate @var{x} to a float, you'll need to bind 372@code{pochhammer_max_index} to 11 or greater and apply @code{float} to @var{x}. 373 374@c CONTINUING PREVIOUS EXAMPLE HERE 375@c ===beg=== 376@c float (x), pochhammer_max_index : 11; 377@c ===end=== 378@example 379(%i3) float (x), pochhammer_max_index : 11; 380(%o3) 3628800.0 381@end example 382 383The default value of @code{pochhammer_max_index} is 100; 384change its value after loading @code{orthopoly}. 385 386Finally, be aware that reference books vary on the definitions of the 387orthogonal polynomials; we've generally used the conventions 388of Abramowitz and Stegun. 389 390Before you suspect a bug in orthopoly, check some special cases 391to determine if your definitions match those used by @code{orthopoly}. 392Definitions often differ by a normalization; occasionally, authors 393use ``shifted'' versions of the functions that makes the family 394orthogonal on an interval other than @math{(-1, 1)}. To define, for example, 395a Legendre polynomial that is orthogonal on @math{(0, 1)}, define 396 397@c ===beg=== 398@c shifted_legendre_p (n, x) := legendre_p (n, 2*x - 1)$ 399@c shifted_legendre_p (2, rat (x)); 400@c legendre_p (2, rat (x)); 401@c ===end=== 402@example 403(%i1) shifted_legendre_p (n, x) := legendre_p (n, 2*x - 1)$ 404 405(%i2) shifted_legendre_p (2, rat (x)); 406 2 407(%o2)/R/ 6 x - 6 x + 1 408(%i3) legendre_p (2, rat (x)); 409 2 410 3 x - 1 411(%o3)/R/ -------- 412 2 413@end example 414 415@anchor{Floating point Evaluation} 416@subsection Floating point Evaluation 417 418Most functions in @code{orthopoly} use a running error analysis to 419estimate the error in floating point evaluation; the 420exceptions are the spherical Bessel functions and the associated Legendre 421polynomials of the second kind. For numerical evaluation, the spherical 422Bessel functions call SLATEC functions. No specialized method is used 423for numerical evaluation of the associated Legendre polynomials of the 424second kind. 425 426The running error analysis ignores errors that are second or higher order 427in the machine epsilon (also known as unit roundoff). It also 428ignores a few other errors. It's possible (although unlikely) 429that the actual error exceeds the estimate. 430 431Intervals have the form @code{interval (@var{c}, @var{r})}, where @var{c} is the 432center of the interval and @var{r} is its radius. The 433center of an interval can be a complex number, and the radius is always a positive real number. 434 435Here is an example. 436 437@c ===beg=== 438@c fpprec : 50$ 439@c y0 : jacobi_p (100, 2, 3, 0.2); 440@c y1 : bfloat (jacobi_p (100, 2, 3, 1/5)); 441@c ===end=== 442 443@example 444(%i1) fpprec : 50$ 445 446(%i2) y0 : jacobi_p (100, 2, 3, 0.2); 447(%o2) interval(0.1841360913516871, 6.8990300925815987E-12) 448(%i3) y1 : bfloat (jacobi_p (100, 2, 3, 1/5)); 449(%o3) 1.8413609135168563091370224958913493690868904463668b-1 450@end example 451 452Let's test that the actual error is smaller than the error estimate 453 454@c CONTINUING PREVIOUS EXAMPLE HERE 455@c ===beg=== 456@c is (abs (part (y0, 1) - y1) < part (y0, 2)); 457@c ===end=== 458@example 459(%i4) is (abs (part (y0, 1) - y1) < part (y0, 2)); 460(%o4) true 461@end example 462 463Indeed, for this example the error estimate is an upper bound for the 464true error. 465 466Maxima does not support arithmetic on intervals. 467 468@c ===beg=== 469@c legendre_p (7, 0.1) + legendre_p (8, 0.1); 470@c ===end=== 471@example 472(%i1) legendre_p (7, 0.1) + legendre_p (8, 0.1); 473(%o1) interval(0.18032072148437508, 3.1477135311021797E-15) 474 + interval(- 0.19949294375000004, 3.3769353084291579E-15) 475@end example 476 477A user could define arithmetic operators that do interval math. To 478define interval addition, we can define 479 480@c ===beg=== 481@c infix ("@+")$ 482@c "@+"(x,y) := interval (part (x, 1) + part (y, 1), part (x, 2) 483@c + part (y, 2))$ 484@c legendre_p (7, 0.1) @+ legendre_p (8, 0.1); 485@c ===end=== 486@example 487(%i1) infix ("@@+")$ 488 489(%i2) "@@+"(x,y) := interval (part (x, 1) + part (y, 1), part (x, 2) 490 + part (y, 2))$ 491 492(%i3) legendre_p (7, 0.1) @@+ legendre_p (8, 0.1); 493(%o3) interval(- 0.019172222265624955, 6.5246488395313372E-15) 494@end example 495 496The special floating point routines get called when the arguments 497are complex. For example, 498 499@c ===beg=== 500@c legendre_p (10, 2 + 3.0*%i); 501@c ===end=== 502@example 503(%i1) legendre_p (10, 2 + 3.0*%i); 504(%o1) interval(- 3.876378825E+7 %i - 6.0787748E+7, 505 1.2089173052721777E-6) 506@end example 507 508Let's compare this to the true value. 509 510@c ===beg=== 511@c float (expand (legendre_p (10, 2 + 3*%i))); 512@c ===end=== 513@example 514(%i1) float (expand (legendre_p (10, 2 + 3*%i))); 515(%o1) - 3.876378825E+7 %i - 6.0787748E+7 516@end example 517 518Additionally, when the arguments are big floats, the special floating point 519routines get called; however, the big floats are converted into double floats 520and the final result is a double. 521 522@c ===beg=== 523@c ultraspherical (150, 0.5b0, 0.9b0); 524@c ===end=== 525@example 526(%i1) ultraspherical (150, 0.5b0, 0.9b0); 527(%o1) interval(- 0.043009481257265, 3.3750051301228864E-14) 528@end example 529 530@subsection Graphics and @code{orthopoly} 531 532To plot expressions that involve the orthogonal polynomials, you 533must do two things: 534@enumerate 535@item 536Set the option variable @code{orthopoly_returns_intervals} to @code{false}, 537@item 538Quote any calls to @code{orthopoly} functions. 539@end enumerate 540If function calls aren't quoted, Maxima evaluates them to polynomials before 541plotting; consequently, the specialized floating point code doesn't get called. 542Here is an example of how to plot an expression that involves 543a Legendre polynomial. 544 545@c ===beg=== 546@c plot2d ('(legendre_p (5, x)), [x, 0, 1]), 547@c orthopoly_returns_intervals : false; 548@c ===end=== 549@example 550(%i1) plot2d ('(legendre_p (5, x)), [x, 0, 1]), 551 orthopoly_returns_intervals : false; 552(%o1) 553@end example 554 555@ifnotinfo 556@image{figures/orthopoly1,8cm} 557@end ifnotinfo 558 559The @i{entire} expression @code{legendre_p (5, x)} is quoted; this is 560different than just quoting the function name using @code{'legendre_p (5, @var{x})}. 561 562@opencatbox 563@category{Plotting} 564@closecatbox 565 566 567@subsection Miscellaneous Functions 568 569The @code{orthopoly} package defines the 570Pochhammer symbol and a unit step function. @code{orthopoly} uses 571the Kronecker delta function and the unit step function in 572@code{gradef} statements. 573 574To convert Pochhammer symbols into quotients of gamma functions, 575use @code{makegamma}. 576 577@c ===beg=== 578@c makegamma (pochhammer (x, n)); 579@c makegamma (pochhammer (1/2, 1/2)); 580@c ===end=== 581@example 582(%i1) makegamma (pochhammer (x, n)); 583 gamma(x + n) 584(%o1) ------------ 585 gamma(x) 586(%i2) makegamma (pochhammer (1/2, 1/2)); 587 1 588(%o2) --------- 589 sqrt(%pi) 590@end example 591 592Derivatives of the Pochhammer symbol are given in terms of the @code{psi} 593function. 594 595@c ===beg=== 596@c diff (pochhammer (x, n), x); 597@c diff (pochhammer (x, n), n); 598@c ===end=== 599@example 600(%i1) diff (pochhammer (x, n), x); 601(%o1) (x) (psi (x + n) - psi (x)) 602 n 0 0 603(%i2) diff (pochhammer (x, n), n); 604(%o2) (x) psi (x + n) 605 n 0 606@end example 607 608You need to be careful with the expression in (%o1); the difference of the 609@code{psi} functions has polynomials when @code{@var{x} = -1, -2, .., -@var{n}}. These polynomials 610cancel with factors in @code{pochhammer (@var{x}, @var{n})} making the derivative a degree 611@code{@var{n} - 1} polynomial when @var{n} is a positive integer. 612 613The Pochhammer symbol is defined for negative orders through its 614representation as a quotient of gamma functions. Consider 615 616@c ===beg=== 617@c q : makegamma (pochhammer (x, n)); 618@c sublis ([x=11/3, n= -6], q); 619@c ===end=== 620@example 621(%i1) q : makegamma (pochhammer (x, n)); 622 gamma(x + n) 623(%o1) ------------ 624 gamma(x) 625(%i2) sublis ([x=11/3, n= -6], q); 626 729 627(%o2) - ---- 628 2240 629@end example 630 631Alternatively, we can get this result directly. 632 633@c ===beg=== 634@c pochhammer (11/3, -6); 635@c ===end=== 636@example 637(%i1) pochhammer (11/3, -6); 638 729 639(%o1) - ---- 640 2240 641@end example 642 643The unit step function is left-continuous; thus 644 645@c ===beg=== 646@c [unit_step (-1/10), unit_step (0), unit_step (1/10)]; 647@c ===end=== 648@example 649(%i1) [unit_step (-1/10), unit_step (0), unit_step (1/10)]; 650(%o1) [0, 0, 1] 651@end example 652 653If you need a unit step function that is neither left or right continuous 654at zero, define your own using @code{signum}; for example, 655 656@c ===beg=== 657@c xunit_step (x) := (1 + signum (x))/2$ 658@c [xunit_step (-1/10), xunit_step (0), xunit_step (1/10)]; 659@c ===end=== 660@example 661(%i1) xunit_step (x) := (1 + signum (x))/2$ 662 663(%i2) [xunit_step (-1/10), xunit_step (0), xunit_step (1/10)]; 664 1 665(%o2) [0, -, 1] 666 2 667@end example 668 669Do not redefine @code{unit_step} itself; some code in @code{orthopoly} 670requires that the unit step function be left-continuous. 671 672@subsection Algorithms 673 674Generally, @code{orthopoly} does symbolic evaluation by using a hypergeometic 675representation of the orthogonal polynomials. The hypergeometic 676functions are evaluated using the (undocumented) functions @code{hypergeo11} 677and @code{hypergeo21}. The exceptions are the half-integer Bessel functions 678and the associated Legendre function of the second kind. The half-integer Bessel functions are 679evaluated using an explicit representation, and the associated Legendre 680function of the second kind is evaluated using recursion. 681 682For floating point evaluation, we again convert most functions into 683a hypergeometic form; we evaluate the hypergeometic functions using 684forward recursion. Again, the exceptions are the half-integer Bessel functions 685and the associated Legendre function of the second kind. Numerically, 686the half-integer Bessel functions are evaluated using the SLATEC code. 687 688 689@node Functions and Variables for orthogonal polynomials, , Introduction to orthogonal polynomials, orthopoly-pkg 690@section Functions and Variables for orthogonal polynomials 691 692@deffn {Function} assoc_legendre_p (@var{n}, @var{m}, @var{x}) 693The associated Legendre function of the first kind of degree @var{n} and 694order @var{m}. 695 696Reference: Abramowitz and Stegun, equations 22.5.37, page 779, 8.6.6 697(second equation), page 334, and 8.2.5, page 333. 698 699@opencatbox 700@category{Package orthopoly} 701@closecatbox 702 703@end deffn 704 705@deffn {Function} assoc_legendre_q (@var{n}, @var{m}, @var{x}) 706The associated Legendre function of the second kind of degree @var{n} 707and order @var{m}. 708 709Reference: Abramowitz and Stegun, equation 8.5.3 and 8.1.8. 710 711@opencatbox 712@category{Package orthopoly} 713@closecatbox 714 715@end deffn 716 717@deffn {Function} chebyshev_t (@var{n}, @var{x}) 718The Chebyshev polynomial of the first kind of degree @var{n}. 719 720Reference: Abramowitz and Stegun, equation 22.5.47, page 779. 721 722@opencatbox 723@category{Package orthopoly} 724@closecatbox 725 726@end deffn 727 728@deffn {Function} chebyshev_u (@var{n}, @var{x}) 729The Chebyshev polynomial of the second kind of degree @var{n}. 730 731Reference: Abramowitz and Stegun, equation 22.5.48, page 779. 732 733@opencatbox 734@category{Package orthopoly} 735@closecatbox 736 737@end deffn 738 739@deffn {Function} gen_laguerre (@var{n}, @var{a}, @var{x}) 740The generalized Laguerre polynomial of degree @var{n}. 741 742Reference: Abramowitz and Stegun, equation 22.5.54, page 780. 743 744@opencatbox 745@category{Package orthopoly} 746@closecatbox 747 748@end deffn 749 750@deffn {Function} hermite (@var{n}, @var{x}) 751The Hermite polynomial of degree @var{n}. 752 753Reference: Abramowitz and Stegun, equation 22.5.55, page 780. 754 755@opencatbox 756@category{Package orthopoly} 757@closecatbox 758 759@end deffn 760 761@deffn {Function} intervalp (@var{e}) 762Return @code{true} if the input is an interval and return false if it isn't. 763 764@opencatbox 765@category{Package orthopoly} 766@category{Predicate functions} 767@closecatbox 768 769@end deffn 770 771@deffn {Function} jacobi_p (@var{n}, @var{a}, @var{b}, @var{x}) 772The Jacobi polynomial. 773 774The Jacobi polynomials are actually defined for all 775@var{a} and @var{b}; however, the Jacobi polynomial 776weight @code{(1 - @var{x})^@var{a} (1 + @var{x})^@var{b}} isn't integrable for @code{@var{a} <= -1} or 777@code{@var{b} <= -1}. 778 779Reference: Abramowitz and Stegun, equation 22.5.42, page 779. 780 781@opencatbox 782@category{Package orthopoly} 783@closecatbox 784 785@end deffn 786 787@deffn {Function} laguerre (@var{n}, @var{x}) 788The Laguerre polynomial of degree @var{n}. 789 790Reference: Abramowitz and Stegun, equations 22.5.16 and 22.5.54, page 780. 791 792@opencatbox 793@category{Package orthopoly} 794@closecatbox 795 796@end deffn 797 798@deffn {Function} legendre_p (@var{n}, @var{x}) 799The Legendre polynomial of the first kind of degree @var{n}. 800 801Reference: Abramowitz and Stegun, equations 22.5.50 and 22.5.51, page 779. 802 803@opencatbox 804@category{Package orthopoly} 805@closecatbox 806 807@end deffn 808 809@deffn {Function} legendre_q (@var{n}, @var{x}) 810The Legendre function of the second kind of degree @var{n}. 811 812Reference: Abramowitz and Stegun, equations 8.5.3 and 8.1.8. 813 814@opencatbox 815@category{Package orthopoly} 816@closecatbox 817 818@end deffn 819 820@deffn {Function} orthopoly_recur (@var{f}, @var{args}) 821Returns a recursion relation for the orthogonal function family 822@var{f} with arguments @var{args}. The recursion is with 823respect to the polynomial degree. 824 825@c ===beg=== 826@c orthopoly_recur (legendre_p, [n, x]); 827@c ===end=== 828@example 829(%i1) orthopoly_recur (legendre_p, [n, x]); 830 (2 n + 1) P (x) x - n P (x) 831 n n - 1 832(%o1) P (x) = ------------------------------- 833 n + 1 n + 1 834@end example 835 836The second argument to @code{orthopoly_recur} must be a list with the 837correct number of arguments for the function @var{f}; if it isn't, 838Maxima signals an error. 839 840@c ===beg=== 841@c orthopoly_recur (jacobi_p, [n, x]); 842@c ===end=== 843@example 844(%i1) orthopoly_recur (jacobi_p, [n, x]); 845 846Function jacobi_p needs 4 arguments, instead it received 2 847 -- an error. Quitting. To debug this try debugmode(true); 848@end example 849 850Additionally, when @var{f} isn't the name of one of the 851families of orthogonal polynomials, an error is signalled. 852 853@c ===beg=== 854@c orthopoly_recur (foo, [n, x]); 855@c ===end=== 856@example 857(%i1) orthopoly_recur (foo, [n, x]); 858 859A recursion relation for foo isn't known to Maxima 860 -- an error. Quitting. To debug this try debugmode(true); 861@end example 862 863@opencatbox 864@category{Package orthopoly} 865@closecatbox 866 867@end deffn 868 869@defvr {Variable} orthopoly_returns_intervals 870Default value: @code{true} 871 872When @code{orthopoly_returns_intervals} is @code{true}, floating point results are returned in 873the form @code{interval (@var{c}, @var{r})}, where @var{c} is the center of an interval 874and @var{r} is its radius. The center can be a complex number; in that 875case, the interval is a disk in the complex plane. 876 877@opencatbox 878@category{Package orthopoly} 879@closecatbox 880 881@end defvr 882 883@deffn {Function} orthopoly_weight (@var{f}, @var{args}) 884 885Returns a three element list; the first element is 886the formula of the weight for the orthogonal polynomial family 887@var{f} with arguments given by the list @var{args}; the 888second and third elements give the lower and upper endpoints 889of the interval of orthogonality. For example, 890 891@c ===beg=== 892@c w : orthopoly_weight (hermite, [n, x]); 893@c integrate (w[1] * hermite (3, x) * hermite (2, x), x, w[2], w[3]); 894@c ===end=== 895@example 896(%i1) w : orthopoly_weight (hermite, [n, x]); 897 2 898 - x 899(%o1) [%e , - inf, inf] 900(%i2) integrate(w[1]*hermite(3, x)*hermite(2, x), x, w[2], w[3]); 901(%o2) 0 902@end example 903 904The main variable of @var{f} must be a symbol; if it isn't, Maxima 905signals an error. 906 907@opencatbox 908@category{Package orthopoly} 909@closecatbox 910 911@end deffn 912 913@deffn {Function} pochhammer (@var{x}, @var{n}) 914The Pochhammer symbol. For nonnegative integers @var{n} with 915@code{@var{n} <= pochhammer_max_index}, the expression @code{pochhammer (@var{x}, @var{n})} 916evaluates to the product @code{@var{x} (@var{x} + 1) (@var{x} + 2) ... (@var{x} + n - 1)} 917when @code{@var{n} > 0} and 918to 1 when @code{@var{n} = 0}. For negative @var{n}, 919@code{pochhammer (@var{x}, @var{n})} is defined as @code{(-1)^@var{n} / pochhammer (1 - @var{x}, -@var{n})}. 920Thus 921 922@c ===beg=== 923@c pochhammer (x, 3); 924@c pochhammer (x, -3); 925@c ===end=== 926@example 927(%i1) pochhammer (x, 3); 928(%o1) x (x + 1) (x + 2) 929(%i2) pochhammer (x, -3); 930 1 931(%o2) - ----------------------- 932 (1 - x) (2 - x) (3 - x) 933@end example 934 935To convert a Pochhammer symbol into a quotient of gamma functions, 936(see Abramowitz and Stegun, equation 6.1.22) use @code{makegamma}; for example 937 938@c ===beg=== 939@c makegamma (pochhammer (x, n)); 940@c ===end=== 941@example 942(%i1) makegamma (pochhammer (x, n)); 943 gamma(x + n) 944(%o1) ------------ 945 gamma(x) 946@end example 947 948When @var{n} exceeds @code{pochhammer_max_index} or when @var{n} 949is symbolic, @code{pochhammer} returns a noun form. 950 951@c ===beg=== 952@c pochhammer (x, n); 953@c ===end=== 954@example 955(%i1) pochhammer (x, n); 956(%o1) (x) 957 n 958@end example 959 960@opencatbox 961@category{Package orthopoly} 962@category{Gamma and factorial functions} 963@closecatbox 964 965@end deffn 966 967@defvr {Variable} pochhammer_max_index 968Default value: 100 969 970@code{pochhammer (@var{n}, @var{x})} expands to a product if and only if 971@code{@var{n} <= pochhammer_max_index}. 972 973Examples: 974 975@c ===beg=== 976@c pochhammer (x, 3), pochhammer_max_index : 3; 977@c pochhammer (x, 4), pochhammer_max_index : 3; 978@c ===end=== 979@example 980(%i1) pochhammer (x, 3), pochhammer_max_index : 3; 981(%o1) x (x + 1) (x + 2) 982(%i2) pochhammer (x, 4), pochhammer_max_index : 3; 983(%o2) (x) 984 4 985@end example 986 987Reference: Abramowitz and Stegun, equation 6.1.16, page 256. 988 989@opencatbox 990@category{Package orthopoly} 991@category{Gamma and factorial functions} 992@closecatbox 993 994@end defvr 995 996@deffn {Function} spherical_bessel_j (@var{n}, @var{x}) 997The spherical Bessel function of the first kind. 998 999Reference: Abramowitz and Stegun, equations 10.1.8, page 437 and 10.1.15, page 439. 1000 1001@opencatbox 1002@category{Package orthopoly} 1003@category{Bessel functions} 1004@closecatbox 1005 1006@end deffn 1007 1008@deffn {Function} spherical_bessel_y (@var{n}, @var{x}) 1009The spherical Bessel function of the second kind. 1010 1011Reference: Abramowitz and Stegun, equations 10.1.9, page 437 and 10.1.15, page 439. 1012 1013@opencatbox 1014@category{Package orthopoly} 1015@category{Bessel functions} 1016@closecatbox 1017 1018@end deffn 1019 1020@deffn {Function} spherical_hankel1 (@var{n}, @var{x}) 1021The spherical Hankel function of the 1022first kind. 1023 1024Reference: Abramowitz and Stegun, equation 10.1.36, page 439. 1025 1026@opencatbox 1027@category{Package orthopoly} 1028@category{Bessel functions} 1029@closecatbox 1030 1031@end deffn 1032 1033@deffn {Function} spherical_hankel2 (@var{n}, @var{x}) 1034The spherical Hankel function of the second kind. 1035 1036Reference: Abramowitz and Stegun, equation 10.1.17, page 439. 1037 1038@opencatbox 1039@category{Package orthopoly} 1040@category{Bessel functions} 1041@closecatbox 1042 1043@end deffn 1044 1045@deffn {Function} spherical_harmonic (@var{n}, @var{m}, @var{x}, @var{y}) 1046The spherical harmonic function. 1047 1048Reference: Merzbacher 9.64. 1049 1050@opencatbox 1051@category{Package orthopoly} 1052@closecatbox 1053 1054@end deffn 1055 1056@anchor{unit_step} 1057@deffn {Function} unit_step (@var{x}) 1058The left-continuous unit step function; thus 1059@code{unit_step (@var{x})} vanishes for @code{@var{x} <= 0} and equals 10601 for @code{@var{x} > 0}. 1061 1062If you want a unit step function that 1063takes on the value 1/2 at zero, use @code{(1 + signum (@var{x}))/2}. 1064 1065@opencatbox 1066@category{Package orthopoly} 1067@category{Mathematical functions} 1068@closecatbox 1069 1070@end deffn 1071 1072@deffn {Function} ultraspherical (@var{n}, @var{a}, @var{x}) 1073The ultraspherical polynomial (also known as the Gegenbauer polynomial). 1074 1075Reference: Abramowitz and Stegun, equation 22.5.46, page 779. 1076 1077@opencatbox 1078@category{Package orthopoly} 1079@closecatbox 1080 1081@end deffn 1082