Changeset 7c5a3186070096ee93e16f2ddf51b2c84e7c5895 for qdbessel.lisp
 Timestamp:
 04/11/12 09:18:31 (2 years ago)
 Children:
 c0e12ddf6b61f571555d46c6f168e6bebabc80b1, 25129063aaa6168763b447db6f033359a19b60e0
 Parents:
 bba9f8940c9f904bf14adc405d795a38ac333c24
 gitcommitter:
 Raymond Toy <rtoy@…> (04/11/12 09:18:31)
 Files:

 1 modified
Legend:
 Unmodified
 Added
 Removed

qdbessel.lisp
re9cd1a r7c5a31 38 38 ;; B[k](p) = 1/2^(k+3/2)*integrate(exp(p*u)*u^(k1/2),u,0,1) 39 39 ;; = 1/2^(k+3/2)/p^(k+1/2)*integrate(t^(k1/2)*exp(t),t,0,p) 40 ;; = 1/2^(k+3/2)/p^(k+1/2) * g(k+1/2, p)40 ;; = 1/2^(k+3/2)/p^(k+1/2) * G(k+1/2, p) 41 41 ;; 42 42 ;; where G(a,z) is the lower incomplete gamma function. … … 110 110 ;; 111 111 ;; Hence I(%i*z, v) = conj(I(%i*z, v)) when both z and v are real. 112 ;; 113 ;; Also note that when v is an integer of the form (2*m+1)/2, then 114 ;; r[2*k+1](2*%i*v) = r[2*k+1](%i*(2*m+1)) 115 ;; = %i*(2*m+1)*product((2*m+1)^2+(2*j1)^2, j, 1, k) 116 ;; so the product is zero when k >= m and the series I(p, q) is 117 ;; finite. 112 118 (defun exparci (p q) 113 119 (let* ((sqrt2 (sqrt (float 2 (realpart p)))) … … 145 151 146 152 (defun exparci2 (p q) 147 (let* ((sqrt2 (sqrt (float 2 (realpart p)))) 148 (exp/p/sqrt2 (/ (exp ( p)) p sqrt2)) 149 (v (* #c(0 2) q)) 153 (let* ((v (* #c(0 2) q)) 150 154 (v2 (expt v 2)) 151 155 (eps (epsilon (realpart p)))) 152 (when *debugexparc*153 (format t "sqrt2 = ~S~%" sqrt2)154 (format t "exp/p/sqrt2 = ~S~%" exp/p/sqrt2))155 156 (do* ((k 0 (1+ k)) 156 157 (bk (bk 0 p) … … 163 164 (sum term (+ sum term))) 164 165 ((< (abs term) (* (abs sum) eps)) 166 (when *debugexparc* 167 (format t "Final k= ~D~%" k) 168 (format t " bk = ~S~%" bk) 169 (format t " ratio = ~S~%" ratio) 170 (format t " term = ~S~%" term) 171 (format t " sum  ~S~%" sum)) 165 172 (* sum #c(0 2) (/ (exp p) q))) 166 173 (when *debugexparc*