recent changes to orthopoly and kron_delta, was: Two Patches for Maxima (on SF)
Subject: recent changes to orthopoly and kron_delta, was: Two Patches for Maxima (on SF)
From: Robert Dodier
Date: Wed, 22 Aug 2012 18:50:56 +0000 (UTC)
I've made some changes recently (committed & pushed) to the orthopoly
package and kron_delta. Details are shown in the attached log. In
summary:
* fixed bugs in definition of kron_delta Lisp function
* make comments in orthopoly agree with code
* return simplifiable kron_delta expressions from orthopoly
* handle assoc_legendre_p in orthopoly_weight
* express legendre_{p,q} recurrences in terms of n + 1, n, and n - 1
About the last one, I suppose it is to some degree a matter of taste.
However, now the legendre_p & legendre_q recurrences are analogous to
other recurrences known to orthopoly_recur, and agree with A&S 22.7.
batch(test_orthopoly) reports no errors, nor does
batch(rtest_pochhammer, test).
Thanks to Kris Katterjohn for suggesting most of these changes.
best,
Robert Dodier
PS.
commit e5a8abaa23caf38de6ed397f8d6a79f5be46ae65
Author: robert_dodier <robert_dodier at users.sourceforge.net>
Date: Wed Aug 22 12:15:05 2012 -0600
Express recurrences for legendre_p and legendre_q in terms of n + 1, n, and n - 1,
as other recurrences are expressed in orthopoly, and as shown in A&S 22.7. Thanks to Kris Katterjohn:
http://sourceforge.net/tracker/?func=detail&aid=3553317&group_id=4933&atid=304933
diff --git a/doc/info/orthopoly.texi b/doc/info/orthopoly.texi
index 12911e5..61cca4a 100644
--- a/doc/info/orthopoly.texi
+++ b/doc/info/orthopoly.texi
@@ -824,10 +824,10 @@ respect to the polynomial degree.
@c ===end===
@example
(%i1) orthopoly_recur (legendre_p, [n, x]);
- (2 n - 1) P (x) x + (1 - n) P (x)
- n - 1 n - 2
-(%o1) P (x) = -----------------------------------------
- n n
+ (2 n + 1) P (x) x - n P (x)
+ n n - 1
+(%o1) P (x) = -------------------------------
+ n + 1 n + 1
@end example
The second argument to @code{orthopoly_recur} must be a list with the
diff --git a/share/orthopoly/orthopoly.lisp b/share/orthopoly/orthopoly.lisp
index d7c643a..9a851b3 100644
--- a/share/orthopoly/orthopoly.lisp
+++ b/share/orthopoly/orthopoly.lisp
@@ -1260,6 +1260,8 @@ Maxima code for evaluating orthogonal polynomials listed in Chapter 22 of Abramo
;; For recursion relations, see A & S 22.7 page 782.
+;; legendre_p(n+1,x) = ((2*n+1)*legendre_p(n,x)*x-n*legendre_p(n-1,x))/(n+1)
+
;; jacobi_p(n+1,a,b,x) = (((2*n+a+b+1)*(a^2-b^2) +
;; x*pochhammer(2*n+a+b,3)) * jacobi_p(n,a,b,x) -
;; 2*(n+a)*(n+b)*(2*n+a+b+2)*jacobi_p(n-1,a,b,x))/(2*(n+1)*(n+a+b+1)*(2*n+a+b))
@@ -1350,16 +1352,16 @@ Maxima code for evaluating orthogonal polynomials listed in Chapter 22 of Abramo
(let* ((n (nth 1 arg))
(x (nth 2 arg))
(z (if (eq fn '$legendre_q)
- `((mtimes) -1 ((%kron_delta) ,n 1)) 0)))
+ `((mtimes) -1 ((%kron_delta) ,n 0)) 0)))
(simplify
- `((mequal) ((,fn) ,n ,x)
+ `((mequal) ((,fn) ((mplus) 1 ,n) ,x)
((mplus)
- ((mtimes) ((mexpt) ,n -1)
+ ((mtimes) ((mexpt) ((mplus) 1 ,n) -1)
((mplus)
- ((mtimes) ((mplus) 1 ((mtimes) -1 ,n))
- ((,fn) ((mplus) -2 ,n) ,x))
- ((mtimes) ((mplus) -1 ((mtimes) 2 ,n))
- ((,fn) ((mplus) -1 ,n) ,x) ,x)))
+ ((mtimes) ((mtimes) -1 ,n)
+ ((,fn) ((mplus) -1 ,n) ,x))
+ ((mtimes) ((mplus) 1 ((mtimes) 2 ,n))
+ ((,fn) ,n ,x) ,x)))
,z)))))
((member fn '($assoc_legendre_p $assoc_legendre_q) :test 'eq)
commit ff248ac1ee83a69eea03a762e365f1635bea30c8
Author: robert_dodier <robert_dodier at users.sourceforge.net>
Date: Wed Aug 22 11:42:00 2012 -0600
Change recurrence relations for ultraspherical and Hermite polynomials to
agree with output of orthopoly_recur and A&S 22.7.
diff --git a/share/orthopoly/orthopoly.lisp b/share/orthopoly/orthopoly.lisp
index 7c63aea..d7c643a 100644
--- a/share/orthopoly/orthopoly.lisp
+++ b/share/orthopoly/orthopoly.lisp
@@ -1265,7 +1265,7 @@ Maxima code for evaluating orthogonal polynomials listed in Chapter 22 of Abramo
;; 2*(n+a)*(n+b)*(2*n+a+b+2)*jacobi_p(n-1,a,b,x))/(2*(n+1)*(n+a+b+1)*(2*n+a+b))
;; ultraspherical(n+1,a,x) = (2*(n+a)*x * ultraspherical(n,a,x) -
-;; (n+2*a+1)*ultraspherical(n-1,a,x))/(n+1)
+;; (n+2*a-1)*ultraspherical(n-1,a,x))/(n+1)
;; chebyshev_t(n+1,x) = 2*x*chebyshev_t(n,x) -chebyshev_t(n-1,x)
@@ -1276,7 +1276,7 @@ Maxima code for evaluating orthogonal polynomials listed in Chapter 22 of Abramo
;; gen_laguerre(n+1,a,x) = (((2*n+a+1) - x)*gen_laguerre(n,a,x)
;; -(n+a)*gen_laguerre(n-1,a,x))/(n+1)
-;; hermite(n+1,x) = 2*x*hermite(n,x) +2*n*hermite(n-1,x)
+;; hermite(n+1,x) = 2*x*hermite(n,x) -2*n*hermite(n-1,x)
;; See G & R 8.733.2; A & S 22.7.11 might be wrong -- or maybe I need
;; reading glasses.
commit 4743ec4737e9dba1b9cbededdf86edc8fa18c2be
Author: robert_dodier <robert_dodier at users.sourceforge.net>
Date: Wed Aug 22 11:18:50 2012 -0600
Apply patch for assoc_legendre_p in orthopoly_weight. Thanks to Kris Katterjohn:
http://sourceforge.net/tracker/?func=detail&aid=3553910&group_id=4933&atid=304933
diff --git a/share/orthopoly/orthopoly.lisp b/share/orthopoly/orthopoly.lisp
index 84607e6..7c63aea 100644
--- a/share/orthopoly/orthopoly.lisp
+++ b/share/orthopoly/orthopoly.lisp
@@ -1484,6 +1484,13 @@ variable ~:M" arg (car (last arg))))
(check-arg-length fn 2 (- (length arg) 1))
`((mlist) 1 -1 1))
+ ; This is for a fixed order. There is also an orthogonality
+ ; condition for fixed degree with weight function 1/(1-x^2).
+ ; See A & S 8.14.11 and 8.14.12.
+ ((eq fn '$assoc_legendre_p)
+ (check-arg-length fn 3 (- (length arg) 1))
+ `((mlist) 1 -1 1))
+
((eq fn '$laguerre)
(check-arg-length fn 2 (- (length arg) 1))
(let ((x (nth 2 arg)))
commit 8bd381f1261024d64a3f4a84906a21699cbd3df0
Author: robert_dodier <robert_dodier at users.sourceforge.net>
Date: Sat Aug 18 15:01:56 2012 -0600
In package orthopoly, return expressions containing %KRON_DELTA instead of $KRON_DELTA
since simplification property is attached to %KRON_DELTA. With this change, batch(test_orthopoly)
reports no errors.
diff --git a/share/orthopoly/orthopoly.lisp b/share/orthopoly/orthopoly.lisp
index a2ff70b..84607e6 100644
--- a/share/orthopoly/orthopoly.lisp
+++ b/share/orthopoly/orthopoly.lisp
@@ -581,7 +581,7 @@ Maxima code for evaluating orthogonal polynomials listed in Chapter 22 of Abramo
'((n x)
((unk) "$first" "$legendre_p")
((mplus)
- ((mtimes) -1 (($kron_delta) 0 n)
+ ((mtimes) -1 ((%kron_delta) 0 n)
((mexpt) ((mplus) -1 ((mexpt) x 2)) -1))
((mtimes)
((mplus)
@@ -997,7 +997,7 @@ Maxima code for evaluating orthogonal polynomials listed in Chapter 22 of Abramo
(defun $spherical_bessel_j (n x)
(cond ((and (eq '$zero (csign ($ratdisrep x)))
(or (integerp n) ($featurep n '$integer)))
- `(($kron_delta) 0 ,n))
+ `((%kron_delta) 0 ,n))
((and (use-float x) (integerp n))
(let ((d 1) (xr) (xi) (z))
@@ -1350,7 +1350,7 @@ Maxima code for evaluating orthogonal polynomials listed in Chapter 22 of Abramo
(let* ((n (nth 1 arg))
(x (nth 2 arg))
(z (if (eq fn '$legendre_q)
- `((mtimes) -1 (($kron_delta) ,n 1)) 0)))
+ `((mtimes) -1 ((%kron_delta) ,n 1)) 0)))
(simplify
`((mequal) ((,fn) ,n ,x)
((mplus)
commit a82e9cbc10664d48a0eb012f4feb0e9029fa561a
Author: robert_dodier <robert_dodier at users.sourceforge.net>
Date: Sat Aug 18 14:01:59 2012 -0600
Correct implementation of Lisp function kron_delta.
(1) In (DEFUN $KRON_DELTA ...), allow for multiple arguments
(2) ... and punt directly to simplifier (no need for detour through TAKE)
(3) ... and simplify with operator %KRON_DELTA (as $KRON_DELTA has no simplification property).
(4) Fix misspelling (knon_delta --> kron_delta) in error message.
diff --git a/src/nset.lisp b/src/nset.lisp
index e6e273b..875152e 100644
--- a/src/nset.lisp
+++ b/src/nset.lisp
@@ -863,7 +863,7 @@
(setf (get '%kron_delta 'verb) '$kron_delta)
(setf (get '$kron_delta 'alias) '%kron_delta)
(setf (get '%kron_delta 'reversealias) '$kron_delta)
-(defun $kron_delta (x) (take '($kron_delta) x))
+(defun $kron_delta (&rest x) (simplifya `((%kron_delta) , at x) t))
(setf (get '%kron_delta 'real-valued) t) ;; conjugate(kron_delta(xxx)) --> kron_delta(xxx)
(setf (get '%kron_delta 'integer-valued) t) ;; featurep(kron_delta(xxx), integer) --> true
@@ -873,7 +873,7 @@
(declare (ignore y))
(setq l (cdr l)) ;; remove (($kron_delta simp)
- (if (and l (null (cdr l))) (wna-err '$knon_delta)) ;; wrong number of arguments error for exactly one argument
+ (if (and l (null (cdr l))) (wna-err '$kron_delta)) ;; wrong number of arguments error for exactly one argument
;; Checking both mnqp and meqp is convenient, but unnecessary. This code misses simplifications that
;; involve three or more arguments. Example: kron_delta(a,b,a+b+1,a-b+5) could (but doesn't) simplify