Changeset 390a7483f5658fe802d5d239070cdaa573adf4a5

Show
Ignore:
Timestamp:
03/12/11 10:35:18 (3 years ago)
Author:
Raymond Toy <toy.raymond@…>
Children:
5bd5df9360268fae90fc41fcdf86b728f8a54e86
Parents:
147fa2c7aa5e99988a7c3fb35800865d22efbc51
git-committer:
Raymond Toy <toy.raymond@…> (03/12/11 10:35:18)
Message:

Add prelimary support for integrals of the 3rd kind.

qd-elliptic.lisp:
o Clean up for unused variable in ELLIPTIC-K
o Add Carlson's Rj functions
o Implement elliptic-pi using Carlson's method.

rt-tests.lisp:
o Add many tests for elliptic-pi. Some tests pass, and some fail. The

failing tests are not enabled because I don't know if the failure is
because the test itself is wrong or if the integral is wrong.

Files:
2 modified

Legend:

Unmodified
Added
Removed
  • qd-elliptic.lisp

    re2d8d6 r390a74  
    409409         (/ (float +pi+ m) 2)) 
    410410        (t 
    411          (let ((precision (float-contagion m))) 
    412            (carlson-rf 0 (- 1 m) 1))))) 
     411         (carlson-rf 0 (- 1 m) 1)))) 
    413412 
    414413;; Elliptic integral of the first kind.  This is computed using 
     
    533532              (* (/ m 3) 
    534533                 (carlson-rd 0 y 1))))))) 
     534 
     535;; Carlson's Rc function. 
     536;; 
     537;; Some interesting identities: 
     538;; 
     539;; log(x)  = (x-1)*rc(((1+x)/2)^2, x), x > 0 
     540;; asin(x) = x * rc(1-x^2, 1), |x|<= 1 
     541;; acos(x) = sqrt(1-x^2)*rc(x^2,1), 0 <= x <=1 
     542;; atan(x) = x * rc(1,1+x^2) 
     543;; asinh(x) = x * rc(1+x^2,1) 
     544;; acosh(x) = sqrt(x^2-1) * rc(x^2,1), x >= 1 
     545;; atanh(x) = x * rc(1,1-x^2), |x|<=1 
     546;; 
     547 
     548(defun carlson-rc (x y) 
     549  "Compute Carlson's Rc function: 
     550 
     551  Rc(x,y) = integrate(1/2*(t+x)^(-1/2)*(t+y)^(-1), t, 0, inf)" 
     552  (let* ((precision (float-contagion x y)) 
     553         (yn (apply-contagion y precision)) 
     554         (x (apply-contagion x precision)) 
     555         xn z w a an pwr4 n epslon lambda sn s) 
     556    (cond ((and (zerop (imagpart yn)) 
     557                (minusp (realpart yn))) 
     558           (setf xn (- x y)) 
     559           (setf yn (- yn)) 
     560           (setf z yn) 
     561           (setf w (sqrt (/ x xn)))) 
     562          (t 
     563           (setf xn x) 
     564           (setf z yn) 
     565           (setf w 1))) 
     566    (setf a (/ (+ xn yn yn) 3)) 
     567    (setf epslon (/ (abs (- a xn)) (errtol x y))) 
     568    (setf an a) 
     569    (setf pwr4 1) 
     570    (setf n 0) 
     571    (loop while (> (* epslon pwr4) (abs an)) 
     572       do 
     573         (setf pwr4 (/ pwr4 4)) 
     574         (setf lambda (+ (* 2 (sqrt xn) (sqrt yn)) yn)) 
     575         (setf an (/ (+ an lambda) 4)) 
     576         (setf xn (/ (+ xn lambda) 4)) 
     577         (setf yn (/ (+ yn lambda) 4)) 
     578         (incf n)) 
     579    ;; c2=3/10,c3=1/7,c4=3/8,c5=9/22,c6=159/208,c7=9/8 
     580    (setf sn (/ (* pwr4 (- z a)) an)) 
     581    (setf s (* sn sn (+ 3/10 
     582                        (* sn (+ 1/7 
     583                                 (* sn (+ 3/8 
     584                                          (* sn (+ 9/22 
     585                                                   (* sn (+ 159/208 
     586                                                            (* sn 9/8)))))))))))) 
     587    (/ (* w (+ 1 s)) 
     588       (sqrt an)))) 
     589 
     590(defun carlson-rj1 (x y z p) 
     591  (let* ((xn x) 
     592         (yn y) 
     593         (zn z) 
     594         (pn p) 
     595         (en (* (- pn xn) 
     596                (- pn yn) 
     597                (- pn zn))) 
     598         (sigma 0) 
     599         (power4 1) 
     600         (k 0) 
     601         (a (/ (+ xn yn zn pn pn) 5)) 
     602         (epslon (/ (max (abs (- a xn)) 
     603                         (abs (- a yn)) 
     604                         (abs (- a zn)) 
     605                         (abs (- a pn))) 
     606                    (errtol x y z p))) 
     607         (an a) 
     608         xnroot ynroot znroot pnroot lam dn) 
     609    (loop while (> (* power4 epslon) (abs an)) 
     610       do 
     611       (setf xnroot (sqrt xn)) 
     612       (setf ynroot (sqrt yn)) 
     613       (setf znroot (sqrt zn)) 
     614       (setf pnroot (sqrt pn)) 
     615       (setf lam (+ (* xnroot ynroot) 
     616                    (* xnroot znroot) 
     617                    (* ynroot znroot))) 
     618       (setf dn (* (+ pnroot xnroot) 
     619                   (+ pnroot ynroot) 
     620                   (+ pnroot znroot))) 
     621       (setf sigma (+ sigma 
     622                      (/ (* power4 
     623                            (carlson-rc 1 (+ 1 (/ en (* dn dn))))) 
     624                         dn))) 
     625       (setf power4 (* power4 1/4)) 
     626       (setf en (/ en 64)) 
     627       (setf xn (* (+ xn lam) 1/4)) 
     628       (setf yn (* (+ yn lam) 1/4)) 
     629       (setf zn (* (+ zn lam) 1/4)) 
     630       (setf pn (* (+ pn lam) 1/4)) 
     631       (setf an (* (+ an lam) 1/4)) 
     632       (incf k)) 
     633    (let* ((xndev (/ (* (- a x) power4) an)) 
     634           (yndev (/ (* (- a y) power4) an)) 
     635           (zndev (/ (* (- a z) power4) an)) 
     636           (pndev (* -0.5 (+ xndev yndev zndev))) 
     637           (ee2 (+ (* xndev yndev) 
     638                   (* xndev zndev) 
     639                   (* yndev zndev) 
     640                   (* -3 pndev pndev))) 
     641           (ee3 (+ (* xndev yndev zndev) 
     642                   (* 2 ee2 pndev) 
     643                   (* 4 pndev pndev pndev))) 
     644           (ee4 (* (+ (* 2 xndev yndev zndev) 
     645                      (* ee2 pndev) 
     646                      (* 3 pndev pndev pndev)) 
     647                   pndev)) 
     648           (ee5 (* xndev yndev zndev pndev pndev)) 
     649           (s (+ 1 
     650                 (* -3/14 ee2) 
     651                 (* 1/6 ee3) 
     652                 (* 9/88 ee2 ee2) 
     653                 (* -3/22 ee4) 
     654                 (* -9/52 ee2 ee3) 
     655                 (* 3/26 ee5) 
     656                 (* -1/16 ee2 ee2 ee2) 
     657                 (* 3/10 ee3 ee3) 
     658                 (* 3/20 ee2 ee4) 
     659                 (* 45/272 ee2 ee2 ee3) 
     660                 (* -9/68 (+ (* ee2 ee5) (* ee3 ee4)))))) 
     661      (+ (* 6 sigma) 
     662         (/ (* power4 s) 
     663            (sqrt (* an an an))))))) 
     664 
     665(defun carlson-rj (x y z p) 
     666  "Compute Carlson's Rj function: 
     667 
     668  Rj(x,y,z,p) = integrate(3/2*(t+x)^(-1/2)*(t+y)^(-1/2)*(t+z)^(-1/2)*(t+p)^(-1), t, 0, inf)" 
     669  (let* ((precision (float-contagion x y z p)) 
     670         (xn (apply-contagion x precision)) 
     671         (yn (apply-contagion y precision)) 
     672         (zn (apply-contagion z precision)) 
     673         (p (apply-contagion p precision)) 
     674         (qn (- p))) 
     675    (cond ((and (and (zerop (imagpart xn)) (>= (realpart xn) 0)) 
     676                (and (zerop (imagpart yn)) (>= (realpart yn) 0)) 
     677                (and (zerop (imagpart zn)) (>= (realpart zn) 0)) 
     678                (and (zerop (imagpart qn)) (> (realpart qn) 0))) 
     679           (destructuring-bind (xn yn zn) 
     680               (sort (list xn yn zn) #'<) 
     681             (let* ((pn (+ yn (* (- zn yn) (/ (- yn xn) (+ yn qn))))) 
     682                    (s (- (* (- pn yn) (carlson-rj1 xn yn zn pn)) 
     683                          (* 3 (carlson-rf xn yn zn))))) 
     684               (setf s (+ s (* 3 (sqrt (/ (* xn yn zn) 
     685                                          (+ (* xn zn) (* pn qn)))) 
     686                               (carlson-rc (+ (* xn zn) (* pn qn)) (* pn qn))))) 
     687               (/ s (+ yn qn))))) 
     688          (t 
     689           (carlson-rj1 x y z p))))) 
     690 
     691;; Elliptic integral of the third kind: 
     692;; 
     693;; (A&S 17.2.14) 
     694;; 
     695;; PI(n; phi|m) = integrate(1/sqrt(1-m*sin(x)^2)/(1-n*sin(x)^2), x, 0, phi) 
     696;; 
     697(defun elliptic-pi (n phi m) 
     698  "Compute elliptic integral of the third kind: 
     699 
     700  PI(n; phi|m) = integrate(1/sqrt(1-m*sin(x)^2)/(1-n*sin(x)^2), x, 0, phi)" 
     701  ;; Note: Carlson's DRJ has n defined as the negative of the n given 
     702  ;; in A&S. 
     703  (let* ((precision (float-contagion n phi m)) 
     704         (n (apply-contagion n precision)) 
     705         (phi (apply-contagion phi precision)) 
     706         (m (apply-contagion m precision)) 
     707         (nn (- n)) 
     708         (sin-phi (sin phi)) 
     709         (cos-phi (cos phi)) 
     710         (k (sqrt m)) 
     711         (k2sin (* (- 1 (* k sin-phi)) 
     712                   (+ 1 (* k sin-phi))))) 
     713    (- (* sin-phi (carlson-rf (expt cos-phi 2) k2sin 1)) 
     714       (* (/ nn 3) (expt sin-phi 3) 
     715          (carlson-rj (expt cos-phi 2) k2sin 1 
     716                      (+ 1 (* nn (expt sin-phi 2)))))))) 
  • rt-tests.lisp

    r147fa2 r390a74  
    930930      (check-accuracy 212 rf true)) 
    931931  nil) 
     932 
     933;; Elliptic integral of the third kind 
     934 
     935;; elliptic-pi(0,phi,m) = elliptic-f(phi, m) 
     936(rt:deftest oct.elliptic-pi.1d 
     937    (loop for k from 0 to 100 
     938       for phi = (random (/ pi 2)) 
     939       for m = (random 1d0) 
     940       for epi = (elliptic-pi 0 phi m) 
     941       for ef = (elliptic-f phi m) 
     942       for result = (check-accuracy 53 epi ef) 
     943       unless (eq nil result) 
     944       append (list (list phi m) result)) 
     945  nil) 
     946 
     947(rt:deftest oct.elliptic-pi.1q 
     948    (loop for k from 0 below 100 
     949       for phi = (random (/ +pi+ 2)) 
     950       for m = (random #q1) 
     951       for epi = (elliptic-pi 0 phi m) 
     952       for ef = (elliptic-f phi m) 
     953       for result = (check-accuracy 53 epi ef) 
     954       unless (eq nil result) 
     955       append (list (list phi m) result)) 
     956  nil) 
     957 
     958;; DLMF 19.6.3 
     959;; 
     960;; PI(n; pi/2 | 0) = pi/(2*sqrt(1-n)) 
     961(rt:deftest oct.elliptic-pi.19.6.3.d 
     962    (loop for k from 0 below 100 
     963       for n = (random 1d0) 
     964       for epi = (elliptic-pi n (/ pi 2) 0) 
     965       for true = (/ pi (* 2 (sqrt (- 1 n)))) 
     966       for result = (check-accuracy 49 epi true) 
     967       unless (eq nil result) 
     968       append (list (list (list k n) result))) 
     969  nil) 
     970 
     971(rt:deftest oct.elliptic-pi.19.6.3.q 
     972    (loop for k from 0 below 100 
     973       for n = (random #q1) 
     974       for epi = (elliptic-pi n (/ (float-pi n) 2) 0) 
     975       for true = (/ (float-pi n) (* 2 (sqrt (- 1 n)))) 
     976       for result = (check-accuracy 210 epi true) 
     977       unless (eq nil result) 
     978       append (list (list (list k n) result))) 
     979  nil) 
     980 
     981#+nil 
     982(rt:deftest oct.elliptic-pi.19.6.2.d 
     983    (loop for k from 0 below 100 
     984       for n = (random 1d0) 
     985       for epi = (elliptic-pi (- n) (/ (float-pi n) 2) n) 
     986       for true = (+ (/ (float-pi n) 4 (sqrt (+ 1 (sqrt n)))) 
     987                     (/ (elliptic-k n) 2)) 
     988       for result = (check-accuracy 53 epi true) 
     989       when result 
     990       append (list (list (list k n) result))) 
     991  nil) 
     992 
     993 
     994#|| 
     995;; elliptic-pi(n, phi, 0) =  
     996;;   atanh(sqrt(1-n)*tan(phi))/sqrt(1-n)  n < 1 
     997;;   atanh(sqrt(n-1)*tan(phi))/sqrt(n-1)  n > 1 
     998;;   tan(phi)                             n = 1 
     999(rt:deftest oct.elliptic-pi.n0.d 
     1000    (loop for k from 0 below 100 
     1001       for phi = (random (/ pi 2)) 
     1002       for n = (random 1d0) 
     1003       for epi = (elliptic-pi n phi 0) 
     1004       for true = (/ (atanh (* (tan phi) (sqrt (- 1 n)))) 
     1005                     (sqrt (- 1 n))) 
     1006       for result = (check-accuracy 53 epi true) 
     1007       unless (eq nil result) 
     1008       append (list (list (list k n phi) result))) 
     1009  nil) 
     1010 
     1011(rt:deftest oct.elliptic-pi.n1.d 
     1012    (loop for k from 0 below 100 
     1013       for phi = (random (/ pi 2)) 
     1014       for epi = (elliptic-pi 0 phi 0) 
     1015       for true = (tan phi) 
     1016       for result = (check-accuracy 53 epi true) 
     1017       unless (eq nil result) 
     1018       append (list (list (list k phi) result))) 
     1019  nil) 
     1020 
     1021(rt:deftest oct.elliptic-pi.n2.d 
     1022    (loop for k from 0 below 100 
     1023       for phi = (random (/ pi 2)) 
     1024       for n = (+ 1d0 (random 100d0)) 
     1025       for epi = (elliptic-pi n phi 0) 
     1026       for true = (/ (atanh (* (tan phi) (sqrt (- n 1)))) 
     1027                     (sqrt (- n 1))) 
     1028       for result = (check-accuracy 52 epi true) 
     1029       ;; Not sure if this formula holds when atanh gives a complex 
     1030       ;; result.  Wolfram doesn't say 
     1031       when (and (not (complexp true)) result) 
     1032       append (list (list (list k n phi) result))) 
     1033  nil) 
     1034 
     1035(rt:deftest oct.elliptic-pi.n0.q 
     1036    (loop for k from 0 below 100 
     1037       for phi = (random (/ +pi+ 2)) 
     1038       for n = (random #q1) 
     1039       for epi = (elliptic-pi n phi 0) 
     1040       for true = (/ (atanh (* (tan phi) (sqrt (- 1 n)))) 
     1041                     (sqrt (- 1 n))) 
     1042       for result = (check-accuracy 212 epi true) 
     1043       unless (eq nil result) 
     1044       append (list (list (list k n phi) result))) 
     1045  nil) 
     1046 
     1047(rt:deftest oct.elliptic-pi.n1.q 
     1048    (loop for k from 0 below 100 
     1049       for phi = (random (/ +pi+ 2)) 
     1050       for epi = (elliptic-pi 0 phi 0) 
     1051       for true = (tan phi) 
     1052       for result = (check-accuracy 212 epi true) 
     1053       unless (eq nil result) 
     1054       append (list (list (list k phi) result))) 
     1055  nil) 
     1056 
     1057(rt:deftest oct.elliptic-pi.n2.q 
     1058    (loop for k from 0 below 100 
     1059       for phi = (random (/ +pi+ 2)) 
     1060       for n = (+ #q1 (random #q1)) 
     1061       for epi = (elliptic-pi n phi 0) 
     1062       for true = (/ (atanh (* (tan phi) (sqrt (- n 1)))) 
     1063                     (sqrt (- n 1))) 
     1064       for result = (check-accuracy 209 epi true) 
     1065       ;; Not sure if this formula holds when atanh gives a complex 
     1066       ;; result.  Wolfram doesn't say 
     1067       when (and (not (complexp true)) result) 
     1068       append (list (list (list k n phi) result))) 
     1069  nil) 
     1070||#