+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Software Libraries (/forum-10.html)
+--- Forum: HP-41C Software Library (/forum-11.html)
+--- Thread: Compact Simpson's Rule (/thread-5324.html)
Compact Simpson's Rule - Namir - 12-13-2015 02:06 AM
This listing calculates the integral for f(X)=1/X (see label E) using Simpson's rule. The code is based on the compact pseudo-code that I posted in a thread in the General area.
Code:
Memory Map
========
R00 = Sum
R01 = A
R02 = B
R03 = h
R04 = 2 or 4
R05 = N
Listing
=======
01 LBL "SIMPS"
02 LBL A
03 "A/^B?"
04 PROMPT
05 STO 02 # store upper limit for integral
06 X<>Y
06 STO 01 # store lower limit for integral
08 "N?"
09 PROMPT
10 STO 05 # store number of divisions. Must be even.
11 RCL 02
12 RCL 01
13 -
14 X<>Y
15 /
16 STO 03 # store h
17 RCL 01
18 XEQ E # Calculate f(A)
19 STO 00 # sum = f(A)
20 RCL 02
21 XEQ E # Calculate f(B)
22 ST- 00 # sum = f(A) - f(B) ... this is part of a math trick
23 RCL 03
24 ST+ 01 # A = A + h
25 4
26 STO 04
27 LBL 00 # start loop
28 RCL 03
29 XEQ E # Calculate f(A)
30 RCL 04
31 * # multiply by 4 or 2
32 ST+ 00 # sum = sum + ...
33 6
34 RCL 4
35 -
36 STO 04 # toggle between 4 and 2
37 RCL 03
38 ST+ 01 # A = A + h
39 DSE 05 # N = N - 1 and test N=0
40 GTO 00 # end of loop
41 RCL 00
42 RCL 03
43 *
44 3
45 /
46 RTN
47 LBL E # Calculate f(x). You can edit code after this label.
49 1/X
50 RTNRE: Compact Simpson's Rule - Thomas Klemm - 12-13-2015 07:13 AM
(12-13-2015 02:06 AM)Namir Wrote:Code:
30 RCL 04
31 * # multiply by 4 or 2
32 ST+ 00 # sum = sum + ...
33 6
34 RCL 4
35 -
36 STO 04 # toggle between 4 and 2
You could use:
Code:
30 6 # toggle between 4 and 2
31 X<> 04
32 ST- 04
33 * # multiply by 4 or 2
34 ST+ 00 # sum = sum + ...Cheers
Thomas
RE: Compact Simpson's Rule - Dieter - 12-13-2015 02:20 PM
(12-13-2015 02:06 AM)Namir Wrote: This listing calculates the integral for f(X)=1/X (see label E) using Simpson's rule. The code is based on the compact pseudo-code that I posted in a thread in the General area.
Here's my version with yet another method of toggling between 2* and 4*f(x): use a flag. This also frees up one register, and even one more can be saved as once h is calculated the value of b is no longer required. So that's two registers less and some bytes of code saved. And one function call as f(b) is not calculated twice.
Code:
01 LBL "SIMP"
02 LBL A
03 " A^B=?"
04 PROMPT
05 STO 02 ' store upper limit for integral
06 X<>Y
07 STO 01 ' store lower limit for integral
08 -
09 STO 03 ' prestore b-a
10 RCL 01
11 XEQ E ' calculate f(a)
12 STO 00 ' sum = f(a)
13 RCL 02
14 XEQ E ' calculate f(b)
15 ST+ 00 ' sum = f(a) + f(b) ... no trick required ;-)
16 " N=?"
17 PROMPT
18 STO 02 ' store number of divisions. Must be even.
19 ST/ 03 ' calculate and store h = (b-a)/n
20 SF 02
21 DSE 02 ' n = n - 1
22 LBL 00 ' start loop
23 RCL 03
24 ST+ 01 ' x := x + h
25 RCL 01
26 XEQ E ' calculate f(x)
27 ST+ X
28 FS? 02
29 ST+ X
30 ST+ 00 ' sum = sum + [2|4]*f(x)
31 FC?C 02 ' toggle flag 2
32 SF 02
33 DSE 02 ' n := n - 1
34 GTO 00 ' if n > 0 do another loop
35 RCL 00
36 RCL 03
37 *
38 3
39 /
40 RTN
41 LBL E ' Calculate f(x).
42 1/x ' Insert your f(x) code here.
43 RTN
Memory Map
----------
R00 = sum
R01 = a, x
R02 = b, n, counter
R03 = hMaybe t's also fun to watch the flashing "2" annunciator. ;-)
Dieter