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

## No comments:

## Post a Comment