Simpson rule for double integral

Question

Compute a quadrature of

$\int_c^d\int_a^b f(x,y)dxdy$

using the Simpson rule and estimate the error.

So the Simpson rule says

$S(f) = (b-a)/6(f(a)+4f((a+b)/2) +f(b))$

So i get

$\int_c^d(b-a)/6(f(a)+4f((a+b)/2) +f(b))dy$

Is that even correct? How do I go on?

Answer 1

Something to get you started. Define an operator: $$ I_x(f) := \int_a^b f(x) \, dx $$ We can see that $I_x(f + g) = I_x(f) + I_x(g)$ and if $\alpha$ does not depend on $x$, $I_x(\alpha f) = \alpha I_x(f)$. In other words the operator is linear.

If it acts on $f(x, y)$ it returns a function of $y$ only.

By analogy define another operator $$ I_y(f) := \int_c^d f(y) \, dy $$

Applying $I_x$ and $I_y$ consequently we get a double integral: $$ I_y(I_x(f)) := \int_c^d \left[ \int_a^b f(x, y) \, dx \right] dy . $$

Now define a "Simpson's rule" operator: $$ S_x(f) := (f(a) + 4 f((a+b)/2) + f(b)) (b - a)/6 $$ It is also linear and if it acts on $f(x, y)$ it "kills" dependency on $x$. The same for $S_y$: $$ S_y(f) := (f(c) + 4 f((c+d)/2) + f(d)) (d - c)/6 $$

By analogy with integral let compute $S_y(S_x(f(x, y)))$ and call it Simpson's rule for double integral. Apply $S_x$ to f(x, y): $$ S_x(f(x, y)) = (f(a, y) + 4 f((a+b)/2, y) + f(b, y)) (b - a)/6 $$ and apply $S_y$ to the result using linearity $$ S_y(S_x(f(x, y)) = (S_y(f(a, y)) + 4 S_y(f((a+b)/2, y)) + S_y(f(b, y))) (b - a)/6 $$ expand all three terms with $S_y$ $$ S_y(S_x(f(x, y)) = \\ \left(16f\left(\frac{b+a}{2},\frac{d+c}{2}\right)+4f\left(\frac{b+a}{2},d\right)+4f\left(\frac{b+a}{2},c\right)+4f\left(b,\frac{d+c}{2}\right) +f\left(b,d\right)+f\left(b,c\right)+4f\left(a,\frac{d+c}{2}\right)+f\left(a,d\right)+f\left(a,c\right)\right) \frac{\left(b-a\right)\left(d-c\right)}{36} $$

Answer 2

A C implementation for applying Simpson's Rule towards solving double integrals can be found here if you are interested. Simpson integration technique for evaluating double integrals

It can be also represented in the following form:

$$ S_x(y_j) = f(x_0, y_j) + f(x_n,y_j) + 4\sum_{i = 1}^{(N_x-2)/2} f(x_{2i-1},y_j) + 2\sum_{i = 1}^{(N_x-2)/2} f(x_{2i},y_j) $$

$$ S = \frac{h_xh_y}{9}[(S_x(y_0) +S_x(y_n) + 4\sum_{j = 1}^{(N_y-2)/2} S_x(y_{2j-1}) + 2\sum_{j = 1}^{(N_y-2)/2} S_x(y_{2j}) $$

Where $$ h_x = \frac{UpperLimit_x - LowerLimit_x}{N_x} $$

and $$ h_y = \frac{UpperLimit_y - LowerLimit_y}{N_y} $$