计算科学 - 在 Matlab 中求解二重积分 - 吾爱随笔录

在 Matlab 中求解二重积分

计算科学 matlab

2021-12-07 13:03:36

我需要解决这个积分：

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

其中定义为： $N$

$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)$

其中和（变量的变化） $\eta = y/b$ $\psi = x/a$

因此积分变为：

$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$

为了解决它，我以这种方式使用了 Matlab 函数trapz()：

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$

但是如果我设置

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

那么 $M = 0.1754$

我究竟做错了什么？

我也尝试使用该quad2d()功能：

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;
M = quad2d(fun,0,1,0,1);

但是这次 $M = 6.8920\times10^{-5}$

3个回答

你的功能 $N$ 是可分的，并且积分的极限是常数，因此可以很容易地将其计算为两个一维积分的乘积。事实上，这些积分是相同的，所以它只是一维积分的平方。

\begin{aligned} M & = ρ h a b \int_{0}^{1} \int_{0}^{1} {[s i n^{2} (\frac{π}{4} η + \frac{3 π}{4}) s i n^{2} (\frac{π}{4} ψ + \frac{3 π}{4})]}^{2} d η d ψ \\ = ρ h a b \int_{0}^{1} s i n^{4} (\frac{π}{4} η + \frac{3 π}{4}) d η \int_{0}^{1} s i n^{4} (\frac{π}{4} ψ + \frac{3 π}{4}) d ψ \\ = ρ h a b {[\int_{0}^{1} s i n^{4} (\frac{π}{4} η + \frac{3 π}{4}) d η]}^{2} \end{aligned}

$\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}$

这个积分可以解析完成（我使用了WolframAlpha），给出

M = ρ h a b {(\frac{3}{8} - \frac{1}{π})}^{2} \approx 6.8920 \times 10^{- 5}

$M = \rho h a b \left( \dfrac{3}{8}-\dfrac{1}{\pi} \right)^2\approx 6.8920\times10^{-5}$ 如果你使用

ρ = 2770

$\rho = 2770$ ,

h = 0.0012

$h = 0.0012$ ,

a = 0.127

$a = 0.127$ ,

b = 0.0508

$b = 0.0508$ .

一般来说，如果函数是可分的，但积分不能解析求解，那么计算两个一维积分的乘积比直接计算二维积分要高效得多。

由于要集成的功能取决于 2 个参数： $\eta$ 和 $\psi$ ，它必须以这种方式集成为 2D 函数：

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));

quad2d这种方法给出了使用该函数获得的相同结果。

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;
M = quad2d(fun,0,1,0,1);

两种方法之间的主要区别在于，在quad2d要集成的函数中，必须将其定义为 2 个参数的函数（在这种情况下，它已定义为 a function_handle）。

符号解决方案直接在Mathematica中：

\[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}$

其它你可能感兴趣的问题