HP 50g Romberg Integration

[peacecalc]

Hello friends,

because of educational reasons I've written a program for the HP 50g in UserRpl.
In some papers of that kind of numerical integration is written, that the Richardson extrapolation can be easily done with a table formed as a triangle.

The only nice feature of my program is, that it doesn't need such a table, it only works with a list with N elements and this list was the input for the DOSUBS command, which seems to me very practical and useful.

The program isn't optimized in performance or memory. Enjoy it.

I included it as a text-file, because I didn't succeed to include it as a well formatted
program text.


  ROMB.txt (Size: 1.99 KB / Downloads: 39)

[Gerson W. Barbosa]

(03-31-2014 03:28 PM)peacecalc Wrote:  I included it as a text-file, because I didn't succeed to include it as a well formatted
program text.

Just instert your code between code tags:

[cоde]
your code
[/cоde]

Would you please post an example? Thanks!

Gerson.

[peacecalc]

Hello Gerson,

first of all, I tried to include with the tags my program, but it doesn't works, the format was spoiled.

Second the example:

The function to integrate is defined in 'Y' (I forgot this to mention in my first post, sorry) (here: << INV >> for 1/x)

Input:
4.
2.
1.

4 is iteration depth
2 upper limit
1 lower limit

output:

6,9412E-1
{ 6,9444E-1 6,9325E-1 6,9315E-1 }
{ 6,9317E-1 6,9315E-1 }
{ 6,9315E-1 }
6,9315E-1

First value: result from pure trapezoid rule in total
Last value: with Richardson extrapolation
the three lists in between are the "table" for the Richardson extrapolation
(this outputs results every DOSUBS call).

[Dieter]

(04-01-2014 05:24 AM)peacecalc Wrote:  first of all, I tried to include with the tags my program, but it doesn't works, the format was spoiled.

Everything works fine. The only exception are the comments that will not line up because you used tab characters. Which is always a bad idea since there is no standard that defines how many blanks will equal one tab. That's why you should always use blanks instead of tabs for such applications.

Without tabs everything looks fine:

Code:

%%HP: T(3)A(R)F(,);        
\<< DUP2 - { }                                 
    \-> N O U H EL                        @@Locale Variables (# Iterations, upper, lower limit integration
                                          @@                  difference of limits, result list) 
    \<< O Y U Y + H 2, / * 'EL' STO+      @@Calc. T11
        IF   N 2, \>=                        
        THEN                              @@Branch if more then one iteration
        H DUP 2, / 0, U 1,                 
        \-> HM HP YK XK K                 @@Local variables (h/2^(J-2), h/2^(J-1), Sum of trapezoid aeras,
        \<<                               @@                 argument for funct., Index for sum of funct.
                                          @@                 values)
            1, N 1, -                     @@Outer loop from 1 to N-1                
            FOR J                                         
            'EL'                          @@load result list on stack 
            EL J GET                      @@load j-element of result list
            HP 'XK' STO+                  @@load first value for XK
             
            1, K                          @@Inner loop from 1 to 2^J
            START                         @@Calc. the sum of function values
            XK Y 'YK' STO+ 
            HM 'XK' STO+
            NEXT                          @@End of inner loop

            YK HM * + 2, /                @@Calc. TJ1 
            STO+                          @@and adding to result list             

            HP 'HM' STO                   @@Initializing locals for next iteration of outer loop 
            'HP' 2, STO/ 
            0, 'YK' STO 
            2, 'K' STO* 
            U 'XK' STO
            NEXT                          @@End of outer loop 

            EL N GET                      @@Reporting the value of trapez rule (can be left)

            4, 3, \-> QZ QN               @@Initializing locals for Richardson extrapolation
           
             \<< 2, N                     @@Loop for Richardson extrapl.                        
                 START                 
                 EL 2,                    @@load result list on stack and parameter 2, (two elements)
                                          @@for DOSUBS command
                 \<< QZ * SWAP - QN /     @@Calc. a new TJL element TJL = (QZ*TJ(L-1)-T(J-1)(L-1))/QN
                 \>>
                 DOSUBS 

                 DUP                      @@load result list of DOSUBS twice on stack (can be left)

                 'EL' STO                 @@store result of DOSUBS in EL for next iteration
                 4, 'QZ' STO*             @@load new values for next iteration
                 QZ 1, - 'QN' STO
                 NEXT                     @@End of extrapolation loop
            \>>
        \>>
    END EL OBJ\-> DROP                    @@Result of numerical integration
    \>>
\>>

Dieter

[peacecalc]

