HP Forums

Full Version: HP 30b numeric integration program
You're currently viewing a stripped down version of our content. View the full version with proper formatting.

Here is a numeric integration program for the HP 30b. This is a translation of the HP-25 Program Library entry # L19 from 65 Notes V4N2P32 originally by Bandes, Davidson, and Derrick.

Usage:

At LBL 99, key in the function to be integrated. When LBL 99 is executed, the value of the X to be evaluated is in the X register (display) at entry to LBL 99. If you need it again to evaluate the function, it is stored in memory 3. Memories 6 and up are available as well.

At the PRGM location in which this program is stored (at the proper entry in the HP 30b program catalog), Key lower limit, press =, key upper limit, press =. Then press = again to start the program.

This program uses Simpson's rule with a default value of 200 for the number of increments.

The program is VERY fast. To solve the integral of e^x from 0 to 1 using 200 as the number of increments takes less than 2 seconds if this program is program 0 in the HP 30b.

The integral of e^x from 0 to 1 is solved to be 1.71828182849. Value of e minus 1 is 1.71828182846, an error of 3 in the last position.

The integral of SIN(PI*X) from 0 to 1 is solved in about 5-6 seconds with 500 as the number of increments with a solved value of 0.636619772363.

Register usage:

Memory 0 = number of increments, default is 200. Must be an even value.
Memory 1 = lower limit of integration
Memory 2 = upper limit of integration
Memory 3 = value of x being evaluated
Memory 4 = a positive or negative 1, depending on the loop.
Memory 5 = value of integrated function (sum)

Version shown here uses 92 bytes with a checksum of 016.

Improvements can be made to this code. The labels chosen correspond to the line numbers of the GTO statements in the original HP 25 program. Suggestions for improvements:

1) I used LBL 99 for the function to make the function occur at the end of the listing. This could be put into the code right after LBL 06, saving bytes and speeding up execution.

2) Lines 33-35 are to make the stack mimic what was in the HP 25 program. I have not found an easier way to do this, but it is clumsy as is.

3) Lines 36-40 and 48-51 are there to accomplish the following: memory 4 is initialized to zero. The first time through the loop, this should have negative 1 placed into it. Thereafter, the sign of the contents of memory 4 should be switched from negative to positive to negative, etc., each time through the loop. I freely admit this is again clumsy, but it does work.

4) Use something other than Simpsons, especially since it will fail if an endpoint is undefined. Integrating LN(X) from 0 to 1 fails, of course. :-)

Listing:

STO 2
SWAP
STO 1
200
STO 0
0
STO 4
STO 5 ...... LINE 10
RCL 1
STO 3
LBL 06
CALL 99
RCL 2
RCL 1
-
RCL 0
/
STO+ 3 ...... LINE 20
3
/
*
STO+ 5
ANS
RCL 2
RCL 3
-
+
0 ...... LINE 30
?<
GT 49
ROLL DOWN
ROLL DOWN
=
RCL 4
GF 48
RCL 4
LBL 47
+/- ...... LINE 40
STO 4
2
+
*
STO+ 5
RCL 3
GTO 06
LBL 48
RCL 4
1 ...... LINE 50
GTO 47
LBL 49
RCL 5
STOP
LBL 99
E^X
RTN

Enjoy and improve.

Edited: 28 Apr 2010, 3:06 p.m.

While the overly sensitive might take a lack of replies to a post like this personally, I'm not (overly) sensitive. :-)

I did some additional comparisons of this integration program with the examples given in John Kennedy's HP 34C review from the PPC Journal V6N5P19-20. He gives 6 examples as shown below:

Using the above program with 500 intervals, the following results are obtained.

Note that if the number of intervals are reduced to 50, the results still appear to be accurate to about 8-10 decimal places and the run times are 1-2 seconds!

1. Result of 1666.6666666633 in about 4 seconds.

2. Result of 32.6666666659 in about 5 seconds

3. Result of -5.4365x10^(-10) in about 4 seconds

4. Result of 3.14159265355 in about 5 seconds

5. Result of 4.40984450786 in about 5 seconds

6. Result of 0.0981747704617 in about 6 seconds.


Perhaps when more of you have an HP 30b, this will generate a bit more excitement? :-)

Gene,

I am not nearly smart enough to make even a stupid remark about this post.

--Sancerre

Quote:
Gene: ...a lack of replies to a post like this ...

Sancerre: I am not nearly smart enough to make even a stupid remark about this post.


Me: Ditto.

This is another numeric integration program for the HP 30b. It is a translation of the HP 41 program found in V7N8P10-11 of the PPC Journal. It uses an iterative 3-point Gaussian Quadrature.

