Improved Simpson's rule

[Dieter]

Recently there were some discussions regarding a compact implementation of Simpson's rule. I would like to add an improved version that calculates a new, usually more exact approximation from two Simpson approximations with n and 2n intervals. This way the user also gets an estimate for the approximation error.

The following code has not seen much testing. In fact I wrote it this morning during breakfast. ;-) So take care and please report all errors you find.

Code:

01  LBL"SIMP"
02  "LIMITS A↑B?"
03  PROMPT
04  STO 02    ' save b temporarily in R02
05  X<>Y
06  STO 01    ' save a
07  -
08  2
09  STO 04    ' n = 2
10  /
11  STO 03    ' h = (b-a)/2
12  RCL 01
13  XEQ 99    ' f(a)
14  X<> 02
15  XEQ 99    ' f(b)
16  ST+ 02    ' save f(a)+f(b) in R02
17  CLX
18  STO 07    ' s2 = 0
19  RCL 01
20  RCL 03
21  +
22  XEQ 99    ' f((a+b)/2)
23  STO 06    ' =: s4
24  4
25  *
26  RCL 02
27  +
28  RCL 03
29  *
30  3
31  /
32  STO 08    ' Simpson approx. for n=2
33  SF 01     ' set "first iteration" flag
34  LBL 01
35  CLX
36  X<> 06    ' new s4 = 0
37  ST+ 07    ' new s2 = old s2 + old s4
38  RCL 03    ' R03 holds 2h
39  2
40  /         ' = h
41  RCL 01
42  +
43  STO 00    ' first x = a + h
44  RCL 04
45  STO 05    ' initialize counter
46  LBL 02
47  RCL 00
48  XEQ 99    ' f(x)
49  ST+ 06    ' s4 = s4 + f(x)
50  RCL 03
51  ST+ 00    ' x = x + 2h
52  DSE 05
53  GTO 02    ' end of loop
54  RCL 08    ' recall previous Simpson approximation
55  RCL 06
56  ST+ X     ' 2*s4
57  RCL 07
58  +         ' 2*s4 + s2
59  ST+ X
60  RCL 02
61  +         ' 4*s4 + 2*s2 + f(a) + f(b) ...
62  RCL 03
63  *
64  6
65  /         ' ... * 2h/6
66  STO 08    ' = new Simpson approximation
67  FS?C 01   ' first iteration?
68  STO 09    ' then save this as previous estimate
69  16
70  *
71  X<>Y
72  -
73  15        ' calculate improved extrapolated value
74  /         ' from new and previous Simpson approximation
75  X<> 09
76  RCL 09    ' return old and new extrapolations in Y and X
77  STOP      ' and stop
78  2
79  ST/ 03    ' h := h/2
80  ST* 04    ' n := 2 n
81  GTO 01    ' start next approximation with 2 n nodes

82  LBL 99    ' f(x) starts here
..  ...       ' If required, R00 can be used to store X
nn  RTN

Insert your function at LBL 99 and start "SIMP". The program returns a sequence of approximations, with the latest one in X and the previous in Y.
So pressing [x<>y] switches between the current and previous approximation.

Example: Integrate f(x) = 1/x  from a=1 to b=2.

Code:

[PRGM]
82  LBL 99
83  1/x
84  RTN
[PRGM]

[FIX] 9

[XEQ]"SIMP"
              LIMITS A↑B?
1 [ENTER] 2
  [R/S]       0,693174603
  [R/S]       0,693147901
  [R/S]       0,693147195
  [R/S]       0,693147181
  [R/S]       0,693147181

In a way this method is similar to the Romberg scheme in that it calculates an improved approximation from a quite good estimate of the error associated with Simpson's rule. On the other hand no matrix has to be stored and the program requires only ten data registers (R00...R09):

Code:

Used registers:
R00  x
R01  a
R02  b, f(a)+f(b)
R03  2h
R04  n/2
R05  loop counter
R06  s4 = sum of f(x) weighted with factor 4
R07  s2 = sum of f(x) weighted with factor 2
R08  current Simpson approximation
R09  improved result calculated from two last Simpson approximations

As usual, all comments and corrections are welcome.

Dieter

Edit: more than two years after posting this I realized that the example does not integrate ln(x) but 1/x. #-)