| blob | 6921f2e18f14956e52a34d3695a6eedbcd193f49 |
1 @c -*- mode: texinfo -*-
2 @menu
3 * Introduction to orthogonal polynomials::
4 * Functions and Variables for orthogonal polynomials::
5 @end menu
7 @node Introduction to orthogonal polynomials, Functions and Variables for orthogonal polynomials, Package orthopoly, Package orthopoly
8 @section Introduction to orthogonal polynomials
10 @code{orthopoly} is a package for symbolic and numerical evaluation of
11 several kinds of orthogonal polynomials, including Chebyshev,
12 Laguerre, Hermite, Jacobi, Legendre, and ultraspherical (Gegenbauer)
13 polynomials. Additionally, @code{orthopoly} includes support for the spherical Bessel,
14 spherical Hankel, and spherical harmonic functions.
16 For the most part, @code{orthopoly} follows the conventions of Abramowitz and Stegun
17 @i{Handbook of Mathematical Functions}, Chapter 22 (10th printing, December 1972);
18 additionally, we use Gradshteyn and Ryzhik,
19 @i{Table of Integrals, Series, and Products} (1980 corrected and
20 enlarged edition), and Eugen Merzbacher @i{Quantum Mechanics} (2nd edition, 1970).
22 @c INSTALLATION INSTRUCTIONS NO LONGER RELEVANT
23 @c BUT MAYBE SOME OF THESE FILES SHOULD BE MENTIONED IN ANOTHER CONTEXT
24 @c This will create a directory @code{orthopoly_x} (again x is the release
25 @c identifier) that contains the source file @code{orthopoly.lisp}, user
26 @c documentation in html and texi formats, a sample maxima initialization file
27 @c @code{orthopoly-init.lisp}, a README file, a testing routine
28 @c @code{test_orthopoly.mac}, and two demonstration files.
30 @c Start Maxima and compile orthopoly. To do this, use the command
31 @c
32 @c (c1) compile_file("orthopoly.lisp");
34 Barton Willis of the University of Nebraska at Kearney (UNK) wrote
35 the @code{orthopoly} package and its documentation. The package
36 is released under the GNU General Public License (GPL).
38 @opencatbox{Categories:}
39 @category{Orthogonal polynomials}
40 @category{Share packages}
41 @category{Package orthopoly}
42 @closecatbox
44 @subsection Getting Started with orthopoly
46 @code{load ("orthopoly")} loads the @code{orthopoly} package.
48 To find the third-order Legendre polynomial,
50 @c ===beg===
51 @c legendre_p (3, x);
52 @c ===end===
53 @example
54 (%i1) legendre_p (3, x);
55 3 2
56 5 (1 - x) 15 (1 - x)
57 (%o1) - ---------- + ----------- - 6 (1 - x) + 1
58 2 2
59 @end example
61 To express this as a sum of powers of @var{x}, apply @var{ratsimp} or @var{rat}
62 to the result.
64 @c CONTINUING PREVIOUS EXAMPLE HERE
65 @c ===beg===
66 @c [ratsimp (%), rat (%)];
67 @c ===end===
68 @example
69 (%i2) [ratsimp (%), rat (%)];
70 3 3
71 5 x - 3 x 5 x - 3 x
72 (%o2)/R/ [----------, ----------]
73 2 2
74 @end example
76 Alternatively, make the second argument to @code{legendre_p} (its ``main'' variable)
77 a canonical rational expression (CRE).
79 @c ===beg===
80 @c legendre_p (3, rat (x));
81 @c ===end===
82 @example
83 (%i1) legendre_p (3, rat (x));
84 3
85 5 x - 3 x
86 (%o1)/R/ ----------
87 2
88 @end example
90 For floating point evaluation, @code{orthopoly} uses a running error analysis
91 to estimate an upper bound for the error. For example,
93 @c ===beg===
94 @c jacobi_p (150, 2, 3, 0.2);
95 @c ===end===
96 @example
97 (%i1) jacobi_p (150, 2, 3, 0.2);
98 (%o1) interval(- 0.062017037936715, 1.533267919277521E-11)
99 @end example
101 Intervals have the form @code{interval (@var{c}, @var{r})}, where @var{c} is the
102 center and @var{r} is the radius of the interval. Since Maxima
103 does not support arithmetic on intervals, in some situations, such
104 as graphics, you want to suppress the error and output only the
105 center of the interval. To do this, set the option
106 variable @code{orthopoly_returns_intervals} to @code{false}.
108 @c ===beg===
109 @c orthopoly_returns_intervals : false;
110 @c jacobi_p (150, 2, 3, 0.2);
111 @c ===end===
112 @example
113 (%i1) orthopoly_returns_intervals : false;
114 (%o1) false
115 (%i2) jacobi_p (150, 2, 3, 0.2);
116 (%o2) - 0.062017037936715
117 @end example
119 Refer to the section @pxref{Floating point Evaluation} for more information.
121 Most functions in @code{orthopoly} have a @code{gradef} property; thus
123 @c ===beg===
124 @c diff (hermite (n, x), x);
125 @c diff (gen_laguerre (n, a, x), x);
126 @c ===end===
127 @example
128 (%i1) diff (hermite (n, x), x);
129 (%o1) 2 n H (x)
130 n - 1
131 (%i2) diff (gen_laguerre (n, a, x), x);
132 (a) (a)
133 n L (x) - (n + a) L (x) unit_step(n)
134 n n - 1
135 (%o2) ------------------------------------------
136 x
137 @end example
139 The unit step function in the second example prevents an error that would
140 otherwise arise by evaluating with @var{n} equal to 0.
142 @c CONTINUING PREVIOUS EXAMPLE HERE
143 @c ===beg===
144 @c ev (%, n = 0);
145 @c ===end===
146 @example
147 (%i3) ev (%, n = 0);
148 (%o3) 0
149 @end example
151 The @code{gradef} property only applies to the ``main'' variable; derivatives with
152 respect other arguments usually result in an error message; for example
154 @c ===beg===
155 @c diff (hermite (n, x), x);
156 @c diff (hermite (n, x), n);
157 @c ===end===
158 @example
159 (%i1) diff (hermite (n, x), x);
160 (%o1) 2 n H (x)
161 n - 1
162 (%i2) diff (hermite (n, x), n);
164 Maxima doesn't know the derivative of hermite with respect the first
165 argument
166 -- an error. Quitting. To debug this try debugmode(true);
167 @end example
169 Generally, functions in @code{orthopoly} map over lists and matrices. For
170 the mapping to fully evaluate, the option variables
171 @code{doallmxops} and @code{listarith} must both be @code{true} (the defaults).
172 To illustrate the mapping over matrices, consider
174 @c ===beg===
175 @c hermite (2, x);
176 @c m : matrix ([0, x], [y, 0]);
177 @c hermite (2, m);
178 @c ===end===
179 @example
180 (%i1) hermite (2, x);
181 2
182 (%o1) - 2 (1 - 2 x )
183 (%i2) m : matrix ([0, x], [y, 0]);
184 [ 0 x ]
185 (%o2) [ ]
186 [ y 0 ]
187 (%i3) hermite (2, m);
188 [ 2 ]
189 [ - 2 - 2 (1 - 2 x ) ]
190 (%o3) [ ]
191 [ 2 ]
192 [ - 2 (1 - 2 y ) - 2 ]
193 @end example
195 In the second example, the @code{i, j} element of the value
196 is @code{hermite (2, m[i,j])}; this is not the same as computing
197 @code{-2 + 4 m . m}, as seen in the next example.
199 @c CONTINUING PREVIOUS EXAMPLE HERE
200 @c ===beg===
201 @c -2 * matrix ([1, 0], [0, 1]) + 4 * m . m;
202 @c ===end===
203 @example
204 (%i4) -2 * matrix ([1, 0], [0, 1]) + 4 * m . m;
205 [ 4 x y - 2 0 ]
206 (%o4) [ ]
207 [ 0 4 x y - 2 ]
208 @end example
210 If you evaluate a function at a point outside its domain, generally
211 @code{orthopoly} returns the function unevaluated. For example,
213 @c ===beg===
214 @c legendre_p (2/3, x);
215 @c ===end===
216 @example
217 (%i1) legendre_p (2/3, x);
218 (%o1) P (x)
219 2/3
220 @end example
222 @code{orthopoly} supports translation into TeX; it also does two-dimensional
223 output on a terminal.
225 @c ===beg===
226 @c spherical_harmonic (l, m, theta, phi);
227 @c tex (%);
228 @c jacobi_p (n, a, a - b, x/2);
229 @c tex (%);
230 @c ===end===
231 @example
232 (%i1) spherical_harmonic (l, m, theta, phi);
233 m
234 (%o1) Y (theta, phi)
235 l
236 (%i2) tex (%);
237 $$Y_@{l@}^@{m@}\left(\vartheta,\varphi\right)$$
238 (%o2) false
239 (%i3) jacobi_p (n, a, a - b, x/2);
240 (a, a - b) x
241 (%o3) P (-)
242 n 2
243 (%i4) tex (%);
244 $$P_@{n@}^@{\left(a,a-b\right)@}\left(@{@{x@}\over@{2@}@}\right)$$
245 (%o4) false
246 @end example
248 @subsection Limitations
250 When an expression involves several orthogonal polynomials with
251 symbolic orders, it's possible that the expression actually
252 vanishes, yet Maxima is unable to simplify it to zero. If you
253 divide by such a quantity, you'll be in trouble. For example,
254 the following expression vanishes for integers @var{n} greater than 1, yet Maxima
255 is unable to simplify it to zero.
257 @c ===beg===
258 @c (2*n - 1) * legendre_p (n - 1, x) * x - n * legendre_p (n, x)
259 @c + (1 - n) * legendre_p (n - 2, x);
260 @c ===end===
261 @example
262 (%i1) (2*n - 1) * legendre_p (n - 1, x) * x - n * legendre_p (n, x)
263 + (1 - n) * legendre_p (n - 2, x);
264 (%o1) (2 n - 1) P (x) x - n P (x) + (1 - n) P (x)
265 n - 1 n n - 2
266 @end example
268 For a specific @var{n}, we can reduce the expression to zero.
270 @c CONTINUING PREVIOUS EXAMPLE HERE
271 @c ===beg===
272 @c ev (% ,n = 10, ratsimp);
273 @c ===end===
274 @example
275 (%i2) ev (% ,n = 10, ratsimp);
276 (%o2) 0
277 @end example
279 Generally, the polynomial form of an orthogonal polynomial is ill-suited
280 for floating point evaluation. Here's an example.
282 @c ACTUALLY NEEDS load("orthopoly"); BEFORE ANYTHING ELSE
283 @c ===beg===
284 @c p : jacobi_p (100, 2, 3, x)$
285 @c subst (0.2, x, p);
286 @c jacobi_p (100, 2, 3, 0.2);
287 @c float(jacobi_p (100, 2, 3, 2/10));
288 @c ===end===
289 @example
290 (%i1) p : jacobi_p (100, 2, 3, x)$
292 (%i2) subst (0.2, x, p);
293 (%o2) 3.4442767023833592E+35
294 (%i3) jacobi_p (100, 2, 3, 0.2);
295 (%o3) interval(0.18413609135169, 6.8990300925815987E-12)
296 (%i4) float(jacobi_p (100, 2, 3, 2/10));
297 (%o4) 0.18413609135169
298 @end example
300 The true value is about 0.184; this calculation suffers from extreme
301 subtractive cancellation error. Expanding the polynomial and then
302 evaluating, gives a better result.
303 @c CONTINUING PREVIOUS EXAMPLE HERE
304 @c ===beg===
305 @c p : expand (p)$
306 @c subst (0.2, x, p);
307 @c ===end===
308 @example
309 (%i5) p : expand(p)$
310 (%i6) subst (0.2, x, p);
311 (%o6) 0.18413609766122982
312 @end example
314 This isn't a general rule; expanding the polynomial does not always
315 result in an expression that is better suited for numerical evaluation.
316 By far, the best way to do numerical evaluation is to make one or more
317 of the function arguments floating point numbers. By doing that,
318 specialized floating point algorithms are used for evaluation.
320 Maxima's @code{float} function is somewhat indiscriminate; if you apply
321 @code{float} to an expression involving an orthogonal polynomial with a
322 symbolic degree or order parameter, these parameters may be
323 converted into floats; after that, the expression will not evaluate
324 fully. Consider
326 @c ===beg===
327 @c assoc_legendre_p (n, 1, x);
328 @c float (%);
329 @c ev (%, n=2, x=0.9);
330 @c ===end===
331 @example
332 (%i1) assoc_legendre_p (n, 1, x);
333 1
334 (%o1) P (x)
335 n
336 (%i2) float (%);
337 1.0
338 (%o2) P (x)
339 n
340 (%i3) ev (%, n=2, x=0.9);
341 1.0
342 (%o3) P (0.9)
343 2
344 @end example
346 The expression in (%o3) will not evaluate to a float; @code{orthopoly} doesn't
347 recognize floating point values where it requires an integer. Similarly,
348 numerical evaluation of the @code{pochhammer} function for orders that
349 exceed @code{pochhammer_max_index} can be troublesome; consider
351 @c ===beg===
352 @c x : pochhammer (1, 10), pochhammer_max_index : 5;
353 @c ===end===
354 @example
355 (%i1) x : pochhammer (1, 10), pochhammer_max_index : 5;
356 (%o1) (1)
357 10
358 @end example
360 Applying @code{float} doesn't evaluate @var{x} to a float
362 @c CONTINUING PREVIOUS EXAMPLE HERE
363 @c ===beg===
364 @c float (x);
365 @c ===end===
366 @example
367 (%i2) float (x);
368 (%o2) (1.0)
369 10.0
370 @end example
372 To evaluate @var{x} to a float, you'll need to bind
373 @code{pochhammer_max_index} to 11 or greater and apply @code{float} to @var{x}.
375 @c CONTINUING PREVIOUS EXAMPLE HERE
376 @c ===beg===
377 @c float (x), pochhammer_max_index : 11;
378 @c ===end===
379 @example
380 (%i3) float (x), pochhammer_max_index : 11;
381 (%o3) 3628800.0
382 @end example
384 The default value of @code{pochhammer_max_index} is 100;
385 change its value after loading @code{orthopoly}.
387 Finally, be aware that reference books vary on the definitions of the
388 orthogonal polynomials; we've generally used the conventions
389 of Abramowitz and Stegun.
391 Before you suspect a bug in orthopoly, check some special cases
392 to determine if your definitions match those used by @code{orthopoly}.
393 Definitions often differ by a normalization; occasionally, authors
394 use ``shifted'' versions of the functions that makes the family
395 orthogonal on an interval other than @math{(-1, 1)}. To define, for example,
396 a Legendre polynomial that is orthogonal on @math{(0, 1)}, define
398 @c ===beg===
399 @c shifted_legendre_p (n, x) := legendre_p (n, 2*x - 1)$
400 @c shifted_legendre_p (2, rat (x));
401 @c legendre_p (2, rat (x));
402 @c ===end===
403 @example
404 (%i1) shifted_legendre_p (n, x) := legendre_p (n, 2*x - 1)$
406 (%i2) shifted_legendre_p (2, rat (x));
407 2
408 (%o2)/R/ 6 x - 6 x + 1
409 (%i3) legendre_p (2, rat (x));
410 2
411 3 x - 1
412 (%o3)/R/ --------
413 2
414 @end example
416 @anchor{Floating point Evaluation}
417 @subsection Floating point Evaluation
419 Most functions in @code{orthopoly} use a running error analysis to
420 estimate the error in floating point evaluation; the
421 exceptions are the spherical Bessel functions and the associated Legendre
422 polynomials of the second kind. For numerical evaluation, the spherical
423 Bessel functions call SLATEC functions. No specialized method is used
424 for numerical evaluation of the associated Legendre polynomials of the
425 second kind.
427 The running error analysis ignores errors that are second or higher order
428 in the machine epsilon (also known as unit roundoff). It also
429 ignores a few other errors. It's possible (although unlikely)
430 that the actual error exceeds the estimate.
432 Intervals have the form @code{interval (@var{c}, @var{r})}, where @var{c} is the
433 center of the interval and @var{r} is its radius. The
434 center of an interval can be a complex number, and the radius is always a positive real number.
436 Here is an example.
438 @c ===beg===
439 @c fpprec : 50$
440 @c y0 : jacobi_p (100, 2, 3, 0.2);
441 @c y1 : bfloat (jacobi_p (100, 2, 3, 1/5));
442 @c ===end===
444 @example
445 (%i1) fpprec : 50$
447 (%i2) y0 : jacobi_p (100, 2, 3, 0.2);
448 (%o2) interval(0.1841360913516871, 6.8990300925815987E-12)
449 (%i3) y1 : bfloat (jacobi_p (100, 2, 3, 1/5));
450 (%o3) 1.8413609135168563091370224958913493690868904463668b-1
451 @end example
453 Let's test that the actual error is smaller than the error estimate
455 @c CONTINUING PREVIOUS EXAMPLE HERE
456 @c ===beg===
457 @c is (abs (part (y0, 1) - y1) < part (y0, 2));
458 @c ===end===
459 @example
460 (%i4) is (abs (part (y0, 1) - y1) < part (y0, 2));
461 (%o4) true
462 @end example
464 Indeed, for this example the error estimate is an upper bound for the
465 true error.
467 Maxima does not support arithmetic on intervals.
469 @c ===beg===
470 @c legendre_p (7, 0.1) + legendre_p (8, 0.1);
471 @c ===end===
472 @example
473 (%i1) legendre_p (7, 0.1) + legendre_p (8, 0.1);
474 (%o1) interval(0.18032072148437508, 3.1477135311021797E-15)
475 + interval(- 0.19949294375000004, 3.3769353084291579E-15)
476 @end example
478 A user could define arithmetic operators that do interval math. To
479 define interval addition, we can define
481 @c ===beg===
482 @c infix ("@+")$
483 @c "@+"(x,y) := interval (part (x, 1) + part (y, 1), part (x, 2)
484 @c + part (y, 2))$
485 @c legendre_p (7, 0.1) @+ legendre_p (8, 0.1);
486 @c ===end===
487 @example
488 (%i1) infix ("@@+")$
490 (%i2) "@@+"(x,y) := interval (part (x, 1) + part (y, 1), part (x, 2)
491 + part (y, 2))$
493 (%i3) legendre_p (7, 0.1) @@+ legendre_p (8, 0.1);
494 (%o3) interval(- 0.019172222265624955, 6.5246488395313372E-15)
495 @end example
497 The special floating point routines get called when the arguments
498 are complex. For example,
500 @c ===beg===
501 @c legendre_p (10, 2 + 3.0*%i);
502 @c ===end===
503 @example
504 (%i1) legendre_p (10, 2 + 3.0*%i);
505 (%o1) interval(- 3.876378825E+7 %i - 6.0787748E+7,
506 1.2089173052721777E-6)
507 @end example
509 Let's compare this to the true value.
511 @c ===beg===
512 @c float (expand (legendre_p (10, 2 + 3*%i)));
513 @c ===end===
514 @example
515 (%i1) float (expand (legendre_p (10, 2 + 3*%i)));
516 (%o1) - 3.876378825E+7 %i - 6.0787748E+7
517 @end example
519 Additionally, when the arguments are big floats, the special floating point
520 routines get called; however, the big floats are converted into double floats
521 and the final result is a double.
523 @c ===beg===
524 @c ultraspherical (150, 0.5b0, 0.9b0);
525 @c ===end===
526 @example
527 (%i1) ultraspherical (150, 0.5b0, 0.9b0);
528 (%o1) interval(- 0.043009481257265, 3.3750051301228864E-14)
529 @end example
531 @subsection Graphics and @code{orthopoly}
533 To plot expressions that involve the orthogonal polynomials, you
534 must do two things:
535 @enumerate
536 @item
537 Set the option variable @code{orthopoly_returns_intervals} to @code{false},
538 @item
539 Quote any calls to @code{orthopoly} functions.
540 @end enumerate
541 If function calls aren't quoted, Maxima evaluates them to polynomials before
542 plotting; consequently, the specialized floating point code doesn't get called.
543 Here is an example of how to plot an expression that involves
544 a Legendre polynomial.
546 @c ===beg===
547 @c plot2d ('(legendre_p (5, x)), [x, 0, 1]),
548 @c orthopoly_returns_intervals : false;
549 @c ===end===
550 @example
551 (%i1) plot2d ('(legendre_p (5, x)), [x, 0, 1]),
552 orthopoly_returns_intervals : false;
553 (%o1)
554 @end example
556 @ifnotinfo
557 @image{figures/orthopoly1,8cm}
558 @end ifnotinfo
560 The @i{entire} expression @code{legendre_p (5, x)} is quoted; this is
561 different than just quoting the function name using @code{'legendre_p (5, @var{x})}.
563 @opencatbox{Categories:}
564 @category{Plotting}
565 @closecatbox
568 @subsection Miscellaneous Functions
570 The @code{orthopoly} package defines the
571 Pochhammer symbol and a unit step function. @code{orthopoly} uses
572 the Kronecker delta function and the unit step function in
573 @code{gradef} statements.
575 To convert Pochhammer symbols into quotients of gamma functions,
576 use @code{makegamma}.
578 @c ===beg===
579 @c makegamma (pochhammer (x, n));
580 @c makegamma (pochhammer (1/2, 1/2));
581 @c ===end===
582 @example
583 (%i1) makegamma (pochhammer (x, n));
584 gamma(x + n)
585 (%o1) ------------
586 gamma(x)
587 (%i2) makegamma (pochhammer (1/2, 1/2));
588 1
589 (%o2) ---------
590 sqrt(%pi)
591 @end example
593 Derivatives of the Pochhammer symbol are given in terms of the @code{psi}
594 function.
596 @c ===beg===
597 @c diff (pochhammer (x, n), x);
598 @c diff (pochhammer (x, n), n);
599 @c ===end===
600 @example
601 (%i1) diff (pochhammer (x, n), x);
602 (%o1) (x) (psi (x + n) - psi (x))
603 n 0 0
604 (%i2) diff (pochhammer (x, n), n);
605 (%o2) (x) psi (x + n)
606 n 0
607 @end example
609 You need to be careful with the expression in (%o1); the difference of the
610 @code{psi} functions has polynomials when @code{@var{x} = -1, -2, .., -@var{n}}. These polynomials
611 cancel with factors in @code{pochhammer (@var{x}, @var{n})} making the derivative a degree
612 @code{@var{n} - 1} polynomial when @var{n} is a positive integer.
614 The Pochhammer symbol is defined for negative orders through its
615 representation as a quotient of gamma functions. Consider
617 @c ===beg===
618 @c q : makegamma (pochhammer (x, n));
619 @c sublis ([x=11/3, n= -6], q);
620 @c ===end===
621 @example
622 (%i1) q : makegamma (pochhammer (x, n));
623 gamma(x + n)
624 (%o1) ------------
625 gamma(x)
626 (%i2) sublis ([x=11/3, n= -6], q);
627 729
628 (%o2) - ----
629 2240
630 @end example
632 Alternatively, we can get this result directly.
634 @c ===beg===
635 @c pochhammer (11/3, -6);
636 @c ===end===
637 @example
638 (%i1) pochhammer (11/3, -6);
639 729
640 (%o1) - ----
641 2240
642 @end example
644 The unit step function is left-continuous; thus
646 @c ===beg===
647 @c [unit_step (-1/10), unit_step (0), unit_step (1/10)];
648 @c ===end===
649 @example
650 (%i1) [unit_step (-1/10), unit_step (0), unit_step (1/10)];
651 (%o1) [0, 0, 1]
652 @end example
654 If you need a unit step function that is neither left or right continuous
655 at zero, define your own using @code{signum}; for example,
657 @c ===beg===
658 @c xunit_step (x) := (1 + signum (x))/2$
659 @c [xunit_step (-1/10), xunit_step (0), xunit_step (1/10)];
660 @c ===end===
661 @example
662 (%i1) xunit_step (x) := (1 + signum (x))/2$
664 (%i2) [xunit_step (-1/10), xunit_step (0), xunit_step (1/10)];
665 1
666 (%o2) [0, -, 1]
667 2
668 @end example
670 Do not redefine @code{unit_step} itself; some code in @code{orthopoly}
671 requires that the unit step function be left-continuous.
673 @subsection Algorithms
675 Generally, @code{orthopoly} does symbolic evaluation by using a hypergeometic
676 representation of the orthogonal polynomials. The hypergeometic
677 functions are evaluated using the (undocumented) functions @code{hypergeo11}
678 and @code{hypergeo21}. The exceptions are the half-integer Bessel functions
679 and the associated Legendre function of the second kind. The half-integer Bessel functions are
680 evaluated using an explicit representation, and the associated Legendre
681 function of the second kind is evaluated using recursion.
683 For floating point evaluation, we again convert most functions into
684 a hypergeometic form; we evaluate the hypergeometic functions using
685 forward recursion. Again, the exceptions are the half-integer Bessel functions
686 and the associated Legendre function of the second kind. Numerically,
687 the half-integer Bessel functions are evaluated using the SLATEC code.
690 @node Functions and Variables for orthogonal polynomials, , Introduction to orthogonal polynomials, Package orthopoly
691 @section Functions and Variables for orthogonal polynomials
693 @anchor{assoc_legendre_p}
694 @deffn {Function} assoc_legendre_p (@var{n}, @var{m}, @var{x})
695 The associated Legendre function of the first kind of degree @math{n} and
696 order @math{m},
697 m4_mathcomma(<<<P_{n}^{m}(z)>>>,<<<assoc_legendre_p(n,m,z)>>>)
698 is a solution of the differential equation:
700 m4_displaymath(
701 <<<(1-z^2){d^2 w\over dz^2} - 2z{dw\over dz} + \left[n(n+1)-{m^2\over 1-z^2}\right] w = 0>>>,
702 <<<(1-z^2)*diff(w,z,2) - 2*z*diff(w,z) + (n*(n+1)-m^2/(1-z^2))*w = 0>>>)
704 This is related to the Legendre polynomial,
705 m4_math(<<<P_n(x)>>>, <<<legendre_p(n,x)>>>)
706 via
708 m4_displaymath(
709 <<<P_n^m(x) = (-1)^m\left(1-x^2\right)^{m/2} {d^m\over dx^m} P_n(x)>>>,
710 <<<assoc_legendre_p(n,m,x) = (-1)^m*(1-x^2)^(m/2) diff(legendre_p(n,x),x,m)>>>)
712 Reference: @urlaands{eqn 22.5.37, 779}, @urlaands{eqn 8.6.6, 334}, and @urlaands{eqn 8.2.5, 333}.
714 Some examples:
715 @c ===beg===
716 @c assoc_legendre_p(2,0,x);
717 @c factor(%);
718 @c factor(assoc_legendre_p(2,1,x));
719 @c (-1)^1*(1-x^2)^(1/2)*diff(legendre_p(2,x),1);
720 @c factor(%);
721 @c ===end===
722 @example
723 (%i1) assoc_legendre_p(2,0,x);
724 2
725 3 (1 - x)
726 (%o1) (- 3 (1 - x)) + ---------- + 1
727 2
728 (%i2) factor(%);
729 2
730 3 x - 1
731 (%o2) --------
732 2
733 (%i3) factor(assoc_legendre_p(2,1,x));
734 2
735 (%o3) - 3 x sqrt(1 - x )
737 (%i4) (-1)^1*(1-x^2)^(1/2)*diff(legendre_p(2,x),x);
738 2
739 (%o4) - (3 - 3 (1 - x)) sqrt(1 - x )
741 (%i5) factor(%);
742 2
743 (%o5) - 3 x sqrt(1 - x )
744 @end example
746 @opencatbox{Categories:}
747 @category{Package orthopoly}
748 @closecatbox
750 @end deffn
752 @anchor{assoc_legendre_q}
753 @deffn {Function} assoc_legendre_q (@var{n}, @var{m}, @var{x})
754 The associated Legendre function of the second kind of degree @math{n}
755 and order @math{m},
756 m4_mathcomma(<<<Q_{n}^{m}(z)>>>,<<<assoc_legendre_q(n,m,z)>>>)
757 is a solution of the differential equation:
759 m4_displaymath(
760 <<<(1-z^2){d^2 w\over dz^2} - 2z{dw\over dz} + \left[n(n+1)-{m^2\over 1-z^2}\right] w = 0>>>,
761 <<<(1-z^2)*diff(w,z,2) - 2*z*diff(w,z) + (n*(n+1)-m^2/(1-z^2))*w = 0>>>)
763 Reference: Abramowitz and Stegun, equation 8.5.3 and 8.1.8.
765 Some examples:
766 @c ===beg===
767 @c assoc_legendre_q(0,0,x);
768 @c assoc_legendre_q(1,0,x);
769 @c assoc_legendre_q(1,1,x);
770 @c ===end===
771 @example
772 (%i1) assoc_legendre_q(0,0,x);
773 x + 1
774 log(- -----)
775 x - 1
776 (%o1) ------------
777 2
778 (%i2) assoc_legendre_q(1,0,x);
779 x + 1
780 log(- -----) x - 2
781 x - 1
782 (%o2)/R/ ------------------
783 2
784 (%i3) assoc_legendre_q(1,1,x);
785 (%o3)/R/
786 x + 1 2 2 2 x + 1 2
787 log(- -----) sqrt(1 - x ) x - 2 sqrt(1 - x ) x - log(- -----) sqrt(1 - x )
788 x - 1 x - 1
789 - ---------------------------------------------------------------------------
790 2
791 2 x - 2
792 @end example
794 @opencatbox{Categories:}
795 @category{Package orthopoly}
796 @closecatbox
798 @end deffn
800 @anchor{chebyshev_t}
801 @deffn {Function} chebyshev_t (@var{n}, @var{x})
802 The Chebyshev polynomial of the first kind of degree @math{n},
803 m4_math(<<<T_n(x).>>>,<<<chebyshev_t(n,x).>>>)
805 Reference: @urlaands{eqn 22.5.47, 779}.
807 The polynomials
808 m4_math(<<<T_n(x)>>>,<<<chebyshev_t(n,x)>>>)
809 can be written in terms of a hypergeometric function:
811 m4_displaymath(
812 <<<T_n(x) = {_{2}}F_{1}\left(-n, n; {1\over 2}; {1-x\over 2}\right)>>>,
813 <<<hypergeometric([-n,n],[1/2],(1-x)/2)>>>)
815 The polynomials can also be defined in terms of the sum
817 m4_displaymath(
818 <<<T_n(x) = {n\over 2} \sum_{r=0}^{\lfloor {n/2}\rfloor} {(-1)^r\over n-r} {n-r\choose k}(2x)^{n-2r}>>>,
819 <<<chebyshev_t(n,x) = n/2*sum((-1)^r/(n-r)*binomial(n-r,r)*(2*x)^(n-2*r), r, 0, floor(n/2))>>>)
821 or the Rodrigues formula
823 m4_displaymath(
824 <<<T_n(x) = {1\over \kappa_n w(x)} {d^n\over dx^n}\left(w(x)(1-x^2)^n\right)>>>,
825 <<<@math{chebyshev_t(n,x) = 1/(k(n)*w(x))*diff(w(x)*(1-x^2)^n, x, n)}>>>
826 )
828 where
830 m4_displaymath(
831 <<<\eqalign{
832 w(x) &= 1/\sqrt{1-x^2} \cr
833 \kappa_n &= (-2)^n\left(1\over 2\right)_n
834 }>>>,
835 <<<@math{w(x) = 1/sqrt(1-x^2)}
837 @math{k_n = (-2)^n*pochhammer(1/2,n)}>>>)
839 Some examples:
840 @c ===beg===
841 @c chebyshev_t(2,x);
842 @c factor(%);
843 @c factor(chebyshev_t(3,x));
844 @c factor(hgfred([-3,3],[1/2],(1-x)/2));
845 @c ===end===
846 @example
847 (%i1) chebyshev_t(2,x);
848 2
849 (%o1) (- 4 (1 - x)) + 2 (1 - x) + 1
850 (%i2) factor(%);
851 2
852 (%o2) 2 x - 1
853 (%i3) factor(chebyshev_t(3,x));
854 2
855 (%o3) x (4 x - 3)
856 (%i4) factor(hgfred([-3,3],[1/2],(1-x)/2));
857 2
858 (%o4) x (4 x - 3)
859 @end example
861 @opencatbox{Categories:}
862 @category{Package orthopoly}
863 @closecatbox
865 @end deffn
867 @anchor{chebyshev_u}
868 @deffn {Function} chebyshev_u (@var{n}, @var{x})
869 The Chebyshev polynomial of the second kind of degree @math{n},
870 m4_mathdot(<<<U_n(x)>>>,<<<chebyshev_u(n,x)>>>)
872 Reference: @urlaands{eqn 22.5.48,779}.
874 The polynomials
875 m4_math(<<<U_n(x)>>>,<<<chebyshev_u(n,x)>>>)
876 can be written in terms of a hypergeometric function:
878 m4_displaymath(
879 <<<U_n(x) = (n+1)\; {_{2}F_{1}}\left(-n, n+2; {3\over 2}; {1-x\over 2}\right)>>>,
880 <<<chebyshev_u(n,x) = (n+1)*hypergeometric([-n,n+1],[3/2],(1-x)/2)>>>)
882 The polynomials can also be defined in terms of the sum
884 m4_displaymath(
885 <<<U_n(x) = \sum_{r=0}^{\lfloor n/2 \rfloor} (-1)^r {n-r \choose r} (2x)^{n-2r}>>>,
886 <<<chebyshev_u(n,x) = sum((-1)^r*binomial(n-r,r)*(2*x)^(n-2*r), r, 0, floor(n/2))>>>)
888 or the Rodrigues formula
890 m4_displaymath(
891 <<<U_n(x) = {1\over \kappa_n w(x)} {d^n\over dx^n}\left(w(x)(1-x^2)^n\right)>>>,
892 <<<@math{chebyshev_u(n,x) = 1/(k(n)*w(x))*diff(w(x)*(1-x^2)^n, x, n)}>>>
893 )
895 where
897 m4_displaymath(
898 <<<\eqalign{
899 w(x) &= \sqrt{1-x^2} \cr
900 \kappa_n &= {(-2)^n\left({3\over 2}\right)_n \over n+1}
901 }>>>,
902 <<<@math{w(x) = sqrt(1-x^2)}
904 @math{k(n) = (-2)^n*pochhammer(3/2,n)/(n+1)}>>>).
906 @c ===beg===
907 @c chebyshev_u(2,x);
908 @c expand(%);
909 @c expand(chebyshev_u(3,x));
910 @c expand(4*hgfred([-3,5],[3/2],(1-x)/2));
911 @c ===end===
912 @example
913 (%i1) chebyshev_u(2,x);
914 2
915 8 (1 - x) 4 (1 - x)
916 (%o1) 3 ((- ---------) + ---------- + 1)
917 3 3
918 (%i2) expand(%);
919 2
920 (%o2) 4 x - 1
921 (%i3) expand(chebyshev_u(3,x));
922 3
923 (%o3) 8 x - 4 x
924 (%i4) expand(4*hgfred([-3,5],[3/2],(1-x)/2));
925 3
926 (%o4) 8 x - 4 x
927 @end example
929 @opencatbox{Categories:}
930 @category{Package orthopoly}
931 @closecatbox
933 @end deffn
935 @anchor{gen_laguerre}
936 @deffn {Function} gen_laguerre (@var{n}, @var{a}, @var{x})
937 The generalized Laguerre polynomial of degree @math{n},
938 m4_mathdot(<<<L_n^{(\alpha)}(x)>>>,<<<gen_laguerre(n,a,x)>>>)
940 These can be defined by
942 m4_displaymath(
943 <<<L_n^{(\alpha)}(x) = {n+\alpha \choose n}\; {_1F_1}(-n; \alpha+1; x)>>>,
944 <<<gen_laguerre(n, a, x) = binomial(n+a,n)*hypergeometric([-n], [a+1], x)>>>)
946 The polynomials can also be defined by the sum
948 m4_displaymath(
949 <<<L_n^{(\alpha)}(x) = \sum_{k=0}^n {(\alpha + k + 1)_{n-k} \over (n-k)! k!} (-x)^k>>>,
950 <<<@math{gen_laguerre(n, a, x) = sum(pochhammer(a+k+1,n-k)/((n-k)!*k!)*(-x)^k, k, 0, n)}>>>)
952 or the Rodrigues formula
954 m4_displaymath(
955 <<<L_n^{(\alpha)}(x) = {1\over \kappa_n w(x)} {d^n\over dx^n}\left(w(x)x^n\right)>>>,
956 <<<@math{gen_laguerre(n, a, x) = 1/(k(n)*w(x))*diff(w(x)*x^n, x, n)}>>>
957 )
959 where
961 m4_displaymath(
962 <<<\eqalign{
963 w(x) &= e^{-x}x^{\alpha} \cr
964 \kappa_n &= n!
965 }>>>,
966 <<<@math{w(x) = %e^(-x)*x^a}
968 @math{k(n) = n!}>>>)
970 Reference: @urlaands{eqn 22.5.54,780}.
972 Some examples:
973 @c ===beg===
974 @c gen_laguerre(1,k,x);
975 @c gen_laguerre(2,k,x);
976 @c binomial(2+k,2)*hgfred([-2],[1+k],x);
977 @c ===end===
978 @example
979 (%i1) gen_laguerre(1,k,x);
980 x
981 (%o1) (k + 1) (1 - -----)
982 k + 1
983 (%i2) gen_laguerre(2,k,x);
984 2
985 x 2 x
986 (k + 1) (k + 2) (--------------- - ----- + 1)
987 (k + 1) (k + 2) k + 1
988 (%o2) ---------------------------------------------
989 2
990 (%i3) binomial(2+k,2)*hgfred([-2],[1+k],x);
991 2
992 x 2 x
993 (k + 1) (k + 2) (--------------- - ----- + 1)
994 (k + 1) (k + 2) k + 1
995 (%o3) ---------------------------------------------
996 2
998 @end example
1000 @opencatbox{Categories:}
1001 @category{Package orthopoly}
1002 @closecatbox
1004 @end deffn
1006 @anchor{hermite}
1007 @deffn {Function} hermite (@var{n}, @var{x})
1008 The Hermite polynomial of degree @math{n},
1009 m4_mathdot(<<<H_n(x)>>>,<<<hermite(n,x)>>>)
1011 These polynomials may be defined by a hypergeometric function
1013 m4_displaymath(
1014 <<<H_n(x) = (2x)^n\; {_2F_0}\left(-{1\over 2} n, -{1\over 2}n+{1\over 2};;-{1\over x^2}\right)>>>,
1015 <<<hermite(n,x) = (2*x)^n * hypergeometric([-n/2, -n/2+1/2],[], -1/x^2)>>>)
1017 or by the series
1019 m4_displaymath(
1020 <<<H_n(x) = n! \sum_{k=0}^{\lfloor n/2 \rfloor} {(-1)^k(2x)^{n-2k} \over k! (n-2k)!}>>>,
1021 <<<hermite(n,x) = n!*sum((-1)^k*(2*x)^(n-2*k)/(k!*(n-2*k)!), k, 0, floor(n/2))>>>)
1023 or the Rodrigues formula
1025 m4_displaymath(
1026 <<<H_n(x) = {1\over \kappa_n w(x)} {d^n\over dx^n}\left(w(x)\right)>>>,
1027 <<<@math{hermite(n,x) = 1/(k(n)*w(x))*diff(w(x), x, n)}>>>
1028 )
1030 where
1032 m4_displaymath(
1033 <<<\eqalign{
1034 w(x) &= e^{-{x^2/2}} \cr
1035 \kappa_n &= (-1)^n
1036 }>>>,
1037 <<<@math{w(x) = %e(-x^2/2)}
1039 @math{k(n) = (-1)^n}>>>)
1041 Reference: @urlaands{eqn 22.5.55,780}.
1043 Some examples:
1044 @c ===beg===
1045 @c hermite(3,x);
1046 @c expand(%);
1047 @c expand(hermite(4,x));
1048 @c expand((2*x)^4*hgfred([-2,-2+1/2],[],-1/x^2));
1049 @c expand(4!*sum((-1)^k*(2*x)^(4-2*k)/(k!*(4-2*k)!),k,0,floor(4/2)));
1050 @c ===end===
1051 @example
1052 (%i1) hermite(3,x);
1053 2
1054 2 x
1055 (%o1) - 12 x (1 - ----)
1056 3
1057 (%i2) expand(%);
1058 3
1059 (%o2) 8 x - 12 x
1060 (%i3) expand(hermite(4,x));
1061 4 2
1062 (%o3) 16 x - 48 x + 12
1063 (%i4) expand((2*x)^4*hgfred([-2,-2+1/2],[],-1/x^2));
1064 4 2
1065 (%o4) 16 x - 48 x + 12
1066 (%i5) expand(4!*sum((-1)^k*(2*x)^(4-2*k)/(k!*(4-2*k)!),k,0,floor(4/2)));
1067 4 2
1068 (%o5) 16 x - 48 x + 12
1069 @end example
1071 @opencatbox{Categories:}
1072 @category{Package orthopoly}
1073 @closecatbox
1075 @end deffn
1077 @anchor{intervalp}
1078 @deffn {Function} intervalp (@var{e})
1079 Return @code{true} if the input is an interval and return false if it isn't.
1081 @opencatbox{Categories:}
1082 @category{Package orthopoly}
1083 @category{Predicate functions}
1084 @closecatbox
1086 @end deffn
1088 @anchor{jacobi_p}
1089 @deffn {Function} jacobi_p (@var{n}, @var{a}, @var{b}, @var{x})
1090 The Jacobi polynomial,
1091 m4_mathdot(<<<P_n^{(a,b)}(x)>>>, <<<jacobi_p(n,a,b,x)>>>)
1093 The Jacobi polynomials are actually defined for all
1094 @math{a} and @math{b}; however, the Jacobi polynomial
1095 weight @math{(1 - x)^a (1 + x)^b} isn't integrable
1096 for
1097 m4_math(<<<a \le -1>>>, <<<a <= -1>>>)
1098 or
1099 m4_mathdot(<<<b \le -1>>>, <<<b <= -1>>>)
1101 Reference: @urlaands{eqn 22.5.42,779}.
1103 The polynomial may be defined in terms of hypergeometric functions:
1105 m4_displaymath(
1106 <<<P_n^{(a,b)}(x) = {n+a\choose n} {_1F_2}\left(-n, n + a + b + 1; a+1; {1-x\over 2}\right)>>>,
1107 <<<jacobi_p(n,a,b,x) = binomial(n+a,n)*hypergeometric([-n,n+a+b+1],[a+1],(1-x)/2)>>>)
1109 or the Rodrigues formula
1111 m4_displaymath(
1112 <<<P_n^{(a, b)}(x) = {1\over \kappa_n w(x)} {d^n\over dx^n}\left(w(x)\left(1-x^2\right)^n\right)>>>,
1113 <<<@math{jacobi_p(n,a,b,x) = 1/(k(n)*w(x))*diff(w(x)*(1-x^2)^n, x, n)}>>>
1114 )
1116 where
1118 m4_displaymath(
1119 <<<\eqalign{
1120 w(x) &= (1-x)^a(1-x)^b \cr
1121 \kappa_n &= (-2)^n n!
1122 }>>>,
1123 <<<@math{w(x) = (1-x)^a*(1-x)^b}
1125 @math{k(n) = (-2)^n*n!}>>>)
1127 Some examples:
1128 @c ===beg===
1129 @c jacobi_p(0,a,b,x);
1130 @c jacobi_p(1,a,b,x);
1131 @c ===end===
1132 @example
1133 (%i1) jacobi_p(0,a,b,x);
1134 (%o1) 1
1135 (%i2) jacobi_p(1,a,b,x);
1136 (b + a + 2) (1 - x)
1137 (%o2) (a + 1) (1 - -------------------)
1138 2 (a + 1)
1139 @end example
1141 @opencatbox{Categories:}
1142 @category{Package orthopoly}
1143 @closecatbox
1145 @end deffn
1147 @anchor{laguerre}
1148 @deffn {Function} laguerre (@var{n}, @var{x})
1149 The Laguerre polynomial,
1150 m4_math(<<<L_n(x)>>>,<<<laguerre(n,x)>>>)
1151 of degree @math{n}.
1153 Reference: @urlaands{eqn 22.5.16, 778} and @urlaands{eqn 22.5.54, 780}.
1155 These are related to the generalized Laguerre polynomial by
1157 m4_displaymath(
1158 <<<L_n(x) = L_n^{(0)}(x)>>>,
1159 <<<laguerre(n,x) = gen_laguerre(n,0,x)>>>)
1161 The polynomials are given by the sum
1163 m4_displaymath(
1164 <<<L_n(x) = \sum_{k=0}^{n} {(-1)^k\over k!}{n \choose k} x^k>>>,
1165 <<<laguerre(n,x) = sum((-1)^k/k!*binomial(n,k)*x^k,k,0,n)>>>)
1167 Some examples:
1168 @c ===beg===
1169 @c laguerre(1,x);
1170 @c laguerre(2,x);
1171 @c gen_laguerre(2,0,x);
1172 @c sum((-1)^k/k!*binomial(2,k)*x^k,k,0,2);
1173 @c ===end===
1174 @example
1175 (%i1) laguerre(1,x);
1176 (%o1) 1 - x
1177 (%i2) laguerre(2,x);
1178 2
1179 x
1180 (%o2) -- - 2 x + 1
1181 2
1182 (%i3) gen_laguerre(2,0,x);
1183 2
1184 x
1185 (%o3) -- - 2 x + 1
1186 2
1187 (%i4) sum((-1)^k/k!*binomial(2,k)*x^k,k,0,2);
1188 2
1189 x
1190 (%o4) -- - 2 x + 1
1191 2
1192 @end example
1194 @opencatbox{Categories:}
1195 @category{Package orthopoly}
1196 @closecatbox
1198 @end deffn
1200 @anchor{legendre_p}
1201 @deffn {Function} legendre_p (@var{n}, @var{x})
1202 The Legendre polynomial of the first kind,
1203 m4_mathcomma(P_n(x),legendre(n,x))
1204 of degree @math{n}.
1206 Reference: @urlaands{eqn 22.5.50, 779} and @urlaands{eqn 22.5.51, 779}.
1208 The Legendre polynomial is related to the Jacobi polynomials by
1210 m4_displaymath(
1211 <<<P_n(x) = P_n^{(0,0)}(x)>>>,
1212 <<<legendre_p(n,x) = jacobi_p(n,0,0,x)>>>)
1214 or the Rodrigues formula
1216 m4_displaymath(
1217 <<<P_n(x) = {1\over \kappa_n w(x)} {d^n\over dx^n}\left(w(x)\left(1-x^2\right)^n\right)>>>,
1218 <<<@math{legendre_p(n,x) = 1/(k(n)*w(x))*diff(w(x)*(1-x^2)^n, x, n)}>>>
1219 )
1221 where
1223 m4_displaymath(
1224 <<<\eqalign{
1225 w(x) &= 1 \cr
1226 \kappa_n &= (-2)^n n!
1227 }>>>,
1228 <<<@math{w(x) = 1}
1230 @math{k(n) = (-2)^n*n!}>>>)
1232 Some examples:
1233 @c ===beg===
1234 @c legendre_p(1,x);
1235 @c legendre_p(2,x);
1236 @c expand(%);
1237 @c expand(legendre_p(3,x));
1238 @c expand(jacobi_p(3,0,0,x));
1239 @c ===end===
1240 @example
1241 (%i1) legendre_p(1,x);
1242 (%o1) x
1243 (%i2) legendre_p(2,x);
1244 2
1245 3 (1 - x)
1246 (%o2) (- 3 (1 - x)) + ---------- + 1
1247 2
1248 (%i3) expand(%);
1249 2
1250 3 x 1
1251 (%o3) ---- - -
1252 2 2
1253 (%i4) expand(legendre_p(3,x));
1254 3
1255 5 x 3 x
1256 (%o4) ---- - ---
1257 2 2
1258 (%i5) expand(jacobi_p(3,0,0,x));
1259 3
1260 5 x 3 x
1261 (%o5) ---- - ---
1262 2 2
1263 @end example
1264 @opencatbox{Categories:}
1265 @category{Package orthopoly}
1266 @closecatbox
1268 @end deffn
1270 @anchor{legendre_q}
1271 @deffn {Function} legendre_q (@var{n}, @var{x})
1272 The Legendre function of the second kind,
1273 m4_math(<<<Q_n(x)>>>, <<<legendre_q(n,x)>>>)
1274 of degree @math{n}.
1276 Reference: Abramowitz and Stegun, equations 8.5.3 and 8.1.8.
1278 These are related to
1279 m4_math(<<<Q_n^m(x)>>>,<<<assoc_legendre_q(n,m,x)>>>)
1280 by
1282 m4_displaymath(
1283 <<<Q_n(x) = Q_n^0(x)>>>,
1284 <<<legendre_q(n,x) = assoc_legendre_q(n,0,x)>>>)
1286 Some examples:
1287 @c ===beg===
1288 @c legendre_q(0,x);
1289 @c legendre_q(1,x);
1290 @c assoc_legendre_q(1,0,x);
1291 @c ===end===
1292 @example
1293 (%i1) legendre_q(0,x);
1294 x + 1
1295 log(- -----)
1296 x - 1
1297 (%o1) ------------
1298 2
1299 (%i2) legendre_q(1,x);
1300 x + 1
1301 log(- -----) x - 2
1302 x - 1
1303 (%o2)/R/ ------------------
1304 2
1305 (%i3) assoc_legendre_q(1,0,x);
1306 x + 1
1307 log(- -----) x - 2
1308 x - 1
1309 (%o3)/R/ ------------------
1310 2
1311 @end example
1313 @opencatbox{Categories:}
1314 @category{Package orthopoly}
1315 @closecatbox
1317 @end deffn
1319 @anchor{orthopoly_recur}
1320 @deffn {Function} orthopoly_recur (@var{f}, @var{args})
1321 Returns a recursion relation for the orthogonal function family
1322 @var{f} with arguments @var{args}. The recursion is with
1323 respect to the polynomial degree.
1325 @c ===beg===
1326 @c orthopoly_recur (legendre_p, [n, x]);
1327 @c ===end===
1328 @example
1329 (%i1) orthopoly_recur (legendre_p, [n, x]);
1330 (2 n + 1) P (x) x - n P (x)
1331 n n - 1
1332 (%o1) P (x) = -------------------------------
1333 n + 1 n + 1
1334 @end example
1336 The second argument to @code{orthopoly_recur} must be a list with the
1337 correct number of arguments for the function @var{f}; if it isn't,
1338 Maxima signals an error.
1340 @c ===beg===
1341 @c orthopoly_recur (jacobi_p, [n, x]);
1342 @c ===end===
1343 @example
1344 (%i1) orthopoly_recur (jacobi_p, [n, x]);
1346 Function jacobi_p needs 4 arguments, instead it received 2
1347 -- an error. Quitting. To debug this try debugmode(true);
1348 @end example
1350 Additionally, when @var{f} isn't the name of one of the
1351 families of orthogonal polynomials, an error is signalled.
1353 @c ===beg===
1354 @c orthopoly_recur (foo, [n, x]);
1355 @c ===end===
1356 @example
1357 (%i1) orthopoly_recur (foo, [n, x]);
1359 A recursion relation for foo isn't known to Maxima
1360 -- an error. Quitting. To debug this try debugmode(true);
1361 @end example
1363 @opencatbox{Categories:}
1364 @category{Package orthopoly}
1365 @closecatbox
1367 @end deffn
1369 @defvr {Variable} orthopoly_returns_intervals
1370 Default value: @code{true}
1372 When @code{orthopoly_returns_intervals} is @code{true}, floating point results are returned in
1373 the form @code{interval (@var{c}, @var{r})}, where @var{c} is the center of an interval
1374 and @var{r} is its radius. The center can be a complex number; in that
1375 case, the interval is a disk in the complex plane.
1377 @opencatbox{Categories:}
1378 @category{Package orthopoly}
1379 @closecatbox
1381 @end defvr
1383 @anchor{orthopoly_weight}
1384 @deffn {Function} orthopoly_weight (@var{f}, @var{args})
1386 Returns a three element list; the first element is
1387 the formula of the weight for the orthogonal polynomial family
1388 @var{f} with arguments given by the list @var{args}; the
1389 second and third elements give the lower and upper endpoints
1390 of the interval of orthogonality. For example,
1392 @c ===beg===
1393 @c w : orthopoly_weight (hermite, [n, x]);
1394 @c integrate (w[1] * hermite (3, x) * hermite (2, x), x, w[2], w[3]);
1395 @c ===end===
1396 @example
1397 (%i1) w : orthopoly_weight (hermite, [n, x]);
1398 2
1399 - x
1400 (%o1) [%e , - inf, inf]
1401 (%i2) integrate(w[1]*hermite(3, x)*hermite(2, x), x, w[2], w[3]);
1402 (%o2) 0
1403 @end example
1405 The main variable of @var{f} must be a symbol; if it isn't, Maxima
1406 signals an error.
1408 @opencatbox{Categories:}
1409 @category{Package orthopoly}
1410 @closecatbox
1412 @end deffn
1414 @anchor{pochhammer}
1415 @deffn {Function} pochhammer (@var{x}, @var{n})
1416 The Pochhammer symbol,
1417 m4_mathdot(<<<(x)_n>>>, <<<pochhammer(x,n)>>>)
1418 (See @urlaands{eqn 6.1.22, 256} and @urldlmf{5.2.iii}).
1420 For nonnegative
1421 integers @var{n} with @code{@var{n} <= pochhammer_max_index}, the
1422 expression
1423 m4_math(<<<(x)_n>>>, <<<pochhammer(x, n)>>>)
1424 evaluates to the
1425 product
1426 m4_math(<<<x(x+1)(x+2)\cdots(x+n-1)>>>, <<<x (x + 1) (x + 2) ... (x +
1427 n - 1)>>>)
1428 when
1429 m4_math(<<<n > 0>>>, n > 0)
1430 and
1431 to 1 when @math{n = 0}.
1432 For negative @math{n},
1433 m4_math(<<<(x)_n>>>, <<<pochhammer (x, n)>>>)
1434 is
1435 defined as
1436 m4_math(<<<(-1)^n/(1-x)_{-n}.>>>, <<<(-1)^n / pochhammer (1 - x, -n).>>>)
1437 Thus
1439 @c ===beg===
1440 @c pochhammer (x, 3);
1441 @c pochhammer (x, -3);
1442 @c ===end===
1443 @example
1444 (%i1) pochhammer (x, 3);
1445 (%o1) x (x + 1) (x + 2)
1446 (%i2) pochhammer (x, -3);
1447 1
1448 (%o2) - -----------------------
1449 (1 - x) (2 - x) (3 - x)
1450 @end example
1452 To convert a Pochhammer symbol into a quotient of gamma functions,
1453 (see @urlaands{eqn 6.1.22, 256}) use @code{makegamma}; for example
1455 @c ===beg===
1456 @c makegamma (pochhammer (x, n));
1457 @c ===end===
1458 @example
1459 (%i1) makegamma (pochhammer (x, n));
1460 gamma(x + n)
1461 (%o1) ------------
1462 gamma(x)
1463 @end example
1465 When @var{n} exceeds @code{pochhammer_max_index} or when @var{n}
1466 is symbolic, @code{pochhammer} returns a noun form.
1468 @c ===beg===
1469 @c pochhammer (x, n);
1470 @c ===end===
1471 @example
1472 (%i1) pochhammer (x, n);
1473 (%o1) (x)
1474 n
1475 @end example
1477 @opencatbox{Categories:}
1478 @category{Package orthopoly}
1479 @category{Gamma and factorial functions}
1480 @closecatbox
1482 @end deffn
1484 @defvr {Variable} pochhammer_max_index
1485 Default value: 100
1487 @code{pochhammer (@var{n}, @var{x})} expands to a product if and only if
1488 @code{@var{n} <= pochhammer_max_index}.
1490 Examples:
1492 @c ===beg===
1493 @c pochhammer (x, 3), pochhammer_max_index : 3;
1494 @c pochhammer (x, 4), pochhammer_max_index : 3;
1495 @c ===end===
1496 @example
1497 (%i1) pochhammer (x, 3), pochhammer_max_index : 3;
1498 (%o1) x (x + 1) (x + 2)
1499 (%i2) pochhammer (x, 4), pochhammer_max_index : 3;
1500 (%o2) (x)
1501 4
1502 @end example
1504 Reference: @urlaands{eqn 6.1.16,256}.
1506 @opencatbox{Categories:}
1507 @category{Package orthopoly}
1508 @category{Gamma and factorial functions}
1509 @closecatbox
1511 @end defvr
1513 @anchor{spherical_bessel_j}
1514 @deffn {Function} spherical_bessel_j (@var{n}, @var{x})
1515 The spherical Bessel function of the first kind,
1516 m4_math(<<<j_n(x).>>>, <<<spherical_bessel_j(n,x).>>>)
1518 Reference: @urlaands{eqn 10.1.8, 437} and @urlaands{eqn 10.1.15, 439}.
1520 It is related to the Bessel function by
1522 m4_displaymath(
1523 <<<j_n(x) = \sqrt{\pi\over 2x} J_{n+1/2}(x)>>>,
1524 <<<spherical_bessel_j(n,x) = sqrt(%pi/(2*x))*bessel_j(n+1/2,x)>>>)
1526 Some examples:
1527 @c ===beg===
1528 @c spherical_bessel_j(1,x);
1529 @c spherical_bessel_j(2,x);
1530 @c expand(%);
1531 @c expand(sqrt(%pi/(2*x))*bessel_j(2+1/2,x)),besselexpand:true;
1532 @c ===end===
1533 @example
1534 (%i1) spherical_bessel_j(1,x);
1535 sin(x)
1536 ------ - cos(x)
1537 x
1538 (%o1) ---------------
1539 x
1540 (%i2) spherical_bessel_j(2,x);
1541 3 3 cos(x)
1542 (- (1 - --) sin(x)) - --------
1543 2 x
1544 x
1545 (%o2) ------------------------------
1546 x
1547 (%i3) expand(%);
1548 sin(x) 3 sin(x) 3 cos(x)
1549 (%o3) (- ------) + -------- - --------
1550 x 3 2
1551 x x
1552 (%i4) expand(sqrt(%pi/(2*x))*bessel_j(2+1/2,x)),besselexpand:true;
1553 sin(x) 3 sin(x) 3 cos(x)
1554 (%o4) (- ------) + -------- - --------
1555 x 3 2
1556 x x
1557 @end example
1558 @opencatbox{Categories:}
1559 @category{Package orthopoly}
1560 @category{Bessel functions}
1561 @closecatbox
1563 @end deffn
1565 @anchor{spherical_bessel_y}
1566 @deffn {Function} spherical_bessel_y (@var{n}, @var{x})
1567 The spherical Bessel function of the second kind,
1568 m4_math(<<<y_n(x).>>>, <<<<spherical_bessel_y(n,x).>>>)
1570 Reference: @urlaands{eqn 10.1.9, 437} and @urlaands{eqn 10.1.15, 439}.
1572 It is related to the Bessel function by
1574 m4_displaymath(
1575 <<<y_n(x) = \sqrt{\pi\over 2x} Y_{n+1/2}(x)>>>,
1576 <<<spherical_bessel_y(n,x) = sqrt(%pi/(2*x))*bessel_y(n+1/2,x)>>>)
1578 @c ===beg===
1579 @c spherical_bessel_y(1,x);
1580 @c spherical_bessel_y(2,x);
1581 @c expand(%);
1582 @c expand(sqrt(%pi/(2*x))*bessel_y(2+1/2,x)),besselexpand:true;
1583 @c ===end===
1584 @example
1585 (%i1) spherical_bessel_y(1,x);
1586 cos(x)
1587 (- sin(x)) - ------
1588 x
1589 (%o1) -------------------
1590 x
1591 (%i2) spherical_bessel_y(2,x);
1592 3 sin(x) 3
1593 -------- - (1 - --) cos(x)
1594 x 2
1595 x
1596 (%o2) - --------------------------
1597 x
1598 (%i3) expand(%);
1599 3 sin(x) cos(x) 3 cos(x)
1600 (%o3) (- --------) + ------ - --------
1601 2 x 3
1602 x x
1603 (%i4) expand(sqrt(%pi/(2*x))*bessel_y(2+1/2,x)),besselexpand:true;
1604 3 sin(x) cos(x) 3 cos(x)
1605 (%o4) (- --------) + ------ - --------
1606 2 x 3
1607 x x
1608 @end example
1609 @opencatbox{Categories:}
1610 @category{Package orthopoly}
1611 @category{Bessel functions}
1612 @closecatbox
1614 @end deffn
1616 @anchor{spherical_hankel1}
1617 @deffn {Function} spherical_hankel1 (@var{n}, @var{x})
1618 The spherical Hankel function of the
1619 first kind,
1620 m4_math(<<<h_n^{(1)}(x).>>>, <<<spherical_hankel1(n,x).>>>)
1622 Reference: @urlaands{eqn 10.1.36,439}.
1624 This is defined by
1626 m4_displaymath(
1627 <<<h_n^{(1)}(x) = j_n(x) + iy_n(x)>>>,
1628 <<<spherical_hankel1(n,x) = spherical_bessel_j(n,x) + %i*spherical_bessel_y(n,x)>>>)
1630 @opencatbox{Categories:}
1631 @category{Package orthopoly}
1632 @category{Bessel functions}
1633 @closecatbox
1635 @end deffn
1637 @anchor{spherical_hankel2}
1638 @deffn {Function} spherical_hankel2 (@var{n}, @var{x})
1639 The spherical Hankel function of the
1640 second kind,
1641 m4_math(<<<h_n^{(2)}(x).>>>, <<<spherical_hankel2(n,x).>>>)
1643 Reference: @urlaands{eqn 10.1.17,439}.
1645 This is defined by
1647 m4_displaymath(
1648 <<<h_n^{(2)}(x) = j_n(x) + iy_n(x)>>>,
1649 <<<spherical_hankel2(n,x) = spherical_bessel_j(n,x) - %i*spherical_bessel_y(n,x)>>>)
1651 @opencatbox{Categories:}
1652 @category{Package orthopoly}
1653 @category{Bessel functions}
1654 @closecatbox
1656 @end deffn
1658 @anchor{spherical_harmonic}
1659 @deffn {Function} spherical_harmonic (@var{n}, @var{m}, @var{theta}, @var{phi})
1660 The spherical harmonic function,
1661 m4_mathdot(<<<Y_n^m(\theta, \phi)>>>, <<<spherical_harmonic(n,m,theta,phi)>>>)
1663 Spherical harmonics satisfy the angular part of Laplace's equation in spherical coordinates.
1665 For integers @math{n} and @math{m} such
1666 that
1667 m4_math(<<<n \geq |m|>>>,<<<n <= abs(m)>>>)
1668 and
1669 for
1670 m4_mathdot(<<<\theta \in [0, \pi]>>>,<<<theta in [0, %pi]>>>)
1671 Maxima’s
1672 spherical harmonic function can be defined by
1674 m4_displaymath(
1675 <<<Y_n^m(\theta, \phi) = (-1)^m \sqrt{{2n+1\over 4\pi} {(n-m)!\over (n+m)!}} P_n^m(\cos\theta) e^{im\phi}>>>,
1676 <<<@math{spherical_harmonic(n, m, theta, phi) := (-1)^m * (((n-m)!*(2*n+1))/(4*%pi*(n+m)!))^(1/2)* exp(%i*m*phi)*assoc_legendre_p(n,m,cos(theta))}>>>)
1679 Further, when
1680 m4_mathcomma(<<<n < |m|>>>, <<<n < abs(m)>>>)
1681 the
1682 spherical harmonic function vanishes.
1684 The factor @math{(-1)^m}, frequently used in Quantum mechanics, is
1685 called the @url{https://en.wikipedia.org/wiki/Spherical_harmonics#Condon%E2%80%93Shortley_phase, Condon-Shortely phase}.
1686 Some references, including @emph{NIST Digital Library of Mathematical Functions} omit
1687 this factor; see @url{http://dlmf.nist.gov/14.30.E1}.
1689 Reference: Merzbacher 9.64.
1691 Some examples:
1692 @c ===beg===
1693 @c spherical_harmonic(1,0,theta,phi);
1694 @c spherical_harmonic(1,1,theta,phi);
1695 @c spherical_harmonic(1,-1,theta,phi);
1696 @c spherical_harmonic(2,0,theta,phi);
1697 @c factor(%);
1698 @c ===end===
1699 @example
1700 (%i1) spherical_harmonic(1,0,theta,phi);
1701 sqrt(3) cos(theta)
1702 (%o1) ------------------
1703 2 sqrt(%pi)
1704 (%i2) spherical_harmonic(1,1,theta,phi);
1705 %i phi
1706 sqrt(3) %e sin(theta)
1707 (%o2) ---------------------------
1708 3/2
1709 2 sqrt(%pi)
1710 (%i3) spherical_harmonic(1,-1,theta,phi);
1711 - %i phi
1712 sqrt(3) %e sin(theta)
1713 (%o3) - -----------------------------
1714 3/2
1715 2 sqrt(%pi)
1716 (%i4) spherical_harmonic(2,0,theta,phi);
1717 2
1718 3 (1 - cos(theta))
1719 sqrt(5) ((- 3 (1 - cos(theta))) + ------------------- + 1)
1720 2
1721 (%o4) ----------------------------------------------------------
1722 2 sqrt(%pi)
1723 (%i5) factor(%);
1724 2
1725 sqrt(5) (3 cos (theta) - 1)
1726 (%o5) ---------------------------
1727 4 sqrt(%pi)
1728 @end example
1730 @opencatbox{Categories:}
1731 @category{Package orthopoly}
1732 @closecatbox
1734 @end deffn
1736 @anchor{unit_step}
1737 @deffn {Function} unit_step (@var{x})
1738 The left-continuous unit step function; thus
1739 @code{unit_step (@var{x})} vanishes for @code{@var{x} <= 0} and equals
1740 1 for @code{@var{x} > 0}.
1742 If you want a unit step function that takes on the value 1/2 at zero,
1743 use @mrefdot{hstep}
1745 @opencatbox{Categories:}
1746 @category{Package orthopoly}
1747 @category{Mathematical functions}
1748 @closecatbox
1750 @end deffn
1752 @anchor{ultraspherical}
1753 @deffn {Function} ultraspherical (@var{n}, @var{a}, @var{x})
1754 The ultraspherical polynomial,
1755 m4_math(<<<C_n^{(a)}(x)>>>, <<<ultraspherical(n,a,x)>>>)
1756 (also known as the Gegenbauer polynomial).
1758 Reference: @urlaands{eqn 22.5.46,779}.
1760 These polynomials can be given in terms of Jacobi polynomials:
1762 m4_displaymath(
1763 <<<C_n^{(\alpha)}(x) = {\Gamma\left(\alpha + {1\over 2}\right) \over \Gamma(2\alpha)}
1764 {\Gamma(n+2\alpha) \over \Gamma\left(n+\alpha + {1\over 2}\right)}
1765 P_n^{(\alpha-1/2, \alpha-1/2)}(x)>>>,
1766 <<<ultraspherical(n,a,x) = gamma(a+1/2)/gamma(2*a)*gamma(n+2*a)/gamma(n+a+1/2)*jacobi_p(n,a-1/2,a-1/2,x)>>>)
1768 or the series
1770 m4_displaymath(
1771 <<<C_n^{(\alpha)}(x) = \sum_{k=0}^{\lfloor n/2 \rfloor} {(-1)^k (\alpha)_{n-k} \over k! (n-2k)!}(2x)^{n-2k}>>>,
1772 <<<@math{ultraspherical(n,a,x) = sum((-1)^k*pochhammer(a,n-k)/k!/(n-2*k)!*(2*x)^(n-2*k),k, 0, floor(n/2))}>>>)
1774 or the Rodrigues formula
1776 m4_displaymath(
1777 <<<C_n^{(\alpha)}(x) = {1\over \kappa_n w(x)} {d^n\over dx^n}\left(w(x)\left(1-x^2\right)^n\right)>>>,
1778 <<<@math{ultraspherical(n,x) = 1/(k(n)*w(x))*diff(w(x)*(1-x^2)^n, x, n)}>>>
1779 )
1781 where
1783 m4_displaymath(
1784 <<<\eqalign{
1785 w(x) &= \left(1-x^2\right)^{\alpha-{1\over 2}} \cr
1786 \kappa_n &= {(-2)^n\left(\alpha + {1\over 2}\right)_n n!\over (2\alpha)_n} \cr
1787 }>>>,
1788 <<<@math{w(x) = (1-x^2)^(a-1/2)}
1790 @math{k(n) = (-2)^n*pochhammer(a+1/2,n)*n!/pochhammer(2*a,n)}>>>)
1792 Some examples:
1793 @c ===beg===
1794 @c ultraspherical(1,a,x);
1795 @c factor(%);
1796 @c factor(ultraspherical(2,a,x));
1797 @c ===end===
1798 @example
1799 (%i1) ultraspherical(1,a,x);
1800 (2 a + 1) (1 - x)
1801 (%o1) 2 a (1 - -----------------)
1802 1
1803 2 (a + -)
1804 2
1805 (%i2) factor(%);
1806 (%o2) 2 a x
1807 (%i3) factor(ultraspherical(2,a,x));
1808 2 2
1809 (%o3) a (2 a x + 2 x - 1)
1810 @end example
1812 @opencatbox{Categories:}
1813 @category{Package orthopoly}
1814 @closecatbox
1816 @end deffn