Hello Dieter,

Quote:The only exception are the comments that will not line up because you used tab characters.

That's the advice I need to avoid such a situation in future, thank you!
And many thanks for deleting all the tabs and filling with spaces in my original text.

Greetings
peaceglue

[Thomas Klemm]

(04-02-2014 05:21 PM)peacecalc Wrote:  And many thanks for deleting all the tabs and filling with spaces in my original text.

That's easy in a decent editor like emacs:
M-x untabify


Cheers
Thomas

[C.Ret]

I just recover from an old notepad the ROMBERG INTEGRATION programm I wrote for my HP-28S a few ago ...

It's regular RPL code, so, it will also run on any "new" RPL systems such as the HP-50g !

Code:


« ->  Fct a b m                           @     Fct - Function or program to integrate  e.g.: « -> x 'sin(x)'» or « SIN » 
                                          @     a b - limits of the integration intervalle
                                          @       m - number of Romberg Iterations (size of table R is m.m)
  « { m m } 0 CON
    1 m FOR i
        1 i FOR j
            IF j 1 == 
            THEN                   
               IF i 1 ==           
               THEN                                    @ R(1,1) = (f(a)+f(b))*(b-a)/2
                  b a - 2 / a Fct EVAL b Fct EVAL + * 
               ELSE                                    @ R(i,1) = R(i-1,1)/2 + h*SUM f(a+ k.h) 
                  DUP i 1 - { 1 } + GET 2 /
                  b a - 2 i 1 - ^ /
                  0
                  1 2 i 1 - ^ FOR k
                        OVER k * a + Fct EVAL +
                  2 STEP * +
               END
            ELSE                                       @ R(i,j) = ( 4^j*R(i,j-1)-R(i-1,j-1) ) / (4^j-1)
               DUP  i 1 - j 1 - 2 ->LIST GET 
               OVER i     j 1 - 2 ->LIST GET 
               4 j 1 - ^ SWAP OVER * ROT -
               SWAP 1 - /
            END
            { i j } SWAP PUT
        NEXT
    NEXT
    { m m } GET
  »
»
'ROMBRG' STO

I just test it with the example \[\int_{-1}^{2} \frac{2}{4t^2+1}dt\] from Wikipedia - Méthode de Romberg :

Code:


« SQ 4 * 1 + 2 SWAP / » 
-1 2 
5  ROMBRG

returns 2.431983278

To see intermediate results simply press [ shift ][ LAST ], the R table is bring back in the stack :

Code:


2.43198 
[[  .77647  0        0        0        0       ]
 [ 1.88824  2.25882  0        0        0       ]
 [ 2.35101  2.50527  2.52170  0        0       ]
 [ 2.42355  2.44773  2.44390  2.44266  0       ]
 [ 2.43077  2.43318  2.43221  2.43202  2.43198 ]]
 { 5 5 }

[peacecalc]

Hello Thomas,

thank you for sharing your code. The code seems to me more clearly and its output is nicer formatted than my program.

Nevertheless I'm tried to make a code, which doesn't use exponential operator "^" and this is replaced by multiplication and division. But this costs more use of the storage operations.

Greetings peacecalc

[Gerson W. Barbosa]

(04-02-2014 10:02 PM)C.Ret Wrote:  I just recover from an old notepad the ROMBERG INTEGRATION programm I wrote for my HP-28S a few ago ...

It's regular RPL code, so, it will also run on any "new" RPL systems such as the HP-50g !

As you said, your code runs with no modification on the HP 50g:

Code:


%%HP: T(3)A(D)F(.);
\<< \-> Fct a b m
  \<< { m m } 0. CON 1. m
    FOR i 1. i
      FOR j
        IF j 1. ==
        THEN
          IF i 1. ==
          THEN b a - 2. / a Fct EVAL b Fct EVAL + *
          ELSE DUP i 1. - { 1. } + GET 2. / b a - 2. i 1. - ^ / 0. 1. 2. i 1. - ^
            FOR k OVER k * a + Fct EVAL + 2.
            STEP * +
          END
        ELSE DUP i 1. - j 1. - 2. \->LIST GET OVER i j 1. - 2. \->LIST GET 4. j 1. - ^ SWAP OVER * ROT - SWAP 1. - /
        END { i j } SWAP PUT
      NEXT
    NEXT { m m } GET
  \>>
\>>

A test:

Code:


« LN NEG LN NEG »
0.000000000001
0.999999999999
20
« ROMBERG » TEVAL

0.577222197843
s:290.1873               (4 m 50.1 s on the emulator --> about three hours on the HP 50g)

