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
Chebyshev−Galerkin方法的计算过程,本文在上面的基础上,继续探讨
L
e
n
g
e
n
d
r
e
−
G
a
r
l
e
r
k
i
n
Lengendre-Garlerkin
Lengendre−Garlerkin方法以及
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
Legendre−Chebyshev−Galerkin方法,可以看到,
C
h
e
b
y
s
h
e
v
−
G
a
l
e
r
k
i
n
Chebyshev-Galerkin
Chebyshev−Galerkin在作正变换以及逆变换的过程中,能够作快速余弦变换,但该方法所推导出的刚度矩阵稀疏程度不够; 最关键的是,刚度矩阵不是对称的!.在本文介绍的
L
e
n
g
e
n
d
r
e
−
G
a
r
l
e
r
k
i
n
Lengendre-Garlerkin
Lengendre−Garlerkin方法中,可以看到刚度矩阵是对角矩阵,但
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
Chebyshev−Legendre−Galerkin方法同时具备了
C
h
e
b
y
s
h
e
v
−
G
a
l
e
r
k
i
n
Chebyshev - Galerkin
Chebyshev−Galerkin方法以及
L
e
g
e
n
d
r
e
−
G
a
l
e
r
k
i
n
Legendre-Galerkin
Legendre−Galerkin方法的优势.
回忆
首先,回忆上篇博文中的重要内容
考虑如下一维齐次
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=0∑N−2aiϕ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
Ritsz−Galerkin离散得到如下变分形式
求
u
N
∈
P
N
0
u_{N} \in P^{0}_{N}
uN∈PN0, 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)ω∀vN∈PN0:
将试函数带入到变分形式中,得到
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=0∑N−2(ϕi′′,ϕj)ωu
i+αi=0∑N−2(ϕi,ϕj)ωu
i=(fN,ϕj)ω, j=0,⋯,N−2
改写方程组为矩阵形式得到:
(
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,ϕN−2(x))ω)T
U
^
=
(
u
^
0
⋯
u
^
N
−
2
)
T
\widehat{U} = (\widehat{u}_{0} \cdots \widehat{u}_{N-2})^{T}
U
=(u
0⋯u
N−2)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)−nLn−1(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[(x2−1)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 even∑n−2(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
Legendre−Gauss−Lobatto型求积公式求积节点以及求积权重
L e g e n d r e − G a u s s − L o b a t t o Legendre-Gauss-Lobatto Legendre−Gauss−Lobatto型求积公式求积节点: 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=1N−1为 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,0≤j≤N
基函数的选取
基函数选取的方式和
C
h
e
b
y
s
h
e
v
−
G
a
l
e
r
k
i
n
Chebyshev-Galerkin
Chebyshev−Galerkin方法基本一致,可以参考上篇博文.
取如下形式的基函数:
ϕ
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)∫−11⎣⎡k=0, k+j even∑j−2(2k+1)Lk(x)+(2j+1)Lj(x)⎦⎤[Li(x)−Li+2(x)]dx
由于矩阵是对称的,不妨设
i
≥
j
i \geq j
i≥j,利用
L
e
g
e
n
d
r
e
Legendre
Legendre多项式的正交性:
- i > j i>j i>j 时, 显然 s i j = 0 s_{ij} = 0 sij=0 利用对称性,立马得到: s j i = 0 s_{ji} = 0 sji=0
-
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
Gauss−Legendre−Lobatto型求积公式的求积节点,求积节点以及权重在上面给出.
设插值多项式:
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=0∑Nf
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=0∑Nω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,0≤i≤N−1i=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=0∑Nf
n(Ln(x),Li(x)−Li+2(x))=2i+12f
i−2i+52f
i+2 (0≤i≤N−2)
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=0∑N−2u
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
Legendre−Galerkin方法中所有关于
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
Chebyshev−Legendre−Galerkin谱方法作改进.
由:
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=0∑N−2u
iϕi(xj)=i=0∑N−2u
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
Legendre−Galerkin谱方法的整个求解流程:
- 计算 L G L LGL LGL求积节点和求积权重(计算 L e g e n d r e Legendre Legendre多项式的零点有一整套算法,可以自行百度).
- 预存储 a k , b k a_{k}, b_{k} ak,bk.
- 用稀疏矩阵存储 S , M S, M S,M.
- 计算向量
(
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=0∑Nω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 0⋮f N⎠⎟⎞=A⎝⎜⎛f0⋮fN⎠⎟⎞
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/γ0⋮1/γn⎠⎟⎞,L=⎝⎜⎛L0(x0)⋮LN(x0)⋯⋯L0(xN)⋮LN(xN)⎠⎟⎞,Ω=(ω0,⋯,ωN)
注意,此处 . ∗ .* .∗ 借用 Matlab 中的点乘运算, 低版本的 Matlab 没有广播机制,需要将 Γ \Gamma Γ和 Ω \Omega Ω先对角化再点乘, Python用户直接用 ∗ * ∗进行运算即可. - 计算右端项
( 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=0∑Nf n(Ln(x),Li(x)−Li+2(x))=2i+12f i−2i+52f i+2 (0≤i≤N−2)
化为矩阵形式得:
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⋱ξN−2⋱0−ηN−2⎠⎟⎟⎞f
求解方程
(
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
0⋮u
N−2⎠⎟⎞
作变换
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=0∑N−2u
iϕi(xj)=i=0∑N−2u
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:N−1,:]−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
Legendre−Chebyshev−Galerkin方法同时结合了
L
e
g
e
n
d
r
e
−
G
a
l
e
r
k
i
n
Legendre-Galerkin
Legendre−Galerkin方法以及
C
h
e
b
y
s
h
e
v
−
G
a
l
e
r
k
i
n
Chebyshev-Galerkin
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
Legendre−Chebyshev−Galerkin方法采取如下的变分形式:
求
u
N
∈
P
N
0
u_{N} \in P^{0}_{N}
uN∈PN0, 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)ω∀vN∈PN0
其中
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
Chebyshev−Gauss−Lobatto求积节点的插值函数.
基函数,刚度矩阵, 质量矩阵
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
Legendre−Chebyshev−Galerkin方法选取基函数与
L
e
g
e
n
d
r
e
−
G
a
l
e
r
k
i
n
Legendre-Galerkin
Legendre−Galerkin方法选取基函数的方式完全一致,即:
ϕ
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
Legendre−Chebyshev变换.
假设插值节点取为
G
a
u
s
s
−
C
h
e
b
y
s
h
e
v
−
L
o
b
a
t
t
o
Gauss-Chebyshev-Lobatto
Gauss−Chebyshev−Lobatto型求积公式的求积节点,求积节点以及权重在上面给出.
设插值多项式:
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=0∑Nf
nLn(x)
L
e
g
e
n
d
r
e
−
C
h
e
b
y
s
h
e
v
Legendre-Chebyshev
Legendre−Chebyshev变换实际上即
G
a
u
s
s
−
C
h
e
b
y
s
h
e
v
−
L
o
b
a
t
t
o
Gauss-Chebyshev-Lobatto
Gauss−Chebyshev−Lobatto求积节点处的值到
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
Legendre−Chebyshev变换可以被分成两部分:
- 计算 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=0∑Nf
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=0∑Nc
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=0∑Nf
nTn(xj)f(xj)=INcf(xj)=n=0∑Nf
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
0⋮f
N⎠⎟⎞=S⎝⎜⎛f
0⋮f
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=0∑Nf
nTn(x)=n=0∑Nf
nLn(x)
对上述
S
S
S的计算在网上能查阅到很多资料以及文章,以及最后的坐标变换有相应的快速算法,使得最终整个变换的复杂度为
O
(
N
l
o
g
N
)
O(NlogN)
O(NlogN),由于过于复杂这里不再给出.
注: 同上一篇最后,代码以及参考资料会在以后补出
标签：Chebyshev,Galerkin,phi,Legendre,widehat,2i 来源： https://blog.csdn.net/liuxinglongqi/article/details/111564923