图

FEDVR 算法

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

FEDVR 基底

   现在以一维 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 的终点坐标, 把式 5 中的 $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) = \leftgroup{ &\frac{1}{\sqrt{a_i}} f_j\qtyRound{\frac{x-b_i}{a_i}} \qquad &&( 1 < j < N) \\ & \frac{1}{\sqrt{a_i+a_{i+1}}} \qtySquare{ f_N\qtyRound{\frac{x-b_i}{a_i}} + f_1\qtyRound{\frac{x-b_{i+1}}{a_{i+1}}}} \quad &&(j = N) } \end{equation}
易证正交归一关系
\begin{equation} \int_{x_1}^{x_{N_e+1}} u_{ij}(x) u_{i'j'}(x) \dd{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}$ 为式 5 中的 $w_i$, 则每个基底在每个节点处的函数值为
\begin{equation} u_{ij}(x_{i'j'}) = \delta_{ii'}\delta_{jj'}\times \leftgroup{ &\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{equation}

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

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

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

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

动能算符

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

\begin{equation}\ali{ K_{ij} &= -\frac12\mel{u_i}{\dv[2]{x}}{u_j}= -\frac12 \int u_i(x) u_j''(x) \dd{x}\\ &= -\frac12 \eval{u_i(x) u'_j(x)}_{x_1}^{x_2} + \frac12 \int u'_i(x) u_j'(x) \dd{x}\\ &= \frac12 \braketTwo{u'_i(x)}{u_j'(x)} }\end{equation}
可以看出 $\mat K$ 是实对称矩阵. 这个积分可以精确地用求和表示, 而且只有同一个 FE 里面的 $u_i(x), u_j(x)$ 才能使积分不为零(bridge function 属于两个 FE), 所以就得到了几乎是 block diagonalized 的矩阵, 只是每一个 block 左上角和右下角的一个矩阵元与其他 block 重叠.

   那如何计算基底函数的导数呢? 可以通过计算勒让德插值多项式的导数得到

\begin{equation}\ali{ &f_i'(t_j) =\frac{1}{\sqrt{w_{0i}}} \times\\ &\leftgroup{ &\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{equation}
代入得
\begin{equation} u'_{ij}(x_{ij'}) = \frac{1}{a_i^{3/2}} f'_j\qtyRound{t_{j'}} \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\qtyRound{t_{j'}} \end{equation}
\begin{equation} u'_{iN}(x_{i+1, j'}) = \frac{1}{a_{i+1}\sqrt{a_i+a_{i+1}}} f'_1\qtyRound{t_{j'}} \end{equation}
注意 $u'_{iN}$ 在 $x_{iN}$ 处的左右导数不相等.

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

\begin{equation} \ali{ K_{(im), (in)} &= \frac12 \int_{x_{i1}}^{x_{iN}} u'_{im}(x) u'_{in}(x) \dd{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{equation}
再看含 bridge function 的矩阵元
\begin{equation} \ali{ K_{(im), (iN)} &= \frac12 \int_{x_{i1}}^{x_{iN}} u'_{im}(x) u'_{iN}(x) \dd{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{equation}
\begin{equation} \ali{ K_{(iN), (iN)} &= \frac12 \int_{x_{i1}}^{x_{iN}} u'_{iN}(x)^2 \dd{x} + \frac12 \int_{x_{i+1,1}}^{x_{i+1,N}} u'_{i+1,1}(x)^2 \dd{x}\\ &= \frac{1}{2(a_i + a_{i+1})} \sum_k w_{0k} \qtySquare{\frac{1}{a_i} f'_N(t_k)^2 + \frac{1}{a_{i+1}} f'_1(t_k)^2} } \end{equation}
\begin{equation} \ali{ K_{(iN), (i+1,n)} &= \frac12 \int_{x_{i+1,1}}^{x_{i+1,N}} u'_{i,N}(x) u'_{i+1,n}(x) \dd{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{equation}
\begin{equation} \ali{ K_{(iN), (i+1,N)} &= \frac12 \int_{x_{i+1,1}}^{x_{i+1,N}} u'_{i,N}(x) u'_{i+1,N}(x) \dd{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{equation}

势能矩阵

\begin{equation} V_{\alpha\beta} = \int_{x_1}^{x_{N_e+1}} u_n^j (x) u_{n'}^{j'} (x) V(x) \dd{x} = \end{equation}


1. 未完成:这么说不严谨! 因为并不是所有基底都是光滑的, 分部积分的常数项也未必能消掉.

致读者: 小时物理百科一直以来坚持所有内容免费且不做广告,这导致我们处于日渐严重的亏损状态。长此以往很可能会最终导致我们不得不选择商业化,例如大量广告,内容付费,会员制,甚至被收购。因此,我们鼓起勇气在此请求广大读者热心捐款,使网站得以健康发展。如果看到这条信息的每位读者能慷慨捐助 10 元,我们几天内就能脱离亏损状态,并保证网站能在接下来的一整年里向所有读者继续免费提供优质内容。感谢您的支持。
—— 小时(项目创始人)

编辑词条 返回目录 返回主页 捐助项目 © 小时物理百科 保留一切权利