一维Helmholtz方程的Chebyshev - Galerkin谱方法以及Python实现

Legendre-Galerkin方法以及Chebyshev-Legendre-Galerkin方法以及Python实现

加权余量法简介

给定算子方程以及边界条件:
$$\{\begin{array}{ll}Tu-f=0& u\in V\\ Bu-g=0& u\in S\end{array}$$
其中$u$为待求函数, $T$ 为区域内部 $V$ 上算子, $B$ 为区域边界上的算子, $f,g$ 分别为定义在区域内部和区域边界的函数.

处理上述问题的加权余量方法如下:
首先,设方程组的近似解为:
$$u_n=\sum _{i=1}^n{a}_i{\phi }_i(x_1,x_2,\cdots ,x_m)$$
其中, $a_i,(i=1,\cdots ,n)$ 待定, ${\phi }_i$为某个容许空间上的基函数(具体什么空间需要具体分析.)
将方程带入原方程后一般无法近似精确求解,即在一般情况下, $T{u}_n\ne f$
因此,引入以下两个残差函数:
$$\{\begin{array}{l}{R}_v=T{u}_n-f\\ R_s=B{u}_n-g\end{array}$$
统称上述两个残差函数为$剩余函数$.
加权余量法的思想是适当选取$W_{vi},W_{si}$使得余量$R_v$以及$R_s$在某种意义上为0,一般的,可以令余量与加权内积的满足正交性条件:
$$\{\begin{array}{l}(R_v,W_{vi}{)}_v={\int }_V{R}_v{W}_{vi}dV={\int }_V\left(T{u}_n-f\right)W_{vi}dV=0\\ (R_s,W_{si}{)}_s={\int }_S{R}_s{W}_{si}dS={\int }_S\left(B{u}_n-g\right)W_{si}dS=0\end{array}$$
如果适当的选取$u_n$,如使得$u_n$能够满足边界条件, 则上面的正交性条件将退化为:
$$(R_v,W_{vi}{)}_v={\int }_V\left(T{u}_n-f\right)W_{vi}dV=0$$
将上式的加权余量法称为$内部法$;
如果$u_n$的选取满足$T{u}_n-f=0$,则上面的方程组将退化为:
$$(R_s,W_{si}{)}_s={\int }_S\left(B{u}_n-g\right)W_{si}dS=0$$
称上述的方法为加权余量的$边界法$.
同理,如果选取的$u_n$既不满足内部方程,又不满足边界方程,则必须同时用两个正交性条件消除剩余.
一般的,会选取合适的基,使得基函数满足边界条件,即考虑内部法.
对于内部法,将近似解带入内积中,选取$n$个权函数$W_{vi},\cdots ,W_{vn}$可得到方程组:
$${\int }_V\left(T{u}_n-f\right)W_{vi}dV=0,(i=1,2,\cdots ,n)$$
再从上述方程组中求解出$a_i$,最终得到逼近解.

由上,总结加权余量方法的步骤如下:

1. 选取适当的基函数.
2. 带入算子方程,求出剩余表达式.
3. 选取适当个数权函数作剩余表达式和权函数的内积,并令其正交消除剩余,从方程中求解出坐标.

注:加权余量法中权函数的选取是非常重要的, 它和所求近似解的精度有着密切的关系, 特别的如果将权函数就取为$W_{vi}={\phi }_i,(i=1,2,\cdots ,n)$, 则将得到方程组:
$${\int }_V\left(T{u}_n-f\right){\phi }_{i}dV=0,(i=1,2,\cdots ,n)$$
称上述方程组为$Galerkin方程组$, 其解称为$Galerkin系数.$

一维$Helmholtz$方程的Galerkin 方法

下面进入主题,考虑如下一维齐次$Dirichlet$型边界条件的$Helmholtz$方程:
$$\{\begin{array}{l}-u^{{}^{\prime \prime }}+\alpha u=finI=(-1,1)\\ u(\pm 1)=0\end{array}$$
定义空间
P N 0 = { ϕ ∈ P N ∣ ϕ ( ± 1 ) = 0 } P^{0}_{N} = \left\lbrace \phi \in P_{N} | \phi(\pm1) = 0\right\rbrace PN0​={ϕ∈PN​∣ϕ(±1)=0}
其中$P_N$ 表示次数不超过N的全体多项式组成的空间, 则显然可知$P_N^0$为$N-2$维线性空间.利用加权余量法的步骤推导上述方程的$Galerkin$方法:

1. 取试函数: $u_N=\sum _{i=0}^{N-2}{a}_i{\varphi }_i(x)$, 其中${\varphi }_i(x)$为空间$P_N^0$的基函数.
2. 将$u_N$代入原方程得到剩余表达式:
   $$R_v=-u_N^{{}^{\prime \prime }}+\alpha u_N-f$$
3. 在权函数$\omega (x)$的意义下得到$Galerkin方程组:$
   $${\int }_{-1}^1{R}_v(x){\varphi }_j(x)\omega (x)dx=0j=0,1,\cdots ,N-2$$
   实际上,利用线性空间的定义, 上述过程一般还可以描述为如下:
   求$u_N\in P_N^0$, s.t.
   $${\int }_{-1}^1{R}_v{\varphi }_j\omega (x)dx=0\forall v_N\in P_N^0$$
   将残量带入,即:
   求$u_N\in P_N^0$, s.t.
   $$-(u_N^{{}^{\prime \prime }},v_N{)}_{\omega }+\alpha (u_N,v_N{)}_{\omega }=(f,v_N{)}_{\omega }\forall v_N\in P_N^0$$

其中$(u,v{)}_{\omega }:={\int }_{-1}^{1}u(x)v(x)\omega (x)dx$为在权重$\omega$意义下定义的内积.

注意,实际上右端项需要计算积分,但是当$f$ 为一般函数时,右端一般无法精确求积,因此在实际过程中需要采用数值积分.在$Galerkin$谱方法中,通常会将$f$用某种意义上的插值多项式$I_{N}f$代替,插值节点一般选为如$Gauss$型求积公式的节点, 由此最终上式子写为:
求$u_N\in P_N^0$, s.t.
$$-(u_N^{{}^{\prime \prime }},v_N{)}_{\omega }+\alpha (u_N,v_N{)}_{\omega }=(I_{N}f,v_N{)}_{\omega }\forall v_N\in P_N^0$$
现在将试函数在基函数下的展开$u_N=\sum _{i=0}^{N-2}{\hat{u}}_i{\varphi }_i(x)dx$带入上面式子, 并且将$v_N$取便所有的基函数即可得到下面方程组:
$$-\sum _{i=0}^{N-2}({\varphi }_i^{{}^{\prime \prime }},{\varphi }_j{)}_{\omega }{\hat{u}}_i+\alpha \sum _{i=0}^{N-2}{\left({\varphi }_i,{\varphi }_j\right)}_{\omega }{\hat{u}}_i=(f_N,{\varphi }_j{)}_{\omega },j=0,\cdots ,N-2$$
将上述方程组改写为矩阵形式得:
$$(S+\alpha M)\hat{U}=f$$
其中:

$(S{)}_{ij}=-({\varphi }_j^{{}^{\prime \prime }}(x),{\varphi }_i(x){)}_{\omega }$ 称为$刚度矩阵$.

$(M{)}_{ij}=({\varphi }_j(x),{\varphi }_i(x){)}_{\omega }$ 称为$质量矩阵$.
$\hat{U}=({\hat{u}}_0\cdots {\hat{u}}_{N-2}{)}^T$
$\mathbf{f}={\left((I_{N}f,{\varphi }_0(x){)}_{\omega },\cdots ,(I_{N}f,{\varphi }_{N-2}(x){)}_{\omega }\right)}^T$
当选取不同的基函数带入上述方程组之后,就会得到不同形式的矩阵,一般希望得到的矩阵需要有良好的性质,(如:希望得到稀疏的矩阵)使得能够快速求解上述方程组.

Chebyshev-Galerkin方法

继续考虑如下一维齐次$Dirichlet$型边界条件的$Helmholtz$方程:
$$\{\begin{array}{l}-u^{{}^{\prime \prime }}+\alpha u=finI=(-1,1)\\ u(\pm 1)=0\end{array}$$
在前面得到上述方程的变分形式:
求$u_N\in P_N^0$, s.t.
$$-(u_N^{{}^{\prime \prime }},v_N{)}_{\omega }+\alpha (u_N,v_N{)}_{\omega }=(I_{N}f,v_N{)}_{\omega }\forall v_N\in P_N^0$$
并选取基函数带入方程,得到以下矩阵方程:
$$(S+\alpha M)\hat{U}=f$$
在本节将会推导将基函数选取为$Chebyshev$基函数, 对应的权函数将会选取为:$\omega (x)=(1-x^2{)}^{-1/2}$
($此时,上面的Galerkin方法称为Chebyshev-Galerkin方法$)时所得到的的刚度矩阵以及质量矩阵,并给出利用$Chebyshev-Galerkin$方法求解上述方程的具体流程.

基函数的选取

由上面可知,要写出具体的方程,首先要选取合适的基函数,直接将基函数选取为$T_n(x)=cos\left(nco{s}^{-1}(x)\right)$肯定是不行的,因为这样的基函数时不在$P_N^0$中的,一般的基函数将会用如下形式选取:
$${\varphi }_k(x)=T_k(x)+a_k{T}_{k+1}(x)+b_k{T}_{k+2}(x)$$
通过将基函数带入边界条件, 求解出 $a_k,b_k$,注意这样的方法对一般的边界条件也适用,可以证明,当边界满足适定性条件的时候(关于适定性条件的定义,可差查阅后面的参考资料1), 对于任意的$k$总能将$a_k,b_k$求解出来.在这里首先考虑齐次$Dirichelet$型边界条件.
现在,为了便于后面的推导,首先给出$Chebyshev$多项式的几个性质:

• 正交性: ${\int }_{-1}^1{T}_k(x)T_j(x)(1-x^2{)}^{-1/2}dx=\frac{c_k\pi }{2}{\delta }_{kj}$
  其中, $c_0=2,c_k=1\forall k>1$
• $T_n(\pm 1)=(\pm 1{)}^n$
• $$\begin{gathered} T_n^{{}^{\prime \prime }}(x)=\sum _{k=0 \\ k+neven}^{n-2}\frac{1}{c_k}n(n^2-k^2)T_k(x) \end{gathered}$$
• Chebyshev-Gauss-Lobatto 型求积公式求积节点以及对应权重:
  $$x_0=1,x_N=-1,{\omega }_0={\omega }_N=\frac{\pi }{2N},x_j=cos\frac{\pi j}{N},{\omega }_j=\frac{\pi }{N},1\le j\le N-1$$
  将基函数${\varphi }_k(x)=T_k(x)+a_k{T}_{k+1}(x)+b_k{T}_{k+2}(x)$带入边界条件, 利用性质$T_n(\pm 1)=(\pm 1{)}^n$得以下方程组:
  $$\{\begin{array}{ll}{a}_k+b_k+1& =0\\ -a_k+b_k+1& =0\end{array}$$
  从上述方程中立刻解得 $a_k=0,b_k=-1$:
  由上可知,在这里选取基函数为:
  ${\varphi }_k(x)=T_k(x)-T_{k+2}(x),0\le k\le N-2$
  基函数取定之后,便可以推导具体的刚度矩阵以及质量矩阵.

刚度矩阵推导

$$\begin{array}{rl}{s}_{ij}& =-{\left({\varphi }_j^{{}^{\prime \prime }},{\varphi }_i\right)}_{\omega }=-{\left(T_j^{{}^{\prime \prime }}-T_{j+2}^{{}^{\prime \prime }},T_i-T_{i+2}\right)}_{\omega }\\ & =-{\left(T_j^{{}^{\prime \prime }},T_i\right)}_{\omega }+{\left(T_j^{{}^{\prime \prime }},T_{i+2}\right)}_{\omega }+{\left(T_{j+2}^{{}^{\prime \prime }},T_i\right)}_{\omega }-{\left(T_{j+2}^{{}^{\prime \prime }},T_{i+2}\right)}_{\omega }\end{array}$$
套用上面的二阶导数与一阶导数的关系,以及正交性条件,容易得到:
$$s_{ij}=\{\begin{array}{ll}2\pi (i+1)(i+2)& j=i\\ 4\pi (i+1)& j=i+2,i+4,i+6,\cdots \\ 0& j<iorj+iodd\end{array}$$

质量矩阵推导

容易知道 ,$m_{ij}=m_{ji}$, 因此,不妨设 $j\ge i$有:
$$\begin{array}{rl}({\varphi }_j,{\varphi }_i{)}_{\omega }& =(T_j-T_{j+2},T_i-T_{i+2}{)}_{\omega }\\ & =(T_j,T_i{)}_{\omega }-(T_j,T_{i+2}{)}_{\omega }-(T_{j+2},T_i{)}_{\omega }+(T_{j+2},T_{i+2}{)}_{\omega }\\ & =\frac{c_i\pi }{2}{\delta }_{ji}-\frac{c_{i+2}\pi }{2}{\delta }_{j,i+2}+\frac{c_{i+2}\pi }{2}{\delta }_{ji}\\ & =\{\begin{array}{ll}\frac{\pi }{2}(c_i+1),& j=i\\ 0,& j=i+1\\ -\frac{\pi }{2},& j=i+2\end{array}\end{array}$$

右端项推导

推导右端项,首先要计算插值多项式 $I_{N}f$, 实际上即 $Chebyshev$变换.
假设插值节点取为$Gauss-Chebyshev-Lobatto$型求积公式的求积节点,求积节点以及权重在上面给出.
则插值多项式设为:
$$I_{N}f=\sum _{n=0}^N{\tilde{f}}_n{T}_n(x)$$
则谱系数${\tilde{f}}_n$ 由$Chebyshev$ 正变换给出如下:
$$\begin{array}{r}{\tilde{f}}_n=\frac{2}{{\tilde{c}}_{n}N}\sum _{k=0}^N\frac{1}{{\tilde{c}}_k}f(x_k)cos(nk\pi /N)\end{array}$$
其中:
$${\tilde{c}}_n=\{\begin{array}{ll}2& n=0,N\\ 1& n=1,2,\cdots N-1\end{array}$$
右端项
$$\begin{array}{rl}(\mathbf{f}{)}_i& =(I_{N}f,{\varphi }_i(x){)}_{\omega }=\sum _{n=0}^N{\tilde{f}}_n\left(T_n(x),T_i(x)-T_{i+2}(x)\right)\\ & =\frac{c_i\pi }{2}{\tilde{f}}_i-\frac{c_{i+2}\pi }{2}{\tilde{f}}_{i+2}(0\le i\le N-2)\end{array}$$
将上面的式子写为向量形式,得:
$$\mathbf{f}=\left(\begin{array}{cccccc}{c}_0& 0& -c_2& & & \\ & c_1& 0& -c_3& & \\ & & \ddots & \ddots & \ddots & \\ & & & c_{N-2}& 0& -c_N\end{array}\right)\left(\begin{array}{c}{\tilde{f}}_0\\ {\tilde{f}}_1\\ ⋮\\ {\tilde{f}}_N\end{array}\right)$$
有了左端项矩阵以及右端项,就可以求解方程:
$$(\mathbf{S}+\alpha \mathbf{M})\hat{U}=\mathbf{f}$$
得到 $\hat{U}=({\hat{u}}_0,{\hat{u}}_1,\cdots ,{\hat{u}}_{N-2})$

Chebyshev逆变换得到解在求积节点处的值.

经过上面的步骤已经可以求得展开系数,在实际过程中,希望最终得到解$u_N=\sum _{i=0}^{N-2}{\hat{u}}_i{\varphi }_i(x)$在各个节点处的的值,这里的节点特指$Chebyshev$求积公式的求积节点,因为选取这些点之后,可以利用$Chebyshev$ 逆变换来计算解在求积节点处的值 $u_N(x_j)$, 并且$Chebyshev$逆变换可以利用$O(NlogN)$的快速算法完成.下面对一般形式的基函数${\varphi }_k(x)=T_k(x)+a_k{T}_{k+1}(x)+b_k{T}_{k+2}(x)$来推导逆变换的过程.
由
$$\begin{gathered} u_N(x_j)=\sum _{i=0}^{N-2}{\hat{u}}_i{\varphi }_i(x_j),j=0,\cdots ,N \\ =\sum _{i=0}^{N-2}{\hat{u}}_i{T}_i(x_j)+\sum _{i=0}^{N-2}{a}_i{\hat{u}}_i{T}_{i+1}(x_j)+\sum _{i=0}^{N-2}{b}_i{\hat{u}}_i{T}_{i+2}(x_j) \end{gathered}$$
注意使用 $T_i(x_j)=cos(ij\pi /N)$则
$$\begin{gathered} u_N(x_j)=\sum _{i=0}^{N-2}{\hat{u}}_{i}cos(ij\pi /N)+\sum _{i=0}^{N-2}{a}_i{\hat{u}}_{i}cos((i+1)j\pi /N)+\sum _{i=0}^{N-2}{b}_i{\hat{u}}_{i}cos((i+2)j\pi /N) \\ =\sum _{i=0}^N({\hat{u}}_i+a_{i-1}{\hat{u}}_{i-1}+b_{i-2}{\hat{u}}_{i-2})cos(ij\pi /N) \end{gathered}$$
这里规定
$$\begin{gathered} {\hat{u}}_{-2}={\hat{u}}_{-1}={\hat{u}}_{N-1}={\hat{u}}_N=0 \\ {a}_{-1}=a_{N-1}=b_{-1}=b_{-2}=0 \end{gathered}$$

Chebyshev-Galerkin方法的计算流程:

最终,给出$Chebyshev-Galerkin$方法的计算流程如下:

• 获取 Chebyshev-Gauss-Lobatto 型求积公式的求积节点
  $$x_0=1,x_N=-1,{\omega }_0={\omega }_N=\frac{\pi }{2N},x_j=cos\frac{\pi j}{N},{\omega }_j=\frac{\pi }{N},1\le j\le N-1$$
  预存储基函数系数: { a k , b k } \left\lbrace a_{k}, b_{k} \right\rbrace {ak​,bk​}(二维数组)
  对于齐次$Dirichlet$型边界条件, 有:$a_k=0,b_k=-1$, 对于$Neumman$型边界条件,有
  稀疏存储刚度矩阵$S$以及质量矩阵$M$
• 计算向量 $\left(\begin{array}{c}f(x_0)\\ f(x_1)\\ ⋮\\ f(x_N)\end{array}\right)$, 对该向量进行$Chebysheb$正变换得到: $\left(\begin{array}{c}{\tilde{f}}_0\\ {\tilde{f}}_1\\ ⋮\\ {\tilde{f}}_N\end{array}\right)$
• 计算右端项:
  $$\mathbf{f}=\left(\begin{array}{cccccc}{c}_0& 0& -c_2& & & \\ & c_1& 0& -c_3& & \\ & & \ddots & \ddots & \ddots & \\ & & & c_{N-2}& 0& -c_N\end{array}\right)\left(\begin{array}{c}{\tilde{f}}_0\\ {\tilde{f}}_1\\ ⋮\\ {\tilde{f}}_N\end{array}\right)$$
  求解代数方程组:
  $$(\mathbf{S}+\alpha \mathbf{M})\hat{U}=\mathbf{f}$$
  得到 $\hat{U}=({\hat{u}}_0,{\hat{u}}_1,\cdots ,{\hat{u}}_{N-2})$
• 作$Chebyshev$逆变换得:
  $$u_N(x_j)=\sum _{i=0}^{N-2}{\hat{u}}_i{\varphi }_i(x_j)$$
  这里利用
  $$\begin{gathered} u_N(x_j)=\sum _{i=0}^N({\hat{u}}_i+a_{i-1}{\hat{u}}_{i-1}+b_{i-2}{\hat{u}}_{i-2})cos(ij\pi /N) \\ {\hat{u}}_{-2}={\hat{u}}_{-1}={\hat{u}}_{N-1}={\hat{u}}_N=0 \\ {a}_{-1}=a_{N-1}=b_{-1}=b_{-2}=0 \end{gathered}$$

具体实现时, 先将输入进行改写,输入有:

1. 第一步中存储 $a_k,b_k$ 的基函数系数矩阵 base_coef: 大小为(N-1, 2)的 二维数组.
2. 解向量 $solution=({\hat{u}}_0,{\hat{u}}_1,\cdots ,{\hat{u}}_{N-2})$,大小为(N-1,)的一维数组.
   首先要改写 解向量数组,让它变成快速切比雪夫逆变换(这个函数预先定义)的输入,用如下公式来改写:
   $$\left[\begin{array}{c}{\tilde{u}}_0\\ {\tilde{u}}_1\\ ⋮\\ {\tilde{u}}_N\end{array}\right]=\left[\begin{array}{cccccccc}0& 0& 1& & & & & \\ & 0& a_0& 1& & & & \\ & & b_0& a_1& 1& & & \\ & & & \ddots & \ddots & \ddots & & \\ & & & & b_{N-3}& a_{N-2}& 1& \\ & & & & & b_{N-2}& 0& 1\end{array}\right]\left[\begin{array}{c}0\\ 0\\ {\hat{u}}_0\\ ⋮\\ 0\\ 0\end{array}\right]$$
   将改写过后的$\left[\begin{array}{c}{\tilde{u}}_0\\ {\tilde{u}}_1\\ ⋮\\ {\tilde{u}}_N\end{array}\right]$送入自定义函数 backward_chebyshev()完成快速切比雪夫逆变换得到最终结果.

注: 由于最近比较忙,代码会在以后补出