And, as you will notice, I am not smart enough to write my own ORIGINAL integration program. I'm having to translate the work of others. :-)

As typed below with LBL 98 containing the SIN function, it is 86 bytes with a checksum of 221. As written, it is very fast. It solves the example problems from the PPC ROM routine IG with varying levels of accuracy, generally 4-8 decimal places - no doubt due to the 3 point nature of this approach compared to a 16-point approach. :-) But, compared to the simpson's approach, it will evaluate integrals where the endpoints are undefined.

Key in your function to integrate at LBL 98. The value of X is in the X register. If you need it more than once, save it into memories 0, 1, or 2. Register usage is 6 - lower limit, 7 - upper limit, 8 - number of iterations, 9 - sums (scratch), 5 - offset (scratch), 4 - midpoints (scratch), 3- dx (scratch)


** Why do these things? Well, here is a BUSINESS calculator with RPN that can be programmed to do all sorts of non-business things...and it is fast. **

Program listing:

STO 7
SWAP
STO 6
STO 4
-
9
9
STO 8
/
STO 3 ... LINE 10
2
/
STO+ 4
0
.
6
SQRT
*
STO 5
0
STO 9 ... LINE 20
LBL 11
RCL 4
RCL 5
-
CALL 98
5
*
STO+ 9
RCL 4
CALL 98 ... LINE 30
8
*
STO+9
RCL 4
RCL 5
+
CALL 98
5
*
STO+ 9 ... LINE 40
RCL 3
STO+ 4
DSE 8
GT 11
RCL 9
1
8
/
STOP ... LINE 50
LBL 98
SIN
RTN

Edited: 30 Apr 2010, 4:55 p.m.

Thanks for posting the programs.
Would you please time the one below? The exact twelve-digit result is 0.03648997398, or 2 + Psi(1/2). (On the 50g Psi is accessible through white left shift MTH NXT SPECI).

/1
|
| (sqrt(x)/(x-1) - 1/ln(x))dx
/0

Quote:
Perhaps when more of you have an HP 30b, this will generate a bit more excitement? :-)

Perhaps I should get myself one.

The simpson's version won't run, of course, since it tries to evaluate LN(0).

The 3-point iterative quadrature program with 50 iterations returns 0.0364973556 in about 1-2 seconds. It's very fast.

Increasing the iterations to 500 takes about 12 seconds, returning 0.0364906580, which is surprisingly good (IMO).

Hey, it IS only a 3-point method. :-)

The code I used for LBL 98 was stack-oriented saving a register. These were my steps:

LBL 98
LN
ANS
SWAP
1/X
SWAP
SQRT
ANS
1
-
/
SWAP
-
RTN

Quote:
The simpson's version won't run, of course, since it tries to evaluate LN(0).

Changing the integration bounds to 1e-38 and 0.9999999 shouldn't affect the 12-digit result, but even Simpson 3/8 with 19500 intervals performs badly on this one:


Program Simpson_3_8;
var a, b, h, s1, s2, x: real;
n: integer;

function f(x: real): real;
begin
f:=Sqrt(x)/(x-1)-1/Ln(x)
end;

begin
ClrScr;
Write('A, B, N: ');
ReadLn(a, b, n);
n:=3*n;
h:=(b-a)/n;
s1:=0;
s2:=0;
x:=a+h;
repeat
s1:=s1+f(x+2*h);
s2:=s2+f(x)+f(x+h);
x:=x+3*h
until x>a+(n-3)*h;
s2:=s2+f(x)+f(x+h);
WriteLn(3*h/8*(f(a)+f(b)+2*s1+3*s2))
end.

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

A, B, N: 1e-38 0.9999999 6500
3.6498295709E-02

Anyway, this is a tough integral for numerical methods:

http://www.cs.berkeley.edu/~wkahan/Math128/INTGTkey.pdf

Quote:
(Gene Wright) did some additional comparisons of this integration program with the examples given in (an) HP-34C review from the PPC Journal V6N5P19-20.

Gene --

Thanks for presenting the table. The author's review of the HP-34C was favorable, yet he chose six integrals that, in a sense, made the HP-34C's INTEG capability 'look bad' by comparison to 16-point Gaussian quadrature. The GQ results, the author points out, were at least as accurate and were obtained in less time. Furthermore, they were produced by an RPN routine programmed on the HP-34C -- not a built-in microcoded function. He speculates that the GQ approach "may not have been used by HP because of internal storage requirements."

