# Solving a double integral in Matlab

I need to solve this integral:

$M = \rho h \int_0^a \int_0^b N^2 \ dx dy$

where $N$ is defined as:

$N = sin\left(\dfrac{\pi}{4}\eta +\dfrac{3\pi}{4}\right)sin\left(\dfrac{\pi}{4}\eta +\dfrac{3\pi}{4}\right)sin\left(\dfrac{\pi}{4}\psi +\dfrac{3\pi}{4}\right)sin\left(\dfrac{\pi}{4}\psi +\dfrac{3\pi}{4}\right)$

where $\eta = y/b$ and $\psi = x/a$ (change of variable)

Therefore the integral becomes:

$M = \rho h a b \int_0^1 \int_0^1 \biggl[sin^2\left(\dfrac{\pi}{4}\eta +\dfrac{3\pi}{4}\right)sin^2\left(\dfrac{\pi}{4}\psi +\dfrac{3\pi}{4}\right) \biggr]^2 \ d\eta d\psi$

To solve it I used the Matlab function trapz() in this way:

eta = 0:0.01:1;
psi = 0:0.01:1;
rho = 2770; h = 0.0012; a = 0.127; b = 0.0508;
fun = rho*h*a*b.*(sin(pi/4.*psi+3*pi/4).*sin(pi/4.*psi+3*pi/4) ...
.*sin(pi/4.*eta+3*pi/4).*sin(pi/4.*eta+3*pi/4)).^2;
M = trapz(fun);


$M = 0.0175$

But if I set

eta = 0:0.001:1;
psi = 0:0.001:1;


then $M = 0.1754$

What am I doing wrong?

I also tried using the quad2d() function:

fun = @(eta,psi)  rho*h*a*b.*(sin(pi/4.*psi+3*pi/4).*sin(pi/4.*psi+3*pi/4) ...
.*sin(pi/4.*eta+3*pi/4).*sin(pi/4.*eta+3*pi/4)).^2;


But this time $M = 6.8920\times10^{-5}$

• What is your actual question? Is there a reason you wan to use trapz instead of integral2 or quad2d? trapz performs 1-D numerical integration, so you'll need to apply it more than once in your case – see the documentation. Jan 25 '15 at 16:28
• The main problem is that I do not know how to do it. First I tried with trapz() but I felt it was not the right method, therefore I tried with quad2d().
– Rhei
Jan 25 '15 at 18:58
• Anyway I do not understand why the value of the integral using trapz() changes of one order of magnitude when I change the $d\eta$ and the $d\psi$
– Rhei
Jan 25 '15 at 18:59
• Sorry, I meant $\Delta \eta$ and $\Delta \psi$
– Rhei
Jan 25 '15 at 19:08
• I tried the code written in the documentation and I obtained the same result given by quad2d . I should look at the online documentation more often, the help doc of my matlab version is not so complete
– Rhei
Jan 26 '15 at 6:38

Your function $N$ is separable and the limits of integration are constants so this can easily be computed as the product of two one dimensional integrals. In fact those integrals are identical so it is just the square of a 1d integral.

\begin{align} M &= \rho h a b \int_0^1 \int_0^1 \left[sin^2\left(\dfrac{\pi}{4}\eta +\dfrac{3\pi}{4}\right)sin^2\left(\dfrac{\pi}{4}\psi +\dfrac{3\pi}{4}\right) \right]^2 \ d\eta d\psi\\ &= \rho h a b \int_0^1 sin^4\left(\dfrac{\pi}{4}\eta +\dfrac{3\pi}{4}\right) d\eta \int_0^1 sin^4\left(\dfrac{\pi}{4}\psi +\dfrac{3\pi}{4}\right) d\psi\\ &= \rho h a b \left[\int_0^1 sin^4\left(\dfrac{\pi}{4}\eta +\dfrac{3\pi}{4}\right) d\eta\right]^2 \end{align}

This integral can be done analytically (I used WolframAlpha), giving $$M = \rho h a b \left( \dfrac{3}{8}-\dfrac{1}{\pi} \right)^2\approx 6.8920\times10^{-5}$$ if you use $\rho = 2770$, $h = 0.0012$, $a = 0.127$, $b = 0.0508$.

In general, if the function is separable but the integral can't be solved analytically it will be much more efficient to compute the product of two one-dimensional integrals than to directly compute the two-dimensional integral.

Since the function to be integrated is dependent on 2 parameters: $\eta$ and $\psi$, it has to be integrated as a 2D function in this way:

x = 0:0.001:1;
y = x;
[eta,csi] = meshgrid(x,y);
fun = rho*h*a*b.*(sin(pi/4.*csi+3*pi/4).*...
...sin(pi/4.*csi+3*pi/4).*sin(pi/4.*eta+3*pi/4).*sin(pi/4.*eta+3*pi/4)).^2;
M = trapz(y,trapz(x,fun,2));


and this approach gives the same result obtained using the quad2d function.

fun = @(eta,psi)  rho*h*a*b.*(sin(pi/4.*psi+3*pi/4).*sin(pi/4.*psi+3*pi/4).*...
...sin(pi/4.*eta+3*pi/4).*sin(pi/4.*eta+3*pi/4)).^2;

The main difference between the 2 methods is that in the quad2d the function to be integrated has to be defined as an function of the 2 parameters (in this case it has been defined as a function_handle).
\[Rho] h a b Integrate[Sin[\[Eta] \[Pi]/4 + 3 \[Pi]/4]^4  Sin[\[Psi] \[Pi]/4 + 3 \[Pi]/4]^4, {\[Eta], 0, 1}, {\[Psi], 0, 1}]

$\frac{(8-3 \pi )^2 a b h \rho }{64 \pi ^2}$