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

  • 时间:
  • 来源:互联网
  • 文章标签:

一维Helmholtz方程的Chebyshev - Galerkin谱方法
在上面这篇博文中,讲述了加权余量法的基本原理,以及系统地推导,阐述了 C h e b y s h e v − G a l e r k i n Chebyshev-Galerkin ChebyshevGalerkin方法的计算过程,本文在上面的基础上,继续探讨 L e n g e n d r e − G a r l e r k i n Lengendre-Garlerkin LengendreGarlerkin方法以及 L e g e n d r e − C h e b y s h e v − G a l e r k i n Legendre-Chebyshev-Galerkin LegendreChebyshevGalerkin方法,可以看到, C h e b y s h e v − G a l e r k i n Chebyshev-Galerkin ChebyshevGalerkin在作正变换以及逆变换的过程中,能够作快速余弦变换,但该方法所推导出的刚度矩阵稀疏程度不够; 最关键的是,刚度矩阵不是对称的!.在本文介绍的 L e n g e n d r e − G a r l e r k i n Lengendre-Garlerkin LengendreGarlerkin方法中,可以看到刚度矩阵是对角矩阵,但 L e g e n d r e Legendre Legendre正变换以及逆变换的过程没有快速计算方法,实际计算复杂度是 O ( N 2 ) O(N^{2}) O(N2);最后给出的 C h e b y s h e v − L e g e n d r e − G a l e r k i n Chebyshev-Legendre-Galerkin ChebyshevLegendreGalerkin方法同时具备了 C h e b y s h e v − G a l e r k i n Chebyshev - Galerkin ChebyshevGalerkin方法以及 L e g e n d r e − G a l e r k i n Legendre-Galerkin LegendreGalerkin方法的优势.

回忆

