FEDVR 网格

                     

贡献者: addis

预备知识 高斯积分(Gauss-Lobatto)

  

未完成:简单介绍什么是 FEDVR,有什么用

1. FEDVR 基底

  1现在以一维 FEDVR 为例,把整个区间划分成 $N_e$ 个有限元,第 $i$ 个有限元的区间为 $[x_i,x_{i+1}]$。每个有限元内进一步加入格点,令 $x_{ij}\ (j = 1\dots N)$ 为第 $i$ 有限元的 $N$ 阶 Gauss-Lobatto 数值积分 的 $N$ 个采样点。这样,$x_{j1}=x_j$, $x_{jN}=x_{j+1}$。

   现在来定义 FEDVR 基底 $u_{ij}(x)$。令 $a_i$ 是第 $i$ 个 FE 长度的一半,$b_i$ 是第 $i$ 个 FE 的中点坐标,把式 12 中的 $f_n(x)$($x$ 超出 $[-1,1]$ 时令 $f_n(x) = 0$)记为 $f_n(t)$,并把 $t \in [-1,1]$ 线性依次映射到每个有限元的区间 $x \in [x_i, x_{i+1}]$ 上,有(下文中 $t$ 和 $x$ 都分别使用该定义)

\begin{equation} x_{ij} = a_i t_j + b_i~, \end{equation}
\begin{equation} x = a_i t + b_i \qquad (x_{i1} \leqslant x \leqslant x_{iN})~. \end{equation}
则基底的定义为
\begin{equation} u_{ij}(x) = \left\{\begin{aligned} &\frac{1}{\sqrt{a_i}} f_j \left(\frac{x-b_i}{a_i} \right) &&( 1 < j < N) \\ & \frac{1}{\sqrt{a_i+a_{i+1}}} \left[f_N \left(\frac{x-b_i}{a_i} \right) + f_1 \left(\frac{x-b_{i+1}}{a_{i+1}} \right) \right] &&(j = N)~. \end{aligned}\right. \end{equation}
每个基底都是由一个 $N-1$ 阶多项式(或两个拼接而成)。若使用求和代替积分(存在误差),基底之间满足正交归一关系
\begin{equation} \int_{x_1}^{x_{N_e+1}} u_{ij}(x) u_{i'j'}(x) \,\mathrm{d}{x} \approx \delta_{ii'} \delta_{jj'}~. \end{equation}
在总区间的两端我们一般采用函数值为零的边界条件,所以不定义 $u_{11}(x)$ 和 $u_{N_e, N}(x)$。这样,我们最终共有 $N_e(N-1)-1$ 个节点 $x_{ij}$ 和对应的基底 $u_{ij}(x)$。令 $w_{0i}$ 为式 12 中的 $w_i$,则每个基底在每个节点处的函数值为
\begin{equation} u_{ij}(x_{i'j'}) = \delta_{ii'}\delta_{jj'}\times \left\{\begin{aligned} &\frac{1}{\sqrt{a_i w_{0j}}} &&( 1 < j < N) \\ & \frac{1}{\sqrt{(a_i+a_{i+1}) w_{0N}}} &&(j = N)~. \end{aligned}\right. \end{equation}

   易证,若要对整个区间积分,所需权重为

\begin{equation} w_{ij} = \frac{1}{u_{ij}^2(x_{ij})} = \begin{cases} a_i w_{0j} &(1 < j < N) \\ (a_i + a_{i+1}) w_{0N} &(j = N)~. \end{cases} \end{equation}
但如果只需要在某个 FE 区间内做积分,只需用 $w_{ij} = a_i w_{0j}$ 即可。

   一个函数 $f(x)$ 用基底展开为

\begin{equation} \left\langle u_{ij}(x) \middle| f(x) \right\rangle = \int_{x_1}^{x_{N_e+1}} u_{ij}(x) f(x) \,\mathrm{d}{x} \approx \sqrt{w_{ij}} f(x_{ij})~. \end{equation}
把所有基底按对应的 $x_{ij}$ 从小到大排序。$f(x)$ 的展开系数就可以记为一个列矢量。

2. 二阶导数矩阵(动能矩阵)

   注意由于同一个 FE 内的基底并不怎么正交($N=6$ 时内积的相对误差甚至可达 10%!更高的 $N=18$ 时也有 3%),主要是求出矩阵后,用数值验证的方法验证是否精确即可(对低阶多项式,精确到机器精度,即使相邻三个 FE 长度都不等)。

   虽然不是所有的 $u_i$ 基底都处处二阶可导,但我们可以假设它们所表示的波函数 $\psi(x) = \sum_k c_k u_k(x)$ 处处二阶可导。求导后,系数为(由于 $u_i, u_j$ 不正交,所以我们在第一个等号中特意避开积分,但由于被积多项式阶数足够低,可以用积分精确代替求和,进而便于使用分部积分。这样求出来的结果和直接计算第一个等号后的求和是完全一样的)

\begin{equation} c_i^{(2)} = \frac{1}{u_i(x_i)} \left. \left( \frac{\mathrm{d}^{2}}{\mathrm{d}{x}^{2}} \sum_k c_k u_k \,\mathrm{d}{x} \right) \right\rvert _{x=x_i} = \int_{-\infty}^{+\infty} {u_i} \frac{\mathrm{d}^{2}}{\mathrm{d}{x}^{2}} \sum_k c_k u_k \,\mathrm{d}{x} ~. \end{equation}
所以为了避免处理有限元的边界,我们可以在积分中直接忽略所有的边界点(忽略有限个点不会影响积分值)。
\begin{equation} c_i^{(2)} = \sum_k c_k \int u_i u''_k \,\mathrm{d}{x} ~. \end{equation}
此时求和中的每个积分范围可以缩小到一个或两个有限元内,此时不可以使用 bridge 的全局权重。
\begin{equation} \begin{aligned} D^{(2)}_{ij} &= \int u_i u''_k \,\mathrm{d}{x} \\ &= \left. u_i(x) u'_j(x) \right\rvert _{x_1}^{x_2} - \int u'_i u'_j \,\mathrm{d}{x} \qquad (u_i \text{ non-bridge})~. \end{aligned} \end{equation}
由直觉可得,第一项可以抵消。可以看出 $ \boldsymbol{\mathbf{D}} ^{(2)}$ 是实对称矩阵。由于被积多项式的阶数只有 $2(N-2)$(小于 $2N-3$),这个积分可以精确地用求和表示,而且只有同一个 FE 里面的 $u_i(x), u_j(x)$ 才能使积分不为零(bridge function 属于两个 FE),所以就得到了几乎是 block diagonalized 的矩阵,只是每一个 block 左上角和右下角的一个矩阵元与其他 block 重叠。

   那如何计算基底函数的导数呢?可以通过计算拉格朗日插值多项式的导数得到

\begin{equation} \begin{aligned}&f_i'(t_j) =\frac{1}{\sqrt{w_{0i}}} \times\\ & \left\{\begin{aligned} &\frac{t_j-t_1}{t_i-t_1} \frac{t_j-t_2}{t_i-t_2} \dots \frac{1}{t_i-t_j} \dots \frac{t_j-t_{i-1}}{t_i-t_{i-1}}\frac{t_j-t_{i+1}}{t_i-t_{i+1}} \dots \frac{t_j-t_N}{t_i-t_N} &&(i \ne j) \\ & \frac{1}{t_i-t_1} + \dots + \frac{1}{t_i-t_{i-1}} + \frac{1}{t_i-t_{i+1}} + \dots \frac{1}{t_i-t_N} &&(i = j) \end{aligned}\right. \\ \\ & (i,j=1,\dots,N)~. \end{aligned} \end{equation}
代入得式 3 求导得
\begin{equation} u'_{ij}(x_{ij'}) = \frac{1}{a_i^{3/2}} f'_j \left(t_{j'} \right) \qquad ( 1 < j < N)~, \end{equation}
\begin{equation} u'_{iN}(x_{i j'}) = \frac{1}{a_i\sqrt{a_i+a_{i+1}}} f'_N \left(t_{j'} \right) ~, \end{equation}
\begin{equation} u'_{iN}(x_{i+1, j'}) = \frac{1}{a_{i+1}\sqrt{a_i+a_{i+1}}} f'_1 \left(t_{j'} \right) ~. \end{equation}
注意 $u'_{iN}$ 在 $x_{iN}$ 处不连续无定义,左右导数分别为 $u'_{iN}(x_{iN})$ 和 $u'_{iN}(x_{i+1, 1})$。

   现在就可以求二阶导数的矩阵元了

