bessel_i quadpack bug 5.26.0
- Subject: bessel_i quadpack bug 5.26.0
- From: Edwin Woollett
- Date: Fri, 2 Mar 2012 14:04:29 -0800
On Mar. 2, 2012, Raymond Toy wrote:
-----------------------------------
>Anyway, I've checked in a couple of changes. Thus, since %iargs defaults
>to true,
>
>(%i2) bessel_j(1,%i*x);
>
>(%o2) %i*bessel_i(1,x)
>(%i3) bessel_i(1,%i*x);
>
>(%o3) %i*bessel_j(1,x)
>(%i4) %iargs:false;
>
>(%o4) false
>(%i5) bessel_j(1,%i*x);
>
>(%o5) bessel_j(1,%i*x)
>(%i6) bessel_i(1,%i*x);
>
>(%o6) bessel_i(1,%i*x)
>(%i7) quad_qags(imagpart(bessel_i(1,%i*x)),x,1,2);
>
>(%o7) [.5413069074167308,6.009713920026927e-15,21,0]
>
>I think in this case you can download src/bessel.lisp and compile and load
>it into your lisp to test it. That's what I did when developing this.
-----------------------------------------
Thanks for the bug fix.
-----------------------------
(%i1) :lisp (compile-file "c:/work2/bessel-new.lisp")
Compiling c:/work2/bessel-new.lisp.
End of Pass 1.
;; Note: Tail-recursive call of BESSEL-J was replaced by iteration.
;; Note: Tail-recursive call of BESSEL-Y was replaced by iteration.
;; Note: Tail-recursive call of BESSEL-I was replaced by iteration.
;; Note: Tail-recursive call of BESSEL-I was replaced by iteration.
End of Pass 2.
OPTIMIZE levels: Safety=2, Space=3, Speed=3
Finished compiling c:/work2/bessel-new.lisp.
#pc:/work2/bessel-new.o
(%i1) load("bessel-new.o");
(%o1) c:/work2/bessel-new.o
(%i2) %iargs;
(%o2) true
(%i3) bessel_j(1,%i*x);
(%o3) %i bessel_i(1, x)
(%i4) bessel_i(1,%i*x);
(%o4) %i bessel_j(1, x)
(%i5) %iargs:false$
(%i6) bessel_j(1,%i*x);
(%o6) bessel_j(1, %i x)
(%i7) bessel_i(1,%i*x);
(%o7) bessel_i(1, %i x)
(%i8) %iargs:true$
(%i9) quad_qags(imagpart(bessel_i(1,%i*x)),x,1,2);
(%o9) [0.54130690741673, 6.0097139200269267E-15, 21, 0]
(%i10) quad_qags(imagpart(bessel_i(1.0,%i*x)),x,1,2);
(%o10) [0.54130690741673, 6.0097139200269267E-15, 21, 0]
(%i1) load(nint);
(%o1) "c:/work2/nint.mac"
(%i2) nint(bessel_i(1,%i*x),x,1,2);
(%o2) 0.541307*%i
(%i3) noutL;
(%o3) [[integrate,0.541307*%i]]
(%i4) nint(bessel_i(1,%i*x),x,1,2,quad);
(%o4) 0.541307*%i
(%i5) noutL;
(%o5) [[qags,0.541307,6.00971392E-15,21,0]]
--------------------------------
Ted