(71B) Simpson’s Rule Approximation for f(x,y) - Printable Version +- HP Forums (https://www.hpmuseum.org/forum) +-- Forum: HP Software Libraries (/forum-10.html) +--- Forum: General Software Library (/forum-13.html) +--- Thread: (71B) Simpson’s Rule Approximation for f(x,y) (/thread-7511.html) |
(71B) Simpson’s Rule Approximation for f(x,y) - Eddie W. Shore - 01-04-2017 04:26 AM integral( integral( f(x,y) dx from A to B) dy from C to D) 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. Code: 5 DESTROY I,J,A,B,C,D 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 http://www.physics.usyd.edu.au/teach_res/mp/doc/math_integration_2D.pdf Retrieved January 30, 2016 Thank you, and wishing you a happy, healthy, and successful 2017! |