From: Robert D. <rob...@us...> - 2005-11-27 05:21:40
|
Update of /cvsroot/maxima/maxima/share/orthopoly In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv18370/orthopoly Added Files: README changes_0.9_to_0.94.txt h_atom.dem maxima-init.lisp orthodisc.lisp orthopoly-init.lisp orthopoly.lisp orthopoly_doc.html orthopoly_doc.texi orthopoly_doc_toc.html test_orthopoly.mac variational_method.dem Log Message: Committing contents of orthopoly-0.94.tar.gz as retrieved from home page of Barton Willis (http://www.unk.edu/facstaff/profiles/willisb/), verbatim, except for excluding variational_method.dem~ (appears to be a backup of variational_method.dem created by a text editor). --- NEW FILE: README --- Copyright (C) 2000, 2001, 2003 Barton Willis, University of Nebraska at Kearney (aka UNK) orthopoly is the new name for specfun. Version 0.94 is a preliminary version. You should expect it to have bugs. To install orthopoly, follow the directions in orthopoly_doc.html. This library is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public License for more details. You should have received a copy of the GNU General Public License along with this library; if not, write to the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. --- NEW FILE: changes_0.9_to_0.94.txt --- 29 August 2003 1. Replaced kron_delta with the one from nset. 2. Made kron_delta and unit_step functions autoload -- actually the simplifying routines for these functions autoload. 3. Added test for (consp (car e)) in intervalp; new function is (defun $intervalp (a) (and (consp a) (consp (car a)) (eq (caar a) '$interval))) 4. Added tex support for unit_step; (C1) unit_step(x); (D1) UNIT_STEP(x) (C2) tex(%); $$\Theta\left(x\right)$$ 4. Made version 0.94 available. 16 July 2003 1. Changed version identifier from 0.93 to 0.94. 2. Fixed bug in summation representation for hypergeo11 -- the arguments to pochhammer were backwards. This bug wasn't detected by anything in test_orthopoly. Yikes. Added some tests for this case in test_orthopoly.mac. All 146 tests pass. 3. At least for now, I now allow hypergeo11 and hypergeo21 to return summations regardless of the main variable. Specifically, I changed use-hypergeo to (defun use-hypergeo (n x) (declare (ignore x)) (or (and (integerp n) (> n -1)) ($featurep n '$integer))) ; (and ($featurep n '$integer) ($ratnump x)))) ------------------------------------------------------------- 9 July 2003 1. Announced version 0.93. ----------------------------------------------------------------- 7 July 2003 1. Added check in orthopoly_weight that checks that all parameters (non-main variables) are free of the integration variable. 2. Added an autoload for orthopoly_weight in orthopoly-init.lisp. 3. Now allow the main variable in orthopoly_weight to be a subscripted variable. 4. Fixed A & S references in documentation. Fixed some spelling errors in documentation. -------------------------------------------------------------------- 6 July 2003 1. Added orthopoly_weight. And that's it folks! No more new functions allowed in orthopoly. Yeah. Added documentation for orthopoly_weight and added checks in test_orthopoly.mac. 2. Changed a (= ...) to a (like ...) in hypergeo21 and hypergeo11 Specifically ((and (integerp n) (> n -1)) (cond ((and (like a (- n)) (use-float b c x)) Now a call such as hypergeo21(a,b,c,x,inf) works okay; however, orthopoly never calls hypergeo21 or hypergeo11 with a symbolic first argument. --------------------------------------------------------------------- 30 June 2003 1. Added tests for kron_delta. 2. Added orthopoly_recur function; added tests in test_orthopoly for the recursion relations; added user documentation for orthopoly_recur; added autoload for orthopoly_recur. 3. Added brief documentation about the algorithms used by orthopoly. 4. Changed "laaguerre" to "laguerre" in gradef for laguerre. -------------------------------------------------------------- 27 June 2003 1. Added note in documentation about loading orthopoly before using any functions; otherwise, a user can be stuck with needing to use upper case function identifiers -- yech. This doesn't apply if the initialization file uses the appropriate setup_autoload statements. 2. Added note in documentation about numerical evaluation of pochhammer for orders exceeding pochhammer_max_index. 3. Added a new function pochhammer-quotient -- it computes the quotient of two pochhammer functions with the same order but different arguments. The new function is used in a few places to compute normalizations. This new scheme prevents some overflow and (I suppose) underflows that could happen for orthogonal polynomials with large order. 4. Added gradef for pochhammer. I expressed the derivatives in terms of the psi function -- when n is a positive integer psi(x+n) - psi(x) is expressible as a finite sum. I didn't use a sum representation because of the trouble of introducing a dummy summation index. The downside of this is Maxima doesn't recognize that psi(x+1)-psi(x) = 1/x, psi(x+2)-psi(x) = .... 5. Bug in pochhammer with negative order; correct relation is pochhammer(x,-n) = (-1)^n / pochhammer(1-x,n). Added tests for pochhammer and its gradient. 6. Extended makegamma1 (in asum.lisp) to be able to convert pochhammer symbols to gamma form; thus pochhammer(x,n) => gamma(x+n) / gamma(x). 7. Re-did documentation on pochhammer. ------------------------------------------------------- 16 June 2003 1. Changed version identifier from 0.92 to 0.93. 2. Replaced hyphens in file names to underscores. 3. Added note in documentation that orthopoly functions map over lists and matrices. 4. Added note in documentation about the shifted functions. 5. Fixed bug in pochhammer; the option variable pochhammer_max_index was being ignored. Added documentation on pochhammer_max_index. 6. Added section in documentation on how to use plot together with orthopoly. ------------------------------------------------------------------- 7 June 2003 1. Fixed a misplaced absolute value in interval-mult. ------------------------------------------------------------------ 14 May 2003 1. New algorithm for spherical bessel functions; they are based on A&S 10.1.8 and 10.1.9. Much faster than old code; spherical_bessel_j now traps for a zero argument and no longer gives a divide by zero for that case. Purloined code from specfun.mac 2. csign doesn't accept CRE arguments; added $ratdisrep to csign in spherical_bessel_j. Do I need $ratdisrep on other csign calls? 3. New algorithm for the associated Legendre polynomials of the second kind. Uses forward degree recursion plus it special-cases asssoc_legendre_q(n,n,x). 4. Changed my generic-autoload function to write the file name to *standard-output* as it loads; nice when you have multiple versions and you're not sure which version gets loaded. --------------------------------------------------------------- 11 May 2003 1. Switched to use bessel.lisp to support the spherical bessel functions. The function $bessel_j has a bug in it; placed a fix in orthopoly-init.lisp. Before, floating point evaluation for spherical_bessel_j and spherical_bessel_y were done via the spherical hankel functions; this approach is floating point "hostile." Now we use slatec code for floating point evaluation -- downside, we don't get an error estimate. For now, I still use my hypergeo21-float function for the spherical hankel functions. ---------------------------------------------------------------- 9 May 2003 1. Added note to user documentation about how apply float to something like assoc_legendre_p(n,2,1/2) is hazardous. 2. Restored the summation representation feature. Added tests for the summation representation. Now I demand that the "main" variable be a constant in order to get a sum representation; this eliminates the problem with the bug / feature ev(sum(x^k,k,0,n),x=0) => 0. 3. Now I define assoc_legendre_q(n,m,x) only for |m| <= n. This solves a problem with the gradef for assoc_legendre_q that gave incorrect values for |m| > n. Note the Wronskian of P^m_n and Q^m_n suggests that we shouldn't allow Q^m_n with |m| > n. ---------------------------------------------------------------- 6 May 2003 1. Replaced code for assoc_legendre_q by a more efficient algorithm. Added more tests for the associated legendre q functions (A & S 8.5.3). 2. Added gradef for legendre_q and a test for it in test_orthopoly. The gradef involves a kron_delta. Removed the claim in the documentation that orthopoly doesn't use kron_delta. 3. Added gradef for assoc_legendre_q and a test for it in test_orthopoly. ---------------------------------------------------------------- 1 May 2003 1. Made a few changes to user documentation in (I hope) increase clarity of installation instructions. 2. Changed file name maxima-init.lisp to orthopoly-init.lisp. 30 April 2003 1. Floating point values returned by assoc_legendre_p have errors that are too large. Fix: missing factor of machine epsilon multiplying dx. ------------------------------------------------------------------ 30 April 2003 1. Fixed miscellaneous grammar problems in the user documentation. --- NEW FILE: h_atom.dem --- if get('orthopoly,'version) = 'false then load("orthopoly")$ showtime : all; /* Let %a0 be the Bohr radius, mu the reduced mass, and %eo the electron charge. Declare these constants to be positive. */ assume(%a0 > 0, mu > 0, %hbar > 0, %eo > 0); /* Merzbacher and orthopoly use different normalizations for the generalized Laguerre polynomials. The function convert_to_mertz supplies a factor that converts from orthopoly normalization to Merzbacher. The hydrogen atom eigenfunctions are given by */ convert_to_mertz(n,l) := ((n+l)!^2 / ((n-l-1)! * (2*l+1)!)) / binomial(n+l, n-l-1); psi(n,l,m) := sqrt((2 / (%a0 * n))^3 * (n-l-1)! / (2*n * ((n+l)!)^3)) * (2 * r / (%a0 * n))^l * exp(-r / (%a0 * n)) * gen_laguerre(n-l-1,2*l+1,2*r / (%a0 * n)) * spherical_harmonic(l,m, theta,phi) * convert_to_mertz(n,l) ; /* Define the L2 inner product using a matchfix operator <<f,g>>. We'll need a conjugate function. */ my_conjugate(e) := subst(-%i, %i,e); matchfix("<<",">>"); "<<"(f,g) := integrate(integrate(integrate(my_conjugate(f)*g*r^2 * sin(theta), theta,0,%pi),phi,0,2*%pi),r,0,inf); /* Find <r> ,<r^2>, and <r^2> -<r>^2 for n=1,2,3, and 4 states. */ makelist(makelist(<<psi(n,l,0), r * psi(n,l,0)>>,l,0,n-1),n,1,4); makelist(makelist(<<psi(n,l,0), r^2 * psi(n,l,0)>>,l,0,n-1),n,1,4); makelist(makelist(<<psi(n,l,0), r^2 * psi(n,l,0)>> - <<psi(n,l,0), r * psi(n,l,0)>>^2,l,0,n-1),n,1,4); /* Find the energies of the n=1,2,3, and 4 states. Let ham be the hamiltonian. The function replace_hbar(e) replaces hbar by %eo * sqrt(mu * %a0). */ ham(e) := (-%hbar^2 / (2*mu)) * (diff(r^2 * diff(e,r),r) / r^2 + diff(sin(theta)*diff(e,theta),theta)/(r^2 * sin(theta)) + diff(e,phi,2)/(r^2 * sin(theta)^2)) - %eo^2 *e / r; replace_hbar(e) := radcan(subst(%hbar=%eo * sqrt(mu * %a0),e)); replace_hbar(makelist(makelist(<<psi(n,l,0), ham(psi(n,l,0))>>,l,0,n-1),n,1,4)); /* Use degenerate perturbational methods to study the n = 3 Stark hamiltonian; this is problem 17.5 in Merzbacher; alpha is the reduced electric field strength. */ n : 3; e : makelist(makelist(psi(n,l,m),m,-l,l),l,0,n-1)$ e : apply(append,e)$ stark_pot : alpha * r * cos(theta); overlap(f,g) := <<f, stark_pot * g>>; m : apply(matrix, outermap(overlap, e,e))$ load("eigen"); energy_shifts : eigenvalues(m); --- NEW FILE: maxima-init.lisp --- ;;---------------------------------------- ;; Optional code starts here ;; The Maxima fixfloat function is buggy. Here is a quick fix. At one ;; time, test_orthopoly needed this fix -- I think it's no longer needed. #+cl (defun fixfloat (x) (cond ((floatp x) (setq x (rationalize x)) (cons (numerator x) (denominator x))) (t x))) ;; fix for apply(rat(f),[x]) bug. (defun error-size (exp) (setq exp (ratdisrep exp)) (if (atom exp) 0 (do ((l (cdr exp) (cdr l)) (n 1 (incf n (+ 1 (error-size (car l)))))) ((or (null l) (> n $error_size)) n) (declare (fixnum n))))) ;; If file fn (with _no_ file extension) is out of date, compile it. We ;; assume that the extension of the source file is "lisp". (defun $make (fn) (setq fn (string-left-trim "&" (string fn))) (let ((fl) (fb) (bin-ext #+gcl "o" #+cmu (c::backend-fasl-file-type c::*target-backend*) #+clisp "fas" #+allegro "fasl" #-(or gcl cmu clisp allegro) "")) (setq fl (concatenate 'string (coerce fn 'string) ".lisp")) (setq fl ($file_search fl)) (setq fb (concatenate 'string (coerce fn 'string) "." bin-ext)) (setq fb ($file_search fb)) (if (or (and fl (not fb)) (and fl fb (> (file-write-date fl) (file-write-date fb)))) ($compile_file fl)))) ;; If you place orthopoly.lisp in a directory that Maxima can't find, ;; you'll need to append its directory to $file_search_lisp. To do this, ;; replace the path "/home/barton/orthopoly-0.93" with the correct path ;; for your machine. ;; The make function recompiles orthopoly automatically on startup if ;; the binary file is older than orthopoly.lisp. Unless you ;; anticipate making changes to orthopoly, you may comment out ;; ($make "orthopoly") (eval-when (load compile eval) (setq $file_search_lisp ($cons "/home/barton/orthopoly-0.94/###.{x86f,lisp}" $file_search_lisp)) ($make "orthopoly")) ;; Load file f with verbose and print bound to true. Useful for loading ;; code with syntax errors. (defun $xload (f) (load ($file_search f) :verbose 't :print 't)) (defprop $entier tex-matchfix tex) (defprop $entier (("\\lfloor ") " \\rfloor") texsym) ;;----------------------------------------------------------- ;; Essential code starts here ;; This fixes a bug in $bessel_j. (defun $bessel_j (arg order) (if (and (numberp order) (numberp ($realpart arg)) (numberp ($imagpart arg))) ;($bessel (complex ($realpart arg) ($imagpart arg)) order) <== bug ($bessel arg order) (subfunmakes '$bessel_j (ncons order) (ncons arg)))) ;; This version of setup_autoload uses file_search. (defun $setup_autoload (filename &rest functions) (let ((file ($file_search filename))) (dolist (func functions) (nonsymchk func '$setup_autoload) (putprop (setq func (dollarify-name func)) file 'autoload) (add2lnc func $props))) '$done) ;; This version of generic-autoload in Maxima 5.9.0 doesn't recognize ;; "fas", "fasl" and "x86f" as valid extensions for compiled Lisp ;; code. #+cl (defun generic-autoload (file &aux type) (setq file (pathname (cdr file))) (setq type (pathname-type file)) (let ((bin-ext #+gcl "o" #+cmu (c::backend-fasl-file-type c::*target-backend*) #+clisp "fas" #+allegro "fasl" #-(or gcl cmu clisp allegro) "")) (if (member type (list bin-ext "lisp" "lsp") :test 'equalp) (load file :verbose 't) ($batchload file)))) ;; Extended version of makegamma1 that converts pochhammer symbols ;; to quotients of gamma functions. (defun makegamma1 (e) (cond ((atom e) e) ((eq (caar e) 'mfactorial) (list '(%gamma) (list '(mplus) 1 (makegamma1 (cadr e))))) ;;------start of new code---------------------------------- ;; Do pochhammer(x,n) ==> gamma(x+n)/gamma(x). ((eq (caar e) '$pochhammer) (let ((x (makegamma1 (nth 1 e))) (n (makegamma1 (nth 2 e)))) `((mtimes) ((%gamma) ((mplus) ,x ,n)) ((mexpt) ((%gamma) ,x) -1)))) ((eq (caar e) '%genfact) (let ((x (makegamma1 (nth 1 e))) (y (makegamma1 (nth 2 e))) (z (makegamma1 (nth 3 e)))) (setq x (add (div x z) 1)) (setq y (simplify `(($entier) ,y))) (setq z (power z y)) `((mtimes) ,z ((%gamma) ,x) ((mexpt) ((%gamma) ((mplus) ,x ((mtimes) -1 ,y))) -1)))) ;;-----that's all folks ... end of new code------------------ ((eq (caar e) '%elliptic_kc) ;; Complete elliptic integral of the first kind (cond ((alike1 (cadr e) '((rat simp) 1 2)) ;; K(1/2) = gamma(1/4)/4/sqrt(pi) '((mtimes simp) ((rat simp) 1 4) ((mexpt simp) $%pi ((rat simp) -1 2)) ((mexpt simp) ((%gamma simp) ((rat simp) 1 4)) 2))) ((or (alike1 (cadr e) '((mtimes simp) ((rat simp) 1 4) ((mplus simp) 2 ((mexpt simp) 3 ((rat simp) 1 2))))) (alike1 (cadr e) '((mplus simp) ((rat simp) 1 2) ((mtimes simp) ((rat simp) 1 4) ((mexpt simp) 3 ((rat simp) 1 2))))) (alike1 (cadr e) ;; 1/(8-4*sqrt(3)) '((mexpt simp) ((mplus simp) 8 ((mtimes simp) -4 ((mexpt simp) 3 ((rat simp) 1 2)))) -1))) ;; K((2+sqrt(3)/4)) '((mtimes simp) ((rat simp) 1 4) ((mexpt simp) 3 ((rat simp) 1 4)) ((mexpt simp) $%pi ((rat simp) -1 2)) ((%gamma simp) ((rat simp) 1 6)) ((%gamma simp) ((rat simp) 1 3)))) ((or (alike1 (cadr e) ;; (2-sqrt(3))/4 '((mtimes simp) ((rat simp) 1 4) ((mplus simp) 2 ((mtimes simp) -1 ((mexpt simp) 3 ((rat simp) 1 2)))))) (alike1 (cadr e) ;; 1/2-sqrt(3)/4 '((mplus simp) ((rat simp) 1 2) ((mtimes simp) ((rat simp) -1 4) ((mexpt simp) 3 ((rat simp) 1 2))))) (alike (cadr e) ;; 1/(4*sqrt(3)+8) '((mexpt simp) ((mplus simp) 8 ((mtimes simp) 4 ((mexpt simp) 3 ((rat simp) 1 2)))) -1))) ;; K((2-sqrt(3))/4) '((mtimes simp) ((rat simp) 1 4) ((mexpt simp) 3 ((rat simp) -1 4)) ((mexpt simp) $%pi ((rat simp) -1 2)) ((%gamma simp) ((rat simp) 1 6)) ((%gamma simp) ((rat simp) 1 3)))) ((or ;; (3-2*sqrt(2))/(3+2*sqrt(2)) (alike1 (cadr e) '((mtimes simp) ((mplus simp) 3 ((mtimes simp) -2 ((mexpt simp) 2 ((rat simp) 1 2)))) ((mexpt simp) ((mplus simp) 3 ((mtimes simp) 2 ((mexpt simp) 2 ((rat simp) 1 2)))) -1))) ;; 17 - 12*sqrt(2) (alike1 (cadr e) '((mplus simp) 17 ((mtimes simp) -12 ((mexpt simp) 2 ((rat simp) 1 2))))) ;; (2*SQRT(2) - 3)/(2*SQRT(2) + 3) (alike1 (cadr e) '((mtimes simp) -1 ((mplus simp) -3 ((mtimes simp) 2 ((mexpt simp) 2 ((rat simp) 1 2)))) ((mexpt simp) ((mplus simp) 3 ((mtimes simp) 2 ((mexpt simp) 2 ((rat simp) 1 2)))) -1)))) '((mtimes simp) ((rat simp) 1 8) ((mexpt simp) 2 ((rat simp) -1 2)) ((mplus simp) 1 ((mexpt simp) 2 ((rat simp) 1 2))) ((mexpt simp) $%pi ((rat simp) -1 2)) ((mexpt simp) ((%gamma simp) ((rat simp) 1 4)) 2))) (t ;; Give up e))) ((eq (caar e) '%elliptic_ec) ;; Complete elliptic integral of the second kind (cond ((alike1 (cadr e) '((rat simp) 1 2)) ;; 2*E(1/2) - K(1/2) = 2*%pi^(3/2)*gamma(1/4)^(-2) '((mplus simp) ((mtimes simp) ((mexpt simp) $%pi ((rat simp) 3 2)) ((mexpt simp) ((%gamma simp irreducible) ((rat simp) 1 4)) -2)) ((mtimes simp) ((rat simp) 1 8) ((mexpt simp) $%pi ((rat simp) -1 2)) ((mexpt simp) ((%gamma simp) ((rat simp) 1 4)) 2)))) ((or (alike1 (cadr e) '((mtimes simp) ((rat simp) 1 4) ((mplus simp) 2 ((mtimes simp) -1 ((mexpt simp) 3 ((rat simp) 1 2)))))) (alike1 (cadr e) '((mplus simp) ((rat simp) 1 2) ((mtimes simp) ((rat simp) -1 4) ((mexpt simp) 3 ((rat simp) 1 2)))))) ;; E((2-sqrt(3))/4) ;; ;; %pi/4/sqrt(3) = K*(E-(sqrt(3)+1)/2/sqrt(3)*K) '((mplus simp) ((mtimes simp) ((mexpt simp) 3 ((rat simp) -1 4)) ((mexpt simp) $%pi ((rat simp) 3 2)) ((mexpt simp) ((%gamma simp) ((rat simp) 1 6)) -1) ((mexpt simp) ((%gamma simp) ((rat simp) 1 3)) -1)) ((mtimes simp) ((rat simp) 1 8) ((mexpt simp) 3 ((rat simp) -3 4)) ((mexpt simp) $%pi ((rat simp) -1 2)) ((%gamma simp) ((rat simp) 1 6)) ((%gamma simp) ((rat simp) 1 3))) ((mtimes simp) ((rat simp) 1 8) ((mexpt simp) 3 ((rat simp) -1 4)) ((mexpt simp) $%pi ((rat simp) -1 2)) ((%gamma simp) ((rat simp) 1 6)) ((%gamma simp) ((rat simp) 1 3))))) ((or (alike1 (cadr e) '((mtimes simp) ((rat simp) 1 4) ((mplus simp) 2 ((mexpt simp) 3 ((rat simp) 1 2))))) (alike1 (cadr e) '((mplus simp) ((rat simp) 1 2) ((mtimes simp) ((rat simp) 1 4) ((mexpt simp) 3 ((rat simp) 1 2)))))) ;; E((2+sqrt(3))/4) ;; ;; %pi*sqrt(3)/4 = K1*(E1-(sqrt(3)-1)/2/sqrt(3)*K1) '((mplus simp) ((mtimes simp) 3 ((mexpt simp) 3 ((rat simp) -3 4)) ((mexpt simp) $%pi ((rat simp) 3 2)) ((mexpt simp) ((%gamma simp) ((rat simp) 1 6)) -1) ((mexpt simp) ((%gamma simp) ((rat simp) 1 3)) -1)) ((mtimes simp) ((rat simp) 3 8) ((mexpt simp) 3 ((rat simp) -3 4)) ((mexpt simp) $%pi ((rat simp) -1 2)) ((%gamma simp) ((rat simp) 1 6)) ((%gamma simp) ((rat simp) 1 3))) ((mtimes simp) ((rat simp) -1 8) ((mexpt simp) 3 ((rat simp) -1 4)) ((mexpt simp) $%pi ((rat simp) -1 2)) ((%gamma simp) ((rat simp) 1 6)) ((%gamma simp) ((rat simp) 1 3))))) (t e))) (t (recur-apply #'makegamma1 e)))) ;; Autoload the nset functions that simplify. (add2lnc '$kron_delta $props) (defprop $kron_delta simp-kron-delta operators) (autof 'simp-kron-delta '|orthopoly|) (add2lnc '$unit_step $props) (defprop $unit_step simp-unit-step operators) (autof 'simp-unit-step '|orthopoly|) ($setup_autoload "orthopoly" '$assoc_legendre_p '$assoc_legendre_q '$chebyshev_t '$chebyshev_u '$gen_laguerre '$hermite '$intervalp '$jacobi_p '$laguerre '$legendre_p '$legendre_q '$orthopoly_recur '$orthopoly_weight '$pochhammer '$spherical_bessel_j '$spherical_bessel_y '$spherical_hankel1 '$spherical_hankel2 '$spherical_harmonic '$ultraspherical) --- NEW FILE: orthodisc.lisp --- (eval-when (compile load eval) ($load "orthopoly")) (defmvar $hyper_f_tol 1.0e-18) (defmvar max-fpprec 1000) (defun float-or-bfloatp (s) (or (floatp s) ($bfloatp s))) (defun $hyper_f (a b x) (if (not ($listp a)) (merror "The first argument to hyper_f must be a list, instead found ~:M" a)) (if (not ($listp a)) (merror "The second argument to hyper_f must be a list, instead found ~:M" b)) (setq a (cdr a)) (setq b (cdr b)) (setq b (cons 1 b)) (let ((n -1) (use-float nil)) (dolist (ai a) (if (and (integerp ai) (< ai 0)) (setq n (max n (- ai))))) (setq use-float (and (or (some 'float-or-bfloatp a) (some 'float-or-bfloatp b) (float-or-bfloatp x)) (every '$numberp a) (every '$numberp b) ($numberp x))) (cond ((= n -1) (setq a (mapcar #'(lambda (s) `(($pochhammer simp) ,s ,$genindex)) a)) (setq b (mapcar #'(lambda (s) `(($pochhammer simp) ,s ,$genindex)) b)) `((%sum) ((mtimes) ,@a ((mexpt) ,x ,$genindex) ((mexpt) ((mtimes) ,@b) -1)) ,$genindex 0 $inf)) (use-float (hyper-f-float a b x)) (t (let ((sum 1) (cf 1)) (dotimes (i n sum) (setq cf (div (mul x cf (apply 'mul (mapcar #'(lambda (s) (add i s)) a))) (apply 'mul (mapcar #'(lambda (s) (add i s)) b)))) (setq sum (add sum cf)))))))) (defun hyper-f-float (aa bb xx &optional (tol $hyper_f_tol)) (let ((n -1) (continue t) (cf) (sum) (err) (save-fpprec $fpprec) ($float2bf t) (new-digits) (bf-digits) (a) (b) (x) (i) (m)) (dolist (ai aa) (if (and (integerp ai) (< ai 0)) (setq n (max n (- ai))))) (setq bf-digits (* 2 (max 19 (- (ceiling (/ (log tol) (log 10.0))))))) (fpprec1 nil bf-digits) (setq $fpprec bf-digits) (setq m (* 2 (+ (length aa) (length bb) 1))) (while continue (if (> $fpprec max-fpprec) (merror "Unable to obtain requested precision")) (print `(fpprec = ,$fpprec)) (setq a (mapcar '$bfloat aa)) (setq b (mapcar '$bfloat bb)) (setq x ($bfloat xx)) (setq sum bigfloatone) (setq err bigfloatzero) (setq i (- n 1)) (while (>= i 0) (setq cf (div (mul x (apply 'mul (mapcar #'(lambda (s) (add i s)) a))) (apply 'mul (mapcar #'(lambda (s) (add i s)) b)))) (setq err (add bigfloatone (mul m (simplify `((mabs) ,cf)) err))) (setq cf (mul sum cf)) (setq sum (add bigfloatone cf)) (setq err (add err (mul (+ 2 m) (simplify `((mabs) ,cf))))) (decf i)) (setq new-digits (- (/ (log (/ ($float tol) ($float err))) (log 10.0)))) (setq err (mul err (bfloat-epsilon))) (setq new-digits (+ 10 (max $fpprec (ceiling new-digits)))) (fpprec1 nil new-digits) (setq $fpprec new-digits) (setq continue (mgrp err tol))) (displa sum) (cond ((or (some '$bfloatp aa) (some '$bfloatp bb) ($bfloatp xx)) (fpprec1 nil bf-digits) (setq $fpprec bf-digits) (setq sum ($bfloat sum))) (t (setq sum ($float sum)))) (fpprec1 nil save-fpprec) (setq $fpprec save-fpprec) sum)) (defun bfloat-epsilon () (power ($bfloat 2) (- fpprec))) (defun $krawtchouk (n a b x) (cond ((use-hypergeo n x) (let ((f) (e)) (multiple-value-setq (f e) ($hypergeo21 (mul -1 n) (mul -1 x) (mul -1 a) (div 1 b) n)) (orthopoly-return-handler 1 f e))) (t `(($krawtchouk simp) ,n ,a ,x)))) (defun $meixner_m (n g mu x) (cond ((use-hypergeo n x) (let ((d) (f) (e)) (setq d ($pochhammer g n)) (multiple-value-setq (f e) ($hypergeo21 (mul -1 n) (mul -1 x) g (sub 1 (div 1 mu)) n)) (orthopoly-return-handler d f e))) (t `(($meixner_m simp) ,n ,g ,mu ,x)))) (defun $charlier_c (n mu x) (cond ((use-hypergeo n x) (let ((f) (e)) (multiple-value-setq (f e) ($hyper_f `((mlist) ,(mul -1 n) ,(mul -1 x)) `((mlist)) (div -1 mu))) (orthopoly-return-handler 1 f e))) (t `(($charlier_c simp) ,n ,mu ,x)))) (defun $hahn_h (n m a b x) (cond ((use-hypergeo n x) (let ((f) (d) (e)) (setq d (mul ($pochhammer (sub m n) n) ($pochhammer (add 1 b) n))) (if (oddp n) (setq d (mul -1 d))) (setq d (div d ($pochhammer 1 n))) (multiple-value-setq (f e) ($hyper_f `((mlist) ,(mul -1 n) ,(mul -1 x) ,(add a b n 1)) `((mlist) ,(add 1 b) ,(sub 1 m)) 1)) (orthopoly-return-handler d f e))) (t `(($hahn_h simp) ,n ,m ,a ,b ,x)))) (defun $hahn_q (n m a b x) (cond ((use-hypergeo n x) (let ((f) (e)) (multiple-value-setq (f e) ($hyper_f `((mlist) ,(mul -1 n) ,(mul -1 x) ,(add a b n 1)) `((mlist) ,(add 1 a) ,(mul -1 m)) 1)) (orthopoly-return-handler 1 f e))) (t `(($hahn_q simp) ,n ,m ,a ,b ,x)))) (defun $discrete_chebyshev_h (n m x) ($hahn_h n m 0 0 x)) (defun $discrete_jacobi_p (n m a b x) ($hahn_h n m a b (div (mul (add x 1) (sub m 1)) 2))) (defun $discrete_legendre_p (n m x) ($discrete_jacobi_p n m 0 0 x)) (defun $orthodisc_weight (fn arg) (if (not ($listp arg)) (merror "The second argument to orthodisc_weight must be a list")) (if (not (or ($symbolp (car (last arg))) ($subvarp (car (last arg))))) (merror "The last element of the second argument to orthopoly_disc must be a symbol or a subscripted symbol, instead found ~:M" (car (last arg)))) (if (not (every #'(lambda (s) ($freeof (car (last arg)) s)) (butlast (cdr arg)))) (merror "Only the last element of ~:M may depend on the summation variable ~:M" arg (car (last arg)))) (cond ((eq fn '$discrete_chebyshev_h) (check-arg-length fn 3 (- (length arg) 1)) (let ((m (nth 2 arg))) `((mlist) 1 ((lambda) ((mlist) k) k) 0 ,($entier m)))) (t (merror "The weight for ~:M isn't known to maxima" fn)))) --- NEW FILE: orthopoly-init.lisp --- ;;---------------------------------------- ;; Optional code starts here ;; The Maxima fixfloat function is buggy. Here is a quick fix. At one ;; time, test_orthopoly needed this fix -- I think it's no longer needed. #+cl (defun fixfloat (x) (cond ((floatp x) (setq x (rationalize x)) (cons (numerator x) (denominator x))) (t x))) ;; fix for apply(rat(f),[x]) bug. (defun error-size (exp) (setq exp (ratdisrep exp)) (if (atom exp) 0 (do ((l (cdr exp) (cdr l)) (n 1 (incf n (+ 1 (error-size (car l)))))) ((or (null l) (> n $error_size)) n) (declare (fixnum n))))) ;; If file fn (with _no_ file extension) is out of date, compile it. We ;; assume that the extension of the source file is "lisp". (defun $make (fn) (setq fn (string-left-trim "&" (string fn))) (let ((fl) (fb) (bin-ext #+gcl "o" #+cmu (c::backend-fasl-file-type c::*target-backend*) #+clisp "fas" #+allegro "fasl" #-(or gcl cmu clisp allegro) "")) (setq fl (concatenate 'string (coerce fn 'string) ".lisp")) (setq fl ($file_search fl)) (setq fb (concatenate 'string (coerce fn 'string) "." bin-ext)) (setq fb ($file_search fb)) (if (or (and fl (not fb)) (and fl fb (> (file-write-date fl) (file-write-date fb)))) ($compile_file fl)))) ;; If you place orthopoly.lisp in a directory that Maxima can't find, ;; you'll need to append its directory to $file_search_lisp. To do this, ;; replace the path "/home/barton/orthopoly-0.93" with the correct path ;; for your machine. ;; The make function recompiles orthopoly automatically on startup if ;; the binary file is older than orthopoly.lisp. Unless you ;; anticipate making changes to orthopoly, you may comment out ;; ($make "orthopoly") (eval-when (load compile eval) (setq $file_search_lisp ($cons "/home/barton/orthopoly-0.93/###.{x86f,lisp}" $file_search_lisp)) ($make "orthopoly")) ;; Load file f with verbose and print bound to true. Useful for loading ;; code with syntax errors. (defun $xload (f) (load ($file_search f) :verbose 't :print 't)) (defprop $entier tex-matchfix tex) (defprop $entier (("\\lfloor ") " \\rfloor") texsym) ;;----------------------------------------------------------- ;; Essential code starts here ;; This fixes a bug in $bessel_j. (defun $bessel_j (arg order) (if (and (numberp order) (numberp ($realpart arg)) (numberp ($imagpart arg))) ;($bessel (complex ($realpart arg) ($imagpart arg)) order) <== bug ($bessel arg order) (subfunmakes '$bessel_j (ncons order) (ncons arg)))) ;; This version of setup_autoload uses file_search. (defun $setup_autoload (filename &rest functions) (let ((file ($file_search filename))) (dolist (func functions) (nonsymchk func '$setup_autoload) (putprop (setq func (dollarify-name func)) file 'autoload) (add2lnc func $props))) '$done) ;; This version of generic-autoload in Maxima 5.9.0 doesn't recognize ;; "fas", "fasl" and "x86f" as valid extensions for compiled Lisp ;; code. #+cl (defun generic-autoload (file &aux type) (setq file (pathname (cdr file))) (setq type (pathname-type file)) (let ((bin-ext #+gcl "o" #+cmu (c::backend-fasl-file-type c::*target-backend*) #+clisp "fas" #+allegro "fasl" #-(or gcl cmu clisp allegro) "")) (if (member type (list bin-ext "lisp" "lsp") :test 'equalp) (load file :verbose 't) ($batchload file)))) ;; Extended version of makegamma1 that converts pochhammer symbols ;; to quotients of gamma functions. (defun makegamma1 (e) (cond ((atom e) e) ((eq (caar e) 'mfactorial) (list '(%gamma) (list '(mplus) 1 (makegamma1 (cadr e))))) ;;------start of new code---------------------------------- ;; Do pochhammer(x,n) ==> gamma(x+n)/gamma(x). ((eq (caar e) '$pochhammer) (let ((x (makegamma1 (nth 1 e))) (n (makegamma1 (nth 2 e)))) `((mtimes) ((%gamma) ((mplus) ,x ,n)) ((mexpt) ((%gamma) ,x) -1)))) ((eq (caar e) '%genfact) (let ((x (makegamma1 (nth 1 e))) (y (makegamma1 (nth 2 e))) (z (makegamma1 (nth 3 e)))) (setq x (add (div x z) 1)) (setq y (simplify `(($entier) ,y))) (setq z (power z y)) `((mtimes) ,z ((%gamma) ,x) ((mexpt) ((%gamma) ((mplus) ,x ((mtimes) -1 ,y))) -1)))) ;;-----that's all folks ... end of new code------------------ ((eq (caar e) '%elliptic_kc) ;; Complete elliptic integral of the first kind (cond ((alike1 (cadr e) '((rat simp) 1 2)) ;; K(1/2) = gamma(1/4)/4/sqrt(pi) '((mtimes simp) ((rat simp) 1 4) ((mexpt simp) $%pi ((rat simp) -1 2)) ((mexpt simp) ((%gamma simp) ((rat simp) 1 4)) 2))) ((or (alike1 (cadr e) '((mtimes simp) ((rat simp) 1 4) ((mplus simp) 2 ((mexpt simp) 3 ((rat simp) 1 2))))) (alike1 (cadr e) '((mplus simp) ((rat simp) 1 2) ((mtimes simp) ((rat simp) 1 4) ((mexpt simp) 3 ((rat simp) 1 2))))) (alike1 (cadr e) ;; 1/(8-4*sqrt(3)) '((mexpt simp) ((mplus simp) 8 ((mtimes simp) -4 ((mexpt simp) 3 ((rat simp) 1 2)))) -1))) ;; K((2+sqrt(3)/4)) '((mtimes simp) ((rat simp) 1 4) ((mexpt simp) 3 ((rat simp) 1 4)) ((mexpt simp) $%pi ((rat simp) -1 2)) ((%gamma simp) ((rat simp) 1 6)) ((%gamma simp) ((rat simp) 1 3)))) ((or (alike1 (cadr e) ;; (2-sqrt(3))/4 '((mtimes simp) ((rat simp) 1 4) ((mplus simp) 2 ((mtimes simp) -1 ((mexpt simp) 3 ((rat simp) 1 2)))))) (alike1 (cadr e) ;; 1/2-sqrt(3)/4 '((mplus simp) ((rat simp) 1 2) ((mtimes simp) ((rat simp) -1 4) ((mexpt simp) 3 ((rat simp) 1 2))))) (alike (cadr e) ;; 1/(4*sqrt(3)+8) '((mexpt simp) ((mplus simp) 8 ((mtimes simp) 4 ((mexpt simp) 3 ((rat simp) 1 2)))) -1))) ;; K((2-sqrt(3))/4) '((mtimes simp) ((rat simp) 1 4) ((mexpt simp) 3 ((rat simp) -1 4)) ((mexpt simp) $%pi ((rat simp) -1 2)) ((%gamma simp) ((rat simp) 1 6)) ((%gamma simp) ((rat simp) 1 3)))) ((or ;; (3-2*sqrt(2))/(3+2*sqrt(2)) (alike1 (cadr e) '((mtimes simp) ((mplus simp) 3 ((mtimes simp) -2 ((mexpt simp) 2 ((rat simp) 1 2)))) ((mexpt simp) ((mplus simp) 3 ((mtimes simp) 2 ((mexpt simp) 2 ((rat simp) 1 2)))) -1))) ;; 17 - 12*sqrt(2) (alike1 (cadr e) '((mplus simp) 17 ((mtimes simp) -12 ((mexpt simp) 2 ((rat simp) 1 2))))) ;; (2*SQRT(2) - 3)/(2*SQRT(2) + 3) (alike1 (cadr e) '((mtimes simp) -1 ((mplus simp) -3 ((mtimes simp) 2 ((mexpt simp) 2 ((rat simp) 1 2)))) ((mexpt simp) ((mplus simp) 3 ((mtimes simp) 2 ((mexpt simp) 2 ((rat simp) 1 2)))) -1)))) '((mtimes simp) ((rat simp) 1 8) ((mexpt simp) 2 ((rat simp) -1 2)) ((mplus simp) 1 ((mexpt simp) 2 ((rat simp) 1 2))) ((mexpt simp) $%pi ((rat simp) -1 2)) ((mexpt simp) ((%gamma simp) ((rat simp) 1 4)) 2))) (t ;; Give up e))) ((eq (caar e) '%elliptic_ec) ;; Complete elliptic integral of the second kind (cond ((alike1 (cadr e) '((rat simp) 1 2)) ;; 2*E(1/2) - K(1/2) = 2*%pi^(3/2)*gamma(1/4)^(-2) '((mplus simp) ((mtimes simp) ((mexpt simp) $%pi ((rat simp) 3 2)) ((mexpt simp) ((%gamma simp irreducible) ((rat simp) 1 4)) -2)) ((mtimes simp) ((rat simp) 1 8) ((mexpt simp) $%pi ((rat simp) -1 2)) ((mexpt simp) ((%gamma simp) ((rat simp) 1 4)) 2)))) ((or (alike1 (cadr e) '((mtimes simp) ((rat simp) 1 4) ((mplus simp) 2 ((mtimes simp) -1 ((mexpt simp) 3 ((rat simp) 1 2)))))) (alike1 (cadr e) '((mplus simp) ((rat simp) 1 2) ((mtimes simp) ((rat simp) -1 4) ((mexpt simp) 3 ((rat simp) 1 2)))))) ;; E((2-sqrt(3))/4) ;; ;; %pi/4/sqrt(3) = K*(E-(sqrt(3)+1)/2/sqrt(3)*K) '((mplus simp) ((mtimes simp) ((mexpt simp) 3 ((rat simp) -1 4)) ((mexpt simp) $%pi ((rat simp) 3 2)) ((mexpt simp) ((%gamma simp) ((rat simp) 1 6)) -1) ((mexpt simp) ((%gamma simp) ((rat simp) 1 3)) -1)) ((mtimes simp) ((rat simp) 1 8) ((mexpt simp) 3 ((rat simp) -3 4)) ((mexpt simp) $%pi ((rat simp) -1 2)) ((%gamma simp) ((rat simp) 1 6)) ((%gamma simp) ((rat simp) 1 3))) ((mtimes simp) ((rat simp) 1 8) ((mexpt simp) 3 ((rat simp) -1 4)) ((mexpt simp) $%pi ((rat simp) -1 2)) ((%gamma simp) ((rat simp) 1 6)) ((%gamma simp) ((rat simp) 1 3))))) ((or (alike1 (cadr e) '((mtimes simp) ((rat simp) 1 4) ((mplus simp) 2 ((mexpt simp) 3 ((rat simp) 1 2))))) (alike1 (cadr e) '((mplus simp) ((rat simp) 1 2) ((mtimes simp) ((rat simp) 1 4) ((mexpt simp) 3 ((rat simp) 1 2)))))) ;; E((2+sqrt(3))/4) ;; ;; %pi*sqrt(3)/4 = K1*(E1-(sqrt(3)-1)/2/sqrt(3)*K1) '((mplus simp) ((mtimes simp) 3 ((mexpt simp) 3 ((rat simp) -3 4)) ((mexpt simp) $%pi ((rat simp) 3 2)) ((mexpt simp) ((%gamma simp) ((rat simp) 1 6)) -1) ((mexpt simp) ((%gamma simp) ((rat simp) 1 3)) -1)) ((mtimes simp) ((rat simp) 3 8) ((mexpt simp) 3 ((rat simp) -3 4)) ((mexpt simp) $%pi ((rat simp) -1 2)) ((%gamma simp) ((rat simp) 1 6)) ((%gamma simp) ((rat simp) 1 3))) ((mtimes simp) ((rat simp) -1 8) ((mexpt simp) 3 ((rat simp) -1 4)) ((mexpt simp) $%pi ((rat simp) -1 2)) ((%gamma simp) ((rat simp) 1 6)) ((%gamma simp) ((rat simp) 1 3))))) (t e))) (t (recur-apply #'makegamma1 e)))) ;; kron_delta and unit_step aren't functions. If you try to autoload ;; them, you'll throw Maxima into an infinite loop! It's okay to add ;; them to props. (add2lnc '$kron_delta $props) (add2lnc '$unit_step $props) ($setup_autoload "orthopoly" '$assoc_legendre_p '$assoc_legendre_q '$chebyshev_t '$chebyshev_u '$gen_laguerre '$hermite '$intervalp '$jacobi_p '$laguerre '$legendre_p '$legendre_q '$orthopoly_recur '$orthopoly_weight '$pochhammer '$spherical_bessel_j '$spherical_bessel_y '$spherical_hankel1 '$spherical_hankel2 '$spherical_harmonic '$ultraspherical) --- NEW FILE: orthopoly.lisp --- ;; 8/8 | 5/5 ;; Copyright (C) 2000, 2001, 2003 Barton Willis ;; Maxima code for evaluating orthogonal polynomials listed in Chapter 22 ;; of Abramowitz and Stegun (A & S). ;; Author: Barton Willis, University of Nebraska at Kearney (aka UNK). ;; License: see README. The user of this code assumes all risk ;; for its use. It has no warranty. If you don't know the meaning ;; of "no warranty," don't use this code. :-) (in-package "MAXIMA") (mfunction-call $put '|$orthopoly| '|&0.94| '$version) ;; A while loop taken from nset.lisp. Someday it should be placed ;; in a file with other Maxima macros and deleted from here and from nset. (defmacro while (cond &rest body) `(do () [...1452 lines suppressed...] (x (nth 3 arg))) (simplify `((mlist) ((mtimes) ((mexpt) ,x ,a) ((mexpt) $%e ((mtimes) -1 ,x))) 0 $inf)))) ((eq fn '$hermite) (check-arg-length fn 2 (- (length arg) 1)) (let ((x (nth 2 arg))) (simplify `((mlist) ((mexpt) $%e ((mtimes) -1 ((mexpt) ,x 2))) ((mtimes ) -1 $inf) $inf)))) (t (merror "a weight for ~:m isn't known to maxima" fn)))) --- NEW FILE: orthopoly_doc.html --- <HTML> <HEAD> <!-- This HTML file has been created by texi2html 1.52 from orthopoly_doc.texi on 30 August 2003 --> <TITLE>Untitled Document</TITLE> </HEAD> <BODY> <H1>Untitled Document</H1> <P> <P><HR><P> <H1><A NAME="SEC1" HREF="orthopoly_doc_toc.html#TOC1">Orthogonal Polynomials</A></H1> [...1034 lines suppressed...] 1 for <EM>x > 0</EM>. If you want a unit step function that takes on the value 1/2 at zero, use <EM> (1 + signum(x))/2</EM>. </DL> </P> <P> <DL> <DT><U>Function:</U> <B>ultraspherical</B> <I>(n,a,x)</I> <DD><A NAME="IDX24"></A> The ultraspherical polynomial (also known the Gegenbauer polynomial). Reference A & S 22.5.46, page 779. </DL> </P> <P><HR><P> This document was generated on 30 August 2003 using the <A HREF="http://wwwcn.cern.ch/dci/texi2html/">texi2html</A> translator version 1.51a.</P> </BODY> </HTML> --- NEW FILE: orthopoly_doc.texi --- @node Orthogonal Polynomials, Special Functions, Command Line, Top @chapter Orthogonal Polynomials @menu * Introduction to Orthogonal Polynomials:: * Definitions for Orthogonal Polynomials:: @end menu @section Introduction to Orthogonal Polynomials Maxima includes support for symbolic and numerical evaluation of most of the orthogonal polynomials listed in Chapter 22 of Abramowitz and Stegun; these functions include the Chebyshev, Laguerre, Hermite, Jacobi, Legendre, and ultraspherical (Gegenbauer) polynomials. Additionally, Maxima includes support for the spherical Bessel, spherical Hankel, and spherical harmonic functions. For the most part, orthopoly follows the conventions of Abramowitz and Stegun (A & S) @emph{Handbook of Mathematical Functions} (10th printing, December 1972); additionally, we use Gradshteyn and Ryzhik (G & R), @emph{Table of Integrals, Series, and Products} (1980 corrected and enlarged edition), and Eugen Merzbacher @emph{Quantum Mechanics} (2 ed, 1970). @subsubsection Installation Download the archive orthopoly_x.tar.gz, where x is the release identifier, from http://www.unk.edu/acad/math/people/willisb. Under Linux, unpack it using gzip -d orthopoly_x.tar.gz tar -xvf orthopoly_x.tar This will create a directory @emph{orthopoly_x} (again x is the release identifier) that contains the source file @emph{orthopoly.lisp}, user documentation in html and texi formats, a sample maxima initialization file @emph{orthopoly-init.lisp}, a README file, a testing routine @emph{test_orthopoly.mac}, and two demonstration files. Start Maxima and compile orthopoly. To do this, use the command (c1) compile_file("orthopoly.lisp"); If Maxima is unable to find orthopoly, use the full pathname. Compiling will create a file orthopoly.x, where the file extension "x" depends on which version of Lisp your Maxima uses; under GCL, the extension is "o", and under cmucl, most likely it is "x86f". Copy the source file "orthopoly.lisp" and its compiled version to a directory that Maxima can find. Since version 0.93 is a preliminary version, you may want to keep orthopoly in the orthopoly_0.93 directory and add this directory to file_search_lisp. The file @emph{orthopoly_init.lisp} shows how to append a directory to file_search_lisp. For a permanent installation, copy the files to Maxima's /share/specfunctions/ directory. Add the contents of "orthopoly-init.lisp" to your maxima-init.lisp file. Some code in orthopoly-init.lisp is marked as optional; you needn't add it to your maxima-init.lisp file. The orthopoly-init.lisp file contains autoload statements for user-level functions in orthopoly and it contains fixes for bugs in the functions $setup_autoload and generic-autoload. Additionally, "orthopoly-init.lisp" contains an extended version of the function @emph{makegamma}. To test orthopoly, load("test-orthopoly.mac"). You may need to give the full pathname for the file. Report errors. @end enumerate @subsubsection Getting Started with Orthopoly Assuming you've placed appropriate autoload statements in your @emph{maxima-init.lisp} file, you needn't manually load @b{orthopoly} to start using it; otherwise, load it with the command @example (C1) load("orthopoly")$ @end example If Maxima isn't able to find the file, use a full pathname. Without autoloading user-level orthopoly functions, you'll need to be careful to load orthopoly @emph{before} using any user-level functions from orthopoly; otherwise, you'll need to use upper-case function names. To find the third order Legendre polynomial, use the command @example (C2) legendre_p(3,x); 3 2 5 (1 - x) 15 (1 - x) (D2) - ---------- + ----------- - 6 (1 - x) + 1 2 2 @end example To express this as a sum of powers of @math{x}, apply ratsimp or rat or to the result @example (C3) [ratsimp(%),rat(%)]; 3 3 5 x - 3 x 5 x - 3 x (D3)/R/ [----------, ----------] 2 2 @end example Alternatively, make the second argument to @math{legendre_p} (its ``main'' variable) a CRE expression @example (C4) legendre_p(3,rat(x)); 3 5 x - 3 x (D4)/R/ ---------- 2 @end example For floating point evaluation, orthopoly uses a running error analysis to @emph{estimate} an upper bound for the error. An example @example (C1) jacobi_p(150,2,3,0.2); (D1) interval(- .0620170379367145, 2.04311850697459e-11) @end example Intervals have the form @math{interval(c, r)}, where @math{c} is the center and @math{r} is the radius of the interval. Since Maxima does not support arithmetic on intervals, in some situations, such as graphics, you want to suppress the error and output only the center of the interval. To do this, set the option variable @emph{orthopoly_returns_intervals} to false @example (C2) orthopoly_returns_intervals : false; (D2) FALSE (C3) jacobi_p(150,2,3,0.2); (D3) - .0620170379367145 @end example Refer to the section @emph{ Floating point Evaluation} for more information. Most functions in orthopoly have a @emph{gradef} property; thus @example (C1) diff(hermite(n,x),x); (D1) 2 n H (x) n - 1 (C2) diff(gen_laguerre(n,a,x),x); (a) (a) n L (x) - (n + a) L (x) UNIT_STEP(n) n n - 1 (D2) ------------------------------------------ x @end example The unit step function in the second example prevents an error that would otherwise arise by evaluating with @math{n = 0}. @example (C3) ev(%,n=0); (D3) 0 @end example The gradef property only applies to the ``main'' variable; derivatives with respect other arguments usually result in an error message; for example @example (C1) diff(hermite(n,x),x); (D1) 2 n H (x) n - 1 (C2) diff(hermite(n,x),n); Maxima doesn't know the derivative of hermite with respect the first argument -- an error. Quitting. To debug this try DEBUGMODE(TRUE);) @end example Generally, functions in orthopoly map over lists and matrices. For the mapping to fully evaluate, the option variables @emph{doallmxops} and @emph{listarith} both must assume their default values (true). To illustrate the mapping over matrices, consider @example (C1) hermite(2,x); 2 (D1) - 2 (1 - 2 x ) (C2) m : matrix([0,x],[y,0]); [ 0 x ] (D2) [ ] [ y 0 ] (C3) hermite(2,m); [ 2 ] [ - 2 - 2 (1 - 2 x ) ] (D3) [ ] [ 2 ] [ - 2 (1 - 2 y ) - 2 ] @end example In the second example, understand that @emph{i,j} element of the value is @emph{hermite(2,m[i,j])}; this is not the same as computing @emph{-2 + 4 m . m} @example (C4) -2 * matrix([1,0],[0,1]) + 4 * m.m; [ 4 x y - 2 0 ] (D4) [ ] [ 0 4 x y - 2 ] @end example If you evaluate a function at a point outside its domain, generally orthopoly will return the function unevaluated; an example @example (C1) legendre_p(2/3,x); (D1) P (x) 2/3 @end example Orthopoly supports translation into TeX; it also does two-dimensional output on a terminal @example (C1) spherical_harmonic(l,m,theta,phi); m (D1) Y (THETA, PHI) l (C2) tex(%); $$Y_{l}^{m}\left(\vartheta,\varphi\right)$$ (D2) FALSE (C3) jacobi_p(n,a,a-b,x/2); (a, a - b) x (D3) P (-) n 2 (C4) tex(%); $$P_{n}^{\left(a,a-b\right)}\left({{x}\over{2}}\right)$$ (D4) FALSE @end example @subsubsection Caveats When an expression involves several orthogonal polynomials with @emph{symbolic} orders, it's possible that the expression actually vanishes, yet Maxima is unable to simplify it to zero. If you divide by such a quantity, you'll be in trouble. For example, the following expression vanishes for integers @emph{n>1}, yet Maxima is unable to simplify it to zero. @example (D5) (2 n - 1) P (x) x - n P (x) + (1 - n) P (x) n - 1 n n - 2 @end example For a specific @emph{n}, we can reduce the expression to zero @example (C6) ev(%,n=10,ratsimp); (D6) 0 @end example @emph{Be careful.} Generally, the polynomial form of an orthogonal polynomial is ill-suited for floating point evaluation. Here's an example @example (C1) p : jacobi_p(150,2,3,x)$ (C2) subst(0.2,x,p); (D2) - 9.470489909945016e+60 @end example The true value is about -0.06; this calculation suffers from extreme subtractive cancellation error. Expanding the polynomial and then evaluating, gives a better result @example (C3) p : expand(p)$ (C4) subst(0.2,x,p); (D4) - 0.062128442689779795 @end example This isn't a general rule; expanding the polynomial does not always result in an expression that is better suited for numerical evaluation. By far, the best way to do numerical evaluation is to make one or more of the function arguments floating point numbers. By doing that, specialized floating point algorithms are used for evaluation. Maxima's float function is somewhat indiscriminant; if you apply float to an an expression involving an orthogonal polynomial with a symbolic degree or order parameter, these parameters may be converted into floats; after that, the expression will not evaluate fully. Consider @example (C1) assoc_legendre_p(n,1,x); 1 (D1) P (x) n (C2) float(%); 1.0 (D2) P (x) n (C3) ev(%,n=2,x=0.9); 1.0 (D3) P (0.9) 2 (C4) @end example The expression in (D3) will not evaluate to a float; orthopoly doesn't recognize floating point values where it requires an integer. Similarly, numerical evaluation of the pochhammer function for orders that exceed @math{pochhammer_max_index} can be troublesome; consider @example (C1) x : pochhammer(1,10), pochhammer_max_index : 5; (D1) (1) 10 @end example Applying @math{float} doesn't evaluate @math{x} to a float @example (C2) float(x); (D2) (1.0) 10.0 @end example To evaluate @math{x} to a float. you'll need to bind @math{pochhammer_max_index} to 11 or greater and apply float to @math{x} @example (C3) float(x), pochhammer_max_index : 11; (D3) 3628800.0 (C4) @end example The default value of @math{pochhammer_max_index} is 100; to change its value, first load orthopoly. Finally, be aware that reference books vary on the definitions of the orthogonal polynomials; we've generally used the conventions of conventions of Abramowitz and Stegun (A & S) @emph{Handbook of Mathematical Functions} (10th printing, December 1972). Before you suspect a bug in orthopoly, check some special cases to determine if your definitions match those used by orthonormal. Definitions often differ by a normalization; occasionally, authors use ``shifted'' versions of the functions that makes the family orthogonal on an interval other than @emph{(-1,1)}. To define, for example, a legendre polynomial that is orthogonal on @emph{(0,1)}, define @example (C1) shifted_legendre_p(n,x) := legendre_p(n,2*x-1)$ (C2) shifted_legendre_p(2,rat(x)); 2 (D2)/R/ 6 x - 6 x + 1 (C3) legendre_p(2,rat(x)); 2 3 x - 1 (D3)/R/ -------- 2 @end example @subsubsection Floating point Evaluation Most functions in orthopoly use a running error analysis to @emph{estimate} the error in floating point evaluation; the exceptions are the spherical Bessel functions and the associated Legendre polynomials of the second kind. For numerical evaluation, the spherical Bessel functions call slatec functions. No specialized method is used for numerical evaluation of the associated Legendre polynomials of the second kind. The running error analysis ignores errors that are second or higher order in the machine epsilon (also known as unit roundoff). It also ignores a few other errors. @emph{It's possible (although unlikely) that the actual error exceeds the estimate.} Intervals have the form @emph{interval(c,r)}, where @emph{c} is the @emph{center} of the interval and @emph{r} is its @emph{radius}. The center of an interval can be a complex number, but it's a bug if the radius isn't a positive real number. Here is an an example @example (C1) fpprec : 50$ (C2) y0 : jacobi_p(150,2,3,0.2); (D2) interval(- 0.06201703793671447, 1.0217141721600247e-11) (C3) y1 : bfloat(jacobi_p(150,2,3,1/5)); (D3) - 6.201703793671628761993584658478664938137943464587B-2 @end example Let's test that the actual error is smaller than the error estimate @example (C4) is(abs(part(y0,1) - y1) < part(y0,2)); (D4) TRUE @end example Indeed, for this example the error estimate is an upper bound for the true error. Maxima does not support arithmetic on intervals @example (C1) legendre_p(7,0.1) + legendre_p(8,0.1); (D1) interval(0.18032072148437508, 2.1508371989568562e-15) + interval(- 0.19949294375000004, 2.301870533171531e-15) @end example A user could define arithmetic operators that do interval math. To define interval addition, we can define @example (C9) infix("@+")$ (C10) "@+"(x,y) := interval(part(x,1) + part(y,1),part(x,2) + part(y,2))$ (C11) legendre_p(7,0.1) @+ legendre_p(8,0.1); (D11) INTERVAL(- 0.019172222265624955, 4.452707732128387e-15) @end example The special floating point routines get called when the arguments are complex. For example @example (C1) legendre_p(10,2 + 3.0*%i); (D1) interval(- 3.8763788250000003e+7 %I - 6.0787748e+7, 8.322322816893654e-7) @end example Let's compare this to the true value @example (C2) float(expand(legendre_p(10,2+3*%i))); (D2) - 3.8763788250000003e+7 %I - 6.0787748e+7 @end example Additionally, when the arguments are big floats, the special floating point routines get called; however, the big floats are converted into double floats and the final result is a double @example (C3) ultraspherical(150,0.5b0,0.9b0); (D3) interval(- 0.0430094812572654, 2.2727803799647722e-14) @end example @subsubsection Graphics and Orthopoly To plot expressions that involve the orthogonal polynomials, you must do two things: @enumerate @item set the option variable @emph{orthopoly_returns_intervals} to false, @item fully quote any calls to orthopoly functions. @end enumerate If function calls aren't quoted, Maxima evaluates them to polynomials before plotting; consequently, the specialized floating point code doesn't get called. Here is an example of how to plot an expression that involves a legendre polynomial. @example (C1) plot2d('(legendre_p(5,x)),[x,0,1]), orthopoly_returns_intervals : false; @end example The @emph{entire} expression @math{legendre_p(5,x)} is quoted; this is different than just quoting the function name using @math{'legendre_p(5,x)}. @subsubsection Miscellaneous Functions The orthopoly package defines the Kronecker delta function, the pochhammer symbol, and a unit step function. Orthopoly uses the Kronecker delta function and the unit step function in gradef statements. To convert pochhammer symbols into quotients of gamma functions, use @math{makegamma} @example (C1) makegamma(pochhammer(x,n)); GAMMA(x + n) (D1) ------------ GAMMA(x) (C2) makegamma(pochhammer(1/2,1/2)); 1 (D2) --------- SQRT(%PI) @end example Derivatives of the pochhammer symbol are given in terms of the @math{psi} function @example (C3) diff(pochhammer(x,n),x); (D3) (x) (PSI (x + n) - PSI (x)) n 0 0 (C4) diff(pochhammer(x,n),n); (D4) (x) PSI (x + n) n 0 @end example You need to be careful with the expression in (D3); the difference of the @math{psi} functions has poles when @math{x = -1,-2,..,-n}. These poles cancel with factors in @math{(x)_n} making the derivative a degree @math{n-1} polynomial when @math{n} is a positive integer. The pochhammer symbol is defined for negative orders through its representation as a quotient of gamma functions. Consider @example (C1) q : makegamma(pochhammer(x,n)); GAMMA(x + n) (D1) ------------ GAMMA(x) (C2) sublis([x=11/3,n=-6],q); 729 (D2) - ---- 2240 @end example Alternatively, we can get this result directly @example (C3) pochhammer(11/3,-6); 729 (D3) - ---- 2240 (C4) @end example The unit step function is @emph{left-continuous}; thus @example (C1) [unit_step(-1/10),unit_step(0),unit_step(1/10)]; (D1) [0, 0, 1] @end example If you need a unit step function that is neither left or right continuous at zero, define your own using signum; for example, @example (C2) xunit_step(x) := (1 + signum(x))/2$ (C3) [xunit_step(-1/10),xunit_step(0),xunit_step(1/10)]; 1 (D3) [0, -, 1] 2 @end example Do not re-define Maxima's unit step function; some code in orthopoly requires that the unit step function is left-continuous. @subsubsection Algorithms Generally, orthopoly does symbolic evaluation by using a hypergeometic representation of the various orthogonal polynomials. The hypergeometic functions are evaluated using the (undocumented) functions @math{hypergeo11} and @math{hypergeo21}. The exceptions are the half-integer Bessel functions and the associated Legendre function of the second kind. The Bessel functions are evaluated using an explicit representation, while the associated Legendre function of the second kind is evaluated using recursion. For floating point evaluation, we again convert most functions into a hypergeometic form; we evaluate the hypergeometic functions using forward recursion. Again, the exceptions are the half-integer Bessel functions and the associated Legendre function of the second kind. Numerically, the half-integer Bessel functions are evaluated using the slatec code, and the associated Legendre functions of the second kind is numerically evaluated using the same algorithm as its symbolic evaluation uses. @subsubsection Author, License, and History Barton Willis of the University of Nebraska at Kearney (aka UNK) wrote and maintains the orthopoly package and its documentation. The package is released under the GNU General Public License (GPL). The first two releases, specfun version 110 and specfun version 111, were released in April 2001 and May 2002. These versions were included in Maxima versions starting with 5.5. An preliminary third release, renamed orthopoly, was announced in May 2003. The third version adds TeX and display support, improved numerical floating point accuracy, and new user documentation. @node Definitions for Orthogonal Polynomials, , Orthogonal Polynomials, Orthogonal Polynomials @section Definitions for Orthogonal Polynomials @defun assoc_legendre_p (n, m, x) The associated Legendre function of the first kind. Reference: A & S equation 22.5.37, page 779, A & S equation 8.6.6 (second equation), page 334, and A & S equation 8.2.5, page 333. @end defun @defun assoc_legendre_q (n, m, x) The associated Legendre function of the second kind. Reference: A & S 8.5.3 and 8.1.8. @end defun @defun chebyshev_t (n, x) The Chebyshev function of the first kind. Reference: A & S 22.5.47, page 779. @end defun @defun chebyshev_u (n, x) The Chebyshev function of the second kind. Reference: A & S, 22.5.48, page 779. @end defun @defun gen_laguerre (n, a, x) The generalized Laguerre polynomial. Reference: A & S 22.5.54, page 780. @end defun @defun hermite (n,x) The Hermite polynomial. Reference: See A&S 22.5.55, page 780. @end defun @defun intervalp(e) Return true if the input is an interval and return false if it isn't. @end defun @defun jacobi_p (n, a, b, x) The Jacobi polynomial. Reference: A & S. 22.5.42, page 779. The Jacobi polynomials are actually defined for all @emph{ a } and @emph{ b }; however, the Jacobi polynomial weight @emph{ (1-x)^a(1+x)^b} isn't integrable for @emph{ a <= -1} or @emph{ b <= -1}. ) @end defun @defun kron_delta (i, j) The Kronecker delta function; @emph{kron_delta(i,j)... [truncated message content] |