\begin{equation} D^{(2)}_{(im), (in)} = \left\{\begin{aligned} &-\int_{x_{i1}}^{x_{i+1,N}} u'_{im}(x) u'_{in}(x) \,\mathrm{d}{x} && (m=n=N)\\ &-\int_{x_{i1}}^{x_{iN}} u'_{im}(x) u'_{in}(x) \,\mathrm{d}{x} && (\text{otherwise})~. \end{aligned}\right. \end{equation}
但注意对 bridge function 的积分不宜直接用全局 weight,还是要在每个 FE 内部分别求和。

   先看不含 bridge function 的矩阵元:

\begin{equation} \begin{aligned} D^{(2)}_{(im), (in)} &= -\frac{1}{a_i^2} \sum_k w_{0k} f'_m(t_k) f'_n(t_k) \qquad (1 < m \leqslant n < N)~. \end{aligned} \end{equation}
再看含 bridge function 的矩阵元:
\begin{equation} \begin{aligned} D^{(2)}_{(im), (iN)} &= \frac{-1}{a_i^{3/2} \sqrt{a_i + a_{i+1}}} \sum_k w_{0k} f'_m(t_k) f'_N(t_k) \qquad (1< m < N)~, \end{aligned} \end{equation}
\begin{equation} \begin{aligned} D^{(2)}_{(iN), (iN)} &= \frac{-1}{(a_i + a_{i+1})} \sum_k w_{0k} \left[\frac{1}{a_i} f'_N(t_k)^2 + \frac{1}{a_{i+1}} f'_1(t_k)^2 \right] ~, \end{aligned} \end{equation}
\begin{equation} \begin{aligned} D^{(2)}_{(iN), (i+1,n)} &= \frac{-1}{a_{i+1}^{3/2} \sqrt{a_i + a_{i+1}}} \sum_k w_{0k} f'_1(t_k) f'_n(t_k) \qquad (1< n < N)~, \end{aligned} \end{equation}
\begin{equation} \begin{aligned} D^{(2)}_{(iN), (i+1,N)} &= \frac{-1}{a_{i+1}\sqrt{(a_i + a_{i+1})(a_{i+1}+a_{i+2})}} \sum_k w_{0k} f'_1(t_k) f'_N(t_k)~. \end{aligned} \end{equation}

3. 一阶导数矩阵

   经数值验证,有

\begin{equation} f'_i(x_i) = 0~. \end{equation}
所以下面可以发现 $ \boldsymbol{\mathbf{D}} $ 的对角元皆为零。对低阶多项式,$ \boldsymbol{\mathbf{D}} $ 对多项式求导结果精确到机器精度(即使相邻两个 FE 长度不等)。

   和二阶导数同理,一阶导数的矩阵元为(积分范围挖去 $x_{iN}$)

\begin{equation} D^{(1)}_{(im), (in)} = \frac{u'_{in}(x_{im})}{u_{im}(x_{im})} = \left\{\begin{aligned} &\int_{x_{i1}}^{x_{i+1,N}} u_{im}(x) u'_{in}(x) \,\mathrm{d}{x} && (m=n=N)\\ &\int_{x_{i1}}^{x_{iN}} u_{im}(x) u'_{in}(x) \,\mathrm{d}{x} && (\text{otherwise})~. \end{aligned}\right. \end{equation}
先看不含 bridge function 的矩阵元:
\begin{equation} \begin{aligned} D_{(im), (in)} &= \frac{\sqrt{w_{0m}}}{a_i} f'_n(t_m) \qquad (1 < m \leqslant n < N)~. \end{aligned} \end{equation}
再看含 bridge function 的矩阵元:
\begin{equation} \begin{aligned} D_{(im), (iN)} &= \sqrt\frac{{w_{0m}}}{{a_i}(a_i + a_{i+1})} f'_N(t_m) \qquad (1< m < N)~, \end{aligned} \end{equation}
\begin{equation} \begin{aligned} D_{(iN), (iN)} &= 0~, \end{aligned} \end{equation}
\begin{equation} \begin{aligned} D_{(iN), (i+1,n)} &= \sqrt\frac{{w_{01}}}{(a_i + a_{i+1}){a_{i+1}}} f'_n(t_1) \qquad (1< n < N)~, \end{aligned} \end{equation}
\begin{equation} \begin{aligned} D_{(iN), (i+1,N)} &= \sqrt{\frac{{w_{01}}}{{(a_i + a_{i+1})(a_{i+1}+a_{i+2})}}} f'_N(t_1)~. \end{aligned} \end{equation}


1. ^ 参考 Aihua Liu 的博士论文 (KSU, 2015)。


致读者: 小时百科一直以来坚持所有内容免费,这导致我们处于严重的亏损状态。 长此以往很可能会最终导致我们不得不选择大量广告以及内容付费等。 因此,我们请求广大读者热心打赏 ,使网站得以健康发展。 如果看到这条信息的每位读者能慷慨打赏 10 元,我们一个星期内就能脱离亏损, 并保证在接下来的一整年里向所有读者继续免费提供优质内容。 但遗憾的是只有不到 1% 的读者愿意捐款, 他们的付出帮助了 99% 的读者免费获取知识, 我们在此表示感谢。

                     

友情链接: 超理论坛 | ©小时科技 保留一切权利