首先,回忆上篇博文中的重要内容
考虑如下一维齐次 D i r i c h l e t Dirichlet Dirichlet型边界条件的 H e l m h o l t z Helmholtz Helmholtz方程( α > 0 \alpha > 0 α>0):
{ − u ′ ′ + α u = f    i n   I = ( − 1 , 1 ) u ( ± 1 ) = 0 \begin{cases} -u ^{''} + \alpha u = f \ \ in \ I = (-1, 1) \\ u(\pm1) = 0 \end{cases} {u+αu=f  in I=(1,1)u(±1)=0
定义空间
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}
取试函数: u N = ∑ i = 0 N − 2 a i ϕ i ( x ) u_{N} = \sum\limits_{i = 0}^{N - 2} a_{i}\phi_{i}(x) uN=i=0N2aiϕi(x), 其中 ϕ i ( x ) \phi_{i}(x) ϕi(x)为空间 P N 0 P^{0}_{N} PN0的基函数.
对上述方程进行 R i t s z − G a l e r k i n Ritsz-Galerkin RitszGalerkin离散得到如下变分形式
u N ∈ P N 0 u_{N} \in P^{0}_{N} uNPN0, s.t.
− ( u N ′ ′ , v N ) ω + α ( u N , v N ) ω = ( I N f , v N ) ω ∀ v N ∈ P N 0 -(u_{N}^{''} , v_{N})_{\omega} + \alpha(u_{N}, v_{N})_{\omega} = (I_{N}f, v_{N})_{\omega} \quad \forall v_{N} \in P_{N}^{0} (uN,vN)ω+α(uN,vN)ω=(INf,vN)ωvNPN0:
将试函数带入到变分形式中,得到 G a l e r k i n Galerkin Galerkin方程组:
− ∑ i = 0 N − 2 ( ϕ i ′ ′ , ϕ j ) ω u ^ i + α ∑ i = 0 N − 2 ( ϕ i , ϕ j ) ω u ^ i = ( f N , ϕ j ) ω ,   j = 0 , ⋯   , N − 2 -\sum\limits_{i = 0}^{N-2}( \phi^{''}_{i}, \phi_{j})_{\omega}\widehat{u}_{i} + \alpha\sum\limits_{i = 0}^{N-2}\left( \phi_{i}, \phi_{j}\right)_{\omega}\widehat{u}_{i} = (f_{N}, \phi_{j})_{\omega} , \ j = 0,\cdots, N-2 i=0N2(ϕi,ϕj)ωu i+αi=0N2(ϕi,ϕj)ωu i=(fN,ϕj)ω, j=0,,N2
改写方程组为矩阵形式得到:
( S + α M ) U ^ = f (S + \alpha M) \widehat{U} = f (S+αM)U =f
其中
( S ) i j = − ( ϕ j ′ ′ ( x ) , ϕ i ( x ) ) ω (S)_{ij} = -(\phi_{j}^{''}(x), \phi_{i}(x))_{\omega} (S)ij=(ϕj(x),ϕi(x))ω 称为 刚度矩阵
( M ) i j = ( ϕ j ( x ) , ϕ i ( x ) ) ω (M)_{ij} = (\phi_{j}(x), \phi_{i}(x))_{\omega} (M)ij=(ϕj(x),ϕi(x))ω 称为 质量矩阵
f = ( ( I N f , ϕ 0 ( x ) ) ω , ⋯   , ( I N f , ϕ N − 2 ( x ) ) ω ) T \mathbf{f} = \left((I_{N}f, \phi_{0}(x))_{\omega}, \cdots, (I_{N}f, \phi_{N-2}(x))_{\omega} \right)^{T} f=((INf,ϕ0(x))ω,,(INf,ϕN2(x))ω)T
U ^ = ( u ^ 0 ⋯ u ^ N − 2 ) T \widehat{U} = (\widehat{u}_{0} \cdots \widehat{u}_{N-2})^{T} U =(u 0u N2)T
计算要点就是通过取不同的基函数来得到不同的代数系统进行求解.

Legendre-Galerkin方法

Legendre多项式及其性质

首先给出下面会用到的 L e g e n d r e Legendre Legendre多项式的一些基本性质,这里考虑区间 [ − 1 , 1 ] [-1, 1] [1,1]上的 L e g e n d r e Legendre Legendre多项式

  • 三项递推关系
    L 0 ( x ) = 1 , L 1 ( x ) = x ( n + 1 ) L n + 1 ( x ) = ( 2 n + 1 ) x L n ( x ) − n L n − 1 ( x ) L_{0}(x) = 1, \quad L_{1}(x) = x \\ (n+1)L_{n+1}(x) = (2n+1)xL_{n}(x) - nL_{n-1}(x) L0(x)=1,L1(x)=x(n+1)Ln+1(x)=(2n+1)xLn(x)nLn1(x)
  • 微分定义
    L n ( x ) = 1 2 n n ! d n d x n [ ( x 2 − 1 ) n ] L_{n}(x) = \frac{1}{2^{n}n!} \frac{d^{n}}{dx^{n}}\left[ (x^{2} - 1)^{n} \right] Ln(x)=2nn!1dxndn[(x21)n]
  • 正交性关系
    ∫ − 1 1 L k ( x ) L j ( x ) d x = 1 k + 1 / 2 δ k j \int_{-1}^{1} L_{k}(x)L_{j}(x) dx = \frac{1}{k + 1/ 2}\delta_{kj} 11Lk(x)Lj(x)dx=k+1/21δkj
    δ k j \delta_{kj} δkj K r o n e c k e r Kronecker Kronecker符号.
  • 二阶导数关系
    L n ′ ′ ( x ) = ∑ k = 0 ,   k + n   e v e n n − 2 ( k + 1 / 2 ) ( n ( n + 1 ) − k ( k + 1 ) ) L k ( x ) L^{\prime\prime}_{n}(x) = \sum\limits_{k = 0 ,\ k+n \ even}^{n-2} (k+1/2)\left( n(n+1) - k(k+1) \right)L_{k}(x) Ln(x)=k=0, k+n evenn2(k+1/2)(n(n+1)k(k+1))Lk(x)
  • L e g e n d r e − G a u s s − L o b a t t o Legendre-Gauss-Lobatto LegendreGaussLobatto型求积公式求积节点以及求积权重
    L e g e n d r e − G a u s s − L o b a t t o Legendre-Gauss-Lobatto LegendreGaussLobatto型求积公式求积节点: x 0 = − 1 , x N = 1 x_{0} = -1, x_{N} = 1 x0=1,xN=1, { x j } j = 1 N − 1 \left\lbrace x_{j} \right\rbrace_{j = 1}^{N-1} {xj}j=1N1 L N ′ ( x ) L^{\prime}_{N}(x) LN(x)的节点.
    权重:
    ω j = 2 N ( N + 1 ) 1 [ L N ( x j ) ] 2 , 0 ≤ j ≤ N \omega_{j} = \dfrac{2}{N(N+1)}\dfrac{1}{[L_{N}(x_{j})]^{2}}, \quad 0 \leq j \leq N ωj=N(N+1)2[LN(xj)]21,0jN

基函数的选取

基函数选取的方式和 C h e b y s h e v − G a l e r k i n Chebyshev-Galerkin ChebyshevGalerkin方法基本一致,可以参考上篇博文.
取如下形式的基函数:
ϕ k ( x ) = L k ( x ) + a k L k + 1 ( x ) + b k L k + 2 ( x ) \phi_{k}(x) = L_{k}(x)+ a_{k}L_{k+1}(x) + b_{k}L_{k+2}(x) ϕk(x)=Lk(x)+akLk+1(x)+bkLk+2(x)
通过将基函数带入边界条件, 求解出 a k , b k a_{k}, b_{k} ak,bk得到满足边界条件的基函数.
当边界条件适定时,满足要求的解是存在唯一的.
特别地
对于齐次 D i r i c h l e t Dirichlet Dirichlet型边界,有: a k = 0 ,   b k = − 1 a_{k} = 0, \ b_{k} = -1 ak=0, bk=1
对于齐次 N e u m a n Neuman Neuman型边界,有: a k = 0 ,   b k = − k ( k + 1 ) ( k + 2 ) ( k + 3 ) a_{k} = 0, \ b_{k} = -\dfrac{k(k+1)}{(k+2)(k+3)} ak=0, bk=(k+2)(k+3)k(k+1)

刚度矩阵的计算

这里仅对齐次 D i r i c h l e t Dirichlet Dirichlet型边界条件给出刚度矩阵的计算;其他边界条件计算方法完全一样.对于齐次 D i r i c h l e t Dirichlet Dirichlet型边界条件,由上可知,基函数将取为:
ϕ k ( x ) = L k ( x ) − L k + 2 ( x ) \phi_{k}(x) = L_{k}(x) - L_{k+2}(x) ϕk(x)=Lk(x)Lk+2(x)
s i j = ( S ) i j s_{ij} = (S)_{ij} sij=(S)ij, 首先我们将说明,在齐次 D i r i c h l e t Dirichlet Dirichlet边界条件下刚度矩阵是对称的:
s i j = − ∫ − 1 1 ϕ j ′ ′ ( x ) ϕ i ( x ) d x = ∫ − 1 1 ϕ j ′ ( x ) ϕ i ′ ( x ) d x = − ∫ − 1 1 ϕ i ′ ′ ( x ) ϕ j ( x ) d x = s j i \begin{aligned} s_{ij} &= -\int_{-1}^{1} \phi^{\prime\prime}_{j}(x)\phi_{i}(x) dx \\ & = \int_{-1}^{1} \phi^{\prime}_{j}(x)\phi^{\prime}_{i}(x) dx \\ & = -\int_{-1}^{1} \phi^{\prime\prime}_{i}(x)\phi_{j}(x) dx = s_{ji} \end{aligned} sij=11ϕj(x)ϕi(x)dx=11ϕj(x)ϕi(x)dx=11ϕi(x)ϕj(x)dx=sji
由:
s i j = − ∫ − 1 1 ( L j ′ ′ ( x ) − L j + 2 ′ ′ ( x ) ) ( L i ( x ) − L i + 2 ( x ) ) d x = ( 2 j + 3 ) ∫ − 1 1 [ ∑ k = 0 ,   k + j   e v e n j − 2 ( 2 k + 1 ) L k ( x ) + ( 2 j + 1 ) L j ( x ) ] [ L i ( x ) − L i + 2 ( x ) ] d x \begin{aligned} s_{ij} &= - \int_{-1}^{1} \left( L^{\prime\prime}_{j}(x) - L^{\prime\prime}_{j+2}(x)\right)\left( L_{i}(x) - L_{i+2}(x)\right) dx \\ & = (2j + 3) \int_{-1}^{1}\left[ \sum\limits_{k = 0, \ k+j \ even}^{j-2} (2k+1)L_{k}(x) + (2j+1)L_{j}(x) \right]\left[ L_{i}(x) - L_{i+2}(x)\right] dx \end{aligned} sij=11(Lj(x)Lj+2(x))(Li(x)Li+2(x))dx=(2j+3)11k=0, k+j evenj2(2k+1)Lk(x)+(2j+1)Lj(x)[Li(x)Li+2(x)]dx
由于矩阵是对称的,不妨设 i ≥ j i \geq j ij,利用 L e g e n d r e Legendre Legendre多项式的正交性:

  1. i > j i>j i>j 时, 显然 s i j = 0 s_{ij} = 0 sij=0 利用对称性,立马得到: s j i = 0 s_{ji} = 0 sji=0
  2. i = j i =j i=j 时,
    s i i = ∫ − 1 1 ( 2 i + 1 ) ( 2 i + 3 ) L i ( x ) L i ( x ) d x = ( 2 i + 1 ) ( 2 i + 3 ) 1 i + 1 / 2 = 4 i + 6 s_{ii} = \int_{-1}^{1}(2i+1)(2i+3)L_{i}(x)L_{i}(x) dx = (2i+1)(2i+3)\dfrac{1}{i+1/2} = 4i+6 sii=11(2i+1)(2i+3)Li(x)Li(x)dx=(2i+1)(2i+3)i+1/21=4i+6
    更一般的情况, 有刚度矩阵 s i i = − b i ( 4 i + 6 ) s_{ii} = -b_{i}(4i + 6) sii=bi(4i+6)

质量矩阵计算

显然,质量矩阵也是对称的
m j i = m i j = ( M ) i j = ∫ − 1 1 ϕ j ( x ) ϕ i ( x ) d x = ∫ − 1 1 ( L j ( x ) − L j + 2 ( x ) ) ( L i ( x ) − L i + 2 ( x ) ) d x \begin{aligned} &m_{ji} = m_{ij} = (M)_{ij} =\int_{-1}^{1} \phi_{j}(x)\phi_{i}(x) dx \\ & = \int_{-1}^{1} (L_{j}(x) - L_{j+2}(x))(L_{i}(x) - L_{i+2}(x)) dx \\ \end{aligned} mji=mij=(M)ij=11ϕj(x)ϕi(x)dx=11(Lj(x)Lj+2(x))(Li(x)Li+2(x))dx
展开,利用正交性容易计算:
m i j = m j i = { 2 2 i + 1 + 2 2 i + 5 , i = j 0 , i = j + 1 − 2 2 i + 1 , i = j + 2 m_{ij} = m_{ji} = \begin{cases} \frac{2}{2i+1} + \frac{2}{2i+5}, & i = j \\ 0, & i =j+1 \\ -\frac{2}{2i+1}, & i =j+2 \end{cases} mij=mji=2i+12+2i+52,0,2i+12,i=ji=j+1i=j+2

右端项计算

推导右端项,首先要计算插值多项式 I N f I_{N}f INf, 实际上即 L e g e n d r e Legendre Legendre变换.
假设插值节点取为 G a u s s − L e g e n d r e − L o b a t t o Gauss-Legendre-Lobatto GaussLegendreLobatto型求积公式的求积节点,求积节点以及权重在上面给出.
设插值多项式:
I N f ( x ) = ∑ n = 0 N f ~ n L n ( x ) I_{N}f(x) = \sum\limits_{n = 0}^{N}\widetilde{f}_{n}L_{n}(x) INf(x)=n=0Nf nLn(x)
其中,谱系数将由 L e g e n d r e Legendre Legendre正变换给出:
f ~ i = 1 γ i ∑ k = 0 N ω k f ( x k ) L i ( x k ) , i = 0 , ⋯   , N \widetilde{f}_{i} = \dfrac{1}{\gamma_{i}} \sum\limits_{k = 0}^{N} \omega_{k}f(x_{k})L_{i}(x_{k}), \quad i = 0, \cdots, N f i=γi1k=0Nωkf(xk)Li(xk),i=0,,N
其中
γ i = { 2 2 i + 1 , 0 ≤ i ≤ N − 1 2 N , i = N \gamma_{i} = \begin{cases} \frac{2}{2i+1}, & 0 \leq i \leq N - 1 \\ \frac{2}{N}, & i = N \end{cases} γi={2i+12,N2,0iN1i=N
则有:
右端项
( f ) i = ( I N f , ϕ i ( x ) ) ω = ∑ n = 0 N f ~ n ( L n ( x ) , L i ( x ) − L i + 2 ( x ) ) = 2 2 i + 1 f ~ i − 2 2 i + 5 f ~ i + 2    ( 0 ≤ i ≤ N − 2 ) \begin{aligned} (\mathbf{f})_{i} & = (I_{N}f, \phi_{i}(x))_{\omega} = \sum_{n = 0}^{N}\widetilde{f}_{n}\left( L_{n}(x), L_{i}(x) - L_{i+2}(x)\right) \\ & = \frac{2}{2i+1}\widetilde{f}_{i} - \frac{2}{2i + 5}\widetilde{f}_{i+2} \ \ (0 \leq i \leq N - 2) \end{aligned} (f)i=(INf,ϕi(x))ω=n=0Nf n(Ln(x),Li(x)Li+2(x))=2i+12f i2i+52f i+2  (0iN2)

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

经过上面的步骤已经可以求得基函数的展开系数,希望最终得到解 u N = ∑ i = 0 N − 2 u ^ i ϕ i ( x ) u_{N} = \sum\limits_{i = 0}^{N - 2} \widehat{u}_{i}\phi_{i}(x) uN=i=0N2u iϕi(x)在各个节点处的的值,这一点可以利用 L e g e n d r e Legendre Legendre逆变换求得,由于 L e g e n d r e Legendre Legendre逆变换计算复杂度为 O ( N 2 ) O(N^{2}) O(N2),因此这里不再做过多赘述,只给出齐次边界条件基函数的下的 L e g e n d r e Legendre Legendre逆变换, L e g e n d r e − G a l e r k i n Legendre-Galerkin LegendreGalerkin方法中所有关于 L e g e n d r e Legendre Legendre变换的部分将由 C h e b y s h e v − L e g e n d r e − G a l e r k i n Chebyshev-Legendre-Galerkin ChebyshevLegendreGalerkin谱方法作改进.
由:
u N ( x j ) = ∑ i = 0 N − 2 u ^ i ϕ i ( x j ) = ∑ i = 0 N − 2 u ^ i ( L i ( x j ) − L i + 2 ( x j ) ) u_{N}(x_{j}) = \sum\limits_{i = 0}^{N - 2} \widehat{u}_{i}\phi_{i}(x_{j}) = \sum\limits_{i = 0}^{N - 2} \widehat{u}_{i} (L_{i}(x_{j}) - L_{i+2}(x_{j})) uN(xj)=i=0N2u iϕi(xj)=i=0N2u i(Li(xj)Li+2(xj))
这个变换是平凡的,在计算中会创建一个有 L e g e n d r e Legendre Legendre多项式组成的数组,这里直接对数组作切片即可.
最后,给出 L e g e n d r e − G a l e r k i n Legendre-Galerkin LegendreGalerkin谱方法的整个求解流程:

  1. 计算 L G L LGL LGL求积节点和求积权重(计算 L e g e n d r e Legendre Legendre多项式的零点有一整套算法,可以自行百度).
  2. 预存储 a k , b k a_{k}, b_{k} ak,bk.
  3. 用稀疏矩阵存储 S , M S, M S,M.
  1. 计算向量 ( f ( x 0 ) , ⋯   , f ( x N ) ) ⊺ (f(x_{0}), \cdots,f(x_{N}))^{\intercal} (f(x0),,f(xN))对其作 L e g e n d r e Legendre Legendre正变换:
    f ~ i = 1 γ i ∑ k = 0 N ω k f ( x k ) L i ( x k ) , i = 0 , ⋯   , N \widetilde{f}_{i} = \dfrac{1}{\gamma_{i}} \sum\limits_{k = 0}^{N} \omega_{k}f(x_{k})L_{i}(x_{k}), \quad i = 0, \cdots, N f i=γi1k=0Nωkf(xk)Li(xk),i=0,,N
    将上式写为矩阵形式得:
    ( f ~ 0 ⋮ f ~ N ) = A ( f 0 ⋮ f N ) \begin{pmatrix} \widetilde{f}_{0} \\ \vdots \\ \widetilde{f}_{N} \end{pmatrix} = A \begin{pmatrix} f_{0} \\ \vdots \\ f_{N} \end{pmatrix} f 0f N=Af0fN
    A A A 可以用以下方式创建: A = Γ . ∗ L . ∗ Ω A = \Gamma .* \mathcal{L} .* \Omega A=Γ.L.Ω
    Γ = ( 1 / γ 0 ⋮ 1 / γ n ) , L = ( L 0 ( x 0 ) ⋯ L 0 ( x N ) ⋮ ⋮ L N ( x 0 ) ⋯ L N ( x N ) ) , Ω = ( ω 0 , ⋯   , ω N ) \Gamma = \begin{pmatrix} 1 / \gamma_{0} \\ \vdots \\ 1/\gamma_{n} \end{pmatrix},\quad \mathcal{L} = \begin{pmatrix} L_{0}(x_{0}) & \cdots & L_{0}(x_{N}) \\ \vdots & & \vdots \\ L_{N}(x_{0}) & \cdots & L_{N}(x_{N}) \end{pmatrix},\quad \Omega = \begin{pmatrix} \omega_{0}, \cdots, \omega_{N} \end{pmatrix} Γ=1/γ01/γn,L=L0(x0)LN(x0)L0(xN)LN(xN),Ω=(ω0,,ωN)
    注意,此处 . ∗ .* . 借用 Matlab 中的点乘运算, 低版本的 Matlab 没有广播机制,需要将 Γ \Gamma Γ Ω \Omega Ω先对角化再点乘, Python用户直接用 ∗ * 进行运算即可.
  2. 计算右端项
    ( f ) i = ( I N f , ϕ i ( x ) ) ω = ∑ n = 0 N f ~ n ( L n ( x ) , L i ( x ) − L i + 2 ( x ) ) = 2 2 i + 1 f ~ i − 2 2 i + 5 f ~ i + 2    ( 0 ≤ i ≤ N − 2 ) \begin{aligned} (\mathbf{f})_{i} & = (I_{N}f, \phi_{i}(x))_{\omega} = \sum_{n = 0}^{N}\widetilde{f}_{n}\left( L_{n}(x), L_{i}(x) - L_{i+2}(x)\right) \\ & = \frac{2}{2i+1}\widetilde{f}_{i} - \frac{2}{2i + 5}\widetilde{f}_{i+2} \ \ (0 \leq i \leq N - 2) \end{aligned} (f)i=(INf,ϕi(x))ω=n=0Nf n(Ln(x),Li(x)Li+2(x))=2i+12f i2i+52f i+2  (0iN2)
    化为矩阵形式得:
    f = ( ξ 0 0 − η 0 ξ 1 0 − η 1 ⋱ ⋱ ⋱ ξ N − 2 0 − η N − 2 ) f ~ \mathbf{f} = \begin{pmatrix} \xi_{0}& 0 & -\eta_{0} & & & \\ & \xi_{1} & 0 & -\eta_{1} & & \\ & & \ddots & \ddots & \ddots & \\ & & & \xi_{N-2} & 0 & -\eta_{N-2} \end{pmatrix} \widetilde{f} f=ξ00ξ1η00η1ξN20ηN2f

求解方程
( s + α M ) U ^ = f (s + \alpha M) \widehat{U} = \mathbf{f} (s+αM)U =f
得到
U ^ = ( u ^ 0 ⋮ u ^ N − 2 ) \widehat{U} = \begin{pmatrix} \widehat{u}_{0} \\ \vdots \\ \widehat{u}_{N-2} \end{pmatrix} U =u 0u N2

作变换
u N ( x j ) = ∑ i = 0 N − 2 u ^ i ϕ i ( x j ) = ∑ i = 0 N − 2 u ^ i ( L i ( x j ) − L i + 2 ( x j ) ) u_{N}(x_{j}) = \sum\limits_{i = 0}^{N - 2} \widehat{u}_{i}\phi_{i}(x_{j}) = \sum\limits_{i = 0}^{N - 2} \widehat{u}_{i} (L_{i}(x_{j}) - L_{i+2}(x_{j})) uN(xj)=i=0N2u iϕi(xj)=i=0N2u i(Li(xj)Li+2(xj))
注意使用上面的 L \mathcal{L} L,将上式改写为矩阵:
U N = ( L [ 0 : N − 1 , : ] − L [ 2 : , : ] ) U ^ U_{N} = \left( \mathcal{L}[0:N-1, :] - \mathcal{L}[2:, :] \right) \widehat{U} UN=(L[0:N1,:]L[2:,:])U

Legendre-Chebyshev-Galerkin方法

L e g e n d r e − C h e b y s h e v − G a l e r k i n Legendre-Chebyshev-Galerkin LegendreChebyshevGalerkin方法同时结合了 L e g e n d r e − G a l e r k i n Legendre-Galerkin LegendreGalerkin方法以及 C h e b y s h e v − G a l e r k i n Chebyshev-Galerkin ChebyshevGalerkin方法的优势.
L e g e n d r e − C h e b y s h e v − G a l e r k i n Legendre-Chebyshev-Galerkin LegendreChebyshevGalerkin方法采取如下的变分形式:
u N ∈ P N 0 u_{N} \in P^{0}_{N} uNPN0, s.t.
− ( u N ′ ′ , v N ) ω + α ( u N , v N ) ω = ( I N C f , v N ) ω ∀ v N ∈ P N 0 -(u_{N}^{''} , v_{N})_{\omega} + \alpha(u_{N}, v_{N})_{\omega} = (I^{C}_{N}f, v_{N})_{\omega} \quad \forall v_{N} \in P_{N}^{0} (uN,vN)ω+α(uN,vN)ω=(INCf,vN)ωvNPN0
其中 I N C f I^{C}_{N}f INCf表示插值节点选取为 C h e b y s h e v − G a u s s − L o b a t t o Chebyshev-Gauss-Lobatto ChebyshevGaussLobatto求积节点的插值函数.

基函数,刚度矩阵, 质量矩阵

L e g e n d r e − C h e b y s h e v − G a l e r k i n Legendre-Chebyshev-Galerkin LegendreChebyshevGalerkin方法选取基函数与 L e g e n d r e − G a l e r k i n Legendre-Galerkin LegendreGalerkin方法选取基函数的方式完全一致,即:
ϕ k ( x ) = L k ( x ) + a k L k + 1 ( x ) + b k L k + 2 ( x ) \phi_{k}(x) = L_{k}(x)+ a_{k}L_{k+1}(x) + b_{k}L_{k+2}(x) ϕk(x)=Lk(x)+akLk+1(x)+bkLk+2(x)
由此,由上;可推导出和上面完全相同的刚度矩阵,质量矩阵:
( S ) i j = ( 4 i + 6 ) δ i j ( M ) i j = { 2 2 i + 1 + 2 2 i + 5 , i = j 0 , i = j + 1 − 2 2 i + 1 , i = j + 2 \begin{aligned} &(S)_{ij} = (4i+6) \delta_{ij} \\ &(M)_{ij} = \begin{cases} \frac{2}{2i+1} + \frac{2}{2i+5}, & i = j \\ 0, & i =j+1 \\ -\frac{2}{2i+1}, & i =j+2 \end{cases} \end{aligned} (S)ij=(4i+6)δij(M)ij=2i+12+2i+52,0,2i+12,i=ji=j+1i=j+2

右端项计算

关键是计算插值多项式 I N C f I^{C}_{N}f INCf, 称为 L e g e n d r e − C h e b y s h e v Legendre-Chebyshev LegendreChebyshev变换.
假设插值节点取为 G a u s s − C h e b y s h e v − L o b a t t o Gauss-Chebyshev-Lobatto GaussChebyshevLobatto型求积公式的求积节点,求积节点以及权重在上面给出.
设插值多项式:
I N C f ( x ) = ∑ n = 0 N f ~ n L n ( x ) I^{C}_{N}f(x) = \sum\limits_{n = 0}^{N}\widetilde{f}_{n}L_{n}(x) INCf(x)=n=0Nf nLn(x)
L e g e n d r e − C h e b y s h e v Legendre-Chebyshev LegendreChebyshev变换实际上即 G a u s s − C h e b y s h e v − L o b a t t o Gauss-Chebyshev-Lobatto GaussChebyshevLobatto求积节点处的值到 L e g e n d r e Legendre Legendre展开系数之间的变换, L e g e n d r e − C h e b y s h e v Legendre-Chebyshev LegendreChebyshev变换可以被分成两部分:

  • 计算 CGL 节点处的值到 C h e b y s h e v Chebyshev Chebyshev 展开系数之间的变换.
  • 计算Chebyshev 展开系数到 L e g e n d r e Legendre Legendre展开系数之间的变换.

第一步容易执行,即令:
I N f ( x j ) = ∑ n = 0 N f ^ n T n ( x j ) , x j = c o s π j N , j = 0 , ⋯   , N I_{N}f(x_{j}) = \sum\limits_{n = 0}^{N} \widehat{f}_{n}T_{n}(x_{j}), \quad x_{j} = cos\dfrac{\pi j}{N} , \quad j = 0,\cdots,N INf(xj)=n=0Nf nTn(xj),xj=cosNπj,j=0,,N
f ^ n \widehat{f}_{n} f n将由 C h e b y s h e v Chebyshev Chebyshev变换:
f ^ n = 2 c ~ n N ∑ k = 0 N 1 c ~ k f ( x k ) c o s ( n k π / N ) \begin{aligned} \widehat{f}_{n} = \dfrac{2}{\widetilde{c}_{n}N}\sum\limits_{k = 0}^{N}\dfrac{1} {\widetilde{c}_{k}}f(x_{k})cos(nk\pi/N) \end{aligned} f n=c nN2k=0Nc k1f(xk)cos(nkπ/N)
给出.
2.
第二步即,设:
f ( x j ) = I N f ( x j ) = ∑ n = 0 N f ^ n T n ( x j ) f ( x j ) = I N c f ( x j ) = ∑ n = 0 N f ~ n L n ( x j ) f(x_{j}) = I_{N}f(x_{j}) = \sum\limits_{n = 0}^{N} \widehat{f}_{n}T_{n}(x_{j}) \\ f(x_{j}) = I_{N}^{c}f(x_{j}) = \sum\limits_{n = 0}^{N} \widetilde{f}_{n}L_{n}(x_{j}) f(xj)=INf(xj)=n=0Nf nTn(xj)f(xj)=INcf(xj)=n=0Nf nLn(xj)
目标,找出坐标变换关系,即求矩阵 S S S,使得:
( f ~ 0 ⋮ f ~ N ) = S ( f ^ 0 ⋮ f ^ N ) \begin{pmatrix} \widetilde{f}_{0} \\ \vdots \\ \widetilde{f}_{N} \end{pmatrix} = S \begin{pmatrix} \widehat{f}_{0} \\ \vdots \\ \widehat{f}_{N} \end{pmatrix} f 0f N=Sf 0f N
由:
∑ n = 0 N f ^ n T n ( x ) = ∑ n = 0 N f ~ n L n ( x ) \sum\limits_{n = 0}^{N} \widehat{f}_{n}T_{n}(x) = \sum\limits_{n = 0}^{N} \widetilde{f}_{n}L_{n}(x) n=0Nf nTn(x)=n=0Nf nLn(x)
对上述 S S S的计算在网上能查阅到很多资料以及文章,以及最后的坐标变换有相应的快速算法,使得最终整个变换的复杂度为 O ( N l o g N ) O(NlogN) O(NlogN),由于过于复杂这里不再给出.

注: 同上一篇最后,代码以及参考资料会在以后补出

本文链接http://www.taodudu.cc/news/show-1944905.html