Maybe so, but I notice that all the examples utilized very well-behaved integrands over narrow intervals, for which even simple methods will tend to produce good results. HP's objective in this inaugural effort was robustness and practicality -- to provide sufficiently-accurate real-world results in science and engineering without the user's having to specify the number of evaluations, or avoid endpoints for which the function was undefined, or risk wrong answers due to an inadequate sampling rate for a periodic function.

I do note that the GQ program uses all 21 storage registers and consumes 40 lines of programming. The user's program defining the integrand must not employ any register, because no unused ones are available. It also cannot exceed 30 lines, lest register R19 ('.9') be automatically deallocated, rendering the GQ program inoperative. (Score a point for the successor HP-15C and its manual memory allocation.)

The exact value for the fifth integral is 6 - 236*e^(-5) = 4.40984450821.

The exact value for the sixth integral is 3*pi/96 = 0.0981747704247.


The HP-34C GQ result and timing for the "Kahan integral" in the paper linked by Gerson:

f(x) = sqrt(x)/(x-1) - 1/ln(x)
0 <= x <= 1

At each endpoint, f(x) is undefined but does approach zero, as shown by limit analysis. From x = 0.000, f(x) rises quickly to a 'smooth' local maximum of 0.117680329841 at x = 0.00473796324935, then gradually and monotonically declines to a local minimum of ...
-0.117680329841 at x = 211.061155894 (which is 1 / 0.00473796324935), as this function has the property f(1/x) = -f(x) for all x (where defined).

The derivative:

             1             (x + 1)
df/dx = ---------- - -----------------
x*(ln x)^2 2*sqrt(x)*(x-1)^2

The exact value of the integral is (2 - ln 4 - gamma), where gamma is the Euler-Mascheroni Constant of 0.577215664901532...

FIX 5            FIX 6           'Exact'        Gaussian

0.03661813695 0.03649023200 0.036489973980 0.0364900935
19 seconds 268 seconds 45 seconds

The 16-point GQ outperforms the built-in Romberg method in this instance, too. Over this short interval of [0, 1], the 16 points seem to provide adequate coverage for high accuracy, and it does not evaluate at endpoints -- a definite 'plus' in this instance.

The highly-accurate GQ result might have been 'dumb luck', as no evaluation point (node) fell between zero and 0.004738.

I do wonder about its performance over large intervals with erratic integrands, but the Gauss - Kronrod method (used by the TI-82 and likely its successors) handles that situation by comparing the difference between estimations and taking more evaluations as necessary, much as HP's modified Romberg method does.

I may add a section to my HP Articles Archive contribution (#556), with a deeper analysis of this integral and computational methods.

Gerson said above:

Quote:
... but even Simpson 3/8 with 19500 intervals performs badly on this one:

3.6498295709E-02

Anyway, this is a tough integral for numerical methods:


In reasonable amounts of time with endpoints of 1E-8 and 1 - 1E-8, the Casio fx-115MS obtains 0.0365 using Simpson's Rule at 512 points, but the successor fx-115ES obtains 0.03648997346 using Gauss-Kronrod (and, I might add, a perfect answer of 0.03648997398 using 1E-10 and 1 - 1E-10 as the endpoints). Then, after a period of inactivity, they shut themselves off and forgot the integrals I'd entered...

Gene said below:

Quote:
Karl, what I'd rather do is make sure you get an HP 30b soon to try your hand at programming it for things like this.

Nah, I'd rather not program things that should be built into a calculator, even though the HP-30b is a business model. Besides, the HP-34C LED 'light show' should be played from time to time...

-- KS


Edited: 11 May 2010, 11:48 a.m. after one or more responses were posted

Karl, what I'd rather do is make sure you get an HP 30b soon to try your hand at programming it for things like this.

In the meantime, I'd be glad to try the routines I've translated for some integrals if you care to suggest some examples.

Or... if you know of an RPN program to do integration that doesn't do lots of X<> IND Z or DSE Y type of stuff, I could try translating that. The 30b stack behavior is slightly different than the HP 41 style "RPN" since it has some RPL tendencies, so sometimes complicated stack-using programs are a pain to translate.

Hello Karl,

Thanks for taking the time to do all these analyses and comparisons among different integrating methods.

Quote:
The 16-point GQ outperforms the built-in Romberg method in this instance, too. Over this short interval of [0, 1], the 16 points seem to provide adequate coverage for high accuracy, and it does not evaluate at endpoints -- a definite 'plus' in this instance.

I had seen JK's program before but I never bothered keying it into the 34C. I will try it later on my double-speed 15C.

Regards,

Gerson.