HP 50g Romberg Integration

[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 }