Changeset 390a74


Ignore:
Timestamp:
03/12/11 18:35:18 (4 years ago)
Author:
Raymond Toy <toy.raymond@…>
Branches:
master
Children:
5bd5df
Parents:
147fa2
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 edited

Legend:

Unmodified
Added
Removed
  • qd-elliptic.lisp

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

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