From: SourceForge.net <noreply@so...>  20080323 17:16:39

Bugs item #1920177, was opened at 20080319 22:07 Message generated for change (Comment added) made by crategus You can respond by visiting: https://sourceforge.net/tracker/?func=detail&atid=104933&aid=1920177&group_id=4933 Please note that this message will contain a full copy of the comment thread, including the initial issue submission, for this request, not just the latest update. Category: Lisp Core Group: Includes proposed fix Status: Open Resolution: None Priority: 5 Private: No Submitted By: Crategus (crategus) Assigned to: Nobody/Anonymous (nobody) Summary: Problems with the bessel functions Initial Comment: There are several problems with the bessel functions. The problems I have found can be divided in the following categories: 1. Inconsistent use of the internal arrays $besselarray, $yarray and $iarray. Example: Restart Maxima and Enter bessel_y(2,1.0). You get an Lisp error. Try first bessel_y(2,1.0) and then repeat bessel_y(2,1.0). Now you get the correct result 1.65...  %i * 3,30... Try bessel_j(2,1.0), you get 0.1149...  %i*2.8142... Next enter bessel_y(2,1.0), you get a Lisp error. The reason for the problems is that the routine bessely uses the global array $YARRAY to store values, but uses also the array $BESSELARRAY to calculate the answer. This is an error. I think the best is to eliminate the use of the global arrays completly. 2. Problematic roundoff errors: Try bessel_j(2,1.0). The result is 0,114903...  %i* 2.8142... e17. The correct result is pure real. The problem is the use of the transformation j[n]=exp(n*%pi*%i)*j[n](x) in the code which produce a small imaginary part. This is a roundoff error for an integer order. In the case of non integer values the imaginary part is correct. So you can not cut off the imaginary part in general. Here is a piece of code which will return an answer which is pure real when the order is an integer (I started to rewrite the code, introduced a function besselj and eliminated the use of the global arrays): (if (integerp order) (if (evenp order) (aref jvals n) ( (aref jvals n)) ) (let ((v (* (cis (* order pi)) (aref jvals n)))) (simplifya `((mplus) ,(realpart v) ((mtimes) $%i ,(imagpart v))) nil ) ) ) 3. Wrong mathematic: Try bessel_j(2.5,1,0). You get 0.04949... Correct is the result 2.8763... * %i Or bessel_j(3,2.0). You get 0,128943... Correct is the result 0.128943... The calculation of the bessel function j[n](x) for real argument x and negativ order n as (realpart (hankel1 order arg)) is not correct. The correct result can be obtained with the formula j[n](x) = 1/2 * (H1[n](x) + H2[n](x)). I tried the following code: (let ((result (* 0.5d0 (+ (hankel1 order arg) (hankel2 order arg))))) (complexify result) ; Problem: you get a small imaginary part ;in the case of a real result (like Problem 2) ; An alternative with a function fpround which round the result ; so that a small imaginary part will vanish ; ; (simplifya ; `((mplus) ; ,(fpround (realpart result) 14) ; ((mtimes) ; ,(fpround (imagpart result) 14) ; $%i ; ) ; ) ; nil ; ) ) This is the definition of the function fpround: (defun fpround (x &optional (digits 1)) (let ((fac (expt 10 digits))) (/ (round x (/ 1 fac)) (float fac)) ) ) I have redesigned a lot of the code for the bessel functions but the work is not finished. Is the work interesting for the project? I use GCL 2.6.6 on Windows XP for programing. Crategus  >Comment By: Crategus (crategus) Date: 20080323 18:16 Message: Logged In: YES user_id=2039760 Originator: YES I have added to this posting the code for a routine besselj which handles the above problems. The global array $BESSELARRAY is not used. Futher, I added numerical test values for the the special cases in the code. I used Maxima 5.14.0 and GCL 2.6.6 ANSI to compile the code. The testsuite runs with no unexpected errors found. 22 of the 29 numerical tests I have added give the result within a presision of 16 digits. 7 results have a precision of about 14 to 15 digits. Perhaps, the code is interesting enough for further investigation of the numerical evaluation of the Bessel functions. Crategus File Added: besselj.lisp  You can respond by visiting: https://sourceforge.net/tracker/?func=detail&atid=104933&aid=1920177&group_id=4933 