16 Point Gaussian Quadrature for the HP-41C/CV/CX (updated)

[Valentin Albillo]

(11-12-2023 05:33 PM)Namir Wrote:  

(11-12-2023 03:24 PM)Valentin Albillo Wrote:  Both values end in 0 ? Really ? It's a 1-in-100 chance. Did you multiply them by 10 as I suggested ? If yes, then they'd be like 6.93... and 9.31... respectively, not the 0.693... and 0.931... values you posted.
Could you please check them out one final time ?

Yes I multiplied both results by 10 ans then moved the decimal to the right, so they correspond to the original answers (now with 10 digits)

Understood, thanks. However, one of your updated coefficients is still in error, namely:

124 6.22535394E-2
125 STO 05

which is missing a digit between the 6^th and 7^th positions. You might want to update it with the correct coefficient.

As a bonus, perhaps you'd like to try your (corrected) program vs. these four integrals and compare your results with the ones provided, which are accurate to all 12 digits shown:

Code:

       / 1
I1 =   |   cos(x)*ln(x) .dx =  -0.946083070367
       / 0

Code:

       / 1
I2 =   |   1/√x .dx =  2.00000000000
       / 0

Code:

       / 1
I3 =   |   Ln(Gamma(x)) .dx =  0.918938533205
       / 0

Code:

       / 1
I4 =   |   1/√(-Ln(x)) .dx =  1.77245385091
       / 0

Regards.
V.

[Namir]

Interesting integrals. Since the lower bounds yield "infinite" values (although the Gaussian quadrature does not evaluate the integrated function at the end points), the resuts of my programs are accurate to 2 decimal places AT BEST. That's what happens when we make the program "walk" on a "mine" field!!!

Namir

[Valentin Albillo]

.
Hi, Namir,

(11-13-2023 12:46 AM)Namir Wrote:  Interesting integrals. Since the lower bounds yield "infinite" values (although the Gaussian quadrature does not evaluate the integrated function at the end points), the resuts of my programs are accurate to 2 decimal places AT BEST. That's what happens when we make the program "walk" on a "mine" field!!!

Yes, but difficult, almost- or outright improper integrals are actually the most common in practice, as Kahan himself pointed out in the Hewlett-Packard Journal August, 1980 issue, p.29, I quote him:

"During the HP-34C's design a suspicion arose that most integrals encountered in practice might be improper or nearly so. Precautions were taken. Now that experience has confirmed the suspicion, we are grateful for those precautions." [...]

Thus it would seem appropriate for us to try and develop our programs so that they can cope with such integrals as best they can. Being effective only when dealing with "tame", easy-peasy integrals doesn't quite cut it.

Also, I see you corrected the erroneous coefficient, attaboy !

Regards.
V.

[Namir]

The four integration tests that you listed are very intgersting for testing HP-15C calculators and emulators. I used the integral of 1/sqrt(x) from 0 to 1. Here are some results!

HP-15 LE gave 1.9998705
HP-15C CE keeps running with NO solution!
JRPN 15C version 2.18 gave 1.9961232
RLM 15 emulator emulator gave 5.878602E44
HP15C emulator Version 4.5.00, Build 6308 © 1997-2023 Torsten Manz gave 1.9932975

I was very suprised that the HP-15C kept going. I stopped the calculator after a few minutes!

I believe that the HP-15C is using a version of Romberg integration. Is this true?

I usually test algorithms with easy cases to see if they work at all, and if my implementatipon is correct.

Namir

[Valentin Albillo]

.
Hi, Namir,

Apart from the 5 wrong coefficients which I pointed out and you already corrected, there's another error in your program, namely (the ► are cosmetic):

Your program begins:

   01  LBL "A-B"
   02  LBL A
   03  XEQ 00 ►
   04  FIX 05
         ...
where the XEQ 00 at line 03 calls this subroutine, which stores the 16 coefficients needed to compute the integral:

   117 ►LBL 00
   118  FS? 01
   119  RTN
   120  2.715245941E-2
   121  STO 03
         ...
   152  SF 01
   153  RTN

Your idea seems to be to use flag 01 in the very first run of your program to signal that the coefficients have already been stored by setting it at line 152 SF 01 so that the next times subroutine 00 is called, flag 01 is checked at line 118 FS? 01, and if set then the subroutine returns immediately without storing the coefficients, as they are assumed to have been already stored.

So far so good. The problem is that if for whatever reason (perhaps as a side effect of previously running another program) flag 01 is already set when you execute this integration program for the first time, then the XEQ 00 will return immediately without storing tha coefficients at all, proceeding to compute the integral without them and thus leaving the user thoroughly perplexed when he gets the thoroughly wrong result.

In short: you're assuming that flag 01 isn't set when you first execute your program, but that's an unsupported assumption and should it be already set, your program will never store the coefficients and thus will fail. A program should never rely on unsupported assumptions, but in this case the solution is to simply clear flag 01 before calling subroutine 00, end of story.

Regards.
V.

[Namir]

(11-15-2023 03:14 AM)Valentin Albillo Wrote:  .
In short: you're assuming that flag 01 isn't set when you first execute your program, but that's an unsupported assumption and should it be already set, your program will never store the coefficients and thus will fail. A program should never rely on unsupported assumptions, but in this case the solution is to simply clear flag 01 before calling subroutine 00, end of story.

Regards.
V.

Clearing flag 01 before the XEQ 00 will make using that flag meaningless. Perhaps adding a CF 01 and RTN after the very first label will do the job!

Code:

LBL "A-B"
CF 01
RTN
LBL A

And ask the use to first do [XEQ][ALPHA]A-B[ALPHA].

I went back and just updated the examples' intructions, asking the user to clear flag 1 if this is the first time he is running the program.