[Maxima-bugs] [ maxima-Bugs-1920177 ] Problems with the bessel functions

[SourceForge.net <no...@so...>]

Bugs item #1920177, was opened at 2008-03-19 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 bessel-y 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 round-off errors:

Try bessel_j(-2,1.0). The result is 0,114903... - %i* 2.8142... e-17. 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 round-off 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 bessel-j 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 (hankel-1 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 (+ (hankel-1 order arg) (hankel-2 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: 2008-03-27 17:21

Message:
Logged In: YES 
user_id=2039760
Originator: YES

I have finished a first redesign of the code for the numerical calculation
of the Bessel functions. 

I have extended the calculation to or changed parts of the calculation to
negative orders and arguments.

Perhaps the ideas and the code are interesting for the project. I have
added the code to this message. A short description of the changes can be
find in the header of the file.

Crategus

File Added: bessel-new.lisp

----------------------------------------------------------------------

Comment By: Crategus (crategus)
Date: 2008-03-23 18:16

Message:
Logged In: YES 
user_id=2039760
Originator: YES

I have added to this posting the code for a routine bessel-j 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: bessel-j.lisp

----------------------------------------------------------------------

You can respond by visiting: 
https://sourceforge.net/tracker/?func=detail&atid=104933&aid=1920177&group_id=4933