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
