Compact Simpson's Rule
12-13-2015, 02:06 AM (This post was last modified: 12-13-2015 02:12 AM by Namir.)
Post: #1 Namir Senior Member Posts: 810 Joined: Dec 2013
Compact Simpson's Rule
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
12-13-2015, 07:13 AM
Post: #2
 Thomas Klemm Senior Member Posts: 1,447 Joined: Dec 2013
RE: Compact Simpson's Rule
(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
12-13-2015, 02:20 PM (This post was last modified: 12-13-2015 02:26 PM by Dieter.)
Post: #3
 Dieter Senior Member Posts: 2,397 Joined: Dec 2013
RE: Compact Simpson's Rule
(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
 « Next Oldest | Next Newest »

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