Commit 78801d6381aaaf4f21967582680f26889582db60

Authored by Raymond Toy
1 parent 5566bc39

* Add comments for {{{integer-bessel-j-exp-arc}}}.

* Simplify {{{sum-an}}} so we stop the sum when the terms no longer
   contribute to the sum.
 * Change {{{big-n}}}.  This still needs work.
Showing 1 changed file with 22 additions and 1 deletions
qd-bessel.lisp
... ... @@ -179,7 +179,13 @@
179 179 (format t " sum - ~S~%" sum)))))
180 180  
181 181  
  182 +;; Not really just for Bessel J for integer orders, but in that case,
  183 +;; this is all that's needed to compute Bessel J. For other values,
  184 +;; this is just part of the computation needed.
182 185 ;;
  186 +;; Compute
  187 +;;
  188 +;; 1/(2*%pi) * (exp(-%i*v*%pi/2) * I(%i*z, v) + exp(%i*v*%pi/2) * I(-%i*z, v))
183 189 (defun integer-bessel-j-exp-arc (v z)
184 190 (let* ((iz (* #c(0 1) z))
185 191 (i+ (exp-arc-i-2 iz v)))
... ... @@ -257,6 +263,7 @@
257 263 ;;
258 264 ;; sum(exp(-k*z)*a[n](k,v), k, 1, N)
259 265 ;;
  266 +#+nil
260 267 (defun sum-an (big-n n v z)
261 268 (let ((sum 0))
262 269 (loop for k from 1 upto big-n
... ... @@ -265,6 +272,20 @@
265 272 (an n k v))))
266 273 sum))
267 274  
  275 +;; Like above, but we just stop when the terms no longer contribute to
  276 +;; the sum.
  277 +(defun sum-an (big-n n v z)
  278 + (let ((eps (epsilon (realpart z))))
  279 + (do* ((k 1 (+ 1 k))
  280 + (term (* (exp (- (* k z)))
  281 + (an n k v))
  282 + (* (exp (- (* k z)))
  283 + (an n k v)))
  284 + (sum term (+ sum term)))
  285 + ((or (<= (abs term) (* eps (abs sum)))
  286 + (> k big-n))
  287 + sum))))
  288 +
268 289 ;; SUM-AB computes the series
269 290 ;;
270 291 ;; sum(alpha[n](z)*a[n](0,v) + beta[n](z)*sum_an(N, n, v, z), n, 0, inf)
... ... @@ -399,7 +420,7 @@
399 420 (integer-bessel-j-exp-arc v z))
400 421 (t
401 422 ;; Need to fine-tune the value of big-n.
402   - (let ((big-n 100)
  423 + (let ((big-n 10)
403 424 (vpi (* v (float-pi (realpart z)))))
404 425 (+ (integer-bessel-j-exp-arc v z)
405 426 (if (= vv v)
... ...