Changeset a5a4c7


Ignore:
Timestamp:
04/07/12 16:28:16 (2 years ago)
Author:
Raymond Toy <toy.raymond@…>
Branches:
master
Children:
015558
Parents:
1d404a
Message:

Add more parts of the exp-arc algorithm. Needs lots of work, but it
seems that bessel_j(n,z) mostly works.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • qd-bessel.lisp

    r1d404a ra5a4c7  
    4040;;         = 1/2^(k+3/2)/p^(k+1/2) * g(k+1/2, p)
    4141;;
    42 ;; where g(a,z) is the lower incomplete gamma function.
    43 ;;
    44 ;; There is the continued fraction expansion for g(a,z) (see
     42;; where G(a,z) is the lower incomplete gamma function.
     43;;
     44;; There is the continued fraction expansion for G(a,z) (see
    4545;; cf-incomplete-gamma in qd-gamma.lisp):
    4646;;
    47 ;;  g(a,z) = z^a*exp(-z)/ CF
     47;;  G(a,z) = z^a*exp(-z)/ CF
    4848;;
    4949;; So
     
    185185       2)))
    186186
     187;; alpha[n](z) = integrate(exp(-z*s)*s^n, s, 0, 1/2)
     188;; beta[n](z)  = integrate(exp(-z*s)*s^n, s, -1/2, 1/2)
     189;;
     190;; The recurrence in [2] is
     191;;
     192;; alpha[n](z) = - exp(-z/2)/2^n/z + n/z*alpha[n-1](z)
     193;; beta[n]z)   = ((-1)^n*exp(z/2)-exp(-z/2))/2^n/z + n/z*beta[n-1](z)
     194;;
     195;; We also note that
     196;;
     197;; alpha[n](z) = G(n+1,z/2)/z^(n+1)
     198;; beta[n](z)  = G(n+1,z/2)/z^(n+1) - G(n+1,-z/2)/z^(n+1)
     199
     200(defun alpha (n z)
     201  (let ((n (float n (realpart z))))
     202    (/ (cf-incomplete-gamma (1+ n) (/ z 2))
     203       (expt z (1+ n)))))
     204
     205(defun beta (n z)
     206  (let ((n (float n (realpart z))))
     207    (/ (- (cf-incomplete-gamma (1+ n) (/ z 2))
     208          (cf-incomplete-gamma (1+ n) (/ z -2)))
     209       (expt z (1+ n)))))
     210
     211;; a[0](k,v) := (k+sqrt(k^2+1))^(-v);
     212;; a[1](k,v) := -v*a[0](k,v)/sqrt(k^2+1);
     213;; a[n](k,v) := 1/(k^2+1)/(n-1)/n*((v^2-(n-2)^2)*a[n-2](k,v)-k*(n-1)*(2*n-3)*a[n-1](k,v));
     214
     215;; Convert this to iteration instead of using this quick-and-dirty
     216;; memoization?
     217(let ((hash (make-hash-table :test 'equal)))
     218  (defun an-clrhash ()
     219    (clrhash hash))
     220  (defun an-dump-hash ()
     221    (maphash #'(lambda (k v)
     222                 (format t "~S -> ~S~%" k v))
     223             hash))
     224  (defun an (n k v)
     225    (or (gethash (list n k v) hash)
     226        (let ((result
     227                (cond ((= n 0)
     228                       (expt (+ k (sqrt (float (1+ (* k k)) (realpart v)))) (- v)))
     229                      ((= n 1)
     230                       (- (/ (* v (an 0 k v))
     231                             (sqrt (float (1+ (* k k)) (realpart v))))))
     232                      (t
     233                       (/ (- (* (- (* v v) (expt (- n 2) 2)) (an (- n 2) k v))
     234                             (* k (- n 1) (+ n n -3) (an (- n 1) k v)))
     235                          (+ 1 (* k k))
     236                          (- n 1)
     237                          n)))))
     238          (setf (gethash (list n k v) hash) result)
     239          result))))
     240
     241;; SUM-AN computes the series
     242;;
     243;; sum(exp(-k*z)*a[n](k,v), k, 1, N)
     244;;
     245(defun sum-an (big-n n v z)
     246  (let ((sum 0))
     247    (loop for k from 1 upto big-n
     248          do
     249             (incf sum (* (exp (- (* k z)))
     250                          (an n k v))))
     251    sum))
     252
     253;; SUM-AB computes the series
     254;;
     255;; sum(alpha[n](z)*a[n](0,v) + beta[n](z)*sum_an(N, n, v, z), n, 0, inf)
     256(defun sum-ab (big-n v z)
     257  (let ((eps (epsilon (realpart z))))
     258    (an-clrhash)
     259    (do* ((n 0 (+ 1 n))
     260          (term (+ (* (alpha n z) (an n 0 v))
     261                   (* (beta n z) (sum-an big-n n v z)))
     262                (+ (* (alpha n z) (an n 0 v))
     263                   (* (beta n z) (sum-an big-n n v z))))
     264          (sum term (+ sum term)))
     265         ((<= (abs term) (* eps (abs sum)))
     266          sum)
     267      (when nil
     268        (format t "n = ~D~%" n)
     269        (format t " term = ~S~%" term)
     270        (format t " sum  = ~S~%" sum)))))
     271
     272;; Convert to iteration instead of this quick-and-dirty memoization?
     273(let ((hash (make-hash-table :test 'equal)))
     274  (defun %big-a-clrhash ()
     275    (clrhash hash))
     276  (defun %big-a-dump-hash ()
     277    (maphash #'(lambda (k v)
     278                 (format t "~S -> ~S~%" k v))
     279             hash))
     280  (defun %big-a (n v)
     281    (or (gethash (list n v) hash)
     282        (let ((result
     283                (cond ((zerop n)
     284                       (expt 2 (- v)))
     285                      (t
     286                       (* (%big-a (- n 1) v)
     287                          (/ (* (+ v n n -2) (+ v n n -1))
     288                             (* 4 n (+ n v))))))))
     289          (setf (gethash (list n v) hash) result)
     290          result))))
     291
     292;; Computes A[n](v) =
     293;; (-1)^n*v*2^(-v)*pochhammer(v+n+1,n-1)/(2^(2*n)*n!)  If v is a
     294;; negative integer -m, use A[n](-m) = (-1)^(m+1)*A[n-m](m) for n >=
     295;; m.
     296(defun big-a (n v)
     297  (let ((m (ftruncate v)))
     298    (cond ((and (= m v) (minusp m))
     299           (if (< n m)
     300               (%big-a n v)
     301               (let ((result (%big-a (+ n m) v)))
     302                 (if (oddp (truncate m))
     303                     result
     304                     (- result)))))
     305          (t
     306           (%big-a n v)))))
     307
     308;; I[n](t, z, v) = exp(-t*z)/t^(2*n+v-1) *
     309;;                  integrate(exp(-t*z*s)*(1+s)^(-2*n-v), s, 0, inf)
     310;;
     311;; Use the substitution u=1+s to get a new integral
     312;;
     313;; integrate(exp(-t*z*s)*(1+s)^(-2*n-v), s, 0, inf)
     314;;   = exp(t*z) * integrate(u^(-v-2*n)*exp(-t*u*z), u, 1, inf)
     315;;   = exp(t*z)*t^(v+2*n-1)*z^(v+2*n-1)*incomplete_gamma_tail(1-v-2*n,t*z)
     316;;
     317;; The continued fraction for incomplete_gamma_tail(a,z) is
     318;;
     319;;   z^a*exp(-z)/CF
     320;;
     321;; So incomplete_gamma_tail(1-v-2*n, t*z) is
     322;;
     323;;   (t*z)^(1-v-2*n)*exp(-t*z)/CF
     324;;
     325;; which finally gives
     326;;
     327;; integrate(exp(-t*z*s)*(1+s)^(-2*n-v), s, 0, inf)
     328;;   = CF
     329;;
     330;; and I[n](t, z, v) = exp(-t*z)/t^(2*n+v-1)/CF
     331(defun big-i (n t z v)
     332  (/ (exp (- (* t z)))
     333     (expt t (+ n n v -1))
     334     (let* ((a (- 1 v n n))
     335            (z-a (- z a)))
     336       (lentz #'(lambda (n)
     337                  (+ n n 1 z-a))
     338              #'(lambda (n)
     339                  (* n (- a n)))))))
     340
     341(defun sum-big-ia (big-n v z)
     342  )
     343
     344(defun bessel-j (v z)
     345  (let ((vv (ftruncate v)))
     346    (cond ((= vv v)
     347           ;; v is an integer
     348           (integer-bessel-j-exp-arc v z))
     349          (t
     350           (let ((big-n 100)
     351                 (vpi (* v (float-pi (realpart z)))))
     352             (+ (integer-bessel-j-exp-arc v z)
     353                (* z
     354                   (/ (sin vpi) vpi)
     355                   (+ (/ -1 z)
     356                      (sum-ab big-n v z)))))))))
     357
    187358
    188359(defun paris-series (v z n)
Note: See TracChangeset for help on using the changeset viewer.