FEDVR 网格

             

预备知识 高斯积分(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 的中点坐标,把式 6 中的 $f_n(x)$($x$ 超出 $[-1,1]$ 时令 $f_n(x) = 0$)记为 $f_n(t)$,并把 $t \in [-1,1]$ 线性依次映射到每个有限元的区间 $x \in [x_i, x_{i+1}]$ 上,有

\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) \qquad &&( 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] \quad &&(j = N) \end{aligned}\right. \end{equation}
易证正交归一关系
\begin{equation} \int_{x_1}^{x_{N_e+1}} u_{ij}(x) u_{i'j'}(x) \,\mathrm{d}{x} = \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}$ 为式 6 中的 $w_i$,则每个基底在每个节点处的函数值为
\begin{equation} u_{ij}(x_{i'j'}) = \delta_{ii'}\delta_{jj'}\times \left\{\begin{aligned} &\frac{1}{\sqrt{a_i w_{0j}}} \qquad &&( 1 < j < N) \\ & \frac{1}{\sqrt{(a_i+a_{i+1}) w_{0N}}} \qquad &&(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. 动能算符

   使用分部积分可以把动能矩阵元表示为两个基底的导函数的内积2

\begin{equation} \begin{aligned} K_{ij} &= -\frac12 \left\langle u_i \middle| \frac{\mathrm{d}^{2}}{\mathrm{d}{x}^{2}} \middle| u_j \right\rangle = -\frac12 \int u_i(x) u_j''(x) \,\mathrm{d}{x} \\ &= -\frac12 \left. u_i(x) u'_j(x) \right\rvert _{x_1}^{x_2} + \frac12 \int u'_i(x) u_j'(x) \,\mathrm{d}{x} \\ &= \frac12 \left\langle u'_i(x) \middle| u_j'(x) \right\rangle \end{aligned} \end{equation}
可以看出 $ \boldsymbol{\mathbf{K}} $ 是实对称矩阵.这个积分可以精确地用求和表示,而且只有同一个 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} \qquad &&(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} \quad &&(i = j) \end{aligned}\right. \end{aligned} \end{equation}
代入得
\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}$ 处的左右导数不相等.

   现在就可以求动能矩阵的矩阵元了(由于是对称矩阵,我们只列出半边),先看不含 bridge function 的矩阵元

\begin{equation} \begin{aligned} K_{(im), (in)} &= \frac12 \int_{x_{i1}}^{x_{iN}} u'_{im}(x) u'_{in}(x) \,\mathrm{d}{x} \\ &= \frac12 \sum_k a_i w_{0k} u'_{im}(x_{ik}) u'_{in}(x_{ik}) \\ &= \frac{1}{2a_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} K_{(im), (iN)} &= \frac12 \int_{x_{i1}}^{x_{iN}} u'_{im}(x) u'_{iN}(x) \,\mathrm{d}{x} \\ &= \frac{1}{2a_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} K_{(iN), (iN)} &= \frac12 \int_{x_{i1}}^{x_{iN}} u'_{iN}(x)^2 \,\mathrm{d}{x} + \frac12 \int_{x_{i+1,1}}^{x_{i+1,N}} u'_{i+1,1}(x)^2 \,\mathrm{d}{x} \\ &= \frac{1}{2(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} K_{(iN), (i+1,n)} &= \frac12 \int_{x_{i+1,1}}^{x_{i+1,N}} u'_{i,N}(x) u'_{i+1,n}(x) \,\mathrm{d}{x} \\ &= \frac{1}{2 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} K_{(iN), (i+1,N)} &= \frac12 \int_{x_{i+1,1}}^{x_{i+1,N}} u'_{i,N}(x) u'_{i+1,N}(x) \,\mathrm{d}{x} \\ &= \frac{1}{2a_{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} V_{\alpha\beta} = \int_{x_1}^{x_{N_e+1}} u_n^j (x) u_{n'}^{j'} (x) V(x) \,\mathrm{d}{x} = \end{equation}
未完成:未完成


1. ^ 参考 Aihua Liu 的博士论文 (KSU, 2015).
2. ^ 未完成:这么说不严谨! 因为并不是所有基底都是光滑的,分部积分的常数项也未必能消掉.

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

         

© 小时科技 保留一切权利