+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Software Libraries (/forum-10.html)
+--- Forum: HP-41C Software Library (/forum-11.html)
+--- Thread: Improved Simpson's rule (/thread-5349.html)
Improved Simpson's rule - Dieter - 12-16-2015 07:12 PM
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 RTNInsert 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,693147181In 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 approximationsAs 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. #-)
RE: Improved Simpson's rule - rprosperi - 12-17-2015 12:52 AM
(12-16-2015 07:12 PM)Dieter Wrote: ...
The following code has not seen much testing. In fact I wrote it this morning during breakfast.
...
Dieter
Wow, meals in my house and yours are quite different.
I was just about to go eat dinner, but now I will feel I'm certainly being lazy, limiting my effort to dining. And maybe reading the daily mail...
