Tuesday, January 3, 2017

HP 71B: Simpson’s Rule Approximation for f(x,y)

HP 71B: Simpson’s Rule Approximation for f(x,y)

Happy New Year everyone!  Hope everyone had a great New Years Day!  J



Today, I will present apply the Simpson’s Rule to functions of two variables like f(x,y) for the HP 71B. Here is the general algorithm:

Start with the function f(x,y) and your integration limits A, B, C, and D.

Then determine Δx and Δy (labeled E and F in the programs below) by:

Δx = (B – A)/(N – 1)
Δy = (D – C)/(N – 1)

Where N is the number of partitions.  Unlike the Simpson’s Rule for one variable, in this case N must be odd.  Generally, the higher N is, the more accurate the approximation is, with the expense of additional computational time.

Next, build a matrix, let’s say [ I ].   This is your Simpson’s Matrix.  The Simpson’s Matrix is built by the expression

[ I ] = [1, 4, 2, 4, 2, 4, 2, 4, 2, …, 4, 1]^T * [1, 4, 2, 4, 2, 4, 2, 4, 2, …, 4, 1]

The length of the vector used to determine [ I ] depends on N.  A way to build it by the routine:

Store 1 into the element 1 of [ I ]  (first element)
Store 1 into the element N of [ I ]  (last element)
For J from 2 to N – 1
If J is divisible by 2, then store 4 in the jth element of [ I ],
Else store 2 in the jth element of [ I ]

For N = 5, the vector would be built is [1, 4, 2, 4, 1]

And

[ I ] = [1,4,2,4,1]^T * [1,4,2,4,1] =

[[1, 4, 2, 4, 1]
[4, 16, 8, 16, 4]
[2, 8, 4, 8, 2]
[4, 16, 8, 16, 4]
[1, 4, 2, 4, 1]]

Build another matrix [ J ].  The elements are determined by the following formulas:

For row j and column k, the element is f(A + Δx*(j – 1),  C + Δy*(k – 1))

Once finished, multiply every element of [ I ] by [ J ].  This is NOT matrix multiplication.   Then sum all of the elements of the results.  In essence:

S = ∑ (j = 1 to N) ∑ (k = 1 to N)   [I](j,k)* [J](j,k)

Determine the final integral approximation as:

Integral = Δx * Δy * 1/9 * S

The program DBLSIMP uses N = 5.  This program works best for f(x,y) where they are polynomials.  On the HP 71B, matrices cannot be typed directly, elements have to be stored and recalled on element at a time.  The program presented does not use modules. 

HP 71B Program DBLSIMP
At least 580 Bytes

Edit f(x,y) at line 10. Use variables X and Y.

5 DESTROY I,J,A,B,C,D
8 DESTROY E,F,K,S,X,Y
10 DEF FNF(X,Y)= [ enter f(X,Y) here ]
12 RADIANS
14 DIM I(5,5)
20 I(1,1) = 1
21 I(1,2) = 4
22 I(1,3) = 2
23 I(1,4) = 4
24 I(1,5) = 1
25 I(2,1) = 4
26 I(2,2) = 16
27 I(2,3) = 8
28 I(2,4) = 16
29 I(2,5) = 4
30 I(3,1) = 2
31 I(3,2) = 8
32 I(3,3) = 4
33 I(3,4) = 8
34 I(3,5) = 2
35 I(4,1) = 4
36 I(4,2) = 16
37 I(4,3) = 8
38 I(4,4) = 16
39 I(4,5) = 4
40 I(5,1) = 1
41 I(5,2) = 4
42 I(5,3) = 2
43 I(5,4) = 4
44 I(5,5) = 1
50 DISP “X: from a to b” @ WAIT 1
52 INPUT “a = “; A
54 INPUT “b = “; B
56 DISP “Y: from c to d” @ WAIT 1
58 INPUT “c = “; C
60 INPUT “d = “; D
62 E = .25 * (B – A)
64 F = .25 * (D – C)
66 S = 0
70 FOR J = 1 TO 5
72 FOR K = 1 TO 5
74 X = A + E * (K – 1)
76 Y = C + F * (J – 1)
78 S = S + FNF(X,Y) * I(J,K)
80 NEXT K
82 NEXT J
90 S = S * E * F/9
95 DISP “INTEGRAL = “ @ WAIT 1
97 DISP S

Examples:

F(X,Y) = 2*Y – 3*X
A = 1, B = 2, C = 2, D = 5
Result:  Integral = 7.5
(Actual answer: 7.5)

F(X,Y) = X^2/Y^2
A = 1, B = 2, C = 2, D = 5
Result:  Integral ≈ 0.70212579101
(Actual answer:  0.7)

F(X,Y) = 0.5*X*EXP(Y)
A = 1, B = 2, C = 2, D = 5
Result:  Integral ≈ 105.942243008
(Actual answer is about 105.768077253)

Source:
Cooper, Ian.  “Doing Physics With Matlab:  Mathematical Routines”  School of Physics, University of Sydney
Retrieved January 30, 2016

Thank you, and wishing you a happy, healthy, and successful 2017!
Eddie


This blog is property of Edward Shore, 2017

No comments:

Post a Comment

HP Prime: Area by Quadratic Splines

HP Prime:  Area by Quadratic Splines Introduction The program QUADSUM calculates the area under the curve described by the set o...