## Thursday, February 4, 2016

### TI-84 Plus: Simpson’s Rule for f(x,y)

TI-84 Plus: Simpson’s Rule for f(x,y)

If you are familiar with calculus and numerical methods, chances are that you are familiar with Simpson’s Rule.  Normally, Simpson’s Rule applies to functions of one variable, say f(x).  Today, we will apply the Simpson’s Rule to functions of two variables like f(x,y).

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:

Δx * Δy * 1/9 * S

The Programs

The following are three programs, one that set N = 5, N = 9, and the last one lets the user set a value for N.

Programs:

DBLSIMP:  N = 5

Disp "2D SIMPSONS"
Disp "INTEGRAL"
[[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]]→[I]
Input "F(X,Y)=",Y₁
Disp "X"
Prompt A,B
Disp "Y"
Prompt C,D
.25(B-A)→E
.25(D-C)→F
{5,5}→dim([J])
0→I
For(J,1,5)
For(K,1,5)
A+E(K-1)→X
C+F(J-1)→Y
Y₁→[J](J,K)
I+Y₁*[I](J,K)→I
End
End
I*E*F/9→I
Disp "INTEGRAL=",I

DBLS9:  N = 9 Version

Disp "2D SIMPSONS"
Disp "INTEGRAL"
[[1,4,2,4,2,4,2,4,1][4,16,8,16,8,16,8,16,4][2,8,4,8,4,8,4,8,2][4,16,8,16,8,16,8,16,4][2,8,4,8,4,8,4,8,2][4,16,8,16,8,16,8,16,4][2,8,4,8,4,8,4,8,2][4,16,8,16,8,16,8,16,4][1,4,2,4,2,4,2,4,1]]→[I]
Input "F(X,Y)=",Y₁
Disp "X"
Prompt A,B
Disp "Y"
Prompt C,D
(B-A)/8→E
(D-C)/8→F
{9,9}→dim([J])
0→I
For(J,1,9)
For(K,1,9)
A+E(K-1)→X
C+F(J-1)→Y
Y₁→[J](J,K)
I+Y₁*[I](J,K)→I
End
End
I*E*F/9→I
Disp "INTEGRAL=",I

DBLSMP: N can be set

Disp "2D SIMPSONS"
Disp "INTEGRAL"
Input "F(X,Y)=",Y₁
Disp "X"
Prompt A,B
Disp "Y"
Prompt C,D

Input "N (ODD)=",N
{1,N}→dim([I])
1→[I](1,1)
1→[I](1,N)
For(J,2,N-1)
If fPart(J/2)=0
Then
4→[I](1,J)
Else
2→[I](1,J)
End
End
[I]^T*[I]→[I]

(B-A)/(N-1)→E
(D-C)/(N-1)→F

{N,N}→dim([J])
0→I
For(J,1,N)
For(K,1,N)
A+E(K-1)→X
C+F(J-1)→Y
Y₁→[J](J,K)
I+Y₁*[I](J,K)→I
End
End
I*E*F/9→I
Disp "INTEGRAL=",I

One thing to note, when you are prompted enter your function, do not forget to use quotes.  Example, for x^2*y, enter “X^2*Y”.

Caution:  Simpson’s Rule is an approximation, 100% accuracy is not guaranteed.

Of the three I recommend DBLSMP so you have control over N.

Examples and Results

∫ ∫ x^2*y^3 dx dy, A = 0, B = 2, C = 1, D = 5
Result is 416 for both N = 5 and N = 9, agrees with the actual answer

∫ ∫ x^2 + y dx dy, A = 0, B = 2, C = 1, D = 1
Result is 5.666666667 for both N = 5 and N = 9, agrees with the actual answer

∫ ∫ sin(x*y) dx dy, A = 0, B = π, C = 0, D = 0.5
N = 5, result is 0.5568449485
N = 9, result is 0.5568006343

To illustrate that Simpson’s Rule may not work out the best when N is low:
∫ x * e^(-y) dx dy,  A = 0, B =5, C = -2, D = 2
N = 5, result is 91.12098525
N = 9, result is 90.70208034

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

Until next time,
Eddie

This blog is property of Edward Shore.  2016