Commit 76ea9b13 authored by Raymond Toy's avatar Raymond Toy
parents 25129063 5daefcc4
......@@ -156,7 +156,8 @@
(do* ((k 0 (1+ k))
(bk (bk 0 p)
(bk k p))
(ratio v
;; Compute g[k](p)/(2*k)!, not r[2*k+1](p)/(2*k)!
(ratio 1
(* ratio (/ (+ v2 (expt (1- (* 2 k)) 2))
(* 2 k (1- (* 2 k))))))
(term (* ratio bk)
......@@ -169,7 +170,7 @@
(format t " ratio = ~S~%" ratio)
(format t " term = ~S~%" term)
(format t " sum - ~S~%" sum))
(* sum #c(0 2) (/ (exp p) q)))
(* sum 4 (exp p)))
(when *debug-exparc*
(format t "k = ~D~%" k)
(format t " bk = ~S~%" bk)
......@@ -182,18 +183,18 @@
(defun integer-bessel-j-exp-arc (v z)
(let* ((iz (* #c(0 1) z))
(i+ (exp-arc-i-2 iz v)))
(cond ((= v (ftruncate v))
(cond ((and (= v (ftruncate v)) (realp z))
;; We can simplify the result
(let ((c (cis (* v (float-pi i+) -1/2))))
(let ((c (exp (* v (float-pi i+) #c(0 -1/2)))))
(/ (+ (* c i+)
(* (conjugate c) (conjugate i+)))
(float-pi i+)
2)))
(t
(let ((i- (exp-arc-i-2 (- iz ) v)))
(/ (+ (* (cis (* v (float-pi i+) -1/2))
(/ (+ (* (exp (* v (float-pi i+) #c(0 -1/2)))
i+)
(* (cis (* v (float-pi i+) 1/2))
(* (exp (* v (float-pi i+) #c(0 1/2)))
i-))
(float-pi i+)
2))))))
......@@ -213,13 +214,13 @@
(defun alpha (n z)
(let ((n (float n (realpart z))))
(/ (cf-incomplete-gamma (1+ n) (/ z 2))
(/ (incomplete-gamma (1+ n) (/ z 2))
(expt z (1+ n)))))
(defun beta (n z)
(let ((n (float n (realpart z))))
(/ (- (cf-incomplete-gamma (1+ n) (/ z 2))
(cf-incomplete-gamma (1+ n) (/ z -2)))
(/ (- (incomplete-gamma (1+ n) (/ z 2))
(incomplete-gamma (1+ n) (/ z -2)))
(expt z (1+ n)))))
;; a[0](k,v) := (k+sqrt(k^2+1))^(-v);
......@@ -375,21 +376,74 @@
(format t " f = ~S~%" f)
(format t " term = ~S~%" term)
(format t " sum = ~S~%" sum))))))
;;
;; TODO:
;; o For |z| <= 1 use the series.
;; o Currently accuracy is not good for large z and half-integer
;; order.
;; o For real v and z, return a real number instead of complex.
;; o Handle the case of Re(z) < 0. (The formulas are for Re(z) > 0:
;; bessel_j(v,z*exp(m*%pi*%i)) = exp(m*v*%pi*%i)*bessel_j(v, z)
;; o The paper suggests using
;; bessel_i(v,z) = exp(-v*%pi*%i/2)*bessel_j(v, %i*z)
;; when Im(z) >> Re(z)
;;
(defun bessel-j (v z)
(let ((vv (ftruncate v)))
(cond ((= vv v)
;; v is an integer
;; Clear the caches for now.
(an-clrhash)
(%big-a-clrhash)
(cond ((and (= vv v) (realp z))
;; v is an integer and z is real
(integer-bessel-j-exp-arc v z))
(t
;; Need to fine-tune the value of big-n.
(let ((big-n 100)
(vpi (* v (float-pi (realpart z)))))
(+ (integer-bessel-j-exp-arc v z)
(* z
(/ (sin vpi) vpi)
(+ (/ -1 z)
(sum-ab big-n v z)
(sum-big-ia big-n v z)))))))))
(if (= vv v)
0
(* z
(/ (sin vpi) vpi)
(+ (/ -1 z)
(sum-ab big-n v z)
(sum-big-ia big-n v z))))))))))
;; Bessel Y
;;
;; bessel_y(v, z) = 1/(2*%pi*%i)*(exp(-%i*v*%pi/2)*I(%i*v,z) - exp(%i*v*%pi/2)*I(-%i*z, v))
;; + z/v/%pi*((1-cos(v*%pi)/z) + S(N,z,v)*cos(v*%pi)-S(N,z,-v))
;;
;; where
;;
;; S(N,z,v) = sum(alpha[n](z)*a[n](0,v) + beta[n](z)*sum(exp(-k*z)*a[n](k,v),k,1,N),n,0,inf)
;; + sum(A[n](v)*I[n](N+1/2,z,v),n,0,inf)
;;
(defun bessel-y (v z)
(flet ((ipart (v z)
(let* ((iz (* #c(0 1) z))
(c+ (exp (* v (float-pi z) 1/2)))
(c- (exp (* v (float-pi z) -1/2)))
(i+ (exp-arc-i-2 iz v))
(i- (exp-arc-i-2 (- iz) v)))
(/ (- (* c- i+) (* c+ i-))
(* #c(0 2) (float-pi z)))))
(s (big-n z v)
(+ (sum-ab big-n v z)
(sum-big-ia big-n v z))))
(let* ((big-n 100)
(vpi (* v (float-pi z)))
(c (cos vpi)))
(+ (ipart v z)
(* (/ z vpi)
(+ (/ (- 1 c)
z)
(* c
(s big-n z v))
(- (s big-n z (- v)))))))))
(defun paris-series (v z n)
(labels ((pochhammer (a k)
......
......@@ -1542,3 +1542,140 @@
(check-accuracy 219.1 e true))
nil)
;; Bessel J tests for negative order
(rt:deftest bessel-j.neg-order.d.1
(let ((b (bessel-j -1d0 2d0))
(true -0.5767248077568734d0))
(check-accuracy 50.2 b true))
nil)
(rt:deftest bessel-j.neg-order.d.2
(let ((b (bessel-j -1d0 1.5d0))
(true -0.5579365079100996d0))
(check-accuracy 50.5 b true))
nil)
(rt:deftest bessel-j.neg-order.d.3
(let ((b (bessel-j -1.5d0 2d0))
(true -0.3956232813587035d0))
(check-accuracy 50.59 b true))
nil)
(rt:deftest bessel-j.neg-order.d.4
(let ((b (bessel-j -1.8d0 1.5d0))
(true -0.251327217627129314d0))
(check-accuracy 49.98 b true))
nil)
(rt:deftest bessel-j.neg-order.d.5
(let ((b (bessel-j -2d0 1.5d0))
(true 0.2320876721442147d0))
(check-accuracy 51.89 b true))
nil)
(rt:deftest bessel-j.neg-order.d.6
(let ((b (bessel-j -2.5d0 1.5d0))
(true 1.315037204805194d0))
(check-accuracy 52.37 b true))
nil)
(rt:deftest bessel-j.neg-order.d.7
(let ((b (bessel-j -2.3d0 1.5d0))
(true 1.012178926325313d0))
(check-accuracy 50.01 b true))
nil)
;; Bessel-J tests for positive order
(rt:deftest bessel-j.pos-order.d.1
(let ((b (bessel-j 1.5d0 1d0))
(true 0.2402978391234270d0))
(check-accuracy 51.83 b true))
nil)
(rt:deftest bessel-j.pos-order.d.2
(let ((b (bessel-j 1.8d0 1d0))
(true 0.1564953153109239d0))
(check-accuracy 51.97 b true))
nil)
(rt:deftest bessel-j.pos-order.d.3
(let ((b (bessel-j 2d0 1d0))
(true 0.1149034849319005d0))
(check-accuracy 51.87 b true))
nil)
(rt:deftest bessel-j.pos-order.d.4
(let ((b (bessel-j 2.5d0 1d0))
(true 0.04949681022847794d0))
(check-accuracy 47.17 b true))
nil)
(rt:deftest bessel-j.pos-order.d.5
(let ((b (bessel-j -2d0 1.5d0))
(true 0.2320876721442147d0))
(check-accuracy 51.89 b true))
nil)
;; Bessel J for half integer order and real args
(rt:deftest bessel-j-1/2.d.1
(loop for k from 0 below 100
;; x in [1,1+pi/2] because we don't want to test the Bessel
;; series and we don't want to test near pi because sin(pi)
;; = 0, where we will lose accuracy.
for x = (+ 1 (random (/ pi 2)))
for b = (bessel-j 0.5d0 x)
for true = (* (/ (sin x) (sqrt x)) (sqrt (/ 2 pi)))
for result = (check-accuracy 48.42 b true)
when result
append (list (list (list k x) result)))
nil)
(rt:deftest bessel-j-1/2.d.1.a
(let* ((x 2.3831631289164497d0)
(b (bessel-j 0.5d0 x))
(true (* (/ (sin x) (sqrt x)) (sqrt (/ 2 pi)))))
(check-accuracy 48.42 b true))
nil)
(rt:deftest bessel-j-1/2.q.1
(loop for k from 0 below 10
;; x in [1,1+pi/2] because we don't want to test the Bessel
;; series and we don't want to test near pi because sin(pi)
;; = 0, where we will lose accuracy.
for x = (+ 1 (random (/ (float-pi #q1) 2)))
for b = (bessel-j #q0.5 x)
for true = (* (/ (sin x) (sqrt x)) (sqrt (/ 2 (float-pi #q1))))
for result = (check-accuracy 169.45 b true)
when result
append (list (list (list k x) result)))
nil)
(rt:deftest bessel-j-1/2.q.1.a
(let* ((x #q1.1288834862545916200627583005758663687705443417892789067029865493882q0)
(b (bessel-j #q0.5 x))
(true (* (/ (sin x) (sqrt x)) (sqrt (/ 2 (float-pi #q1))))))
(check-accuracy 182.92 b true))
nil)
(rt:deftest bessel-j-1/2.q.1.b
(let* ((x #q1.1288834862545916200627583005758663687705443417892789067029865493882q0)
(b (bessel-j #q0.5 x))
(true (* (/ (sin x) (sqrt x)) (sqrt (/ 2 (float-pi #q1))))))
(check-accuracy 173.28 b true))
nil)
(rt:deftest bessel-j-1/2.q.1.c
(let* ((x #q1.0360263937639582798798376485114581552570020473846457752365459851056q0)
(b (bessel-j #q0.5 x))
(true (* (/ (sin x) (sqrt x)) (sqrt (/ 2 (float-pi #q1))))))
(check-accuracy 169.45 b true))
nil)
;; Bessel J for complex args
(rt:deftest bessel-j-complex.pos-order.d.1
(let ((b (bessel-j 0d0 #c(1d0 1)))
(true #c(0.9376084768060293d0 -0.4965299476091221d0)))
(check-accuracy 50.73 b true))
nil)
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment