Changeset f203389e5e78d1f001f68447ac2a9dd86dcfbbf6

Show
Ignore:
Timestamp:
04/10/12 00:45:13 (2 years ago)
Author:
Raymond Toy <toy.raymond@…>
Children:
e9cd1a46fcf2a5c0101b3473648cb242b55987e8
Parents:
f8943af6bff60e23d679089db5207b4834aa83ff
git-committer:
Raymond Toy <toy.raymond@…> (04/10/12 00:45:13)
Message:
  • Implement sum-big-ia.
  • Add series for Bessel J. (Not working yet.)
Files:
1 modified

Legend:

Unmodified
Added
Removed
  • qd-bessel.lisp

    rc7fc98 rf20338  
    325325 
    326326(defun sum-big-ia (big-n v z) 
    327   ) 
    328  
     327  (let ((big-n-1/2 (+ big-n 1/2)) 
     328        (eps (epsilon z))) 
     329    (do* ((n 0 (1+ n)) 
     330          (term (* (big-a 0 v) 
     331                   (big-i 0 big-n-1/2 z v)) 
     332                (* (big-a n v) 
     333                   (big-i n big-n-1/2 z v))) 
     334          (sum term (+ sum term))) 
     335         ((<= (abs term) (* eps (abs sum))) 
     336          sum) 
     337      #+nil 
     338      (progn 
     339        (format t "n = ~D~%" n) 
     340        (format t " term = ~S~%" term) 
     341        (format t " sum  = ~S~%" sum))))) 
     342 
     343;; Series for bessel J: 
     344;; 
     345;; (z/2)^v*sum((-1)^k/Gamma(k+v+1)/k!*(z^2//4)^k, k, 0, inf) 
     346(defun s-bessel-j (v z) 
     347  (with-floating-point-contagion (v z) 
     348    (let ((z2/4 (* z z 1/4)) 
     349          (eps (epsilon z))) 
     350      (do* ((k 0 (+ 1 k)) 
     351            (f (gamma (+ v 1)) 
     352               (* f (* k (+ v k)))) 
     353            (term (/ f) 
     354                  (/ (* (- term) z2/4) f)) 
     355            (sum term (+ sum term))) 
     356           ((<= (abs term) (* eps (abs sum))) 
     357            (* sum (expt (* z 1/2) v))) 
     358        (format t "k = ~D~%" k) 
     359        (format t " term = ~S~%" term) 
     360        (format t " sum  = ~S~%" sum))))) 
     361   
    329362(defun bessel-j (v z) 
    330363  (let ((vv (ftruncate v))) 
     
    339372                   (/ (sin vpi) vpi) 
    340373                   (+ (/ -1 z) 
    341                       (sum-ab big-n v z))))))))) 
     374                      (sum-ab big-n v z) 
     375                      (sum-big-ia big-n v z))))))))) 
    342376 
    343377(defun paris-series (v z n)