From: SourceForge.net <no...@so...> - 2008-03-19 21:07:22
|
Bugs item #1920177, was opened at 2008-03-19 22:07 Message generated for change (Tracker Item Submitted) made by Item Submitter 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: Yes 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 ---------------------------------------------------------------------- You can respond by visiting: https://sourceforge.net/tracker/?func=detail&atid=104933&aid=1920177&group_id=4933 |
From: SourceForge.net <no...@so...> - 2008-03-22 12:52:48
|
Bugs item #1920177, was opened at 2008-03-19 22:07 Message generated for change (Settings changed) 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 ---------------------------------------------------------------------- You can respond by visiting: https://sourceforge.net/tracker/?func=detail&atid=104933&aid=1920177&group_id=4933 |
From: SourceForge.net <no...@so...> - 2008-03-23 17:16:39
|
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-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 |
From: SourceForge.net <no...@so...> - 2008-03-27 16:22:00
|
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 |
From: SourceForge.net <no...@so...> - 2008-03-27 16:50:51
|
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:50 Message: Logged In: YES user_id=2039760 Originator: YES I have tested the Bessel functions with the numerical data attached to this message. For the tests I used a function TEST_BESSEL(value, result, digits) which allow to compare the value with the result within the specified digits. Crategus File Added: mytest_bessel.mac ---------------------------------------------------------------------- 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 |
From: SourceForge.net <no...@so...> - 2008-03-28 16:37:43
|
Bugs item #1920177, was opened at 2008-03-19 17:07 Message generated for change (Comment added) made by rtoy 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: Raymond Toy (rtoy) Date: 2008-03-28 12:37 Message: Logged In: YES user_id=28849 Originator: NO First, thank you very much for the fixes to the Bessel routines. I'm sure we all want the routines to be correct! Second, about $besselarray. I see the documentation for that is gone, but I think the intent for besselarray was to let the user have access to all the computed values, since some algorithms end up computing all the intermediate orders to get to the desired order. We need to think if this should go away or not. Third, it would certainly help me a lot if you actually produced a diff between your new code and the existing code. Right now, I have to go look at every single function in different places to figure out what's changed. Trying to minimize the changes would help. I know this is a lot of work for you too. Just identifying which function you actually changed would help a lot too. About the specific changes you mention in the comments. No problem with changing bessel-x-simp to simp-bessel-x. (I think most of the code uses simpbesselx, but that's too messy for my tastes.) Extending the range is very, very good. Not sure about the special tests for pure real or imaginary answers. Need to look at that some more. Definitions for Hankel functions are good. May want to name them bessel_hankel_1 to be consistent with other Bessel functions. Certainly want to add some properties like derivatives, but that can wait. All in all, I would very much like to take in your very nice changes. ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-27 12:50 Message: Logged In: YES user_id=2039760 Originator: YES I have tested the Bessel functions with the numerical data attached to this message. For the tests I used a function TEST_BESSEL(value, result, digits) which allow to compare the value with the result within the specified digits. Crategus File Added: mytest_bessel.mac ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-27 12: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 13: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 |
From: SourceForge.net <no...@so...> - 2008-03-29 00:35:05
|
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-29 01:35 Message: Logged In: YES user_id=2039760 Originator: YES Thank you very much for your answer. I have redesigned the following routines: simp-bessel-j (old name bessel-j-simp) bessel-j (old name $bessel) simp-bessel-y (old name bessel-y-simp) bessel-y simp-bessel-i (old name bessel-i-simp) bessel-i simp-bessel-k (old name bessel-k-simp) bessel-k The new routines are: $hankel_1 simp-hankel-1 $hankel_2 simp-hankel-2 If have further introduced the nouns %hankel_1 und %hankel_2. I would prefer to call the functions $hankel_1 and $hankel_2 and not bessel_hankel. The reason is that these functions are well known as Hankel functions. Because I have collected the routines in a new file it would be a good idea to do the changes in the orginal file bessel.lisp. Than it easier to produce a diff for you. I do this the next time. But first, I would like to have a second look at the global arrays $besselarray, $yarray and $iarray. I think it is not so easy to manage these global arrays in a consistent and predictable way for the user. One reason is, that the code for the numerical calculations use the other functions too (i.e. to calculate bessel-y we need bessel-j or for the calculation of bessel-i we need bessel-j and bessel-k). So $besselarray can be destroyed when calculating bessel-i. It might be possible to prevent this effect by introducing a global flag which signals the internal call of one of the routines. A second point is, that we often dont't use the numerical slatec routines directly. In this case it is necessary to do a lot of extra calculations to obtain the values for the arrays. Perhaps it is useful to decide to fill the global arrays only for the case when we can use the values of the slatec routines directly. This is possible for positive order and argument or positiv order and complex argument of the Bessel functions. Further, an additional flag could be introduced to switch the filling of the arrays on and off. It is also to decide what we do when the user calculate a new value for a Bessel function and the calculation don't fill the global array. In this case, the user may be confused by the fact that the global arrays still hold old values. But for which order and argument? For this case it may be helpful to invalidate the array every time we calculate a new value or to store additionally the order and argument for which the values in the global array are calculated. All these problems with the global arrays arise because of the extension of the routines to negative arguments and orders. In my opinion it would be the best to use the concept of the global arrays no longer. The code becomes much more complicated and the benefit for the user might be small. Crategus ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-03-28 17:37 Message: Logged In: YES user_id=28849 Originator: NO First, thank you very much for the fixes to the Bessel routines. I'm sure we all want the routines to be correct! Second, about $besselarray. I see the documentation for that is gone, but I think the intent for besselarray was to let the user have access to all the computed values, since some algorithms end up computing all the intermediate orders to get to the desired order. We need to think if this should go away or not. Third, it would certainly help me a lot if you actually produced a diff between your new code and the existing code. Right now, I have to go look at every single function in different places to figure out what's changed. Trying to minimize the changes would help. I know this is a lot of work for you too. Just identifying which function you actually changed would help a lot too. About the specific changes you mention in the comments. No problem with changing bessel-x-simp to simp-bessel-x. (I think most of the code uses simpbesselx, but that's too messy for my tastes.) Extending the range is very, very good. Not sure about the special tests for pure real or imaginary answers. Need to look at that some more. Definitions for Hankel functions are good. May want to name them bessel_hankel_1 to be consistent with other Bessel functions. Certainly want to add some properties like derivatives, but that can wait. All in all, I would very much like to take in your very nice changes. ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-27 17:50 Message: Logged In: YES user_id=2039760 Originator: YES I have tested the Bessel functions with the numerical data attached to this message. For the tests I used a function TEST_BESSEL(value, result, digits) which allow to compare the value with the result within the specified digits. Crategus File Added: mytest_bessel.mac ---------------------------------------------------------------------- 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 |
From: SourceForge.net <no...@so...> - 2008-03-29 21:00:21
|
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-29 22:00 Message: Logged In: YES user_id=2039760 Originator: YES I have put all changes of the code for the Bessel functions in the orginal file bessel.lisp which is distributed with Maxima 5.14.0 and attached this file to this message. I have tested the code with the testsuite and the values of the mytest_bessel.mac. Crategus File Added: bessel-changed.lisp ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-29 01:35 Message: Logged In: YES user_id=2039760 Originator: YES Thank you very much for your answer. I have redesigned the following routines: simp-bessel-j (old name bessel-j-simp) bessel-j (old name $bessel) simp-bessel-y (old name bessel-y-simp) bessel-y simp-bessel-i (old name bessel-i-simp) bessel-i simp-bessel-k (old name bessel-k-simp) bessel-k The new routines are: $hankel_1 simp-hankel-1 $hankel_2 simp-hankel-2 If have further introduced the nouns %hankel_1 und %hankel_2. I would prefer to call the functions $hankel_1 and $hankel_2 and not bessel_hankel. The reason is that these functions are well known as Hankel functions. Because I have collected the routines in a new file it would be a good idea to do the changes in the orginal file bessel.lisp. Than it easier to produce a diff for you. I do this the next time. But first, I would like to have a second look at the global arrays $besselarray, $yarray and $iarray. I think it is not so easy to manage these global arrays in a consistent and predictable way for the user. One reason is, that the code for the numerical calculations use the other functions too (i.e. to calculate bessel-y we need bessel-j or for the calculation of bessel-i we need bessel-j and bessel-k). So $besselarray can be destroyed when calculating bessel-i. It might be possible to prevent this effect by introducing a global flag which signals the internal call of one of the routines. A second point is, that we often dont't use the numerical slatec routines directly. In this case it is necessary to do a lot of extra calculations to obtain the values for the arrays. Perhaps it is useful to decide to fill the global arrays only for the case when we can use the values of the slatec routines directly. This is possible for positive order and argument or positiv order and complex argument of the Bessel functions. Further, an additional flag could be introduced to switch the filling of the arrays on and off. It is also to decide what we do when the user calculate a new value for a Bessel function and the calculation don't fill the global array. In this case, the user may be confused by the fact that the global arrays still hold old values. But for which order and argument? For this case it may be helpful to invalidate the array every time we calculate a new value or to store additionally the order and argument for which the values in the global array are calculated. All these problems with the global arrays arise because of the extension of the routines to negative arguments and orders. In my opinion it would be the best to use the concept of the global arrays no longer. The code becomes much more complicated and the benefit for the user might be small. Crategus ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-03-28 17:37 Message: Logged In: YES user_id=28849 Originator: NO First, thank you very much for the fixes to the Bessel routines. I'm sure we all want the routines to be correct! Second, about $besselarray. I see the documentation for that is gone, but I think the intent for besselarray was to let the user have access to all the computed values, since some algorithms end up computing all the intermediate orders to get to the desired order. We need to think if this should go away or not. Third, it would certainly help me a lot if you actually produced a diff between your new code and the existing code. Right now, I have to go look at every single function in different places to figure out what's changed. Trying to minimize the changes would help. I know this is a lot of work for you too. Just identifying which function you actually changed would help a lot too. About the specific changes you mention in the comments. No problem with changing bessel-x-simp to simp-bessel-x. (I think most of the code uses simpbesselx, but that's too messy for my tastes.) Extending the range is very, very good. Not sure about the special tests for pure real or imaginary answers. Need to look at that some more. Definitions for Hankel functions are good. May want to name them bessel_hankel_1 to be consistent with other Bessel functions. Certainly want to add some properties like derivatives, but that can wait. All in all, I would very much like to take in your very nice changes. ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-27 17:50 Message: Logged In: YES user_id=2039760 Originator: YES I have tested the Bessel functions with the numerical data attached to this message. For the tests I used a function TEST_BESSEL(value, result, digits) which allow to compare the value with the result within the specified digits. Crategus File Added: mytest_bessel.mac ---------------------------------------------------------------------- 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 |
From: SourceForge.net <no...@so...> - 2008-03-29 21:05:44
|
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-29 22:05 Message: Logged In: YES user_id=2039760 Originator: YES Next I have added a diff between the orginal and the changed file. Question: Is it possible to attach at once more than one file to a message? Crategus File Added: diff.txt ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-29 22:00 Message: Logged In: YES user_id=2039760 Originator: YES I have put all changes of the code for the Bessel functions in the orginal file bessel.lisp which is distributed with Maxima 5.14.0 and attached this file to this message. I have tested the code with the testsuite and the values of the mytest_bessel.mac. Crategus File Added: bessel-changed.lisp ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-29 01:35 Message: Logged In: YES user_id=2039760 Originator: YES Thank you very much for your answer. I have redesigned the following routines: simp-bessel-j (old name bessel-j-simp) bessel-j (old name $bessel) simp-bessel-y (old name bessel-y-simp) bessel-y simp-bessel-i (old name bessel-i-simp) bessel-i simp-bessel-k (old name bessel-k-simp) bessel-k The new routines are: $hankel_1 simp-hankel-1 $hankel_2 simp-hankel-2 If have further introduced the nouns %hankel_1 und %hankel_2. I would prefer to call the functions $hankel_1 and $hankel_2 and not bessel_hankel. The reason is that these functions are well known as Hankel functions. Because I have collected the routines in a new file it would be a good idea to do the changes in the orginal file bessel.lisp. Than it easier to produce a diff for you. I do this the next time. But first, I would like to have a second look at the global arrays $besselarray, $yarray and $iarray. I think it is not so easy to manage these global arrays in a consistent and predictable way for the user. One reason is, that the code for the numerical calculations use the other functions too (i.e. to calculate bessel-y we need bessel-j or for the calculation of bessel-i we need bessel-j and bessel-k). So $besselarray can be destroyed when calculating bessel-i. It might be possible to prevent this effect by introducing a global flag which signals the internal call of one of the routines. A second point is, that we often dont't use the numerical slatec routines directly. In this case it is necessary to do a lot of extra calculations to obtain the values for the arrays. Perhaps it is useful to decide to fill the global arrays only for the case when we can use the values of the slatec routines directly. This is possible for positive order and argument or positiv order and complex argument of the Bessel functions. Further, an additional flag could be introduced to switch the filling of the arrays on and off. It is also to decide what we do when the user calculate a new value for a Bessel function and the calculation don't fill the global array. In this case, the user may be confused by the fact that the global arrays still hold old values. But for which order and argument? For this case it may be helpful to invalidate the array every time we calculate a new value or to store additionally the order and argument for which the values in the global array are calculated. All these problems with the global arrays arise because of the extension of the routines to negative arguments and orders. In my opinion it would be the best to use the concept of the global arrays no longer. The code becomes much more complicated and the benefit for the user might be small. Crategus ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-03-28 17:37 Message: Logged In: YES user_id=28849 Originator: NO First, thank you very much for the fixes to the Bessel routines. I'm sure we all want the routines to be correct! Second, about $besselarray. I see the documentation for that is gone, but I think the intent for besselarray was to let the user have access to all the computed values, since some algorithms end up computing all the intermediate orders to get to the desired order. We need to think if this should go away or not. Third, it would certainly help me a lot if you actually produced a diff between your new code and the existing code. Right now, I have to go look at every single function in different places to figure out what's changed. Trying to minimize the changes would help. I know this is a lot of work for you too. Just identifying which function you actually changed would help a lot too. About the specific changes you mention in the comments. No problem with changing bessel-x-simp to simp-bessel-x. (I think most of the code uses simpbesselx, but that's too messy for my tastes.) Extending the range is very, very good. Not sure about the special tests for pure real or imaginary answers. Need to look at that some more. Definitions for Hankel functions are good. May want to name them bessel_hankel_1 to be consistent with other Bessel functions. Certainly want to add some properties like derivatives, but that can wait. All in all, I would very much like to take in your very nice changes. ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-27 17:50 Message: Logged In: YES user_id=2039760 Originator: YES I have tested the Bessel functions with the numerical data attached to this message. For the tests I used a function TEST_BESSEL(value, result, digits) which allow to compare the value with the result within the specified digits. Crategus File Added: mytest_bessel.mac ---------------------------------------------------------------------- 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 |
From: SourceForge.net <no...@so...> - 2008-04-01 16:41:37
|
Bugs item #1920177, was opened at 2008-03-19 17:07 Message generated for change (Comment added) made by rtoy 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: Raymond Toy (rtoy) Date: 2008-04-01 12:41 Message: Logged In: YES user_id=28849 Originator: NO Thank you very much for the updated files. (The diff should have used diff -u or diff -c, but no matter. The bessel-changed.lisp file is perfect.) Given the imminent release of Maxima, I will wait on these changes until after the release. Then I'll try to incorporate your changes. I noticed that you changed the derivatives of some of the Bessel functions. Was that intentional? Were they wrong? (Some of the derivatives need to be in a certain form for some simplifications to be done.) ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-29 17:05 Message: Logged In: YES user_id=2039760 Originator: YES Next I have added a diff between the orginal and the changed file. Question: Is it possible to attach at once more than one file to a message? Crategus File Added: diff.txt ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-29 17:00 Message: Logged In: YES user_id=2039760 Originator: YES I have put all changes of the code for the Bessel functions in the orginal file bessel.lisp which is distributed with Maxima 5.14.0 and attached this file to this message. I have tested the code with the testsuite and the values of the mytest_bessel.mac. Crategus File Added: bessel-changed.lisp ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-28 20:35 Message: Logged In: YES user_id=2039760 Originator: YES Thank you very much for your answer. I have redesigned the following routines: simp-bessel-j (old name bessel-j-simp) bessel-j (old name $bessel) simp-bessel-y (old name bessel-y-simp) bessel-y simp-bessel-i (old name bessel-i-simp) bessel-i simp-bessel-k (old name bessel-k-simp) bessel-k The new routines are: $hankel_1 simp-hankel-1 $hankel_2 simp-hankel-2 If have further introduced the nouns %hankel_1 und %hankel_2. I would prefer to call the functions $hankel_1 and $hankel_2 and not bessel_hankel. The reason is that these functions are well known as Hankel functions. Because I have collected the routines in a new file it would be a good idea to do the changes in the orginal file bessel.lisp. Than it easier to produce a diff for you. I do this the next time. But first, I would like to have a second look at the global arrays $besselarray, $yarray and $iarray. I think it is not so easy to manage these global arrays in a consistent and predictable way for the user. One reason is, that the code for the numerical calculations use the other functions too (i.e. to calculate bessel-y we need bessel-j or for the calculation of bessel-i we need bessel-j and bessel-k). So $besselarray can be destroyed when calculating bessel-i. It might be possible to prevent this effect by introducing a global flag which signals the internal call of one of the routines. A second point is, that we often dont't use the numerical slatec routines directly. In this case it is necessary to do a lot of extra calculations to obtain the values for the arrays. Perhaps it is useful to decide to fill the global arrays only for the case when we can use the values of the slatec routines directly. This is possible for positive order and argument or positiv order and complex argument of the Bessel functions. Further, an additional flag could be introduced to switch the filling of the arrays on and off. It is also to decide what we do when the user calculate a new value for a Bessel function and the calculation don't fill the global array. In this case, the user may be confused by the fact that the global arrays still hold old values. But for which order and argument? For this case it may be helpful to invalidate the array every time we calculate a new value or to store additionally the order and argument for which the values in the global array are calculated. All these problems with the global arrays arise because of the extension of the routines to negative arguments and orders. In my opinion it would be the best to use the concept of the global arrays no longer. The code becomes much more complicated and the benefit for the user might be small. Crategus ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-03-28 12:37 Message: Logged In: YES user_id=28849 Originator: NO First, thank you very much for the fixes to the Bessel routines. I'm sure we all want the routines to be correct! Second, about $besselarray. I see the documentation for that is gone, but I think the intent for besselarray was to let the user have access to all the computed values, since some algorithms end up computing all the intermediate orders to get to the desired order. We need to think if this should go away or not. Third, it would certainly help me a lot if you actually produced a diff between your new code and the existing code. Right now, I have to go look at every single function in different places to figure out what's changed. Trying to minimize the changes would help. I know this is a lot of work for you too. Just identifying which function you actually changed would help a lot too. About the specific changes you mention in the comments. No problem with changing bessel-x-simp to simp-bessel-x. (I think most of the code uses simpbesselx, but that's too messy for my tastes.) Extending the range is very, very good. Not sure about the special tests for pure real or imaginary answers. Need to look at that some more. Definitions for Hankel functions are good. May want to name them bessel_hankel_1 to be consistent with other Bessel functions. Certainly want to add some properties like derivatives, but that can wait. All in all, I would very much like to take in your very nice changes. ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-27 12:50 Message: Logged In: YES user_id=2039760 Originator: YES I have tested the Bessel functions with the numerical data attached to this message. For the tests I used a function TEST_BESSEL(value, result, digits) which allow to compare the value with the result within the specified digits. Crategus File Added: mytest_bessel.mac ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-27 12: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 13: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 |
From: SourceForge.net <no...@so...> - 2008-04-01 21:02:41
|
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-04-01 23:02 Message: Logged In: YES user_id=2039760 Originator: YES Again, thank you very much for your interest. I have not changed any of the derivatives. I have only changed the routines which I have listed in this thread. But I have not used the last CVS file. I have used the file which was send with the Maxima 5.14.0 release. There was one change for the derviative for bessel_k since the last release. On Sunday I have installed a tool to download the CVS files. Now I can get the really actual files. The next time I would like to optimize the routines further. Additionally, I am preparing extensive numerical tests for the functions. At the moment I have no idea how to extent the calculation to complex order. In my textbooks and in the book of Abramowitz and Stegun I can't find a relationship which help to calculate the Bessel functions for complex order. Crategus ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-01 18:41 Message: Logged In: YES user_id=28849 Originator: NO Thank you very much for the updated files. (The diff should have used diff -u or diff -c, but no matter. The bessel-changed.lisp file is perfect.) Given the imminent release of Maxima, I will wait on these changes until after the release. Then I'll try to incorporate your changes. I noticed that you changed the derivatives of some of the Bessel functions. Was that intentional? Were they wrong? (Some of the derivatives need to be in a certain form for some simplifications to be done.) ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-29 22:05 Message: Logged In: YES user_id=2039760 Originator: YES Next I have added a diff between the orginal and the changed file. Question: Is it possible to attach at once more than one file to a message? Crategus File Added: diff.txt ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-29 22:00 Message: Logged In: YES user_id=2039760 Originator: YES I have put all changes of the code for the Bessel functions in the orginal file bessel.lisp which is distributed with Maxima 5.14.0 and attached this file to this message. I have tested the code with the testsuite and the values of the mytest_bessel.mac. Crategus File Added: bessel-changed.lisp ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-29 01:35 Message: Logged In: YES user_id=2039760 Originator: YES Thank you very much for your answer. I have redesigned the following routines: simp-bessel-j (old name bessel-j-simp) bessel-j (old name $bessel) simp-bessel-y (old name bessel-y-simp) bessel-y simp-bessel-i (old name bessel-i-simp) bessel-i simp-bessel-k (old name bessel-k-simp) bessel-k The new routines are: $hankel_1 simp-hankel-1 $hankel_2 simp-hankel-2 If have further introduced the nouns %hankel_1 und %hankel_2. I would prefer to call the functions $hankel_1 and $hankel_2 and not bessel_hankel. The reason is that these functions are well known as Hankel functions. Because I have collected the routines in a new file it would be a good idea to do the changes in the orginal file bessel.lisp. Than it easier to produce a diff for you. I do this the next time. But first, I would like to have a second look at the global arrays $besselarray, $yarray and $iarray. I think it is not so easy to manage these global arrays in a consistent and predictable way for the user. One reason is, that the code for the numerical calculations use the other functions too (i.e. to calculate bessel-y we need bessel-j or for the calculation of bessel-i we need bessel-j and bessel-k). So $besselarray can be destroyed when calculating bessel-i. It might be possible to prevent this effect by introducing a global flag which signals the internal call of one of the routines. A second point is, that we often dont't use the numerical slatec routines directly. In this case it is necessary to do a lot of extra calculations to obtain the values for the arrays. Perhaps it is useful to decide to fill the global arrays only for the case when we can use the values of the slatec routines directly. This is possible for positive order and argument or positiv order and complex argument of the Bessel functions. Further, an additional flag could be introduced to switch the filling of the arrays on and off. It is also to decide what we do when the user calculate a new value for a Bessel function and the calculation don't fill the global array. In this case, the user may be confused by the fact that the global arrays still hold old values. But for which order and argument? For this case it may be helpful to invalidate the array every time we calculate a new value or to store additionally the order and argument for which the values in the global array are calculated. All these problems with the global arrays arise because of the extension of the routines to negative arguments and orders. In my opinion it would be the best to use the concept of the global arrays no longer. The code becomes much more complicated and the benefit for the user might be small. Crategus ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-03-28 17:37 Message: Logged In: YES user_id=28849 Originator: NO First, thank you very much for the fixes to the Bessel routines. I'm sure we all want the routines to be correct! Second, about $besselarray. I see the documentation for that is gone, but I think the intent for besselarray was to let the user have access to all the computed values, since some algorithms end up computing all the intermediate orders to get to the desired order. We need to think if this should go away or not. Third, it would certainly help me a lot if you actually produced a diff between your new code and the existing code. Right now, I have to go look at every single function in different places to figure out what's changed. Trying to minimize the changes would help. I know this is a lot of work for you too. Just identifying which function you actually changed would help a lot too. About the specific changes you mention in the comments. No problem with changing bessel-x-simp to simp-bessel-x. (I think most of the code uses simpbesselx, but that's too messy for my tastes.) Extending the range is very, very good. Not sure about the special tests for pure real or imaginary answers. Need to look at that some more. Definitions for Hankel functions are good. May want to name them bessel_hankel_1 to be consistent with other Bessel functions. Certainly want to add some properties like derivatives, but that can wait. All in all, I would very much like to take in your very nice changes. ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-27 17:50 Message: Logged In: YES user_id=2039760 Originator: YES I have tested the Bessel functions with the numerical data attached to this message. For the tests I used a function TEST_BESSEL(value, result, digits) which allow to compare the value with the result within the specified digits. Crategus File Added: mytest_bessel.mac ---------------------------------------------------------------------- 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 |
From: SourceForge.net <no...@so...> - 2008-04-10 16:57:52
|
Bugs item #1920177, was opened at 2008-03-19 17:07 Message generated for change (Comment added) made by rtoy 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: Raymond Toy (rtoy) Date: 2008-04-10 12:57 Message: Logged In: YES user_id=28849 Originator: NO I have now integrated your changes into bessel.lisp. The testsuite passes and your bessel tests pass too. I changed test_bessel(a,b,d) to mean abs(a-b)<10^(-d), which I hope is close enough. (Perhaps the real and imaginary parts should be tested separately?) I'll check in these changes after some code clean up. ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-04-01 17:02 Message: Logged In: YES user_id=2039760 Originator: YES Again, thank you very much for your interest. I have not changed any of the derivatives. I have only changed the routines which I have listed in this thread. But I have not used the last CVS file. I have used the file which was send with the Maxima 5.14.0 release. There was one change for the derviative for bessel_k since the last release. On Sunday I have installed a tool to download the CVS files. Now I can get the really actual files. The next time I would like to optimize the routines further. Additionally, I am preparing extensive numerical tests for the functions. At the moment I have no idea how to extent the calculation to complex order. In my textbooks and in the book of Abramowitz and Stegun I can't find a relationship which help to calculate the Bessel functions for complex order. Crategus ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-01 12:41 Message: Logged In: YES user_id=28849 Originator: NO Thank you very much for the updated files. (The diff should have used diff -u or diff -c, but no matter. The bessel-changed.lisp file is perfect.) Given the imminent release of Maxima, I will wait on these changes until after the release. Then I'll try to incorporate your changes. I noticed that you changed the derivatives of some of the Bessel functions. Was that intentional? Were they wrong? (Some of the derivatives need to be in a certain form for some simplifications to be done.) ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-29 17:05 Message: Logged In: YES user_id=2039760 Originator: YES Next I have added a diff between the orginal and the changed file. Question: Is it possible to attach at once more than one file to a message? Crategus File Added: diff.txt ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-29 17:00 Message: Logged In: YES user_id=2039760 Originator: YES I have put all changes of the code for the Bessel functions in the orginal file bessel.lisp which is distributed with Maxima 5.14.0 and attached this file to this message. I have tested the code with the testsuite and the values of the mytest_bessel.mac. Crategus File Added: bessel-changed.lisp ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-28 20:35 Message: Logged In: YES user_id=2039760 Originator: YES Thank you very much for your answer. I have redesigned the following routines: simp-bessel-j (old name bessel-j-simp) bessel-j (old name $bessel) simp-bessel-y (old name bessel-y-simp) bessel-y simp-bessel-i (old name bessel-i-simp) bessel-i simp-bessel-k (old name bessel-k-simp) bessel-k The new routines are: $hankel_1 simp-hankel-1 $hankel_2 simp-hankel-2 If have further introduced the nouns %hankel_1 und %hankel_2. I would prefer to call the functions $hankel_1 and $hankel_2 and not bessel_hankel. The reason is that these functions are well known as Hankel functions. Because I have collected the routines in a new file it would be a good idea to do the changes in the orginal file bessel.lisp. Than it easier to produce a diff for you. I do this the next time. But first, I would like to have a second look at the global arrays $besselarray, $yarray and $iarray. I think it is not so easy to manage these global arrays in a consistent and predictable way for the user. One reason is, that the code for the numerical calculations use the other functions too (i.e. to calculate bessel-y we need bessel-j or for the calculation of bessel-i we need bessel-j and bessel-k). So $besselarray can be destroyed when calculating bessel-i. It might be possible to prevent this effect by introducing a global flag which signals the internal call of one of the routines. A second point is, that we often dont't use the numerical slatec routines directly. In this case it is necessary to do a lot of extra calculations to obtain the values for the arrays. Perhaps it is useful to decide to fill the global arrays only for the case when we can use the values of the slatec routines directly. This is possible for positive order and argument or positiv order and complex argument of the Bessel functions. Further, an additional flag could be introduced to switch the filling of the arrays on and off. It is also to decide what we do when the user calculate a new value for a Bessel function and the calculation don't fill the global array. In this case, the user may be confused by the fact that the global arrays still hold old values. But for which order and argument? For this case it may be helpful to invalidate the array every time we calculate a new value or to store additionally the order and argument for which the values in the global array are calculated. All these problems with the global arrays arise because of the extension of the routines to negative arguments and orders. In my opinion it would be the best to use the concept of the global arrays no longer. The code becomes much more complicated and the benefit for the user might be small. Crategus ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-03-28 12:37 Message: Logged In: YES user_id=28849 Originator: NO First, thank you very much for the fixes to the Bessel routines. I'm sure we all want the routines to be correct! Second, about $besselarray. I see the documentation for that is gone, but I think the intent for besselarray was to let the user have access to all the computed values, since some algorithms end up computing all the intermediate orders to get to the desired order. We need to think if this should go away or not. Third, it would certainly help me a lot if you actually produced a diff between your new code and the existing code. Right now, I have to go look at every single function in different places to figure out what's changed. Trying to minimize the changes would help. I know this is a lot of work for you too. Just identifying which function you actually changed would help a lot too. About the specific changes you mention in the comments. No problem with changing bessel-x-simp to simp-bessel-x. (I think most of the code uses simpbesselx, but that's too messy for my tastes.) Extending the range is very, very good. Not sure about the special tests for pure real or imaginary answers. Need to look at that some more. Definitions for Hankel functions are good. May want to name them bessel_hankel_1 to be consistent with other Bessel functions. Certainly want to add some properties like derivatives, but that can wait. All in all, I would very much like to take in your very nice changes. ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-27 12:50 Message: Logged In: YES user_id=2039760 Originator: YES I have tested the Bessel functions with the numerical data attached to this message. For the tests I used a function TEST_BESSEL(value, result, digits) which allow to compare the value with the result within the specified digits. Crategus File Added: mytest_bessel.mac ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-27 12: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 13: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 |
From: SourceForge.net <no...@so...> - 2008-04-11 14:17:51
|
Bugs item #1920177, was opened at 2008-03-19 17:07 Message generated for change (Comment added) made by rtoy 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: Raymond Toy (rtoy) Date: 2008-04-11 10:17 Message: Logged In: YES user_id=28849 Originator: NO Changes checked in to bessel.lisp, rev 1.50. ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-10 12:57 Message: Logged In: YES user_id=28849 Originator: NO I have now integrated your changes into bessel.lisp. The testsuite passes and your bessel tests pass too. I changed test_bessel(a,b,d) to mean abs(a-b)<10^(-d), which I hope is close enough. (Perhaps the real and imaginary parts should be tested separately?) I'll check in these changes after some code clean up. ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-04-01 17:02 Message: Logged In: YES user_id=2039760 Originator: YES Again, thank you very much for your interest. I have not changed any of the derivatives. I have only changed the routines which I have listed in this thread. But I have not used the last CVS file. I have used the file which was send with the Maxima 5.14.0 release. There was one change for the derviative for bessel_k since the last release. On Sunday I have installed a tool to download the CVS files. Now I can get the really actual files. The next time I would like to optimize the routines further. Additionally, I am preparing extensive numerical tests for the functions. At the moment I have no idea how to extent the calculation to complex order. In my textbooks and in the book of Abramowitz and Stegun I can't find a relationship which help to calculate the Bessel functions for complex order. Crategus ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-01 12:41 Message: Logged In: YES user_id=28849 Originator: NO Thank you very much for the updated files. (The diff should have used diff -u or diff -c, but no matter. The bessel-changed.lisp file is perfect.) Given the imminent release of Maxima, I will wait on these changes until after the release. Then I'll try to incorporate your changes. I noticed that you changed the derivatives of some of the Bessel functions. Was that intentional? Were they wrong? (Some of the derivatives need to be in a certain form for some simplifications to be done.) ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-29 17:05 Message: Logged In: YES user_id=2039760 Originator: YES Next I have added a diff between the orginal and the changed file. Question: Is it possible to attach at once more than one file to a message? Crategus File Added: diff.txt ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-29 17:00 Message: Logged In: YES user_id=2039760 Originator: YES I have put all changes of the code for the Bessel functions in the orginal file bessel.lisp which is distributed with Maxima 5.14.0 and attached this file to this message. I have tested the code with the testsuite and the values of the mytest_bessel.mac. Crategus File Added: bessel-changed.lisp ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-28 20:35 Message: Logged In: YES user_id=2039760 Originator: YES Thank you very much for your answer. I have redesigned the following routines: simp-bessel-j (old name bessel-j-simp) bessel-j (old name $bessel) simp-bessel-y (old name bessel-y-simp) bessel-y simp-bessel-i (old name bessel-i-simp) bessel-i simp-bessel-k (old name bessel-k-simp) bessel-k The new routines are: $hankel_1 simp-hankel-1 $hankel_2 simp-hankel-2 If have further introduced the nouns %hankel_1 und %hankel_2. I would prefer to call the functions $hankel_1 and $hankel_2 and not bessel_hankel. The reason is that these functions are well known as Hankel functions. Because I have collected the routines in a new file it would be a good idea to do the changes in the orginal file bessel.lisp. Than it easier to produce a diff for you. I do this the next time. But first, I would like to have a second look at the global arrays $besselarray, $yarray and $iarray. I think it is not so easy to manage these global arrays in a consistent and predictable way for the user. One reason is, that the code for the numerical calculations use the other functions too (i.e. to calculate bessel-y we need bessel-j or for the calculation of bessel-i we need bessel-j and bessel-k). So $besselarray can be destroyed when calculating bessel-i. It might be possible to prevent this effect by introducing a global flag which signals the internal call of one of the routines. A second point is, that we often dont't use the numerical slatec routines directly. In this case it is necessary to do a lot of extra calculations to obtain the values for the arrays. Perhaps it is useful to decide to fill the global arrays only for the case when we can use the values of the slatec routines directly. This is possible for positive order and argument or positiv order and complex argument of the Bessel functions. Further, an additional flag could be introduced to switch the filling of the arrays on and off. It is also to decide what we do when the user calculate a new value for a Bessel function and the calculation don't fill the global array. In this case, the user may be confused by the fact that the global arrays still hold old values. But for which order and argument? For this case it may be helpful to invalidate the array every time we calculate a new value or to store additionally the order and argument for which the values in the global array are calculated. All these problems with the global arrays arise because of the extension of the routines to negative arguments and orders. In my opinion it would be the best to use the concept of the global arrays no longer. The code becomes much more complicated and the benefit for the user might be small. Crategus ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-03-28 12:37 Message: Logged In: YES user_id=28849 Originator: NO First, thank you very much for the fixes to the Bessel routines. I'm sure we all want the routines to be correct! Second, about $besselarray. I see the documentation for that is gone, but I think the intent for besselarray was to let the user have access to all the computed values, since some algorithms end up computing all the intermediate orders to get to the desired order. We need to think if this should go away or not. Third, it would certainly help me a lot if you actually produced a diff between your new code and the existing code. Right now, I have to go look at every single function in different places to figure out what's changed. Trying to minimize the changes would help. I know this is a lot of work for you too. Just identifying which function you actually changed would help a lot too. About the specific changes you mention in the comments. No problem with changing bessel-x-simp to simp-bessel-x. (I think most of the code uses simpbesselx, but that's too messy for my tastes.) Extending the range is very, very good. Not sure about the special tests for pure real or imaginary answers. Need to look at that some more. Definitions for Hankel functions are good. May want to name them bessel_hankel_1 to be consistent with other Bessel functions. Certainly want to add some properties like derivatives, but that can wait. All in all, I would very much like to take in your very nice changes. ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-27 12:50 Message: Logged In: YES user_id=2039760 Originator: YES I have tested the Bessel functions with the numerical data attached to this message. For the tests I used a function TEST_BESSEL(value, result, digits) which allow to compare the value with the result within the specified digits. Crategus File Added: mytest_bessel.mac ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-27 12: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 13: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 |
From: SourceForge.net <no...@so...> - 2008-04-12 21:15:44
|
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-04-12 23:15 Message: Logged In: YES user_id=2039760 Originator: YES Thank you very much for your work and your interest in the changes of the code. I have seen that you have deleted the use of the global $besselarray too. I have done further numerical calculations to test the code. I think a good test to show that the different parts of the code work well is to calculate the Wronskians which are defined in A&S. I have attached a file which does the calculations in test mode. (Excuse my programing style, but I have no experience programing in Maxima. Excuse my Englisch too.) As a side effect of the extension of the Bessel functions it is now possible to get nice plots to show the symmetry of the Bessel functions. A hint: You can kill the comment in German I had added at line 465. The question asks for the sign. The sign is correct! A second hint: There a still some parts of the code which don't use the symmetries of the Bessel functions completly. I will work further on this subject. At last: The documentation has to be updated too. I would like to do this work, but my English is not good enough. Crategus File Added: test_bessel_wronskian.mac ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-11 16:17 Message: Logged In: YES user_id=28849 Originator: NO Changes checked in to bessel.lisp, rev 1.50. ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-10 18:57 Message: Logged In: YES user_id=28849 Originator: NO I have now integrated your changes into bessel.lisp. The testsuite passes and your bessel tests pass too. I changed test_bessel(a,b,d) to mean abs(a-b)<10^(-d), which I hope is close enough. (Perhaps the real and imaginary parts should be tested separately?) I'll check in these changes after some code clean up. ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-04-01 23:02 Message: Logged In: YES user_id=2039760 Originator: YES Again, thank you very much for your interest. I have not changed any of the derivatives. I have only changed the routines which I have listed in this thread. But I have not used the last CVS file. I have used the file which was send with the Maxima 5.14.0 release. There was one change for the derviative for bessel_k since the last release. On Sunday I have installed a tool to download the CVS files. Now I can get the really actual files. The next time I would like to optimize the routines further. Additionally, I am preparing extensive numerical tests for the functions. At the moment I have no idea how to extent the calculation to complex order. In my textbooks and in the book of Abramowitz and Stegun I can't find a relationship which help to calculate the Bessel functions for complex order. Crategus ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-01 18:41 Message: Logged In: YES user_id=28849 Originator: NO Thank you very much for the updated files. (The diff should have used diff -u or diff -c, but no matter. The bessel-changed.lisp file is perfect.) Given the imminent release of Maxima, I will wait on these changes until after the release. Then I'll try to incorporate your changes. I noticed that you changed the derivatives of some of the Bessel functions. Was that intentional? Were they wrong? (Some of the derivatives need to be in a certain form for some simplifications to be done.) ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-29 22:05 Message: Logged In: YES user_id=2039760 Originator: YES Next I have added a diff between the orginal and the changed file. Question: Is it possible to attach at once more than one file to a message? Crategus File Added: diff.txt ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-29 22:00 Message: Logged In: YES user_id=2039760 Originator: YES I have put all changes of the code for the Bessel functions in the orginal file bessel.lisp which is distributed with Maxima 5.14.0 and attached this file to this message. I have tested the code with the testsuite and the values of the mytest_bessel.mac. Crategus File Added: bessel-changed.lisp ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-29 01:35 Message: Logged In: YES user_id=2039760 Originator: YES Thank you very much for your answer. I have redesigned the following routines: simp-bessel-j (old name bessel-j-simp) bessel-j (old name $bessel) simp-bessel-y (old name bessel-y-simp) bessel-y simp-bessel-i (old name bessel-i-simp) bessel-i simp-bessel-k (old name bessel-k-simp) bessel-k The new routines are: $hankel_1 simp-hankel-1 $hankel_2 simp-hankel-2 If have further introduced the nouns %hankel_1 und %hankel_2. I would prefer to call the functions $hankel_1 and $hankel_2 and not bessel_hankel. The reason is that these functions are well known as Hankel functions. Because I have collected the routines in a new file it would be a good idea to do the changes in the orginal file bessel.lisp. Than it easier to produce a diff for you. I do this the next time. But first, I would like to have a second look at the global arrays $besselarray, $yarray and $iarray. I think it is not so easy to manage these global arrays in a consistent and predictable way for the user. One reason is, that the code for the numerical calculations use the other functions too (i.e. to calculate bessel-y we need bessel-j or for the calculation of bessel-i we need bessel-j and bessel-k). So $besselarray can be destroyed when calculating bessel-i. It might be possible to prevent this effect by introducing a global flag which signals the internal call of one of the routines. A second point is, that we often dont't use the numerical slatec routines directly. In this case it is necessary to do a lot of extra calculations to obtain the values for the arrays. Perhaps it is useful to decide to fill the global arrays only for the case when we can use the values of the slatec routines directly. This is possible for positive order and argument or positiv order and complex argument of the Bessel functions. Further, an additional flag could be introduced to switch the filling of the arrays on and off. It is also to decide what we do when the user calculate a new value for a Bessel function and the calculation don't fill the global array. In this case, the user may be confused by the fact that the global arrays still hold old values. But for which order and argument? For this case it may be helpful to invalidate the array every time we calculate a new value or to store additionally the order and argument for which the values in the global array are calculated. All these problems with the global arrays arise because of the extension of the routines to negative arguments and orders. In my opinion it would be the best to use the concept of the global arrays no longer. The code becomes much more complicated and the benefit for the user might be small. Crategus ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-03-28 17:37 Message: Logged In: YES user_id=28849 Originator: NO First, thank you very much for the fixes to the Bessel routines. I'm sure we all want the routines to be correct! Second, about $besselarray. I see the documentation for that is gone, but I think the intent for besselarray was to let the user have access to all the computed values, since some algorithms end up computing all the intermediate orders to get to the desired order. We need to think if this should go away or not. Third, it would certainly help me a lot if you actually produced a diff between your new code and the existing code. Right now, I have to go look at every single function in different places to figure out what's changed. Trying to minimize the changes would help. I know this is a lot of work for you too. Just identifying which function you actually changed would help a lot too. About the specific changes you mention in the comments. No problem with changing bessel-x-simp to simp-bessel-x. (I think most of the code uses simpbesselx, but that's too messy for my tastes.) Extending the range is very, very good. Not sure about the special tests for pure real or imaginary answers. Need to look at that some more. Definitions for Hankel functions are good. May want to name them bessel_hankel_1 to be consistent with other Bessel functions. Certainly want to add some properties like derivatives, but that can wait. All in all, I would very much like to take in your very nice changes. ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-27 17:50 Message: Logged In: YES user_id=2039760 Originator: YES I have tested the Bessel functions with the numerical data attached to this message. For the tests I used a function TEST_BESSEL(value, result, digits) which allow to compare the value with the result within the specified digits. Crategus File Added: mytest_bessel.mac ---------------------------------------------------------------------- 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 |
From: SourceForge.net <no...@so...> - 2008-04-15 16:42:49
|
Bugs item #1920177, was opened at 2008-03-19 17:07 Message generated for change (Comment added) made by rtoy 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: Raymond Toy (rtoy) Date: 2008-04-15 12:42 Message: Logged In: YES user_id=28849 Originator: NO Thank you for the fixes! I'll fix the comment in German and I'll update the documentation too. I'll add your new tests as well. Once these are done, I'll close this bug report as fixed. If you find more issues, please open a new bug report. Thanks! ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-04-12 17:15 Message: Logged In: YES user_id=2039760 Originator: YES Thank you very much for your work and your interest in the changes of the code. I have seen that you have deleted the use of the global $besselarray too. I have done further numerical calculations to test the code. I think a good test to show that the different parts of the code work well is to calculate the Wronskians which are defined in A&S. I have attached a file which does the calculations in test mode. (Excuse my programing style, but I have no experience programing in Maxima. Excuse my Englisch too.) As a side effect of the extension of the Bessel functions it is now possible to get nice plots to show the symmetry of the Bessel functions. A hint: You can kill the comment in German I had added at line 465. The question asks for the sign. The sign is correct! A second hint: There a still some parts of the code which don't use the symmetries of the Bessel functions completly. I will work further on this subject. At last: The documentation has to be updated too. I would like to do this work, but my English is not good enough. Crategus File Added: test_bessel_wronskian.mac ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-11 10:17 Message: Logged In: YES user_id=28849 Originator: NO Changes checked in to bessel.lisp, rev 1.50. ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-10 12:57 Message: Logged In: YES user_id=28849 Originator: NO I have now integrated your changes into bessel.lisp. The testsuite passes and your bessel tests pass too. I changed test_bessel(a,b,d) to mean abs(a-b)<10^(-d), which I hope is close enough. (Perhaps the real and imaginary parts should be tested separately?) I'll check in these changes after some code clean up. ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-04-01 17:02 Message: Logged In: YES user_id=2039760 Originator: YES Again, thank you very much for your interest. I have not changed any of the derivatives. I have only changed the routines which I have listed in this thread. But I have not used the last CVS file. I have used the file which was send with the Maxima 5.14.0 release. There was one change for the derviative for bessel_k since the last release. On Sunday I have installed a tool to download the CVS files. Now I can get the really actual files. The next time I would like to optimize the routines further. Additionally, I am preparing extensive numerical tests for the functions. At the moment I have no idea how to extent the calculation to complex order. In my textbooks and in the book of Abramowitz and Stegun I can't find a relationship which help to calculate the Bessel functions for complex order. Crategus ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-01 12:41 Message: Logged In: YES user_id=28849 Originator: NO Thank you very much for the updated files. (The diff should have used diff -u or diff -c, but no matter. The bessel-changed.lisp file is perfect.) Given the imminent release of Maxima, I will wait on these changes until after the release. Then I'll try to incorporate your changes. I noticed that you changed the derivatives of some of the Bessel functions. Was that intentional? Were they wrong? (Some of the derivatives need to be in a certain form for some simplifications to be done.) ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-29 17:05 Message: Logged In: YES user_id=2039760 Originator: YES Next I have added a diff between the orginal and the changed file. Question: Is it possible to attach at once more than one file to a message? Crategus File Added: diff.txt ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-29 17:00 Message: Logged In: YES user_id=2039760 Originator: YES I have put all changes of the code for the Bessel functions in the orginal file bessel.lisp which is distributed with Maxima 5.14.0 and attached this file to this message. I have tested the code with the testsuite and the values of the mytest_bessel.mac. Crategus File Added: bessel-changed.lisp ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-28 20:35 Message: Logged In: YES user_id=2039760 Originator: YES Thank you very much for your answer. I have redesigned the following routines: simp-bessel-j (old name bessel-j-simp) bessel-j (old name $bessel) simp-bessel-y (old name bessel-y-simp) bessel-y simp-bessel-i (old name bessel-i-simp) bessel-i simp-bessel-k (old name bessel-k-simp) bessel-k The new routines are: $hankel_1 simp-hankel-1 $hankel_2 simp-hankel-2 If have further introduced the nouns %hankel_1 und %hankel_2. I would prefer to call the functions $hankel_1 and $hankel_2 and not bessel_hankel. The reason is that these functions are well known as Hankel functions. Because I have collected the routines in a new file it would be a good idea to do the changes in the orginal file bessel.lisp. Than it easier to produce a diff for you. I do this the next time. But first, I would like to have a second look at the global arrays $besselarray, $yarray and $iarray. I think it is not so easy to manage these global arrays in a consistent and predictable way for the user. One reason is, that the code for the numerical calculations use the other functions too (i.e. to calculate bessel-y we need bessel-j or for the calculation of bessel-i we need bessel-j and bessel-k). So $besselarray can be destroyed when calculating bessel-i. It might be possible to prevent this effect by introducing a global flag which signals the internal call of one of the routines. A second point is, that we often dont't use the numerical slatec routines directly. In this case it is necessary to do a lot of extra calculations to obtain the values for the arrays. Perhaps it is useful to decide to fill the global arrays only for the case when we can use the values of the slatec routines directly. This is possible for positive order and argument or positiv order and complex argument of the Bessel functions. Further, an additional flag could be introduced to switch the filling of the arrays on and off. It is also to decide what we do when the user calculate a new value for a Bessel function and the calculation don't fill the global array. In this case, the user may be confused by the fact that the global arrays still hold old values. But for which order and argument? For this case it may be helpful to invalidate the array every time we calculate a new value or to store additionally the order and argument for which the values in the global array are calculated. All these problems with the global arrays arise because of the extension of the routines to negative arguments and orders. In my opinion it would be the best to use the concept of the global arrays no longer. The code becomes much more complicated and the benefit for the user might be small. Crategus ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-03-28 12:37 Message: Logged In: YES user_id=28849 Originator: NO First, thank you very much for the fixes to the Bessel routines. I'm sure we all want the routines to be correct! Second, about $besselarray. I see the documentation for that is gone, but I think the intent for besselarray was to let the user have access to all the computed values, since some algorithms end up computing all the intermediate orders to get to the desired order. We need to think if this should go away or not. Third, it would certainly help me a lot if you actually produced a diff between your new code and the existing code. Right now, I have to go look at every single function in different places to figure out what's changed. Trying to minimize the changes would help. I know this is a lot of work for you too. Just identifying which function you actually changed would help a lot too. About the specific changes you mention in the comments. No problem with changing bessel-x-simp to simp-bessel-x. (I think most of the code uses simpbesselx, but that's too messy for my tastes.) Extending the range is very, very good. Not sure about the special tests for pure real or imaginary answers. Need to look at that some more. Definitions for Hankel functions are good. May want to name them bessel_hankel_1 to be consistent with other Bessel functions. Certainly want to add some properties like derivatives, but that can wait. All in all, I would very much like to take in your very nice changes. ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-27 12:50 Message: Logged In: YES user_id=2039760 Originator: YES I have tested the Bessel functions with the numerical data attached to this message. For the tests I used a function TEST_BESSEL(value, result, digits) which allow to compare the value with the result within the specified digits. Crategus File Added: mytest_bessel.mac ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-27 12: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 13: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 |
From: SourceForge.net <no...@so...> - 2008-04-15 16:42:56
|
Bugs item #1920177, was opened at 2008-03-19 17:07 Message generated for change (Comment added) made by rtoy 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: Raymond Toy (rtoy) Date: 2008-04-15 12:42 Message: Logged In: YES user_id=28849 Originator: NO Thank you for the fixes! I'll fix the comment in German and I'll update the documentation too. I'll add your new tests as well. Once these are done, I'll close this bug report as fixed. If you find more issues, please open a new bug report. Thanks! ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-15 12:42 Message: Logged In: YES user_id=28849 Originator: NO Thank you for the fixes! I'll fix the comment in German and I'll update the documentation too. I'll add your new tests as well. Once these are done, I'll close this bug report as fixed. If you find more issues, please open a new bug report. Thanks! ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-04-12 17:15 Message: Logged In: YES user_id=2039760 Originator: YES Thank you very much for your work and your interest in the changes of the code. I have seen that you have deleted the use of the global $besselarray too. I have done further numerical calculations to test the code. I think a good test to show that the different parts of the code work well is to calculate the Wronskians which are defined in A&S. I have attached a file which does the calculations in test mode. (Excuse my programing style, but I have no experience programing in Maxima. Excuse my Englisch too.) As a side effect of the extension of the Bessel functions it is now possible to get nice plots to show the symmetry of the Bessel functions. A hint: You can kill the comment in German I had added at line 465. The question asks for the sign. The sign is correct! A second hint: There a still some parts of the code which don't use the symmetries of the Bessel functions completly. I will work further on this subject. At last: The documentation has to be updated too. I would like to do this work, but my English is not good enough. Crategus File Added: test_bessel_wronskian.mac ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-11 10:17 Message: Logged In: YES user_id=28849 Originator: NO Changes checked in to bessel.lisp, rev 1.50. ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-10 12:57 Message: Logged In: YES user_id=28849 Originator: NO I have now integrated your changes into bessel.lisp. The testsuite passes and your bessel tests pass too. I changed test_bessel(a,b,d) to mean abs(a-b)<10^(-d), which I hope is close enough. (Perhaps the real and imaginary parts should be tested separately?) I'll check in these changes after some code clean up. ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-04-01 17:02 Message: Logged In: YES user_id=2039760 Originator: YES Again, thank you very much for your interest. I have not changed any of the derivatives. I have only changed the routines which I have listed in this thread. But I have not used the last CVS file. I have used the file which was send with the Maxima 5.14.0 release. There was one change for the derviative for bessel_k since the last release. On Sunday I have installed a tool to download the CVS files. Now I can get the really actual files. The next time I would like to optimize the routines further. Additionally, I am preparing extensive numerical tests for the functions. At the moment I have no idea how to extent the calculation to complex order. In my textbooks and in the book of Abramowitz and Stegun I can't find a relationship which help to calculate the Bessel functions for complex order. Crategus ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-01 12:41 Message: Logged In: YES user_id=28849 Originator: NO Thank you very much for the updated files. (The diff should have used diff -u or diff -c, but no matter. The bessel-changed.lisp file is perfect.) Given the imminent release of Maxima, I will wait on these changes until after the release. Then I'll try to incorporate your changes. I noticed that you changed the derivatives of some of the Bessel functions. Was that intentional? Were they wrong? (Some of the derivatives need to be in a certain form for some simplifications to be done.) ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-29 17:05 Message: Logged In: YES user_id=2039760 Originator: YES Next I have added a diff between the orginal and the changed file. Question: Is it possible to attach at once more than one file to a message? Crategus File Added: diff.txt ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-29 17:00 Message: Logged In: YES user_id=2039760 Originator: YES I have put all changes of the code for the Bessel functions in the orginal file bessel.lisp which is distributed with Maxima 5.14.0 and attached this file to this message. I have tested the code with the testsuite and the values of the mytest_bessel.mac. Crategus File Added: bessel-changed.lisp ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-28 20:35 Message: Logged In: YES user_id=2039760 Originator: YES Thank you very much for your answer. I have redesigned the following routines: simp-bessel-j (old name bessel-j-simp) bessel-j (old name $bessel) simp-bessel-y (old name bessel-y-simp) bessel-y simp-bessel-i (old name bessel-i-simp) bessel-i simp-bessel-k (old name bessel-k-simp) bessel-k The new routines are: $hankel_1 simp-hankel-1 $hankel_2 simp-hankel-2 If have further introduced the nouns %hankel_1 und %hankel_2. I would prefer to call the functions $hankel_1 and $hankel_2 and not bessel_hankel. The reason is that these functions are well known as Hankel functions. Because I have collected the routines in a new file it would be a good idea to do the changes in the orginal file bessel.lisp. Than it easier to produce a diff for you. I do this the next time. But first, I would like to have a second look at the global arrays $besselarray, $yarray and $iarray. I think it is not so easy to manage these global arrays in a consistent and predictable way for the user. One reason is, that the code for the numerical calculations use the other functions too (i.e. to calculate bessel-y we need bessel-j or for the calculation of bessel-i we need bessel-j and bessel-k). So $besselarray can be destroyed when calculating bessel-i. It might be possible to prevent this effect by introducing a global flag which signals the internal call of one of the routines. A second point is, that we often dont't use the numerical slatec routines directly. In this case it is necessary to do a lot of extra calculations to obtain the values for the arrays. Perhaps it is useful to decide to fill the global arrays only for the case when we can use the values of the slatec routines directly. This is possible for positive order and argument or positiv order and complex argument of the Bessel functions. Further, an additional flag could be introduced to switch the filling of the arrays on and off. It is also to decide what we do when the user calculate a new value for a Bessel function and the calculation don't fill the global array. In this case, the user may be confused by the fact that the global arrays still hold old values. But for which order and argument? For this case it may be helpful to invalidate the array every time we calculate a new value or to store additionally the order and argument for which the values in the global array are calculated. All these problems with the global arrays arise because of the extension of the routines to negative arguments and orders. In my opinion it would be the best to use the concept of the global arrays no longer. The code becomes much more complicated and the benefit for the user might be small. Crategus ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-03-28 12:37 Message: Logged In: YES user_id=28849 Originator: NO First, thank you very much for the fixes to the Bessel routines. I'm sure we all want the routines to be correct! Second, about $besselarray. I see the documentation for that is gone, but I think the intent for besselarray was to let the user have access to all the computed values, since some algorithms end up computing all the intermediate orders to get to the desired order. We need to think if this should go away or not. Third, it would certainly help me a lot if you actually produced a diff between your new code and the existing code. Right now, I have to go look at every single function in different places to figure out what's changed. Trying to minimize the changes would help. I know this is a lot of work for you too. Just identifying which function you actually changed would help a lot too. About the specific changes you mention in the comments. No problem with changing bessel-x-simp to simp-bessel-x. (I think most of the code uses simpbesselx, but that's too messy for my tastes.) Extending the range is very, very good. Not sure about the special tests for pure real or imaginary answers. Need to look at that some more. Definitions for Hankel functions are good. May want to name them bessel_hankel_1 to be consistent with other Bessel functions. Certainly want to add some properties like derivatives, but that can wait. All in all, I would very much like to take in your very nice changes. ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-27 12:50 Message: Logged In: YES user_id=2039760 Originator: YES I have tested the Bessel functions with the numerical data attached to this message. For the tests I used a function TEST_BESSEL(value, result, digits) which allow to compare the value with the result within the specified digits. Crategus File Added: mytest_bessel.mac ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-27 12: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 13: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 |
From: SourceForge.net <no...@so...> - 2008-04-15 16:43:13
|
Bugs item #1920177, was opened at 2008-03-19 17:07 Message generated for change (Comment added) made by rtoy 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: Raymond Toy (rtoy) Date: 2008-04-15 12:43 Message: Logged In: YES user_id=28849 Originator: NO Thank you for the fixes! I'll fix the comment in German and I'll update the documentation too. I'll add your new tests as well. Once these are done, I'll close this bug report as fixed. If you find more issues, please open a new bug report. Thanks! ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-15 12:42 Message: Logged In: YES user_id=28849 Originator: NO Thank you for the fixes! I'll fix the comment in German and I'll update the documentation too. I'll add your new tests as well. Once these are done, I'll close this bug report as fixed. If you find more issues, please open a new bug report. Thanks! ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-15 12:42 Message: Logged In: YES user_id=28849 Originator: NO Thank you for the fixes! I'll fix the comment in German and I'll update the documentation too. I'll add your new tests as well. Once these are done, I'll close this bug report as fixed. If you find more issues, please open a new bug report. Thanks! ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-04-12 17:15 Message: Logged In: YES user_id=2039760 Originator: YES Thank you very much for your work and your interest in the changes of the code. I have seen that you have deleted the use of the global $besselarray too. I have done further numerical calculations to test the code. I think a good test to show that the different parts of the code work well is to calculate the Wronskians which are defined in A&S. I have attached a file which does the calculations in test mode. (Excuse my programing style, but I have no experience programing in Maxima. Excuse my Englisch too.) As a side effect of the extension of the Bessel functions it is now possible to get nice plots to show the symmetry of the Bessel functions. A hint: You can kill the comment in German I had added at line 465. The question asks for the sign. The sign is correct! A second hint: There a still some parts of the code which don't use the symmetries of the Bessel functions completly. I will work further on this subject. At last: The documentation has to be updated too. I would like to do this work, but my English is not good enough. Crategus File Added: test_bessel_wronskian.mac ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-11 10:17 Message: Logged In: YES user_id=28849 Originator: NO Changes checked in to bessel.lisp, rev 1.50. ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-10 12:57 Message: Logged In: YES user_id=28849 Originator: NO I have now integrated your changes into bessel.lisp. The testsuite passes and your bessel tests pass too. I changed test_bessel(a,b,d) to mean abs(a-b)<10^(-d), which I hope is close enough. (Perhaps the real and imaginary parts should be tested separately?) I'll check in these changes after some code clean up. ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-04-01 17:02 Message: Logged In: YES user_id=2039760 Originator: YES Again, thank you very much for your interest. I have not changed any of the derivatives. I have only changed the routines which I have listed in this thread. But I have not used the last CVS file. I have used the file which was send with the Maxima 5.14.0 release. There was one change for the derviative for bessel_k since the last release. On Sunday I have installed a tool to download the CVS files. Now I can get the really actual files. The next time I would like to optimize the routines further. Additionally, I am preparing extensive numerical tests for the functions. At the moment I have no idea how to extent the calculation to complex order. In my textbooks and in the book of Abramowitz and Stegun I can't find a relationship which help to calculate the Bessel functions for complex order. Crategus ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-01 12:41 Message: Logged In: YES user_id=28849 Originator: NO Thank you very much for the updated files. (The diff should have used diff -u or diff -c, but no matter. The bessel-changed.lisp file is perfect.) Given the imminent release of Maxima, I will wait on these changes until after the release. Then I'll try to incorporate your changes. I noticed that you changed the derivatives of some of the Bessel functions. Was that intentional? Were they wrong? (Some of the derivatives need to be in a certain form for some simplifications to be done.) ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-29 17:05 Message: Logged In: YES user_id=2039760 Originator: YES Next I have added a diff between the orginal and the changed file. Question: Is it possible to attach at once more than one file to a message? Crategus File Added: diff.txt ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-29 17:00 Message: Logged In: YES user_id=2039760 Originator: YES I have put all changes of the code for the Bessel functions in the orginal file bessel.lisp which is distributed with Maxima 5.14.0 and attached this file to this message. I have tested the code with the testsuite and the values of the mytest_bessel.mac. Crategus File Added: bessel-changed.lisp ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-28 20:35 Message: Logged In: YES user_id=2039760 Originator: YES Thank you very much for your answer. I have redesigned the following routines: simp-bessel-j (old name bessel-j-simp) bessel-j (old name $bessel) simp-bessel-y (old name bessel-y-simp) bessel-y simp-bessel-i (old name bessel-i-simp) bessel-i simp-bessel-k (old name bessel-k-simp) bessel-k The new routines are: $hankel_1 simp-hankel-1 $hankel_2 simp-hankel-2 If have further introduced the nouns %hankel_1 und %hankel_2. I would prefer to call the functions $hankel_1 and $hankel_2 and not bessel_hankel. The reason is that these functions are well known as Hankel functions. Because I have collected the routines in a new file it would be a good idea to do the changes in the orginal file bessel.lisp. Than it easier to produce a diff for you. I do this the next time. But first, I would like to have a second look at the global arrays $besselarray, $yarray and $iarray. I think it is not so easy to manage these global arrays in a consistent and predictable way for the user. One reason is, that the code for the numerical calculations use the other functions too (i.e. to calculate bessel-y we need bessel-j or for the calculation of bessel-i we need bessel-j and bessel-k). So $besselarray can be destroyed when calculating bessel-i. It might be possible to prevent this effect by introducing a global flag which signals the internal call of one of the routines. A second point is, that we often dont't use the numerical slatec routines directly. In this case it is necessary to do a lot of extra calculations to obtain the values for the arrays. Perhaps it is useful to decide to fill the global arrays only for the case when we can use the values of the slatec routines directly. This is possible for positive order and argument or positiv order and complex argument of the Bessel functions. Further, an additional flag could be introduced to switch the filling of the arrays on and off. It is also to decide what we do when the user calculate a new value for a Bessel function and the calculation don't fill the global array. In this case, the user may be confused by the fact that the global arrays still hold old values. But for which order and argument? For this case it may be helpful to invalidate the array every time we calculate a new value or to store additionally the order and argument for which the values in the global array are calculated. All these problems with the global arrays arise because of the extension of the routines to negative arguments and orders. In my opinion it would be the best to use the concept of the global arrays no longer. The code becomes much more complicated and the benefit for the user might be small. Crategus ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-03-28 12:37 Message: Logged In: YES user_id=28849 Originator: NO First, thank you very much for the fixes to the Bessel routines. I'm sure we all want the routines to be correct! Second, about $besselarray. I see the documentation for that is gone, but I think the intent for besselarray was to let the user have access to all the computed values, since some algorithms end up computing all the intermediate orders to get to the desired order. We need to think if this should go away or not. Third, it would certainly help me a lot if you actually produced a diff between your new code and the existing code. Right now, I have to go look at every single function in different places to figure out what's changed. Trying to minimize the changes would help. I know this is a lot of work for you too. Just identifying which function you actually changed would help a lot too. About the specific changes you mention in the comments. No problem with changing bessel-x-simp to simp-bessel-x. (I think most of the code uses simpbesselx, but that's too messy for my tastes.) Extending the range is very, very good. Not sure about the special tests for pure real or imaginary answers. Need to look at that some more. Definitions for Hankel functions are good. May want to name them bessel_hankel_1 to be consistent with other Bessel functions. Certainly want to add some properties like derivatives, but that can wait. All in all, I would very much like to take in your very nice changes. ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-27 12:50 Message: Logged In: YES user_id=2039760 Originator: YES I have tested the Bessel functions with the numerical data attached to this message. For the tests I used a function TEST_BESSEL(value, result, digits) which allow to compare the value with the result within the specified digits. Crategus File Added: mytest_bessel.mac ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-27 12: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 13: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 |
From: SourceForge.net <no...@so...> - 2008-04-16 13:39:54
|
Bugs item #1920177, was opened at 2008-03-19 17:07 Message generated for change (Comment added) made by rtoy 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: Raymond Toy (rtoy) Date: 2008-04-16 09:39 Message: Logged In: YES user_id=28849 Originator: NO What changes should be made to the documentation? I've already removed the part about besselarray. Also, I've run the test_bessel_wronskian tests. These need to be modified to fit in with how regression tests are done. However, I get the following messages: W_II: MAX EPS 7.44e-5 for order = -0.50 and arg = -5.00 W_IK: MAX EPS 1.17e-4 for order = -0.50 and arg = -5.00 Are these expected? ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-15 12:43 Message: Logged In: YES user_id=28849 Originator: NO Thank you for the fixes! I'll fix the comment in German and I'll update the documentation too. I'll add your new tests as well. Once these are done, I'll close this bug report as fixed. If you find more issues, please open a new bug report. Thanks! ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-15 12:42 Message: Logged In: YES user_id=28849 Originator: NO Thank you for the fixes! I'll fix the comment in German and I'll update the documentation too. I'll add your new tests as well. Once these are done, I'll close this bug report as fixed. If you find more issues, please open a new bug report. Thanks! ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-15 12:42 Message: Logged In: YES user_id=28849 Originator: NO Thank you for the fixes! I'll fix the comment in German and I'll update the documentation too. I'll add your new tests as well. Once these are done, I'll close this bug report as fixed. If you find more issues, please open a new bug report. Thanks! ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-04-12 17:15 Message: Logged In: YES user_id=2039760 Originator: YES Thank you very much for your work and your interest in the changes of the code. I have seen that you have deleted the use of the global $besselarray too. I have done further numerical calculations to test the code. I think a good test to show that the different parts of the code work well is to calculate the Wronskians which are defined in A&S. I have attached a file which does the calculations in test mode. (Excuse my programing style, but I have no experience programing in Maxima. Excuse my Englisch too.) As a side effect of the extension of the Bessel functions it is now possible to get nice plots to show the symmetry of the Bessel functions. A hint: You can kill the comment in German I had added at line 465. The question asks for the sign. The sign is correct! A second hint: There a still some parts of the code which don't use the symmetries of the Bessel functions completly. I will work further on this subject. At last: The documentation has to be updated too. I would like to do this work, but my English is not good enough. Crategus File Added: test_bessel_wronskian.mac ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-11 10:17 Message: Logged In: YES user_id=28849 Originator: NO Changes checked in to bessel.lisp, rev 1.50. ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-10 12:57 Message: Logged In: YES user_id=28849 Originator: NO I have now integrated your changes into bessel.lisp. The testsuite passes and your bessel tests pass too. I changed test_bessel(a,b,d) to mean abs(a-b)<10^(-d), which I hope is close enough. (Perhaps the real and imaginary parts should be tested separately?) I'll check in these changes after some code clean up. ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-04-01 17:02 Message: Logged In: YES user_id=2039760 Originator: YES Again, thank you very much for your interest. I have not changed any of the derivatives. I have only changed the routines which I have listed in this thread. But I have not used the last CVS file. I have used the file which was send with the Maxima 5.14.0 release. There was one change for the derviative for bessel_k since the last release. On Sunday I have installed a tool to download the CVS files. Now I can get the really actual files. The next time I would like to optimize the routines further. Additionally, I am preparing extensive numerical tests for the functions. At the moment I have no idea how to extent the calculation to complex order. In my textbooks and in the book of Abramowitz and Stegun I can't find a relationship which help to calculate the Bessel functions for complex order. Crategus ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-01 12:41 Message: Logged In: YES user_id=28849 Originator: NO Thank you very much for the updated files. (The diff should have used diff -u or diff -c, but no matter. The bessel-changed.lisp file is perfect.) Given the imminent release of Maxima, I will wait on these changes until after the release. Then I'll try to incorporate your changes. I noticed that you changed the derivatives of some of the Bessel functions. Was that intentional? Were they wrong? (Some of the derivatives need to be in a certain form for some simplifications to be done.) ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-29 17:05 Message: Logged In: YES user_id=2039760 Originator: YES Next I have added a diff between the orginal and the changed file. Question: Is it possible to attach at once more than one file to a message? Crategus File Added: diff.txt ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-29 17:00 Message: Logged In: YES user_id=2039760 Originator: YES I have put all changes of the code for the Bessel functions in the orginal file bessel.lisp which is distributed with Maxima 5.14.0 and attached this file to this message. I have tested the code with the testsuite and the values of the mytest_bessel.mac. Crategus File Added: bessel-changed.lisp ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-28 20:35 Message: Logged In: YES user_id=2039760 Originator: YES Thank you very much for your answer. I have redesigned the following routines: simp-bessel-j (old name bessel-j-simp) bessel-j (old name $bessel) simp-bessel-y (old name bessel-y-simp) bessel-y simp-bessel-i (old name bessel-i-simp) bessel-i simp-bessel-k (old name bessel-k-simp) bessel-k The new routines are: $hankel_1 simp-hankel-1 $hankel_2 simp-hankel-2 If have further introduced the nouns %hankel_1 und %hankel_2. I would prefer to call the functions $hankel_1 and $hankel_2 and not bessel_hankel. The reason is that these functions are well known as Hankel functions. Because I have collected the routines in a new file it would be a good idea to do the changes in the orginal file bessel.lisp. Than it easier to produce a diff for you. I do this the next time. But first, I would like to have a second look at the global arrays $besselarray, $yarray and $iarray. I think it is not so easy to manage these global arrays in a consistent and predictable way for the user. One reason is, that the code for the numerical calculations use the other functions too (i.e. to calculate bessel-y we need bessel-j or for the calculation of bessel-i we need bessel-j and bessel-k). So $besselarray can be destroyed when calculating bessel-i. It might be possible to prevent this effect by introducing a global flag which signals the internal call of one of the routines. A second point is, that we often dont't use the numerical slatec routines directly. In this case it is necessary to do a lot of extra calculations to obtain the values for the arrays. Perhaps it is useful to decide to fill the global arrays only for the case when we can use the values of the slatec routines directly. This is possible for positive order and argument or positiv order and complex argument of the Bessel functions. Further, an additional flag could be introduced to switch the filling of the arrays on and off. It is also to decide what we do when the user calculate a new value for a Bessel function and the calculation don't fill the global array. In this case, the user may be confused by the fact that the global arrays still hold old values. But for which order and argument? For this case it may be helpful to invalidate the array every time we calculate a new value or to store additionally the order and argument for which the values in the global array are calculated. All these problems with the global arrays arise because of the extension of the routines to negative arguments and orders. In my opinion it would be the best to use the concept of the global arrays no longer. The code becomes much more complicated and the benefit for the user might be small. Crategus ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-03-28 12:37 Message: Logged In: YES user_id=28849 Originator: NO First, thank you very much for the fixes to the Bessel routines. I'm sure we all want the routines to be correct! Second, about $besselarray. I see the documentation for that is gone, but I think the intent for besselarray was to let the user have access to all the computed values, since some algorithms end up computing all the intermediate orders to get to the desired order. We need to think if this should go away or not. Third, it would certainly help me a lot if you actually produced a diff between your new code and the existing code. Right now, I have to go look at every single function in different places to figure out what's changed. Trying to minimize the changes would help. I know this is a lot of work for you too. Just identifying which function you actually changed would help a lot too. About the specific changes you mention in the comments. No problem with changing bessel-x-simp to simp-bessel-x. (I think most of the code uses simpbesselx, but that's too messy for my tastes.) Extending the range is very, very good. Not sure about the special tests for pure real or imaginary answers. Need to look at that some more. Definitions for Hankel functions are good. May want to name them bessel_hankel_1 to be consistent with other Bessel functions. Certainly want to add some properties like derivatives, but that can wait. All in all, I would very much like to take in your very nice changes. ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-27 12:50 Message: Logged In: YES user_id=2039760 Originator: YES I have tested the Bessel functions with the numerical data attached to this message. For the tests I used a function TEST_BESSEL(value, result, digits) which allow to compare the value with the result within the specified digits. Crategus File Added: mytest_bessel.mac ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-27 12: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 13: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 |
From: SourceForge.net <no...@so...> - 2008-04-17 20:14:16
|
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-04-17 22:14 Message: Logged In: YES user_id=2039760 Originator: YES First my last build info: Maxima version: 5.14.0cvs Maxima build date: 5:26 4/18/2008 host type: @host@ lisp-implementation-type: GNU Common Lisp (GCL) lisp-implementation-version: GCL 2.6.8 I have build Maxima from the new CVS files today. I run all tests and had no problems. For the tests W_II and W_IK I had adjusted eps to 1e-10. With this value for eps all tests pass in the range of -5.0 to +5.0 for the argument. On my system I get the following results for the tests with real arg (%i23) abs(w_ii(-0.50,-5.00)-(-2.0*sin(-0.5*%pi)/((-5.0)*%pi))),numer; (%o23) 9.4075686847027995E-14 (%i24) abs(w_ik(-0.50,-5.00)-(1/(-5.0))),numer; (%o24) 4.7815106425210348E-13 And the following results with complex argument: (%i28) abs(w_ii(-0.50,-5.00*%i)-(-2.0*sin(-0.5*%pi)/((-5.0*%i)*%pi))),numer; (%o28) 1.0825261750342018E-15 (%i29) abs(w_ik(-0.50,-5.00*%i)-(1/(-5.0*%i))),numer; (%o29) 1.6378141611741369E-15 So I get the values for the Wronskians with much higher accuracy. Crategus ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-16 15:39 Message: Logged In: YES user_id=28849 Originator: NO What changes should be made to the documentation? I've already removed the part about besselarray. Also, I've run the test_bessel_wronskian tests. These need to be modified to fit in with how regression tests are done. However, I get the following messages: W_II: MAX EPS 7.44e-5 for order = -0.50 and arg = -5.00 W_IK: MAX EPS 1.17e-4 for order = -0.50 and arg = -5.00 Are these expected? ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-15 18:43 Message: Logged In: YES user_id=28849 Originator: NO Thank you for the fixes! I'll fix the comment in German and I'll update the documentation too. I'll add your new tests as well. Once these are done, I'll close this bug report as fixed. If you find more issues, please open a new bug report. Thanks! ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-15 18:42 Message: Logged In: YES user_id=28849 Originator: NO Thank you for the fixes! I'll fix the comment in German and I'll update the documentation too. I'll add your new tests as well. Once these are done, I'll close this bug report as fixed. If you find more issues, please open a new bug report. Thanks! ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-15 18:42 Message: Logged In: YES user_id=28849 Originator: NO Thank you for the fixes! I'll fix the comment in German and I'll update the documentation too. I'll add your new tests as well. Once these are done, I'll close this bug report as fixed. If you find more issues, please open a new bug report. Thanks! ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-04-12 23:15 Message: Logged In: YES user_id=2039760 Originator: YES Thank you very much for your work and your interest in the changes of the code. I have seen that you have deleted the use of the global $besselarray too. I have done further numerical calculations to test the code. I think a good test to show that the different parts of the code work well is to calculate the Wronskians which are defined in A&S. I have attached a file which does the calculations in test mode. (Excuse my programing style, but I have no experience programing in Maxima. Excuse my Englisch too.) As a side effect of the extension of the Bessel functions it is now possible to get nice plots to show the symmetry of the Bessel functions. A hint: You can kill the comment in German I had added at line 465. The question asks for the sign. The sign is correct! A second hint: There a still some parts of the code which don't use the symmetries of the Bessel functions completly. I will work further on this subject. At last: The documentation has to be updated too. I would like to do this work, but my English is not good enough. Crategus File Added: test_bessel_wronskian.mac ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-11 16:17 Message: Logged In: YES user_id=28849 Originator: NO Changes checked in to bessel.lisp, rev 1.50. ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-10 18:57 Message: Logged In: YES user_id=28849 Originator: NO I have now integrated your changes into bessel.lisp. The testsuite passes and your bessel tests pass too. I changed test_bessel(a,b,d) to mean abs(a-b)<10^(-d), which I hope is close enough. (Perhaps the real and imaginary parts should be tested separately?) I'll check in these changes after some code clean up. ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-04-01 23:02 Message: Logged In: YES user_id=2039760 Originator: YES Again, thank you very much for your interest. I have not changed any of the derivatives. I have only changed the routines which I have listed in this thread. But I have not used the last CVS file. I have used the file which was send with the Maxima 5.14.0 release. There was one change for the derviative for bessel_k since the last release. On Sunday I have installed a tool to download the CVS files. Now I can get the really actual files. The next time I would like to optimize the routines further. Additionally, I am preparing extensive numerical tests for the functions. At the moment I have no idea how to extent the calculation to complex order. In my textbooks and in the book of Abramowitz and Stegun I can't find a relationship which help to calculate the Bessel functions for complex order. Crategus ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-01 18:41 Message: Logged In: YES user_id=28849 Originator: NO Thank you very much for the updated files. (The diff should have used diff -u or diff -c, but no matter. The bessel-changed.lisp file is perfect.) Given the imminent release of Maxima, I will wait on these changes until after the release. Then I'll try to incorporate your changes. I noticed that you changed the derivatives of some of the Bessel functions. Was that intentional? Were they wrong? (Some of the derivatives need to be in a certain form for some simplifications to be done.) ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-29 22:05 Message: Logged In: YES user_id=2039760 Originator: YES Next I have added a diff between the orginal and the changed file. Question: Is it possible to attach at once more than one file to a message? Crategus File Added: diff.txt ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-29 22:00 Message: Logged In: YES user_id=2039760 Originator: YES I have put all changes of the code for the Bessel functions in the orginal file bessel.lisp which is distributed with Maxima 5.14.0 and attached this file to this message. I have tested the code with the testsuite and the values of the mytest_bessel.mac. Crategus File Added: bessel-changed.lisp ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-29 01:35 Message: Logged In: YES user_id=2039760 Originator: YES Thank you very much for your answer. I have redesigned the following routines: simp-bessel-j (old name bessel-j-simp) bessel-j (old name $bessel) simp-bessel-y (old name bessel-y-simp) bessel-y simp-bessel-i (old name bessel-i-simp) bessel-i simp-bessel-k (old name bessel-k-simp) bessel-k The new routines are: $hankel_1 simp-hankel-1 $hankel_2 simp-hankel-2 If have further introduced the nouns %hankel_1 und %hankel_2. I would prefer to call the functions $hankel_1 and $hankel_2 and not bessel_hankel. The reason is that these functions are well known as Hankel functions. Because I have collected the routines in a new file it would be a good idea to do the changes in the orginal file bessel.lisp. Than it easier to produce a diff for you. I do this the next time. But first, I would like to have a second look at the global arrays $besselarray, $yarray and $iarray. I think it is not so easy to manage these global arrays in a consistent and predictable way for the user. One reason is, that the code for the numerical calculations use the other functions too (i.e. to calculate bessel-y we need bessel-j or for the calculation of bessel-i we need bessel-j and bessel-k). So $besselarray can be destroyed when calculating bessel-i. It might be possible to prevent this effect by introducing a global flag which signals the internal call of one of the routines. A second point is, that we often dont't use the numerical slatec routines directly. In this case it is necessary to do a lot of extra calculations to obtain the values for the arrays. Perhaps it is useful to decide to fill the global arrays only for the case when we can use the values of the slatec routines directly. This is possible for positive order and argument or positiv order and complex argument of the Bessel functions. Further, an additional flag could be introduced to switch the filling of the arrays on and off. It is also to decide what we do when the user calculate a new value for a Bessel function and the calculation don't fill the global array. In this case, the user may be confused by the fact that the global arrays still hold old values. But for which order and argument? For this case it may be helpful to invalidate the array every time we calculate a new value or to store additionally the order and argument for which the values in the global array are calculated. All these problems with the global arrays arise because of the extension of the routines to negative arguments and orders. In my opinion it would be the best to use the concept of the global arrays no longer. The code becomes much more complicated and the benefit for the user might be small. Crategus ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-03-28 17:37 Message: Logged In: YES user_id=28849 Originator: NO First, thank you very much for the fixes to the Bessel routines. I'm sure we all want the routines to be correct! Second, about $besselarray. I see the documentation for that is gone, but I think the intent for besselarray was to let the user have access to all the computed values, since some algorithms end up computing all the intermediate orders to get to the desired order. We need to think if this should go away or not. Third, it would certainly help me a lot if you actually produced a diff between your new code and the existing code. Right now, I have to go look at every single function in different places to figure out what's changed. Trying to minimize the changes would help. I know this is a lot of work for you too. Just identifying which function you actually changed would help a lot too. About the specific changes you mention in the comments. No problem with changing bessel-x-simp to simp-bessel-x. (I think most of the code uses simpbesselx, but that's too messy for my tastes.) Extending the range is very, very good. Not sure about the special tests for pure real or imaginary answers. Need to look at that some more. Definitions for Hankel functions are good. May want to name them bessel_hankel_1 to be consistent with other Bessel functions. Certainly want to add some properties like derivatives, but that can wait. All in all, I would very much like to take in your very nice changes. ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-27 17:50 Message: Logged In: YES user_id=2039760 Originator: YES I have tested the Bessel functions with the numerical data attached to this message. For the tests I used a function TEST_BESSEL(value, result, digits) which allow to compare the value with the result within the specified digits. Crategus File Added: mytest_bessel.mac ---------------------------------------------------------------------- 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 |
From: SourceForge.net <no...@so...> - 2008-04-18 02:04:08
|
Bugs item #1920177, was opened at 2008-03-19 17:07 Message generated for change (Comment added) made by rtoy 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: Raymond Toy (rtoy) Date: 2008-04-17 22:04 Message: Logged In: YES user_id=28849 Originator: NO Ok. Something is wrong. Here are some samples with CVS sources and cmucl on linux and sparc: besselexpand:true; bessel_i(-1/2,-5) -> -sqrt(2)*%i*cosh(5)/sqrt(5)/sqrt(%pi) -> -26.47995176430595*%i bessel_i(-0.5,-5.0) -> -26.47995216244259*%i-5.787377599375532e-7 But, curiously, gcl gives -26.47995174630617*%i+1.77635683940025e-15. Need to investigate why the difference... ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-04-17 16:14 Message: Logged In: YES user_id=2039760 Originator: YES First my last build info: Maxima version: 5.14.0cvs Maxima build date: 5:26 4/18/2008 host type: @host@ lisp-implementation-type: GNU Common Lisp (GCL) lisp-implementation-version: GCL 2.6.8 I have build Maxima from the new CVS files today. I run all tests and had no problems. For the tests W_II and W_IK I had adjusted eps to 1e-10. With this value for eps all tests pass in the range of -5.0 to +5.0 for the argument. On my system I get the following results for the tests with real arg (%i23) abs(w_ii(-0.50,-5.00)-(-2.0*sin(-0.5*%pi)/((-5.0)*%pi))),numer; (%o23) 9.4075686847027995E-14 (%i24) abs(w_ik(-0.50,-5.00)-(1/(-5.0))),numer; (%o24) 4.7815106425210348E-13 And the following results with complex argument: (%i28) abs(w_ii(-0.50,-5.00*%i)-(-2.0*sin(-0.5*%pi)/((-5.0*%i)*%pi))),numer; (%o28) 1.0825261750342018E-15 (%i29) abs(w_ik(-0.50,-5.00*%i)-(1/(-5.0*%i))),numer; (%o29) 1.6378141611741369E-15 So I get the values for the Wronskians with much higher accuracy. Crategus ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-16 09:39 Message: Logged In: YES user_id=28849 Originator: NO What changes should be made to the documentation? I've already removed the part about besselarray. Also, I've run the test_bessel_wronskian tests. These need to be modified to fit in with how regression tests are done. However, I get the following messages: W_II: MAX EPS 7.44e-5 for order = -0.50 and arg = -5.00 W_IK: MAX EPS 1.17e-4 for order = -0.50 and arg = -5.00 Are these expected? ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-15 12:43 Message: Logged In: YES user_id=28849 Originator: NO Thank you for the fixes! I'll fix the comment in German and I'll update the documentation too. I'll add your new tests as well. Once these are done, I'll close this bug report as fixed. If you find more issues, please open a new bug report. Thanks! ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-15 12:42 Message: Logged In: YES user_id=28849 Originator: NO Thank you for the fixes! I'll fix the comment in German and I'll update the documentation too. I'll add your new tests as well. Once these are done, I'll close this bug report as fixed. If you find more issues, please open a new bug report. Thanks! ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-15 12:42 Message: Logged In: YES user_id=28849 Originator: NO Thank you for the fixes! I'll fix the comment in German and I'll update the documentation too. I'll add your new tests as well. Once these are done, I'll close this bug report as fixed. If you find more issues, please open a new bug report. Thanks! ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-04-12 17:15 Message: Logged In: YES user_id=2039760 Originator: YES Thank you very much for your work and your interest in the changes of the code. I have seen that you have deleted the use of the global $besselarray too. I have done further numerical calculations to test the code. I think a good test to show that the different parts of the code work well is to calculate the Wronskians which are defined in A&S. I have attached a file which does the calculations in test mode. (Excuse my programing style, but I have no experience programing in Maxima. Excuse my Englisch too.) As a side effect of the extension of the Bessel functions it is now possible to get nice plots to show the symmetry of the Bessel functions. A hint: You can kill the comment in German I had added at line 465. The question asks for the sign. The sign is correct! A second hint: There a still some parts of the code which don't use the symmetries of the Bessel functions completly. I will work further on this subject. At last: The documentation has to be updated too. I would like to do this work, but my English is not good enough. Crategus File Added: test_bessel_wronskian.mac ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-11 10:17 Message: Logged In: YES user_id=28849 Originator: NO Changes checked in to bessel.lisp, rev 1.50. ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-10 12:57 Message: Logged In: YES user_id=28849 Originator: NO I have now integrated your changes into bessel.lisp. The testsuite passes and your bessel tests pass too. I changed test_bessel(a,b,d) to mean abs(a-b)<10^(-d), which I hope is close enough. (Perhaps the real and imaginary parts should be tested separately?) I'll check in these changes after some code clean up. ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-04-01 17:02 Message: Logged In: YES user_id=2039760 Originator: YES Again, thank you very much for your interest. I have not changed any of the derivatives. I have only changed the routines which I have listed in this thread. But I have not used the last CVS file. I have used the file which was send with the Maxima 5.14.0 release. There was one change for the derviative for bessel_k since the last release. On Sunday I have installed a tool to download the CVS files. Now I can get the really actual files. The next time I would like to optimize the routines further. Additionally, I am preparing extensive numerical tests for the functions. At the moment I have no idea how to extent the calculation to complex order. In my textbooks and in the book of Abramowitz and Stegun I can't find a relationship which help to calculate the Bessel functions for complex order. Crategus ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-01 12:41 Message: Logged In: YES user_id=28849 Originator: NO Thank you very much for the updated files. (The diff should have used diff -u or diff -c, but no matter. The bessel-changed.lisp file is perfect.) Given the imminent release of Maxima, I will wait on these changes until after the release. Then I'll try to incorporate your changes. I noticed that you changed the derivatives of some of the Bessel functions. Was that intentional? Were they wrong? (Some of the derivatives need to be in a certain form for some simplifications to be done.) ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-29 17:05 Message: Logged In: YES user_id=2039760 Originator: YES Next I have added a diff between the orginal and the changed file. Question: Is it possible to attach at once more than one file to a message? Crategus File Added: diff.txt ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-29 17:00 Message: Logged In: YES user_id=2039760 Originator: YES I have put all changes of the code for the Bessel functions in the orginal file bessel.lisp which is distributed with Maxima 5.14.0 and attached this file to this message. I have tested the code with the testsuite and the values of the mytest_bessel.mac. Crategus File Added: bessel-changed.lisp ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-28 20:35 Message: Logged In: YES user_id=2039760 Originator: YES Thank you very much for your answer. I have redesigned the following routines: simp-bessel-j (old name bessel-j-simp) bessel-j (old name $bessel) simp-bessel-y (old name bessel-y-simp) bessel-y simp-bessel-i (old name bessel-i-simp) bessel-i simp-bessel-k (old name bessel-k-simp) bessel-k The new routines are: $hankel_1 simp-hankel-1 $hankel_2 simp-hankel-2 If have further introduced the nouns %hankel_1 und %hankel_2. I would prefer to call the functions $hankel_1 and $hankel_2 and not bessel_hankel. The reason is that these functions are well known as Hankel functions. Because I have collected the routines in a new file it would be a good idea to do the changes in the orginal file bessel.lisp. Than it easier to produce a diff for you. I do this the next time. But first, I would like to have a second look at the global arrays $besselarray, $yarray and $iarray. I think it is not so easy to manage these global arrays in a consistent and predictable way for the user. One reason is, that the code for the numerical calculations use the other functions too (i.e. to calculate bessel-y we need bessel-j or for the calculation of bessel-i we need bessel-j and bessel-k). So $besselarray can be destroyed when calculating bessel-i. It might be possible to prevent this effect by introducing a global flag which signals the internal call of one of the routines. A second point is, that we often dont't use the numerical slatec routines directly. In this case it is necessary to do a lot of extra calculations to obtain the values for the arrays. Perhaps it is useful to decide to fill the global arrays only for the case when we can use the values of the slatec routines directly. This is possible for positive order and argument or positiv order and complex argument of the Bessel functions. Further, an additional flag could be introduced to switch the filling of the arrays on and off. It is also to decide what we do when the user calculate a new value for a Bessel function and the calculation don't fill the global array. In this case, the user may be confused by the fact that the global arrays still hold old values. But for which order and argument? For this case it may be helpful to invalidate the array every time we calculate a new value or to store additionally the order and argument for which the values in the global array are calculated. All these problems with the global arrays arise because of the extension of the routines to negative arguments and orders. In my opinion it would be the best to use the concept of the global arrays no longer. The code becomes much more complicated and the benefit for the user might be small. Crategus ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-03-28 12:37 Message: Logged In: YES user_id=28849 Originator: NO First, thank you very much for the fixes to the Bessel routines. I'm sure we all want the routines to be correct! Second, about $besselarray. I see the documentation for that is gone, but I think the intent for besselarray was to let the user have access to all the computed values, since some algorithms end up computing all the intermediate orders to get to the desired order. We need to think if this should go away or not. Third, it would certainly help me a lot if you actually produced a diff between your new code and the existing code. Right now, I have to go look at every single function in different places to figure out what's changed. Trying to minimize the changes would help. I know this is a lot of work for you too. Just identifying which function you actually changed would help a lot too. About the specific changes you mention in the comments. No problem with changing bessel-x-simp to simp-bessel-x. (I think most of the code uses simpbesselx, but that's too messy for my tastes.) Extending the range is very, very good. Not sure about the special tests for pure real or imaginary answers. Need to look at that some more. Definitions for Hankel functions are good. May want to name them bessel_hankel_1 to be consistent with other Bessel functions. Certainly want to add some properties like derivatives, but that can wait. All in all, I would very much like to take in your very nice changes. ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-27 12:50 Message: Logged In: YES user_id=2039760 Originator: YES I have tested the Bessel functions with the numerical data attached to this message. For the tests I used a function TEST_BESSEL(value, result, digits) which allow to compare the value with the result within the specified digits. Crategus File Added: mytest_bessel.mac ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-27 12: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 13: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 |
From: SourceForge.net <no...@so...> - 2008-04-18 02:44:35
|
Bugs item #1920177, was opened at 2008-03-19 17:07 Message generated for change (Comment added) made by rtoy 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: Raymond Toy (rtoy) Date: 2008-04-17 22:44 Message: Logged In: YES user_id=28849 Originator: NO Found it. In a few places, the argument was an integer which was not converted to a double for the computations. This resulted in the answers only having single-float accuracy. And the use of truncate should have been changed to floor because the second value of truncate of a negative number is negative. We wanted a positive value, which is what floor gives. In gcl, a single-float type is a double precision number. (A short-float in gcl is a single precision number.) The wronskian test passes now. ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-17 22:04 Message: Logged In: YES user_id=28849 Originator: NO Ok. Something is wrong. Here are some samples with CVS sources and cmucl on linux and sparc: besselexpand:true; bessel_i(-1/2,-5) -> -sqrt(2)*%i*cosh(5)/sqrt(5)/sqrt(%pi) -> -26.47995176430595*%i bessel_i(-0.5,-5.0) -> -26.47995216244259*%i-5.787377599375532e-7 But, curiously, gcl gives -26.47995174630617*%i+1.77635683940025e-15. Need to investigate why the difference... ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-04-17 16:14 Message: Logged In: YES user_id=2039760 Originator: YES First my last build info: Maxima version: 5.14.0cvs Maxima build date: 5:26 4/18/2008 host type: @host@ lisp-implementation-type: GNU Common Lisp (GCL) lisp-implementation-version: GCL 2.6.8 I have build Maxima from the new CVS files today. I run all tests and had no problems. For the tests W_II and W_IK I had adjusted eps to 1e-10. With this value for eps all tests pass in the range of -5.0 to +5.0 for the argument. On my system I get the following results for the tests with real arg (%i23) abs(w_ii(-0.50,-5.00)-(-2.0*sin(-0.5*%pi)/((-5.0)*%pi))),numer; (%o23) 9.4075686847027995E-14 (%i24) abs(w_ik(-0.50,-5.00)-(1/(-5.0))),numer; (%o24) 4.7815106425210348E-13 And the following results with complex argument: (%i28) abs(w_ii(-0.50,-5.00*%i)-(-2.0*sin(-0.5*%pi)/((-5.0*%i)*%pi))),numer; (%o28) 1.0825261750342018E-15 (%i29) abs(w_ik(-0.50,-5.00*%i)-(1/(-5.0*%i))),numer; (%o29) 1.6378141611741369E-15 So I get the values for the Wronskians with much higher accuracy. Crategus ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-16 09:39 Message: Logged In: YES user_id=28849 Originator: NO What changes should be made to the documentation? I've already removed the part about besselarray. Also, I've run the test_bessel_wronskian tests. These need to be modified to fit in with how regression tests are done. However, I get the following messages: W_II: MAX EPS 7.44e-5 for order = -0.50 and arg = -5.00 W_IK: MAX EPS 1.17e-4 for order = -0.50 and arg = -5.00 Are these expected? ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-15 12:43 Message: Logged In: YES user_id=28849 Originator: NO Thank you for the fixes! I'll fix the comment in German and I'll update the documentation too. I'll add your new tests as well. Once these are done, I'll close this bug report as fixed. If you find more issues, please open a new bug report. Thanks! ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-15 12:42 Message: Logged In: YES user_id=28849 Originator: NO Thank you for the fixes! I'll fix the comment in German and I'll update the documentation too. I'll add your new tests as well. Once these are done, I'll close this bug report as fixed. If you find more issues, please open a new bug report. Thanks! ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-15 12:42 Message: Logged In: YES user_id=28849 Originator: NO Thank you for the fixes! I'll fix the comment in German and I'll update the documentation too. I'll add your new tests as well. Once these are done, I'll close this bug report as fixed. If you find more issues, please open a new bug report. Thanks! ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-04-12 17:15 Message: Logged In: YES user_id=2039760 Originator: YES Thank you very much for your work and your interest in the changes of the code. I have seen that you have deleted the use of the global $besselarray too. I have done further numerical calculations to test the code. I think a good test to show that the different parts of the code work well is to calculate the Wronskians which are defined in A&S. I have attached a file which does the calculations in test mode. (Excuse my programing style, but I have no experience programing in Maxima. Excuse my Englisch too.) As a side effect of the extension of the Bessel functions it is now possible to get nice plots to show the symmetry of the Bessel functions. A hint: You can kill the comment in German I had added at line 465. The question asks for the sign. The sign is correct! A second hint: There a still some parts of the code which don't use the symmetries of the Bessel functions completly. I will work further on this subject. At last: The documentation has to be updated too. I would like to do this work, but my English is not good enough. Crategus File Added: test_bessel_wronskian.mac ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-11 10:17 Message: Logged In: YES user_id=28849 Originator: NO Changes checked in to bessel.lisp, rev 1.50. ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-10 12:57 Message: Logged In: YES user_id=28849 Originator: NO I have now integrated your changes into bessel.lisp. The testsuite passes and your bessel tests pass too. I changed test_bessel(a,b,d) to mean abs(a-b)<10^(-d), which I hope is close enough. (Perhaps the real and imaginary parts should be tested separately?) I'll check in these changes after some code clean up. ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-04-01 17:02 Message: Logged In: YES user_id=2039760 Originator: YES Again, thank you very much for your interest. I have not changed any of the derivatives. I have only changed the routines which I have listed in this thread. But I have not used the last CVS file. I have used the file which was send with the Maxima 5.14.0 release. There was one change for the derviative for bessel_k since the last release. On Sunday I have installed a tool to download the CVS files. Now I can get the really actual files. The next time I would like to optimize the routines further. Additionally, I am preparing extensive numerical tests for the functions. At the moment I have no idea how to extent the calculation to complex order. In my textbooks and in the book of Abramowitz and Stegun I can't find a relationship which help to calculate the Bessel functions for complex order. Crategus ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-01 12:41 Message: Logged In: YES user_id=28849 Originator: NO Thank you very much for the updated files. (The diff should have used diff -u or diff -c, but no matter. The bessel-changed.lisp file is perfect.) Given the imminent release of Maxima, I will wait on these changes until after the release. Then I'll try to incorporate your changes. I noticed that you changed the derivatives of some of the Bessel functions. Was that intentional? Were they wrong? (Some of the derivatives need to be in a certain form for some simplifications to be done.) ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-29 17:05 Message: Logged In: YES user_id=2039760 Originator: YES Next I have added a diff between the orginal and the changed file. Question: Is it possible to attach at once more than one file to a message? Crategus File Added: diff.txt ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-29 17:00 Message: Logged In: YES user_id=2039760 Originator: YES I have put all changes of the code for the Bessel functions in the orginal file bessel.lisp which is distributed with Maxima 5.14.0 and attached this file to this message. I have tested the code with the testsuite and the values of the mytest_bessel.mac. Crategus File Added: bessel-changed.lisp ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-28 20:35 Message: Logged In: YES user_id=2039760 Originator: YES Thank you very much for your answer. I have redesigned the following routines: simp-bessel-j (old name bessel-j-simp) bessel-j (old name $bessel) simp-bessel-y (old name bessel-y-simp) bessel-y simp-bessel-i (old name bessel-i-simp) bessel-i simp-bessel-k (old name bessel-k-simp) bessel-k The new routines are: $hankel_1 simp-hankel-1 $hankel_2 simp-hankel-2 If have further introduced the nouns %hankel_1 und %hankel_2. I would prefer to call the functions $hankel_1 and $hankel_2 and not bessel_hankel. The reason is that these functions are well known as Hankel functions. Because I have collected the routines in a new file it would be a good idea to do the changes in the orginal file bessel.lisp. Than it easier to produce a diff for you. I do this the next time. But first, I would like to have a second look at the global arrays $besselarray, $yarray and $iarray. I think it is not so easy to manage these global arrays in a consistent and predictable way for the user. One reason is, that the code for the numerical calculations use the other functions too (i.e. to calculate bessel-y we need bessel-j or for the calculation of bessel-i we need bessel-j and bessel-k). So $besselarray can be destroyed when calculating bessel-i. It might be possible to prevent this effect by introducing a global flag which signals the internal call of one of the routines. A second point is, that we often dont't use the numerical slatec routines directly. In this case it is necessary to do a lot of extra calculations to obtain the values for the arrays. Perhaps it is useful to decide to fill the global arrays only for the case when we can use the values of the slatec routines directly. This is possible for positive order and argument or positiv order and complex argument of the Bessel functions. Further, an additional flag could be introduced to switch the filling of the arrays on and off. It is also to decide what we do when the user calculate a new value for a Bessel function and the calculation don't fill the global array. In this case, the user may be confused by the fact that the global arrays still hold old values. But for which order and argument? For this case it may be helpful to invalidate the array every time we calculate a new value or to store additionally the order and argument for which the values in the global array are calculated. All these problems with the global arrays arise because of the extension of the routines to negative arguments and orders. In my opinion it would be the best to use the concept of the global arrays no longer. The code becomes much more complicated and the benefit for the user might be small. Crategus ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-03-28 12:37 Message: Logged In: YES user_id=28849 Originator: NO First, thank you very much for the fixes to the Bessel routines. I'm sure we all want the routines to be correct! Second, about $besselarray. I see the documentation for that is gone, but I think the intent for besselarray was to let the user have access to all the computed values, since some algorithms end up computing all the intermediate orders to get to the desired order. We need to think if this should go away or not. Third, it would certainly help me a lot if you actually produced a diff between your new code and the existing code. Right now, I have to go look at every single function in different places to figure out what's changed. Trying to minimize the changes would help. I know this is a lot of work for you too. Just identifying which function you actually changed would help a lot too. About the specific changes you mention in the comments. No problem with changing bessel-x-simp to simp-bessel-x. (I think most of the code uses simpbesselx, but that's too messy for my tastes.) Extending the range is very, very good. Not sure about the special tests for pure real or imaginary answers. Need to look at that some more. Definitions for Hankel functions are good. May want to name them bessel_hankel_1 to be consistent with other Bessel functions. Certainly want to add some properties like derivatives, but that can wait. All in all, I would very much like to take in your very nice changes. ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-27 12:50 Message: Logged In: YES user_id=2039760 Originator: YES I have tested the Bessel functions with the numerical data attached to this message. For the tests I used a function TEST_BESSEL(value, result, digits) which allow to compare the value with the result within the specified digits. Crategus File Added: mytest_bessel.mac ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-27 12: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 13: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 |
From: SourceForge.net <no...@so...> - 2008-04-18 16:28:26
|
Bugs item #1920177, was opened at 2008-03-19 17:07 Message generated for change (Comment added) made by rtoy 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: Closed >Resolution: Fixed 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: Raymond Toy (rtoy) Date: 2008-04-18 12:28 Message: Logged In: YES user_id=28849 Originator: NO New tests were added, which pass on clisp, cmucl, and gcl. Closing this bug as fixed. ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-17 22:44 Message: Logged In: YES user_id=28849 Originator: NO Found it. In a few places, the argument was an integer which was not converted to a double for the computations. This resulted in the answers only having single-float accuracy. And the use of truncate should have been changed to floor because the second value of truncate of a negative number is negative. We wanted a positive value, which is what floor gives. In gcl, a single-float type is a double precision number. (A short-float in gcl is a single precision number.) The wronskian test passes now. ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-17 22:04 Message: Logged In: YES user_id=28849 Originator: NO Ok. Something is wrong. Here are some samples with CVS sources and cmucl on linux and sparc: besselexpand:true; bessel_i(-1/2,-5) -> -sqrt(2)*%i*cosh(5)/sqrt(5)/sqrt(%pi) -> -26.47995176430595*%i bessel_i(-0.5,-5.0) -> -26.47995216244259*%i-5.787377599375532e-7 But, curiously, gcl gives -26.47995174630617*%i+1.77635683940025e-15. Need to investigate why the difference... ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-04-17 16:14 Message: Logged In: YES user_id=2039760 Originator: YES First my last build info: Maxima version: 5.14.0cvs Maxima build date: 5:26 4/18/2008 host type: @host@ lisp-implementation-type: GNU Common Lisp (GCL) lisp-implementation-version: GCL 2.6.8 I have build Maxima from the new CVS files today. I run all tests and had no problems. For the tests W_II and W_IK I had adjusted eps to 1e-10. With this value for eps all tests pass in the range of -5.0 to +5.0 for the argument. On my system I get the following results for the tests with real arg (%i23) abs(w_ii(-0.50,-5.00)-(-2.0*sin(-0.5*%pi)/((-5.0)*%pi))),numer; (%o23) 9.4075686847027995E-14 (%i24) abs(w_ik(-0.50,-5.00)-(1/(-5.0))),numer; (%o24) 4.7815106425210348E-13 And the following results with complex argument: (%i28) abs(w_ii(-0.50,-5.00*%i)-(-2.0*sin(-0.5*%pi)/((-5.0*%i)*%pi))),numer; (%o28) 1.0825261750342018E-15 (%i29) abs(w_ik(-0.50,-5.00*%i)-(1/(-5.0*%i))),numer; (%o29) 1.6378141611741369E-15 So I get the values for the Wronskians with much higher accuracy. Crategus ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-16 09:39 Message: Logged In: YES user_id=28849 Originator: NO What changes should be made to the documentation? I've already removed the part about besselarray. Also, I've run the test_bessel_wronskian tests. These need to be modified to fit in with how regression tests are done. However, I get the following messages: W_II: MAX EPS 7.44e-5 for order = -0.50 and arg = -5.00 W_IK: MAX EPS 1.17e-4 for order = -0.50 and arg = -5.00 Are these expected? ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-15 12:43 Message: Logged In: YES user_id=28849 Originator: NO Thank you for the fixes! I'll fix the comment in German and I'll update the documentation too. I'll add your new tests as well. Once these are done, I'll close this bug report as fixed. If you find more issues, please open a new bug report. Thanks! ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-15 12:42 Message: Logged In: YES user_id=28849 Originator: NO Thank you for the fixes! I'll fix the comment in German and I'll update the documentation too. I'll add your new tests as well. Once these are done, I'll close this bug report as fixed. If you find more issues, please open a new bug report. Thanks! ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-15 12:42 Message: Logged In: YES user_id=28849 Originator: NO Thank you for the fixes! I'll fix the comment in German and I'll update the documentation too. I'll add your new tests as well. Once these are done, I'll close this bug report as fixed. If you find more issues, please open a new bug report. Thanks! ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-04-12 17:15 Message: Logged In: YES user_id=2039760 Originator: YES Thank you very much for your work and your interest in the changes of the code. I have seen that you have deleted the use of the global $besselarray too. I have done further numerical calculations to test the code. I think a good test to show that the different parts of the code work well is to calculate the Wronskians which are defined in A&S. I have attached a file which does the calculations in test mode. (Excuse my programing style, but I have no experience programing in Maxima. Excuse my Englisch too.) As a side effect of the extension of the Bessel functions it is now possible to get nice plots to show the symmetry of the Bessel functions. A hint: You can kill the comment in German I had added at line 465. The question asks for the sign. The sign is correct! A second hint: There a still some parts of the code which don't use the symmetries of the Bessel functions completly. I will work further on this subject. At last: The documentation has to be updated too. I would like to do this work, but my English is not good enough. Crategus File Added: test_bessel_wronskian.mac ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-11 10:17 Message: Logged In: YES user_id=28849 Originator: NO Changes checked in to bessel.lisp, rev 1.50. ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-10 12:57 Message: Logged In: YES user_id=28849 Originator: NO I have now integrated your changes into bessel.lisp. The testsuite passes and your bessel tests pass too. I changed test_bessel(a,b,d) to mean abs(a-b)<10^(-d), which I hope is close enough. (Perhaps the real and imaginary parts should be tested separately?) I'll check in these changes after some code clean up. ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-04-01 17:02 Message: Logged In: YES user_id=2039760 Originator: YES Again, thank you very much for your interest. I have not changed any of the derivatives. I have only changed the routines which I have listed in this thread. But I have not used the last CVS file. I have used the file which was send with the Maxima 5.14.0 release. There was one change for the derviative for bessel_k since the last release. On Sunday I have installed a tool to download the CVS files. Now I can get the really actual files. The next time I would like to optimize the routines further. Additionally, I am preparing extensive numerical tests for the functions. At the moment I have no idea how to extent the calculation to complex order. In my textbooks and in the book of Abramowitz and Stegun I can't find a relationship which help to calculate the Bessel functions for complex order. Crategus ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-04-01 12:41 Message: Logged In: YES user_id=28849 Originator: NO Thank you very much for the updated files. (The diff should have used diff -u or diff -c, but no matter. The bessel-changed.lisp file is perfect.) Given the imminent release of Maxima, I will wait on these changes until after the release. Then I'll try to incorporate your changes. I noticed that you changed the derivatives of some of the Bessel functions. Was that intentional? Were they wrong? (Some of the derivatives need to be in a certain form for some simplifications to be done.) ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-29 17:05 Message: Logged In: YES user_id=2039760 Originator: YES Next I have added a diff between the orginal and the changed file. Question: Is it possible to attach at once more than one file to a message? Crategus File Added: diff.txt ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-29 17:00 Message: Logged In: YES user_id=2039760 Originator: YES I have put all changes of the code for the Bessel functions in the orginal file bessel.lisp which is distributed with Maxima 5.14.0 and attached this file to this message. I have tested the code with the testsuite and the values of the mytest_bessel.mac. Crategus File Added: bessel-changed.lisp ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-28 20:35 Message: Logged In: YES user_id=2039760 Originator: YES Thank you very much for your answer. I have redesigned the following routines: simp-bessel-j (old name bessel-j-simp) bessel-j (old name $bessel) simp-bessel-y (old name bessel-y-simp) bessel-y simp-bessel-i (old name bessel-i-simp) bessel-i simp-bessel-k (old name bessel-k-simp) bessel-k The new routines are: $hankel_1 simp-hankel-1 $hankel_2 simp-hankel-2 If have further introduced the nouns %hankel_1 und %hankel_2. I would prefer to call the functions $hankel_1 and $hankel_2 and not bessel_hankel. The reason is that these functions are well known as Hankel functions. Because I have collected the routines in a new file it would be a good idea to do the changes in the orginal file bessel.lisp. Than it easier to produce a diff for you. I do this the next time. But first, I would like to have a second look at the global arrays $besselarray, $yarray and $iarray. I think it is not so easy to manage these global arrays in a consistent and predictable way for the user. One reason is, that the code for the numerical calculations use the other functions too (i.e. to calculate bessel-y we need bessel-j or for the calculation of bessel-i we need bessel-j and bessel-k). So $besselarray can be destroyed when calculating bessel-i. It might be possible to prevent this effect by introducing a global flag which signals the internal call of one of the routines. A second point is, that we often dont't use the numerical slatec routines directly. In this case it is necessary to do a lot of extra calculations to obtain the values for the arrays. Perhaps it is useful to decide to fill the global arrays only for the case when we can use the values of the slatec routines directly. This is possible for positive order and argument or positiv order and complex argument of the Bessel functions. Further, an additional flag could be introduced to switch the filling of the arrays on and off. It is also to decide what we do when the user calculate a new value for a Bessel function and the calculation don't fill the global array. In this case, the user may be confused by the fact that the global arrays still hold old values. But for which order and argument? For this case it may be helpful to invalidate the array every time we calculate a new value or to store additionally the order and argument for which the values in the global array are calculated. All these problems with the global arrays arise because of the extension of the routines to negative arguments and orders. In my opinion it would be the best to use the concept of the global arrays no longer. The code becomes much more complicated and the benefit for the user might be small. Crategus ---------------------------------------------------------------------- Comment By: Raymond Toy (rtoy) Date: 2008-03-28 12:37 Message: Logged In: YES user_id=28849 Originator: NO First, thank you very much for the fixes to the Bessel routines. I'm sure we all want the routines to be correct! Second, about $besselarray. I see the documentation for that is gone, but I think the intent for besselarray was to let the user have access to all the computed values, since some algorithms end up computing all the intermediate orders to get to the desired order. We need to think if this should go away or not. Third, it would certainly help me a lot if you actually produced a diff between your new code and the existing code. Right now, I have to go look at every single function in different places to figure out what's changed. Trying to minimize the changes would help. I know this is a lot of work for you too. Just identifying which function you actually changed would help a lot too. About the specific changes you mention in the comments. No problem with changing bessel-x-simp to simp-bessel-x. (I think most of the code uses simpbesselx, but that's too messy for my tastes.) Extending the range is very, very good. Not sure about the special tests for pure real or imaginary answers. Need to look at that some more. Definitions for Hankel functions are good. May want to name them bessel_hankel_1 to be consistent with other Bessel functions. Certainly want to add some properties like derivatives, but that can wait. All in all, I would very much like to take in your very nice changes. ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-27 12:50 Message: Logged In: YES user_id=2039760 Originator: YES I have tested the Bessel functions with the numerical data attached to this message. For the tests I used a function TEST_BESSEL(value, result, digits) which allow to compare the value with the result within the specified digits. Crategus File Added: mytest_bessel.mac ---------------------------------------------------------------------- Comment By: Crategus (crategus) Date: 2008-03-27 12: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 13: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 |