Compact Simpson's Rule

[Namir]

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 RTN

[Thomas Klemm]

(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

[Dieter]

(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 = h

Maybe t's also fun to watch the flashing "2" annunciator. ;-)

Dieter