source: qd-bessel.lisp @ 4c1ed0

Last change on this file since 4c1ed0 was 4c1ed0, checked in by Raymond Toy <toy.raymond@…>, 3 years ago

First cut at Bessel functions. Needs lots of work.

  • Property mode set to 100644
File size: 5.9 KB
Line 
1;;;; -*- Mode: lisp -*-
2;;;;
3;;;; Copyright (c) 2011 Raymond Toy
4;;;; Permission is hereby granted, free of charge, to any person
5;;;; obtaining a copy of this software and associated documentation
6;;;; files (the "Software"), to deal in the Software without
7;;;; restriction, including without limitation the rights to use,
8;;;; copy, modify, merge, publish, distribute, sublicense, and/or sell
9;;;; copies of the Software, and to permit persons to whom the
10;;;; Software is furnished to do so, subject to the following
11;;;; conditions:
12;;;;
13;;;; The above copyright notice and this permission notice shall be
14;;;; included in all copies or substantial portions of the Software.
15;;;;
16;;;; THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
17;;;; EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
18;;;; OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
19;;;; NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
20;;;; HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
21;;;; WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
22;;;; FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
23;;;; OTHER DEALINGS IN THE SOFTWARE.
24
25(in-package #:oct)
26
27;;; References:
28;;;
29;;; [1] Borwein, Borwein, Crandall, "Effective Laguerre Asymptotics",
30;;; http://people.reed.edu/~crandall/papers/Laguerre-f.pdf
31;;;
32;;; [2] Borwein, Borwein, Chan, "The Evaluation of Bessel Functions
33;;; via Exp-Arc Integrals", http://web.cs.dal.ca/~jborwein/bessel.pdf
34;;;
35
36(defvar *debug-exparc* nil)
37
38;; B[k](p) = 1/2^(k+3/2)*integrate(exp(-p*u)*u^(k-1/2),u,0,1)
39;;         = 1/2^(k+3/2)/p^(k+1/2)*integrate(t^(k-1/2)*exp(-t),t,0,p)
40;;         = 1/2^(k+3/2)/p^(k+1/2) * g(k+1/2, p)
41;;
42;; where g(a,z) is the lower incomplete gamma function.
43;;
44;; There is the continued fraction expansion for g(a,z) (see
45;; cf-incomplete-gamma in qd-gamma.lisp):
46;;
47;;  g(a,z) = z^a*exp(-z)/ CF
48;;
49;; So
50;;
51;;  B[k](p) = 1/2^(k+3/2)/p^(k+1/2)*p^(k+1/2)*exp(-p)/CF
52;;          = exp(-p)/2^(k+3/2)/CF
53;;
54(defun bk (k p)
55  (/ (exp (- p))
56     (* (sqrt (float 2 (realpart p))) (ash 1 (+ k 1)))
57     (let ((a (float (+ k 1/2) (realpart p))))
58       (lentz #'(lambda (n)
59                  (+ n a))
60              #'(lambda (n)
61                  (if (evenp n)
62                      (* (ash n -1) p)
63                      (- (* (+ a (ash n -1)) p))))))))
64
65;; exp-arc I function, as given in the Laguerre paper
66;;
67;; I(p, q) = 4*exp(p) * sum(g[k](-2*%i*q)/(2*k)!*B[k](p), k, 0, inf)
68;;
69;; where g[k](p) = product(p^2+(2*j-1)^2, j, 1, k) and B[k](p) as above.
70;;
71;; For computation, note that g[k](p) = g[k-1](p) * (p^2 + (2*k-1)^2)
72;; and (2*k)! = (2*k-2)! * (2*k-1) * (2*k).  Then, let
73;;
74;;  R[k](p) = g[k](p)/(2*k)!
75;;
76;; Then
77;;
78;;  R[k](p) = g[k](p)/(2*k)!
79;;          = g[k-1](p)/(2*k-2)! * (p^2 + (2*k-1)^2)/((2*k-1)*(2*k)
80;;          = R[k-1](p) * (p^2 + (2*k-1)^2)/((2*k-1)*(2*k)
81;;
82;; In the exp-arc paper, the function is defined (equivalently) as
83;;
84;; I(p, q) = 2*%i*exp(p)/q * sum(r[2*k+1](-2*%i*q)/(2*k)!*B[k](p), k, 0, inf)
85;;
86;; where r[2*k+1](p) = p*product(p^2 + (2*j-1)^2, j, 1, k)
87;;
88;; Let's note some properties of I(p, q).
89;;
90;; I(-%i*z, v) = 2*%i*exp(-%i*z)/q * sum(r[2*k+1](-2*%i*v)/(2*k)!*B[k](-%i*z))
91;;
92;; Note thate B[k](-%i*z) = 1/2^(k+3/2)*integrate(exp(%i*z*u)*u^(k-1/2),u,0,1)
93;;                        = conj(B[k](%i*z).
94;;
95;; Hence I(-%i*z, v) = conj(I(%i*z, v)) when both z and v are real.
96(defun exp-arc-i (p q)
97  (let* ((sqrt2 (sqrt (float 2 (realpart p))))
98         (exp/p/sqrt2 (/ (exp (- p)) p sqrt2))
99         (v (* #c(0 -2) q))
100         (v2 (expt v 2))
101         (eps (epsilon (realpart p))))
102    (when *debug-exparc*
103      (format t "sqrt2 = ~S~%" sqrt2)
104      (format t "exp/p/sqrt2 = ~S~%" exp/p/sqrt2))
105    (do* ((k 0 (1+ k))
106          (bk (/ (incomplete-gamma 1/2 p)
107                 2 sqrt2 (sqrt p))
108              (- (/ (* bk (- k 1/2)) 2 p)
109                 (/ exp/p/sqrt2 (ash 1 (+ k 1)))))
110          ;; ratio[k] = r[2*k+1](v)/(2*k)!.
111          ;; r[1] = v and r[2*k+1](v) = r[2*k-1](v)*(v^2 + (2*k-1)^2)
112          ;; ratio[0] = v
113          ;; and ratio[k] = r[2*k-1](v)*(v^2+(2*k-1)^2) / ((2*k-2)! * (2*k-1) * 2*k)
114          ;;              = ratio[k]*(v^2+(2*k-1)^2)/((2*k-1) * 2 * k)
115          (ratio v
116                 (* ratio (/ (+ v2 (expt (1- (* 2 k)) 2))
117                             (* 2 k (1- (* 2 k))))))
118          (term (* ratio bk)
119                (* ratio bk))
120          (sum term (+ sum term)))
121         ((< (abs term) (* (abs sum) eps))
122          (* sum #c(0 2) (/ (exp p) q)))
123      (when *debug-exparc*
124        (format t "k      = ~D~%" k)
125        (format t " bk    = ~S~%" bk)
126        (format t " ratio = ~S~%" ratio)
127        (format t " term  = ~S~%" term)
128        (format t " sum   - ~S~%" sum)))))
129
130(defun exp-arc-i-2 (p q)
131  (let* ((sqrt2 (sqrt (float 2 (realpart p))))
132         (exp/p/sqrt2 (/ (exp (- p)) p sqrt2))
133         (v (* #c(0 -2) q))
134         (v2 (expt v 2))
135         (eps (epsilon (realpart p))))
136    (when *debug-exparc*
137      (format t "sqrt2 = ~S~%" sqrt2)
138      (format t "exp/p/sqrt2 = ~S~%" exp/p/sqrt2))
139    (do* ((k 0 (1+ k))
140          (bk (bk 0 p)
141              (bk k p))
142          (ratio v
143                 (* ratio (/ (+ v2 (expt (1- (* 2 k)) 2))
144                             (* 2 k (1- (* 2 k))))))
145          (term (* ratio bk)
146                (* ratio bk))
147          (sum term (+ sum term)))
148         ((< (abs term) (* (abs sum) eps))
149          (* sum #c(0 2) (/ (exp p) q)))
150      (when *debug-exparc*
151        (format t "k      = ~D~%" k)
152        (format t " bk    = ~S~%" bk)
153        (format t " ratio = ~S~%" ratio)
154        (format t " term  = ~S~%" term)
155        (format t " sum   - ~S~%" sum)))))
156
157
158;; This currently only works for v an integer.
159;;
160(defun bessel-j-exp-arc (v z)
161  (let* ((iz (* #c(0 1) z))
162         (i+ (exp-arc-i-2 iz v))
163         (i- (exp-arc-i-2 (- iz ) v)))
164    (/ (+ (* (cis (* v (float-pi i+) -1/2))
165             i+)
166          (* (cis (* v (float-pi i+) 1/2))
167             i-))
168       (float-pi i+)
169       2)))
170
171(defun paris-series (v z n)
172  (labels ((pochhammer (a k)
173             (/ (gamma (+ a k))
174                (gamma a)))
175           (a (v k)
176             (* (/ (pochhammer (+ 1/2 v) k)
177                   (gamma (float (1+ k) z)))
178                (pochhammer (- 1/2 v) k))))
179    (* (loop for k from 0 below n
180             sum (* (/ (a v k)
181                       (expt (* 2 z) k))
182                    (/ (cf-incomplete-gamma (+ k v 1/2) (* 2 z))
183                       (gamma (+ k v 1/2)))))
184       (/ (exp z)
185          (sqrt (* 2 (float-pi z) z))))))
186
Note: See TracBrowser for help on using the repository browser.