QAG48 Beta 0.1 Aug 26, 1997 QAG48 is a replacement for the internal numerical integration of the HP48's. Bytes: 3383.5 #4258h Library #1548 To install: Download the library to your calculator. Put the library on stack. Type 0 (or the port number of your choice) STO. Turn the calculator off & on. Delete the old copy of the library to save space. The library uses some unsupported but stable entries, they all appear in Entries.doc. The library has been tested on a HP48S rev ?, but not extensively. The library should not give any problems, although there might be some interference with other libraries. What it is: QAG48 is a Numerical Integration library for the HP48's. It has been tested on a HP48G rev R, and a HP48S rev ?. QAG48 is based on QAGS from the QUADPACK library of FORTRAN programs. This library uses an adaptive scheme that estimates errors for intervals and bisects only those whose errors are not within the limit set by the user. For those who care; QAG48 (& QAGS) uses a Gauss Kronrod method with 10/21 points. The 10 point rule is exact for up to 21st order polynomials and is identical to the 10 point Gauss-Legendre rule. The extra 11 points are added in between the first 10 points & the endpoints - which are reused with new weights for all 21 points to be exact for up to 31st order polynomials. The difference between these two rules is the basis for estimating the error of the 21 point rule. Furthermore, QAG48 (& QAGS) also uses the Epsilon Extrapolation method to speed up convergence of the results from the Gauss Kronrod method. The end result of all of this is a very efficient and accurate method for calculating Definite Integrals numerically. It also means that we must store lists of results, errors, integration ranges, and flags for each of the intervals integrated over. This can result with up to 1.1 bytes stored per function evaluation. I have tested some difficult integrals to about 10,000 function evaluations, thus resulting in about 11,000 bytes stored. For the most part this should not cause any problems. Syntax: Enter the integral with the integral symbol. Then press or type QAG48. The library currently only works for a single integration. You can use the HP48's internal numerical integration for the inner integration and QAG48 for the outer integration of a double integral. e.g. 'INTEGRAL(0,pi,sin(X),X)' (where INTEGRAL is the integral sign) then QAG48. The command will run, returning the final result to stack level 1, and the approximated error to the global variable IERR. (same as the HP48's). At any time you can press CANCEL and the command will return the latest result with the latest approximated error in IERR. The program uses a global variable named: 'intpar' and is defined: intpar: {risk level, time limit, relative error, display} risk level: 0= low accuracy, low confidence - Default 1= high accuracy, high confidence The "High accuracy" affects the Epsilon extrapolation only. It will generally run the Epsilon routine another time, and thus take longer. time limit: currently disabled, it has a bug that I haven't been able to fix yet. relative error: 1E-11<= relative error <= .1 ; Default=1E-6 the error bound is approximated by (relative error)*INTEGRAL(a,b,ABS(integrand),var)where the integral is INTEGRAL(a,b,integrand,var) This is the same way that the HP48's do it, except that the relative error is given by the FIX/ENG/SCI modes. display: 0= QAG48 display OFF 1= QAG48 display ON - Default When QAG48 display is on, it shows the latest result, error and the number of function evaluations. You can then decide if the result is accurate enough & CANCEL. QAG48 will then return the latest result with "BEST" tagged onto it, and with the latest error in IERR. 'intpar' is not created automatically (I dislike programs that do this), but does need to be within the directory tree up to HOME. Otherwise QAG48 will use the default values. Performance: The command must do 21 function evaluations to get a result and thus will take over 2 seconds to return the first result. The HP48's internal routine can give results in under a second, particularly with 3 FIX (or a relative error of 1E-3). I have tested QAG48 with 35 integrals. With a relative error of 1E-3, QAG48 was faster than the internal routine for only 7 of them, which includes one that the internal routine gave the wrong answer for. With a relative error of 1E-11, QAG48 was faster than the internal routine for 27 of the 35 test integrals. The internal routine failed to give results to the accuracy desired or totally incorrect results, for 13 of the integrals. QAG48 failed to give results to the accuracy desired for only 2 of the integrals, and even then still gave 5 and 10 places accuracy. One of the integrals ran for 80 minutes on the internal routine and gave only 3 places accuracy, but QAG48 gave the answer in 23 seconds. Comparison to the TI-85: The TI-85 uses a different method of which I've been unable to determine. The TI-85 outperforms the HP48G's about half the time. QAG48 is not faster than the TI-85 for low accuracy (i.e. relative error=1E-3), but does better against the TI-85 at high accuracy (1E-11). One awkward feature of the TI-85 is that the error limit appears to be an absolute limit, thus you may have to adjust it after getting some results. Comparison to the TI-92: The TI-92 has a very powerful integration package. Most of it's "numerical" integration is apparently done by first solving the symbolic integral (lookup tables?) and then applying the limits. Thus it's fast for simple high accuracy integrals. Of the 8 integrals I tested, that the TI-92 used a numerical method (unknown) QAG48 was able to outperform the TI-92 on 4 of them and generally not by much. On Integrands: The type of functions been integrated affects the speed of the results. 1.For "smooth" functions both QAG48 and the HP48's internal routines can give quick and accurate results. 2.For Rapidly Oscillatory functions; for example functions including terms like COS(50*X) over a range of 0 to Pi.; typically the simple Rectangular rule or the Trapezoidal rule can give the best results. For the HP48 user, the internal routine usually outperforms QAG48, except when the coefficient is a multiple of 512 for certain ranges - the internal routine gives the wrong answer for INTEGRAL(0,PI,COS(512*X),X). 3. Singularities on or near (outside) the endpoints: Only for certain types of singularities can the internal routine give decent results. QAG48, however, can and does make short work on these types of integrands, including infinite singularities unless the integrand is oscillatory. 5. Singularities within the interval of integration: These types of integrands can be difficult to get results in a reasonable time. For finite singularities, QAG48 can give results eventually (perhaps up to 3-4 minutes). For infinite singularities, QAG48 can give ballpark figures - indeed virtually no numerical integration routines can give even decent results on such an integrand. The solution is to break up the integrating interval into intervals spanning between the singularities. This can considerably help for finite singularities also. Future improvements: 1. Recode some routines into ML to speed up things. It's all SysRPL right now, and only about 40% of the time is actually used for the function evaluations - plenty of room for improvements. 2. Find and fix bug for time limit stuff. My goal here is to allow the user to set a time limit for QAG48 to run, and then end with "BEST" tagged onto the result to indicate that the accuracy desired was not met. I'm trying to use ABORT, but I get a garbage result. 3. Add tests for detecting divergent integrals, detecting roundoff error, and detecting slowly convergent or divergent integrals. 4. Add support for calculating some improper integrals, such as with one or both integral bounds being infinite Send comments, suggestions, etc. to: John Edry, JEEjohn@aol.com