在matlab中数值求解双积分答案

【问题标题】：Solving double integral numerically in matlab在matlab中数值求解双积分
【发布时间】：2013-10-21 19:26:23
【问题描述】：

在论文"The fractional Laplacian operator on bounded domains as a special case of the nonlocal diffusion operator"。作者将有界域上的分数拉普拉斯方程求解为非局部扩散方程。

我正在尝试在matlab中实现一维问题的有限元逼近（请参阅上述论文的第14页）。

我正在使用 $\phi_k$ 的以下定义，因为论文中提到 $\phi$ 是$hat\;function$ \begin{equation} \phi_{k}(x)=\begin{cases} {x-x_{k-1} \over x_k\,-x_{k-1}} & \mbox{ if } x \in [x_{k-1},x_k], \\ {x_{k+1}\,-x \over x_{k+1}\,-x_k} & \mbox{ if } x \in [x_k,x_{k+1}], \\ 0 & \mbox{ otherwise},\end{cases} \end{equation}

$\Omega=(-1,1)$ 和 $\Omega_I=(-1-\lambda,-1) \cup (1,1+\lambda)$ 使得 $\Omega\cup\Omega_I= (-1-\lambda,1+\lambda)$

对于整数 K,N，我们将 $\overline{\Omega\cup\Omega_I}=[-1-\lambda,1+\lambda]$ 的划分定义为，

\开始{方程} -1-\lambda=x_{-K}<...>

最后，对于某些系数 $U_j$，我们必须解出的方程 $\tilde{u_N}=\sum_{i=-K}^{K+N}U_j\phi_j(x)$ 是：

其中 $i=1,...,N-1$。

我需要指针来简化和解决matlab中的LHS双积分。在论文（第15页）中写到我应该使用四点高斯正交进行内积分，使用quadgk.m函数进行外积分，但是由于内部积分的限制是 x 我怎么能在它上面应用四点高斯正交？？任何帮助将不胜感激。谢谢。

你可以找到原问题here。（由于SO不支持Latex）

【问题讨论】：

建议 math.stackexchange.com 代替（是的，我知道，这是编程，但你可能有更多的运气。）
我建议你把这个问题分成两部分：1）要求简化 mathoverflow，2）要求在 SO 上实现。在 SO 上，除了双积分（我为你放了一张图片）之外的所有数学都省略了。在 mathoverflow 上，省略在 MATLAB 中实现的请求。

标签： algorithm matlab math

【解决方案1】：

要首次尝试解决问题，请查看dblquad 和/或quad2d。

最后，您将需要自定义正交方法，因此您应该执行以下操作：

% The integrand is of course a function of both x and y
integrand = @(x,y) (phi_j(y) - phi_j(x))*(phi_i(y) - phi_i(x))/abs(y-x)^(2*s+1);

% The inner integral is a function of x, and integrates over y
inner = @(x) quadgk(@(y)integrand(x,y), x-lambda, x+lambda);

% The inner integral is integrated over x to yield the value of the double integral 
dblIntegral = quadgk(inner, -(1+lambda), 1+lambda)

我已经使用了两次quadgk，但您可以使用任何其他（自定义）求积方法替换。

顺便问一下——作者提出（非自适应）4 点高斯方法的原因是什么？这样，您就无法估计（和/或控制）内部积分中的误差...

【讨论】：

作者一定建议在分区上使用4点高斯求积以获得良好的精度。我还给她写了关于求解积分的信，她回答说：求解系统矩阵中的积分是解决数值非局部问题的难点。这个想法是为内部和外部积分实现 2 个不同的求积规则。
内积分：这个从 x-epsilon 到 x+epsilon。由于我们想避免 y=x 处的奇点，我们将积分拆分为 (x-epsilon,x) 和 (x,x+epsilon)。然后，采取例如(x-epsilon,x) 并根据分区中的元素对其进行拆分。你会有类似 (x-epsilon, x_j) (x_j, x_{j+1}) (x_{j+1}, x_{j+2}) (x_{j+2}, x) 对于 x_j , x_{j+1} 和 x_{j+2} 在 (x-epsilon, x) 中。在这些区间中，您可以使用高斯求积法则（除非您想得到非常准确的解，否则您应该使用自适应求积法则）。
外积分：这里需要更准确的规则。您必须对分区的每个元素进行积分，即 (x_i, x_{i+1}) 并且应该使用自适应正交规则（如 Matlab 中的 quadgk）。被积函数是内积分的结果。
@VikasGautam 我明白了。正如作者所说，我怀疑在内积分上的自适应求积不会花费太多性能，但会为您提供更多信息（误差估计）。但老实说，我对问题的细节完全不熟悉，听起来作者确实已经投入了很多的摆弄和调整；我不会反对的。因此：您必须将我上面放置的被积函数分成两部分，并使用自定义求积方法来评估每个部分（这构成了第二个细分）。
当我运行此代码 integrand = @(x,y) (phi(j,y,X) - phi(j,x,X))*(phi(i,y,X) - phi(i,x,X))/abs(y-x)^(2*s+1) % The inner integral is a function of x, and integrates over y inner = @(x) quadgk(@(y)integrand(x,y), x-lambda, x+lambda) % The inner integral is integrated over x to yield the value of the double integral dblIntegral = quadgk(inner, -(1+lambda), 1+lambda) 时，出现以下错误：（请参阅下一条评论）

【解决方案2】：

您可以进行 4 点一维高斯求积。您似乎认为这意味着二维积分。不是这样 - 这是假设一维的高阶正交。

如果您要解决一维有限元问题，那么在二维域上进行积分毫无意义。

我没有读过论文，但这是我从 FEA 中学到的。

【讨论】：