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
Actual answer is about 0.556797718751
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
Actual answer is about 90.6715102
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