Changeset a5a4c7
 Timestamp:
 04/07/12 16:28:16 (3 years ago)
 Branches:
 master
 Children:
 015558
 Parents:
 1d404a
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

qdbessel.lisp
r1d404a ra5a4c7 40 40 ;; = 1/2^(k+3/2)/p^(k+1/2) * g(k+1/2, p) 41 41 ;; 42 ;; where g(a,z) is the lower incomplete gamma function.43 ;; 44 ;; There is the continued fraction expansion for g(a,z) (see42 ;; where G(a,z) is the lower incomplete gamma function. 43 ;; 44 ;; There is the continued fraction expansion for G(a,z) (see 45 45 ;; cfincompletegamma in qdgamma.lisp): 46 46 ;; 47 ;; g(a,z) = z^a*exp(z)/ CF47 ;; G(a,z) = z^a*exp(z)/ CF 48 48 ;; 49 49 ;; So … … 185 185 2))) 186 186 187 ;; alpha[n](z) = integrate(exp(z*s)*s^n, s, 0, 1/2) 188 ;; beta[n](z) = integrate(exp(z*s)*s^n, s, 1/2, 1/2) 189 ;; 190 ;; The recurrence in [2] is 191 ;; 192 ;; alpha[n](z) =  exp(z/2)/2^n/z + n/z*alpha[n1](z) 193 ;; beta[n]z) = ((1)^n*exp(z/2)exp(z/2))/2^n/z + n/z*beta[n1](z) 194 ;; 195 ;; We also note that 196 ;; 197 ;; alpha[n](z) = G(n+1,z/2)/z^(n+1) 198 ;; beta[n](z) = G(n+1,z/2)/z^(n+1)  G(n+1,z/2)/z^(n+1) 199 200 (defun alpha (n z) 201 (let ((n (float n (realpart z)))) 202 (/ (cfincompletegamma (1+ n) (/ z 2)) 203 (expt z (1+ n))))) 204 205 (defun beta (n z) 206 (let ((n (float n (realpart z)))) 207 (/ ( (cfincompletegamma (1+ n) (/ z 2)) 208 (cfincompletegamma (1+ n) (/ z 2))) 209 (expt z (1+ n))))) 210 211 ;; a[0](k,v) := (k+sqrt(k^2+1))^(v); 212 ;; a[1](k,v) := v*a[0](k,v)/sqrt(k^2+1); 213 ;; a[n](k,v) := 1/(k^2+1)/(n1)/n*((v^2(n2)^2)*a[n2](k,v)k*(n1)*(2*n3)*a[n1](k,v)); 214 215 ;; Convert this to iteration instead of using this quickanddirty 216 ;; memoization? 217 (let ((hash (makehashtable :test 'equal))) 218 (defun anclrhash () 219 (clrhash hash)) 220 (defun andumphash () 221 (maphash #'(lambda (k v) 222 (format t "~S > ~S~%" k v)) 223 hash)) 224 (defun an (n k v) 225 (or (gethash (list n k v) hash) 226 (let ((result 227 (cond ((= n 0) 228 (expt (+ k (sqrt (float (1+ (* k k)) (realpart v)))) ( v))) 229 ((= n 1) 230 ( (/ (* v (an 0 k v)) 231 (sqrt (float (1+ (* k k)) (realpart v)))))) 232 (t 233 (/ ( (* ( (* v v) (expt ( n 2) 2)) (an ( n 2) k v)) 234 (* k ( n 1) (+ n n 3) (an ( n 1) k v))) 235 (+ 1 (* k k)) 236 ( n 1) 237 n))))) 238 (setf (gethash (list n k v) hash) result) 239 result)))) 240 241 ;; SUMAN computes the series 242 ;; 243 ;; sum(exp(k*z)*a[n](k,v), k, 1, N) 244 ;; 245 (defun suman (bign n v z) 246 (let ((sum 0)) 247 (loop for k from 1 upto bign 248 do 249 (incf sum (* (exp ( (* k z))) 250 (an n k v)))) 251 sum)) 252 253 ;; SUMAB computes the series 254 ;; 255 ;; sum(alpha[n](z)*a[n](0,v) + beta[n](z)*sum_an(N, n, v, z), n, 0, inf) 256 (defun sumab (bign v z) 257 (let ((eps (epsilon (realpart z)))) 258 (anclrhash) 259 (do* ((n 0 (+ 1 n)) 260 (term (+ (* (alpha n z) (an n 0 v)) 261 (* (beta n z) (suman bign n v z))) 262 (+ (* (alpha n z) (an n 0 v)) 263 (* (beta n z) (suman bign n v z)))) 264 (sum term (+ sum term))) 265 ((<= (abs term) (* eps (abs sum))) 266 sum) 267 (when nil 268 (format t "n = ~D~%" n) 269 (format t " term = ~S~%" term) 270 (format t " sum = ~S~%" sum))))) 271 272 ;; Convert to iteration instead of this quickanddirty memoization? 273 (let ((hash (makehashtable :test 'equal))) 274 (defun %bigaclrhash () 275 (clrhash hash)) 276 (defun %bigadumphash () 277 (maphash #'(lambda (k v) 278 (format t "~S > ~S~%" k v)) 279 hash)) 280 (defun %biga (n v) 281 (or (gethash (list n v) hash) 282 (let ((result 283 (cond ((zerop n) 284 (expt 2 ( v))) 285 (t 286 (* (%biga ( n 1) v) 287 (/ (* (+ v n n 2) (+ v n n 1)) 288 (* 4 n (+ n v)))))))) 289 (setf (gethash (list n v) hash) result) 290 result)))) 291 292 ;; Computes A[n](v) = 293 ;; (1)^n*v*2^(v)*pochhammer(v+n+1,n1)/(2^(2*n)*n!) If v is a 294 ;; negative integer m, use A[n](m) = (1)^(m+1)*A[nm](m) for n >= 295 ;; m. 296 (defun biga (n v) 297 (let ((m (ftruncate v))) 298 (cond ((and (= m v) (minusp m)) 299 (if (< n m) 300 (%biga n v) 301 (let ((result (%biga (+ n m) v))) 302 (if (oddp (truncate m)) 303 result 304 ( result))))) 305 (t 306 (%biga n v))))) 307 308 ;; I[n](t, z, v) = exp(t*z)/t^(2*n+v1) * 309 ;; integrate(exp(t*z*s)*(1+s)^(2*nv), s, 0, inf) 310 ;; 311 ;; Use the substitution u=1+s to get a new integral 312 ;; 313 ;; integrate(exp(t*z*s)*(1+s)^(2*nv), s, 0, inf) 314 ;; = exp(t*z) * integrate(u^(v2*n)*exp(t*u*z), u, 1, inf) 315 ;; = exp(t*z)*t^(v+2*n1)*z^(v+2*n1)*incomplete_gamma_tail(1v2*n,t*z) 316 ;; 317 ;; The continued fraction for incomplete_gamma_tail(a,z) is 318 ;; 319 ;; z^a*exp(z)/CF 320 ;; 321 ;; So incomplete_gamma_tail(1v2*n, t*z) is 322 ;; 323 ;; (t*z)^(1v2*n)*exp(t*z)/CF 324 ;; 325 ;; which finally gives 326 ;; 327 ;; integrate(exp(t*z*s)*(1+s)^(2*nv), s, 0, inf) 328 ;; = CF 329 ;; 330 ;; and I[n](t, z, v) = exp(t*z)/t^(2*n+v1)/CF 331 (defun bigi (n t z v) 332 (/ (exp ( (* t z))) 333 (expt t (+ n n v 1)) 334 (let* ((a ( 1 v n n)) 335 (za ( z a))) 336 (lentz #'(lambda (n) 337 (+ n n 1 za)) 338 #'(lambda (n) 339 (* n ( a n))))))) 340 341 (defun sumbigia (bign v z) 342 ) 343 344 (defun besselj (v z) 345 (let ((vv (ftruncate v))) 346 (cond ((= vv v) 347 ;; v is an integer 348 (integerbesseljexparc v z)) 349 (t 350 (let ((bign 100) 351 (vpi (* v (floatpi (realpart z))))) 352 (+ (integerbesseljexparc v z) 353 (* z 354 (/ (sin vpi) vpi) 355 (+ (/ 1 z) 356 (sumab bign v z))))))))) 357 187 358 188 359 (defun parisseries (v z n)
Note: See TracChangeset
for help on using the changeset viewer.