# How to approximate a double integral numerically using the Simpson's Rule

.I'm trying to approximate the integral :

$$\int_{0}^{1} \int_{0}^{1} \ln(e^x+e^y+x+y)\, dx\,dy$$ with an error $$\leq 10^{-4}$$. I tried to use the Simpson's rule, and trapezoidal rule.

But I had difficulty with this, since in these rules I have to find $$n$$ where $$h=\frac{b-a}{n}$$. In these two rules we have a formula for the error, that includes a derivative of $$f$$ in a point (that belongs to $$[a,b]$$). But I could manage to how calculate $$n$$.

It is the first time for me to calculate this, so I apologize if I was not accurate.

• It is $$\int_{0}^{1}f(y)\,dy$$ where $$f(y) = \int_{0}^{1}\ln(e^{x} + e^{y} + x + y)\,dx.$$ You can differentiate $f(y)$ by interchanging the order of derivative and integration, though there is still the problem of calculating $f(y)$ accurately. Dec 7, 2020 at 23:34
• Can you show me please how to calculate it, I'm not getting it @Mason Dec 8, 2020 at 10:22
• Yes these methods are defined like that @moo Dec 8, 2020 at 13:56
• I am trying a different method, assuming that my N=M=1 , and substitue in Sympson's Rule which also contain the error , because the first approch with finding N was too tiring with all derivatives...and in addition the error formula is connected with too constants not only one, N and M Dec 9, 2020 at 11:31

The double integral Composite Simpson's error term is given by

$$E = \dfrac{-(d-c)(b-a)}{180}\left[h^4 ~\mbox{max} \left(\dfrac{\partial^4f}{\partial x^4}\right) + k^4 ~\mbox{max} \left(\dfrac{\partial^4f}{\partial y^4}\right) \right]$$

For this particular problem, we will assume that we use an equal step size, that is $$h = k = \frac{b-a}{n} = \frac{1-0}{n} = \frac{1}{n}$$, which makes the error

$$E = \dfrac{-(d-c)(b-a)}{180 ~n^4}\left[~\mbox{max} \left(\dfrac{\partial^4f}{\partial x^4}\right) + ~\mbox{max} \left(\dfrac{\partial^4f}{\partial y^4}\right) \right]$$

Our function is $$f(x, y) = \ln(e^x+e^y+x+y)$$, so we have

$$\dfrac{\partial^4f}{\partial x^4} = -\dfrac{6 \left(e^x+1\right)^4}{\left(x+e^x+e^y+y\right)^4}+\dfrac{12 e^x \left(e^x+1\right)^2}{\left(x+e^x+e^y+y\right)^3}-\dfrac{4 e^x \left(e^x+1\right)}{\left(x+e^x+e^y+y\right)^2}+\dfrac{e^x}{x+e^x+e^y+y}-\dfrac{3 e^{2 x}}{\left(x+e^x+e^y+y\right)^2}$$

$$\dfrac{\partial^4f}{\partial y^4} = -\dfrac{6 \left(e^y+1\right)^4}{\left(x+e^x+e^y+y\right)^4}+\dfrac{12 e^y \left(e^y+1\right)^2}{\left(x+e^x+e^y+y\right)^3}-\dfrac{4 e^y \left(e^y+1\right)}{\left(x+e^x+e^y+y\right)^2}+\dfrac{e^y}{x+e^x+e^y+y}-\dfrac{3 e^{2 y}}{\left(x+e^x+e^y+y\right)^2}$$

Using a numerical solver, the max of $$\dfrac{\partial^4f}{\partial x^4}$$ and $$\dfrac{\partial^4f}{\partial y^4}$$ over the region $$0 \le x \le 1, 0 \le y \le 1$$ is $$-0.0186659$$.

We want an $$E \le 10^{-4}$$, and have

$$E = -\dfrac{1}{180~ n^4}\left(-0.0186659 -0.0186659\right) = \frac{0.000207399}{n^4} \le 10^{-4}$$

This results in $$n \ge 1.2$$, and because we need an integer $$n$$, we choose $$n = 2$$, but this is not very satisfying as maybe we need to choose $$n = 3$$ or more.

Let's check this result, using Composite Simpson with $$n = 2$$, we get an error of $$|1.4562429599654196 - 1.4573986597429938| = 0.0011557 = 1.1157 \times 10^{-3}$$

That was pretty close, but if we had chosen $$n = 3$$, we have $$E = 4.63677 \times 10^{-4}$$, which aligns with the ask much better.

I think what this question is asking is for you to use $$n = 2, 3, \ldots$$ until you get a result that has a smaller error than required.

Update

This is the 2D-Simpson Formula I used above for the case $$n = m = 2$$

$$h=\frac{b-a}{2 m}$$

$$k=\frac{b-a}{2 n}$$

$$X_{\text{i\_}}=a+h i$$

$$Y_{\text{j\_}}=a+j k$$

$$\int_{a}^{b} \int_{c}^{d} f(x, y)\, dx\,dy = \frac{1}{9} h k \left(4 \sum _{j=1}^n f\left(a,Y_{2 j-1}\right)+2 \sum _{j=1}^{n-1} f\left(a,Y_{2 j}\right)+4 \sum _{j=1}^n f\left(b,Y_{2 j-1}\right)+2 \sum _{j=1}^{n-1} f\left(b,Y_{2 j}\right)+4 \sum _{i=1}^m f\left(X_{2 i-1},c\right)+2 \sum _{i=1}^{m-1} f\left(X_{2 i},c\right)+4 \sum _{i=1}^m f\left(X_{2 i-1},d\right)+2 \sum _{i=1}^{m-1} f\left(X_{2 i},d\right)+16 \sum _{j=1}^n \left(\sum _{i=1}^m f\left(X_{2 i-1},Y_{2 j-1}\right)\right)+8 \sum _{j=1}^{n-1} \left(\sum _{i=1}^m f\left(X_{2 i-1},Y_{2 j}\right)\right)+8 \sum _{j=1}^n \left(\sum _{i=1}^{m-1} f\left(X_{2 i},Y_{2 j-1}\right)\right)+4 \sum _{j=1}^{n-1} \left(\sum _{i=1}^{m-1} f\left(X_{2 i},Y_{2 j}\right)\right)+f(a,c)+f(a,d)+f(b,c)+f(b,d)\right).$$

• I would note that my k and h are defined as you mentioned but divided by 2, will it change anything? Dec 9, 2020 at 17:52
• Indeed! Wishing us all well! So far so good re: staying healthy! Jan 1, 2021 at 19:09