DG方法：一维ODE
作者：互联网
DG方法：一维ODE
考虑一维ODE的边值问题：
{
u
x
=
f
(
x
)
u
(
0
)
=
a
,
x
∈
[
0
,
1
]
\left\{\begin{matrix}u_x=f(x)\\ u(0)=a\end{matrix}\right. \space,\space x\in [0,1]
{ux=f(x)u(0)=a , x∈[0,1]
有限差分方法
对 x ∈ [ 0 , 1 ] x\in[0,1] x∈[0,1]做uniform划分： x j = j Δ x = j 1 N x_j=j\Delta x=j \frac{1}{N} xj=jΔx=jN1
如果是一阶FD：对于第 j j j个网格点，需要利用 u j − 1 , u j , u j + 1 u_{j-1}, u_j, u_{j+1} uj−1,uj,uj+1等信息逼近 u x u_x ux，其中 u j u_j uj表示准确解 u ( x j ) u(x_j) u(xj)的近似，则
u x ∣ x j ≈ u j − u j − 1 Δ x u_x|_{x_j} \approx\frac{u_j-u_{j-1}}{\Delta x} ux∣xj≈Δxuj−uj−1
要选择用 x j − 1 x_{j-1} xj−1是因为边界条件定义在左侧，信息从迎风的左边来。代入原式有
u j − u j − 1 Δ x = f ( x j ) ⇒ u j = u j − 1 + Δ x f ( x j ) \frac{u_j-u_{j-1}}{\Delta x}=f(x_j) \Rightarrow u_j=u_{j-1}+\Delta x f(x_j) Δxuj−uj−1=f(xj)⇒uj=uj−1+Δxf(xj)
所以求解过程可以从左向右扫描（sweep，最左边是 u 0 = a u_0=a u0=a），非常易于计算。
但只有一阶精度：
∣ u ( x j ) − u j ∣ ≈ O ( Δ x ) |u(x_j)-u_j|\approx O(\Delta x) ∣u(xj)−uj∣≈O(Δx)
当从 N N N个网格点加密到 2 N 2N 2N时，工作量增加一倍，但误差仅能减少一半。
而如果使用更高阶精度的差分方法，同理，将 u x ( x j ) u_x(x_j) ux(xj)表示为 u j − 2 , u j − 1 , u j u_{j-2}, u_{j-1}, u_{j} uj−2,uj−1,uj的线性组合（可通过Taylor展开推导）。这对于 j = 2 , 3 , . . . , N j=2,3,...,N j=2,3,...,N是可行的，但 j = 1 j=1 j=1处显然没有这么多stencil可供构造，所以在边界上使用一阶格式，但这又会导致计算降阶！因为 u 1 u_1 u1只有一阶精度，那么根据 u 1 u_1 u1计算出的 u 2 u_2 u2也达不到二阶，如此类推。
Discrete Galerkin
格式导出
DG希望既易于计算，又能有高阶精度。它首先定义一系列半网格点 x 2 k + 1 2 ( k = 0 , 1 , . . . , N ) x_{\frac{2k+1}{2}}(k=0,1,...,N) x22k+1(k=0,1,...,N)，将区间划分为很多个单元，单元 I j = [ x j − 1 2 , x j + 1 2 ] I_j=[x_{j-\frac{1}{2}}, x_{j+\frac{1}{2}}] Ij=[xj−21,xj+21]。注意仅要求 0 = x 1 2 < x 3 2 < . . . < x N − 1 2 < x N + 1 2 = 1 0=x_{\frac{1}{2}}<x_{\frac{3}{2}}<...<x_{N-\frac{1}{2}}<x_{N+\frac{1}{2}}=1 0=x21<x23<...<xN−21<xN+21=1，而不需要每个单元都是均匀长度的。每个单元的中点 x j = 1 2 ( x j − 1 2 + x j + 1 2 ) x_j=\frac{1}{2}(x_{j-\frac{1}{2}}+x_{j+\frac{1}{2}}) xj=21(xj−21+xj+21)，长度 Δ x j = x j − 1 2 − x j + 1 2 \Delta x_j=x_{j-\frac{1}{2}}-x_{j+\frac{1}{2}} Δxj=xj−21−xj+21。记 h = max j Δ x j h=\underset{j}{\max} \Delta x_j h=jmaxΔxj。
定义解空间（solution space）：
V
h
k
:
=
{
v
:
v
∣
I
j
∈
P
k
(
I
j
)
,
j
=
1
,
2
,
.
.
.
,
N
}
\bm{V}_h^k:=\{v: v|_{I_j} \in \bm{P}^k(I_j),j=1,2,...,N\}
Vhk:={v:v∣Ij∈Pk(Ij),j=1,2,...,N}
其中 P k ( I j ) \bm{P}^k(I_j) Pk(Ij)表示第 j j j个单元上的次数小于等于 k k k的多项式空间。当 k = 1 k=1 k=1时即为分片线性函数（piecewise linear functions）。注意在cell边界处，函数值时可以间断的，而区别于传统有限元要求内部单元的边界值连续。
由于解在单元边界处可以完全不连续，因此 x j + 1 2 x_{j+\frac{1}{2}} xj+21处存在一个左值 u j + 1 2 − u_{j+\frac{1}{2}}^- uj+21−和一个右值 u j + 1 2 + u_{j+\frac{1}{2}}^+ uj+21+，分别属于左右两个单元。
不同单元的多项式空间 P k \bm{P}^k Pk的阶数 k k k可以不同（ P \bm{P} P adaptivity）。
定义测试函数（test function） v ( x ) v(x) v(x)，原问题成立，则如下积分式一定成立：
∫ I j u x v = ∫ I j f ( x ) v , j = 1 , 2 , . . . , N \int _{I_j} u_xv =\int_{I_j}f(x)v \space, \space j=1,2,...,N ∫Ijuxv=∫Ijf(x)v , j=1,2,...,N
假定 v ( x ) v(x) v(x)是至少可微的，那么分部积分：
∫ I j f ( x ) v = ∫ I j u x v = − ∫ I j u v x + u ( x j + 1 2 ) v ( x j + 1 2 ) − u ( x j − 1 2 ) v ( x j − 1 2 ) \int_{I_j}f(x)v=\int _{I_j} u_xv =-\int _{I_j} uv_x + u(x_{j+\frac{1}{2}})v(x_{j+\frac{1}{2}}) - u(x_{j-\frac{1}{2}})v(x_{j-\frac{1}{2}}) ∫Ijf(x)v=∫Ijuxv=−∫Ijuvx+u(xj+21)v(xj+21)−u(xj−21)v(xj−21)
所以问题转换为：
find u h ( x ) ∈ V h k u_h(x) \in \bm{V}_h^k uh(x)∈Vhk来逼近 u ( x ) u(x) u(x)，代入上式有
− ∫ I j u v x + ( u h v ) ∣ x j + 1 2 − ( u h v ) ∣ x j − 1 2 = ∫ I j f ( x ) v -\int _{I_j} uv_x +(u_hv)|_{x_{j+\frac{1}{2}}} - (u_hv)|_{x_{j-\frac{1}{2}}}=\int_{I_j}f(x)v −∫Ijuvx+(uhv)∣xj+21−(uhv)∣xj−21=∫Ijf(x)v
但问题是 v ∈ v\in v∈？，如果任意的可微函数都可以，那么自由度太多会造成over-determined！注意到
u h k ∈ V h k : = { v : v ∣ I j ∈ P k ( I j ) , j = 1 , 2 , . . . , N } u_h^k \in \bm{V}_h^k:=\{v: v|_{I_j} \in \bm{P}^k(I_j),j=1,2,...,N\} uhk∈Vhk:={v:v∣Ij∈Pk(Ij),j=1,2,...,N}
在每个单元 I j I_j Ij上确定 P k ( I j ) \bm{P}^k(I_j) Pk(Ij)需要 k + 1 k+1 k+1个参数，将 P k ( I j ) \bm{P}^k(I_j) Pk(Ij)这个线性空间的基取为 φ j 0 ( x ) , φ j 1 ( x ) , . . . , φ j k ( x ) \varphi_j^0(x),\varphi_j^1(x),...,\varphi_j^k(x) φj0(x),φj1(x),...,φjk(x)，则线性表出为
u h ( x ) ∣ I j = ∑ l = 0 k u j l φ j l ( x ) u_h(x)|_{I_j}=\sum_{l=0}^k u_j^l\varphi_j^l(x) uh(x)∣Ij=l=0∑kujlφjl(x)
因为共有 N N N个单元，则 dim V h k = N ( k + 1 ) \dim V_h^k=N(k+1) dimVhk=N(k+1)。需要求解的值是线性组合的系数 u j l , l = 0 , 1 , . . . , k , j = 1 , 2 , . . . , N u_j^l\space,\space l=0,1,...,k,\space j=1,2,...,N ujl , l=0,1,...,k, j=1,2,...,N，因此需要有 N ( k + 1 ) N(k+1) N(k+1)个约束条件。假定 v ∈ W h v\in \bm{W}_h v∈Wh，那么应该要有 dim W h = dim V h k = N ( k + 1 ) \dim \bm{W}_h=\dim \bm{V}_h^k=N(k+1) dimWh=dimVhk=N(k+1)。
这两者是可以取成不一样的线性空间的，而如果取成一样的（ W h = V h k \bm{W}_h=\bm{V}_h^k Wh=Vhk，都来自分片线性多项式的空间），则成为Galerkin方法。
于是问题转换为：find u h ∈ V h k , v ∈ V h k u_h \in \bm{V}_h^k\space, \space v \in \bm{V}_h^k uh∈Vhk , v∈Vhk ……
但仍不对，因为单元边界上有两个值（左值 u j + 1 2 − u_{j+\frac{1}{2}}^- uj+21−和右值 u j + 1 2 + u_{j+\frac{1}{2}}^+ uj+21+）。需要额外加条件。而我们希望当 k = 0 k=0 k=0时，DG方法能还原回普通的一阶（迎风）有限差分方法： u j − u j − 1 Δ x = f ( x j ) \frac{u_j - u_{j-1}}{\Delta x} = f(x_j) Δxuj−uj−1=f(xj)。
此时解空间和测试函数空间都为分片常值函数。由于测试函数是任意的，取 v = { 1 , I j 0 , o t h e r v=\left\{\begin{matrix}1\space,\space \space \space I_j\\ 0, \space\rm{other} \end{matrix}\right. v={1 , Ij0, other，则 v x = 0 v_x=0 vx=0， ∫ I j u h v x = 0 \int_{I_j}u_hv_x=0 ∫Ijuhvx=0。于是有
u h ( x j + 1 2 ) v ( x j + 1 2 ) − u h ( x j − 1 2 ) v ( x j − 1 2 ) = ∫ I j f ( x ) u_h(x_{j+\frac{1}{2}})v(x_{j+\frac{1}{2}})-u_h(x_{j-\frac{1}{2}})v(x_{j-\frac{1}{2}})=\int_{I_j}f(x) uh(xj+21)v(xj+21)−uh(xj−21)v(xj−21)=∫Ijf(x)
希望能还原出： u j − u j − 1 = Δ x f ( x j ) u_j-u_{j-1}=\Delta x f(x_j) uj−uj−1=Δxf(xj)。而在 k = 0 k=0 k=0时， u j u_j uj就是第 j j j个单元的函数近似解的值 u h ( I j ) u_h(I_j) uh(Ij)。
注意到 v ( x ) v(x) v(x)的取值方式，因此取 v ( x j + 1 2 ) = v ( x j + 1 2 − ) , u h ( x j + 1 2 ) = u h ( x j + 1 2 − ) = u j v(x_{j+\frac{1}{2}})=v(x_{j+\frac{1}{2}}^-), u_h(x_{j+\frac{1}{2}})=u_h(x_{j+\frac{1}{2}}^-)=u_j v(xj+21)=v(xj+21−),uh(xj+21)=uh(xj+21−)=uj，和 v ( x j − 1 2 ) = v ( x j − 1 2 + ) , u h ( x j − 1 2 ) = u h ( x j − 1 2 − ) = u j − 1 v(x_{j-\frac{1}{2}})=v(x_{j-\frac{1}{2}}^+), u_h(x_{j-\frac{1}{2}})=u_h(x_{j-\frac{1}{2}}^-)=u_{j-1} v(xj−21)=v(xj−21+),uh(xj−21)=uh(xj−21−)=uj−1。
所以完整的DG格式为：
find u h ( x ) ∈ V h k u_h(x) \in \bm{V}_h^k uh(x)∈Vhk，使得对于任意的测试函数 v ∈ V h k v\in\bm{V}_h^k v∈Vhk和任意的单元 I j I_j Ij，都有
− ∫ I j u h ( x ) v x d x + u j + 1 2 − v j + 1 2 − − u j − 1 2 − v j − 1 2 + = ∫ I j f ( x ) v ( x ) d x (1) -\int_{I_j}u_h(x)v_x {\rm d}x + u_{j+\frac{1}{2}}^{-} v_{j+\frac{1}{2}}^{-} - u_{j-\frac{1}{2}}^{-} v_{j-\frac{1}{2}}^{+}=\int_{I_j}f(x)v(x){\rm d}x \tag{1} −∫Ijuh(x)vxdx+uj+21−vj+21−−uj−21−vj−21+=∫Ijf(x)v(x)dx(1)
此时的格式与计算中如何选取 P k \bm{P}^k Pk的基函数的具体形式无关！基函数可以有多种选择，例如Legendre多项式函数，Lagrange插值基函数，自然幂函数（ 1 , x , x 2 , . . . , x k 1,x,x^2,...,x^k 1,x,x2,...,xk），和 1 , x − x j Δ x j , ( x − x j Δ x j ) 2 , . . . , ( x − x j Δ x j ) k 1,\frac{x-x_j}{\Delta x_j}, \left(\frac{x-x_j}{\Delta x_j}\right)^2, ..., \left(\frac{x-x_j}{\Delta x_j}\right)^k 1,Δxjx−xj,(Δxjx−xj)2,...,(Δxjx−xj)k等等。
求解方式
一般地，假定选取的基函数为
φ
j
0
(
x
)
,
φ
j
1
(
x
)
,
.
.
.
,
φ
j
k
(
x
)
\varphi_j^0(x),\varphi_j^1(x),...,\varphi_j^k(x)
φj0(x),φj1(x),...,φjk(x)，线性表出
u
h
(
x
)
∣
I
j
=
∑
l
=
0
k
u
j
l
φ
j
l
(
x
)
u_h(x)|_{I_j}=\sum_{l=0}^k u_j^l \varphi_j^l(x)
uh(x)∣Ij=∑l=0kujlφjl(x)。需要求解向量
u
j
=
[
u
j
0
u
j
1
.
.
.
u
j
k
]
\bm{u}_j=\left[ \begin{matrix} u_j^0\\ u_j^1\\ ...\\u_j^k \end{matrix} \right]
uj=⎣⎢⎢⎡uj0uj1...ujk⎦⎥⎥⎤
由于测试函数是任意的，只要在线性空间 W h = V h k \bm{W}_h=\bm{V}_h^k Wh=Vhk内就行。所以不妨取 v = φ j m ( x ) , m = 0 , 1 , . . . , k v=\varphi_j^m(x), \space m=0,1,...,k v=φjm(x), m=0,1,...,k，则上一节中的格式写为
− ∫ I j ∑ l = 0 k u j l φ j l ( x ) ( φ j m ) x d x + ∑ l = 0 k u j l φ j l ( x j + 1 2 ) φ j m ( x j + 1 2 ) − u j − 1 2 − φ j m ( x j − 1 2 ) = ∫ I j f ( x ) φ j m ( x ) d x -\int_{I_j}\sum_{l=0}^k u_j^l \varphi_j^l(x) (\varphi_j^m)_x {\rm d}x + \sum_{l=0}^k u_j^l \varphi_j^l(x_{j+\frac{1}{2}}) \varphi_j^m(x_{j+\frac{1}{2}}) - u_{j-\frac{1}{2}}^{-} \varphi_j^m(x_{j-\frac{1}{2}})=\int_{I_j}f(x) \varphi_j^m(x) {\rm d}x −∫Ijl=0∑kujlφjl(x)(φjm)xdx+l=0∑kujlφjl(xj+21)φjm(xj+21)−uj−21−φjm(xj−21)=∫Ijf(x)φjm(x)dx
上式要对于 m = 0 , 1 , . . . , k m=0,1,...,k m=0,1,...,k均成立。注意其中的 u j − 1 2 − u_{j-\frac{1}{2}}^{-} uj−21−是隔壁单元算出来的结果值，传给本单元的。上式提取出 ∑ l = 0 k u j l \sum_{l=0}^k u_j^l ∑l=0kujl，有
∑ l = 0 k u j l [ − ∫ I j φ j l ( x ) ( φ j m ) x d x + φ j l ( x j + 1 2 ) φ j m ( x j + 1 2 ) ] = u j − 1 2 − φ j m ( x j − 1 2 ) + ∫ I j f ( x ) φ j m ( x ) d x (2) \sum_{l=0}^k u_j^l\left[ - \int_{I_j} \varphi_j^l(x) (\varphi_j^m)_x {\rm d}x + \varphi_j^l(x_{j+\frac{1}{2}}) \varphi_j^m(x_{j+\frac{1}{2}}) \right] = u_{j-\frac{1}{2}}^{-} \varphi_j^m(x_{j-\frac{1}{2}}) + \int_{I_j}f(x) \varphi_j^m(x) {\rm d}x \tag{2} l=0∑kujl[−∫Ijφjl(x)(φjm)xdx+φjl(xj+21)φjm(xj+21)]=uj−21−φjm(xj−21)+∫Ijf(x)φjm(x)dx(2)
记：
a
m
l
=
−
∫
I
j
φ
j
l
(
x
)
(
φ
j
m
)
x
d
x
+
φ
j
l
(
x
j
+
1
2
)
φ
j
m
(
x
j
+
1
2
)
(2-1)
a_{ml}=- \int_{I_j} \varphi_j^l(x) (\varphi_j^m)_x {\rm d}x + \varphi_j^l(x_{j+\frac{1}{2}}) \varphi_j^m(x_{j+\frac{1}{2}}) \tag{2-1}
aml=−∫Ijφjl(x)(φjm)xdx+φjl(xj+21)φjm(xj+21)(2-1)
b m = u j − 1 2 − φ j m ( x j − 1 2 ) + ∫ I j f ( x ) φ j m ( x ) d x (2-2) b_{m} = u_{j-\frac{1}{2}}^{-} \varphi_j^m(x_{j-\frac{1}{2}}) + \int_{I_j}f(x) \varphi_j^m(x) {\rm d}x \tag{2-2} bm=uj−21−φjm(xj−21)+∫Ijf(x)φjm(x)dx(2-2)
则在每个单元内都组成一个线性方程组： A j u j = b j \bm{A}_j\bm{u}_j=\bm{b}_j Ajuj=bj，其中 A j = ( a m l ) \bm{A}_j = (a_{ml}) Aj=(aml)是一个 ( k + 1 ) × ( k + 1 ) (k+1)\times(k+1) (k+1)×(k+1)的矩阵。
一旦 u j − 1 2 − u_{j-\frac{1}{2}}^{-} uj−21−给定，就可以求解该线性方程组。由此可见，DG的好处在于可以从左（给定了边界条件的一侧）向右逐次求解各个单元。
分析单元矩阵
A
j
\bm{A}_j
Aj的性质。由于不同单元
I
j
I_j
Ij中的基函数
φ
j
m
(
x
)
\varphi_j^m(x)
φjm(x)映射到参考单元（reference element）后是一样的，
所以元素
a
m
l
a_{ml}
aml通过坐标参数变换（一维中的
ξ
=
2
x
−
x
j
Δ
x
j
\xi=2 \frac{x-x_j}{\Delta x_j}
ξ=2Δxjx−xj能将physical element中的坐标
x
∈
[
x
j
−
1
2
,
x
j
+
1
2
]
x\in[x_{j-\frac{1}{2}}, x_{j+\frac{1}{2}}]
x∈[xj−21,xj+21]转成reference element的坐标
ξ
∈
[
−
1
,
1
]
\xi\in[-1,1]
ξ∈[−1,1]）映射到参考单元后，它的数值是一样的。这意味着，
A
j
\bm{A}_j
Aj与实际的单元
I
j
I_j
Ij的具体性质无关（independent of
j
j
j），整体只要计算一次
A
j
\bm{A}_j
Aj即可。下面简记为
A
\bm{A}
A。
解的存在唯一性
A u j = b j \bm{A}\bm{u}_j=\bm{b}_j Auj=bj的解存在唯一等价于 A u j = 0 \bm{A}\bm{u}_j=\bm{0} Auj=0有且仅有零解。
这对应到原方程，相当于已知 f = 0 f=0 f=0和边界条件 a = 0 a=0 a=0时，能否推出解 u = 0 u=0 u=0。
-
先证两端点为0。由测试函数的任意性，取 v ∣ I 1 = u h ∣ I 1 ∈ P k ( I 1 ) v|_{I_1}=u_h|_{I_1}\in\bm{P}^k(I_1) v∣I1=uh∣I1∈Pk(I1)。代入式（1）即为：
− ∫ I j u h ( u h ) x d x + ( u h ) x 3 2 − ( u h ) x 3 2 − = 0 -\int_{I_j}u_h(u_h)_x {\rm d}x + (u_h)_{x_\frac{3}{2}}^{-} (u_h)_{x_\frac{3}{2}}^{-}=0 −∫Ijuh(uh)xdx+(uh)x23−(uh)x23−=0
推出 ( u h 2 ) 2 ∣ x 1 2 x 3 2 + [ ( u h ) x 3 2 − ] 2 = 0 \left( \frac{u_h}{2} \right)^2|_{x_\frac{1}{2}}^{x_\frac{3}{2}}+\left[(u_h)_{x_\frac{3}{2}}^{-}\right]^2=0 (2uh)2∣x21x23+[(uh)x23−]2=0，显然有 ( u h ) x 3 2 − = ( u h ) x 1 2 + = a = 0 (u_h)_{x_\frac{3}{2}}^{-}=(u_h)_{x_\frac{1}{2}}^{+}=a=0 (uh)x23−=(uh)x21+=a=0。 -
再证在整个单元内为0。由上一步知道有 ∫ I 1 u h v x d x = 0 \int_{I_1}u_hv_x {\rm d}x =0 ∫I1uhvxdx=0。而由条件知道近似解 u h u_h uh多项式有 ( x − x 1 2 ) (x-x_\frac{1}{2}) (x−x21)的因子，可写为 u h ( x ) = ( x − x 1 2 ) u ^ ( x ) u_h(x)=(x-x_\frac{1}{2})\hat{u}(x) uh(x)=(x−x21)u^(x)，其中 u ^ ( x ) ∈ P k − 1 ( I 1 ) \hat{u}(x)\in \bm{P}^{k-1}(I_1) u^(x)∈Pk−1(I1)。由测试函数的任意性，取 v ( x ) = ∫ x 1 2 x u ^ ( ξ ) d ξ v(x)=\int_{x_\frac{1}{2}}^x\hat{u}(\xi) {\rm d} \xi v(x)=∫x21xu^(ξ)dξ。则因为 ( x − x 1 2 ) (x-x_\frac{1}{2}) (x−x21)恒大于0，积分 ∫ I 1 ( x − x 1 2 ) u ^ ( x ) v x d x = 0 \int_{I_1}(x-x_\frac{1}{2})\hat{u}(x)v_x {\rm d}x=0 ∫I1(x−x21)u^(x)vxdx=0当且仅当 u ^ ( x ) = 0 \hat{u}(x)=0 u^(x)=0，故整个单元 I 1 I_1 I1均为0.
-
以此类推逐个单元。
误差量度
L
2
L_2
L2范数误差：
L
2
:
=
∣
∣
u
−
u
h
∣
∣
L
2
=
(
∫
0
1
(
u
(
x
)
−
u
h
(
x
)
)
2
d
x
)
1
2
L_2:=||u-u_h||_{L_2}=\left( \int_0^1 \left( u(x) - u_h(x)\right)^2 {\rm d}x \right)^\frac{1}{2}
L2:=∣∣u−uh∣∣L2=(∫01(u(x)−uh(x))2dx)21
如何说明一个格式的阶数 r r r呢？假定在 h h h宽度下计算得到 e h ≈ C h r e_h\approx Ch^r eh≈Chr，加密一倍得到 e 2 h ≈ C ( 2 h ) r e_{2h} \approx C(2h)^r e2h≈C(2h)r，则
e 2 h e h = 2 r ⇒ r = log 2 ( e 2 h e h ) \frac{e_{2h}}{e_h}=2^r \space \Rightarrow \space r = \log_2\left(\frac{e_{2h}}{e_h}\right) ehe2h=2r ⇒ r=log2(ehe2h)
标签：xj,frac,21,DG,uj,Ij,一维,bm,ODE 来源： https://blog.csdn.net/weixin_43614211/article/details/123168090