The buit-in integrator, OTOH:

Code:


« '∫(0.000000000001,0.999999999999,-LN(-LN(X)),X)' EVAL » TEVAL
0.577215664728
s: 52.8994               (52.9 s on the emulator)


« '∫(0.,1.,-LN(-LN(X)),X)' EVAL » TEVAL
0.577215664748
s:1934.7719              (32 m 15.5 s on the HP 50g)

The actual result is 0.577215664902 ( « 1. Psi NEG» ). I understand your program and peacecalc's are not meant to replace the built-in integrator, but I've read the legacy algorithm is not the best one around (it's only there because it's been tried and tested since the HP-34C, the Saturn code since the HP-71B, I think). Does anyone know of alternative programs that have better performance on this example, that is, faster and more accurate? Thanks!

Gerson.

[HP67]

(04-02-2014 05:32 PM)Thomas Klemm Wrote:  That's easy in a decent editor like emacs:
M-x untabify

Thanks! I have been using Emacs on and off since the 1980s and there are still plenty of tricks for me to learn. This is a particularly good one!

[Thomas Klemm]

(04-04-2014 06:22 AM)HP67 Wrote:  This is a particularly good one!

Quote:Oh yeah! Good ol' C-x M-c M-butterfly ...

[peacecalc]

Hello Gershon,

Quote:but I've read the legacy algorithm is not the best one around (it's only there because it's been tried and tested since the HP-34C, the Saturn code since the HP-71B, I think).

That maybe is true, but it is not fair, because when you use the Romberg method with the iteration depth of 20, the program have to calculate 2^20 - 1 function values, this are at least over a million, no wonder about the time you have to wait...
I don't know wether a side effect occures: so many additions can cause a summation error???

Did you try the lower iteration depths and how the results behavier is?

And your test function is a hard piece of cake because of its limits at 0 and 1, both are divergent.

Greetings
peacecalc

[Gerson W. Barbosa]

(04-04-2014 06:25 PM)peacecalc Wrote:  

Quote:but I've read the legacy algorithm is not the best one around (it's only there because it's been tried and tested since the HP-34C, the Saturn code since the HP-71B, I think).

That maybe is true, but it is not fair, because when you use the Romberg method with the iteration depth of 20, the program have to calculate 2^20 - 1 function values, this are at least over a million, no wonder about the time you have to wait...
I don't know wether a side effect occures: so many additions can cause a summation error???

Did you try the lower iteration depths and how the results behavier is?

Yes, I tried successively 6, 8, 10 and noticed the exponential growth of the running times while the results would be improved only by a little. Without analyzing the code, this story soon crossed my mind :-)

Wheat and chessboard problem

(04-04-2014 06:25 PM)peacecalc Wrote:  And your test function is a hard piece of cake because of its limits at 0 and 1, both are divergent.

I'll try integrating it from 0 to 1/e then from 1/e to 1 and see if the results get better.

In 1982 I wrote what was perhaps the shortest numerical integration program ever, using midpoint rectangles: it fit in 30 steps, on an algebraic calculator. Totally useless though, as only the few built-in functions could be integrated, to only a few decimal places:

http://www.hpmuseum.org/cgi-sys/cgiwrap/...ead=216249 (Message #11)

Regards,

Gerson.

[Thomas Klemm]

(04-03-2014 09:51 PM)Gerson W. Barbosa Wrote:  The actual result is 0.577215664902 ( « 1. Psi NEG» ).

Did you try the substitution \(u=log(x)\) leading to \(\int_{0}^{\infty}-\frac{\log(u)}{e^u}du\)? Or have you started with this and tried to avoid the calculation of an improper integral?

Cheers
Thomas

[Gerson W. Barbosa]

(04-04-2014 08:23 PM)Thomas Klemm Wrote:  

(04-03-2014 09:51 PM)Gerson W. Barbosa Wrote:  The actual result is 0.577215664902 ( « 1. Psi NEG» ).

Did you try the substitution \(u=log(x)\) leading to \(\int_{0}^{\infty}-\frac{\log(u)}{e^u}du\)? Or have you started with this and tried to avoid the calculation of an improper integral?

Cheers
Thomas

Hi Thomas,

No, I just used the integral in the form it was presented by Valentin, at the end of an old thread I started:

http://www.hpmuseum.org/cgi-sys/cgiwrap/...ead=109056

There are a few mistakes in the post, hopefully mostly linguistic ones. I intended to revised it, but I haven't done it yet. Also, at least one link is not working anymore.

Cheers,

Gerson.