Improved Simpson's rule
12-16-2015, 07:12 PM (This post was last modified: 01-12-2018 11:03 AM by Dieter.)
Post: #1
 Dieter Senior Member Posts: 2,397 Joined: Dec 2013
Improved Simpson's rule
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. #-)
 « Next Oldest | Next Newest »

 Messages In This Thread Improved Simpson's rule - Dieter - 12-16-2015 07:12 PM RE: Improved Simpson's rule - rprosperi - 12-17-2015, 12:52 AM

User(s) browsing this thread: 1 Guest(s)