-----Urspr?ngliche Nachricht-----
Von: maxima-bounces at math.utexas.edu [mailto:maxima-bounces at math.utexas.edu] Im
Auftrag von Barton Willis
Gesendet: Sonntag, 29. M?rz 2009 17:02
An: maxima at math.utexas.edu
Betreff: [Maxima] Clozure CL bugs
> (defun bessel-j (order arg)
> <cut>
>
> ;; We have a real arg and order > 0 and order not 0 or 1
> ;; for this case we can call the function dbesj
> (multiple-value-bind (n alpha)
> (floor (float order))
> (let ((jvals (make-array (1+ n) :element-type 'double-float)))
> (slatec:dbesj (abs (float arg)) alpha (1+ n) jvals 0)
> After the call to slatec:dbesj, each member of the array jvals is 0.0. And
> that
> seems correct to me, but it's not. I don't understand how calling
> slatec:dbesj
> modifies jvals. How is this a Clozure CL bug? The complex case works:
> (%i2) bessel_j(3, 2.0 + 0.0001 * %i);
> (%o2) 1.5941915455754587E-5 %i + 0.128943249067055
>
> I don't see all that many differences between the real and complex cases in
> bessel-j.
> So I'm really stuck. Advice?
Hello Barton,
the implementation for real and complex arguments is quite different and uses
complete different algorithm.
Therefore, it is the best to concentrate only on the real case.
slatec:dbesj does not return the result of the numerical evaluation, but writes
the result into the array jvals, which is passed to the routine. The array
contains at the end not only the result for the order n, but for all orders
0,1,...,n
Here is a sample session in Lisp mode:
We define the array jvals:
MAXIMA> (setq jvals (make-array 4 :element-type 'double-float))
#(0.0 0.0 0.0 0.0)
We call dbesj the first value is the argument, the second value is the
fractional part of the order (for our order 3, the fractional part is 0), the
third argument is the integral part of the order plus 1, next the array jvals,
and the last argument might be a variable to get back an error status:
MAXIMA> (slatec:dbesj 2.0 0 4 jvals 0)
NIL
NIL
NIL
NIL
0
dbesj does not return the result of the numerical evaluation, but writes the
result into the array jvals. The entries into the array correspond to the
results of bessel_j(0,2.0), bessel_j(1,2.0), bessel_j(2,2.0) and
bessel_j(3,2.0):
MAXIMA> jvals
#(0.22389077914123579 0.57672480775687363 0.35283402861563784
0.12894324947440208)
You can compare the entries in the array with a call to $bessel_j:
MAXIMA> ($bessel_j 3 2.0)
0.12894324947440208
MAXIMA> ($bessel_j 2 2.0)
0.35283402861563773
MAXIMA> ($bessel_j 1 2.0)
0.5767248077568734
MAXIMA> ($bessel_j 0 2.0)
0.22389077914123562
So the best would be at first to check if the array jvals after the calculation
contains the expected values.
Dieter Kaiser