DG方法：一维ODE

• 有限差分方法
• Discrete Galerkin
• • 格式导出
  • 求解方式
  • 解的存在唯一性
  • 误差量度

考虑一维ODE的边值问题：
$$\{\begin{array}{c}{u}_x=f(x)\\ u(0)=a\end{array},x\in [0,1]$$

有限差分方法

对$x\in [0,1]$做uniform划分：$x_j=j\Delta x=j\frac{1}{N}$

如果是一阶FD：对于第$j$个网格点，需要利用$u_{j-1},u_j,u_{j+1}$等信息逼近$u_x$，其中$u_j$表示准确解$u(x_j)$的近似，则

$$u_x{|}_{x_j}\approx \frac{u_j-u_{j-1}}{\Delta x}$$

要选择用$x_{j-1}$是因为边界条件定义在左侧，信息从迎风的左边来。代入原式有

$$\frac{u_j-u_{j-1}}{\Delta x}=f(x_j)\Rightarrow u_j=u_{j-1}+\Delta xf(x_j)$$

所以求解过程可以从左向右扫描（sweep，最左边是$u_0=a$），非常易于计算。

但只有一阶精度：

$$|u(x_j)-u_j|\approx O(\Delta x)$$

当从$N$个网格点加密到$2N$时，工作量增加一倍，但误差仅能减少一半。

而如果使用更高阶精度的差分方法，同理，将$u_x(x_j)$表示为$u_{j-2},u_{j-1},u_j$的线性组合（可通过Taylor展开推导）。这对于$j=2,3,...,N$是可行的，但$j=1$处显然没有这么多stencil可供构造，所以在边界上使用一阶格式，但这又会导致计算降阶！因为$u_1$只有一阶精度，那么根据$u_1$计算出的$u_2$也达不到二阶，如此类推。

Discrete Galerkin

格式导出

DG希望既易于计算，又能有高阶精度。它首先定义一系列半网格点$x_{\frac{2k+1}{2}}(k=0,1,...,N)$，将区间划分为很多个单元，单元$I_j=[x_{j-\frac{1}{2}},x_{j+\frac{1}{2}}]$。注意仅要求$0=x_{\frac{1}{2}}<x_{\frac{3}{2}}<...<x_{N-\frac{1}{2}}<x_{N+\frac{1}{2}}=1$，而不需要每个单元都是均匀长度的。每个单元的中点$x_j=\frac{1}{2}(x_{j-\frac{1}{2}}+x_{j+\frac{1}{2}})$，长度$\Delta x_j=x_{j-\frac{1}{2}}-x_{j+\frac{1}{2}}$。记$h=\underset{j}{\max }\Delta x_j$。

定义解空间（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}

其中${\mathbf{P}}^k(I_j)$表示第$j$个单元上的次数小于等于$k$的多项式空间。当$k=1$时即为分片线性函数（piecewise linear functions）。注意在cell边界处，函数值时可以间断的，而区别于传统有限元要求内部单元的边界值连续。

由于解在单元边界处可以完全不连续，因此$x_{j+\frac{1}{2}}$处存在一个左值$u_{j+\frac{1}{2}}^{-}$和一个右值$u_{j+\frac{1}{2}}^{+}$，分别属于左右两个单元。

不同单元的多项式空间${\mathbf{P}}^k$的阶数$k$可以不同（$\mathbf{P}$ adaptivity）。

定义测试函数（test function）$v(x)$，原问题成立，则如下积分式一定成立：

$${\int }_{I_j}{u}_{x}v={\int }_{I_j}f(x)v,j=1,2,...,N$$

假定$v(x)$是至少可微的，那么分部积分：

$${\int }_{I_j}f(x)v={\int }_{I_j}{u}_{x}v=-{\int }_{I_j}u{v}_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}})$$

所以问题转换为：

find $u_h(x)\in {\mathbf{V}}_h^k$来逼近$u(x)$，代入上式有

$$-{\int }_{I_j}u{v}_x+(u_{h}v){|}_{x_{j+\frac{1}{2}}}-(u_{h}v){|}_{x_{j-\frac{1}{2}}}={\int }_{I_j}f(x)v$$

但问题是$v\in$？，如果任意的可微函数都可以，那么自由度太多会造成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$上确定${\mathbf{P}}^k(I_j)$需要$k+1$个参数，将${\mathbf{P}}^k(I_j)$这个线性空间的基取为${\phi }_j^0(x),{\phi }_j^1(x),...,{\phi }_j^k(x)$，则线性表出为

$$u_h(x){|}_{I_j}=\sum _{l=0}^k{u}_j^l{\phi }_j^l(x)$$

因为共有$N$个单元，则$\dim V_h^k=N(k+1)$。需要求解的值是线性组合的系数$u_j^l,l=0,1,...,k,j=1,2,...,N$，因此需要有$N(k+1)$个约束条件。假定$v\in {\mathbf{W}}_h$，那么应该要有$\dim {\mathbf{W}}_h=\dim {\mathbf{V}}_h^k=N(k+1)$。

这两者是可以取成不一样的线性空间的，而如果取成一样的（${\mathbf{W}}_h={\mathbf{V}}_h^k$，都来自分片线性多项式的空间），则成为Galerkin方法。

于是问题转换为：find $u_h\in {\mathbf{V}}_h^k,v\in {\mathbf{V}}_h^k$ ……

但仍不对，因为单元边界上有两个值（左值$u_{j+\frac{1}{2}}^{-}$和右值$u_{j+\frac{1}{2}}^{+}$）。需要额外加条件。而我们希望当$k=0$时，DG方法能还原回普通的一阶（迎风）有限差分方法：$\frac{u_j-u_{j-1}}{\Delta x}=f(x_j)$。

此时解空间和测试函数空间都为分片常值函数。由于测试函数是任意的，取$v=\{\begin{array}{c}1,I_j\\ 0,\mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\end{array}$，则$v_x=0$，${\int }_{I_j}{u}_h{v}_x=0$。于是有

$$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)$$

希望能还原出：$u_j-u_{j-1}=\Delta xf(x_j)$。而在$k=0$时，$u_j$就是第$j$个单元的函数近似解的值$u_h(I_j)$。

注意到$v(x)$的取值方式，因此取$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(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}$。

所以完整的DG格式为：

find $u_h(x)\in {\mathbf{V}}_h^k$，使得对于任意的测试函数$v\in {\mathbf{V}}_h^k$和任意的单元$I_j$，都有

$$-{\int }_{I_j}{u}_h(x)v_x\mathrm{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)\mathrm{d}x\text{\hspace{1em}(1)}$$

此时的格式与计算中如何选取${\mathbf{P}}^k$的基函数的具体形式无关！基函数可以有多种选择，例如Legendre多项式函数，Lagrange插值基函数，自然幂函数（$1,x,x^2,...,x^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$等等。

求解方式

一般地，假定选取的基函数为${\phi }_j^0(x),{\phi }_j^1(x),...,{\phi }_j^k(x)$，线性表出$u_h(x){|}_{I_j}={\sum }_{l=0}^k{u}_j^l{\phi }_j^l(x)$。需要求解向量
$${\mathbf{u}}_j=\left[\begin{array}{c}{u}_j^0\\ u_j^1\\ ...\\ u_j^k\end{array}\right]$$

由于测试函数是任意的，只要在线性空间${\mathbf{W}}_h={\mathbf{V}}_h^k$内就行。所以不妨取$v={\phi }_j^m(x),m=0,1,...,k$，则上一节中的格式写为

$$-{\int }_{I_j}\sum _{l=0}^k{u}_j^l{\phi }_j^l(x)({\phi }_j^m{)}_x\mathrm{d}x+\sum _{l=0}^k{u}_j^l{\phi }_j^l(x_{j+\frac{1}{2}}){\phi }_j^m(x_{j+\frac{1}{2}})-u_{j-\frac{1}{2}}^{-}{\phi }_j^m(x_{j-\frac{1}{2}})={\int }_{I_j}f(x){\phi }_j^m(x)\mathrm{d}x$$

上式要对于$m=0,1,...,k$均成立。注意其中的$u_{j-\frac{1}{2}}^{-}$是隔壁单元算出来的结果值，传给本单元的。上式提取出${\sum }_{l=0}^k{u}_j^l$，有

$$\sum _{l=0}^k{u}_j^l\left[-{\int }_{I_j}{\phi }_j^l(x)({\phi }_j^m{)}_x\mathrm{d}x+{\phi }_j^l(x_{j+\frac{1}{2}}){\phi }_j^m(x_{j+\frac{1}{2}})\right]=u_{j-\frac{1}{2}}^{-}{\phi }_j^m(x_{j-\frac{1}{2}})+{\int }_{I_j}f(x){\phi }_j^m(x)\mathrm{d}x\text{\hspace{1em}(2)}$$

记：
$$a_{ml}=-{\int }_{I_j}{\phi }_j^l(x)({\phi }_j^m{)}_x\mathrm{d}x+{\phi }_j^l(x_{j+\frac{1}{2}}){\phi }_j^m(x_{j+\frac{1}{2}})\text{\hspace{1em}(2-1)}$$

$$b_m=u_{j-\frac{1}{2}}^{-}{\phi }_j^m(x_{j-\frac{1}{2}})+{\int }_{I_j}f(x){\phi }_j^m(x)\mathrm{d}x\text{\hspace{1em}(2-2)}$$

则在每个单元内都组成一个线性方程组：${\mathbf{A}}_j{\mathbf{u}}_j={\mathbf{b}}_j$，其中${\mathbf{A}}_j=(a_{ml})$是一个$(k+1)\times (k+1)$的矩阵。

一旦$u_{j-\frac{1}{2}}^{-}$给定，就可以求解该线性方程组。由此可见，DG的好处在于可以从左（给定了边界条件的一侧）向右逐次求解各个单元。

分析单元矩阵${\mathbf{A}}_j$的性质。由于不同单元$I_j$中的基函数${\phi }_j^m(x)$映射到参考单元（reference element）后是一样的，
所以元素$a_{ml}$通过坐标参数变换（一维中的$\xi =2\frac{x-x_j}{\Delta x_j}$能将physical element中的坐标$x\in [x_{j-\frac{1}{2}},x_{j+\frac{1}{2}}]$转成reference element的坐标$\xi \in [-1,1]$）映射到参考单元后，它的数值是一样的。这意味着，${\mathbf{A}}_j$与实际的单元$I_j$的具体性质无关（independent of $j$），整体只要计算一次${\mathbf{A}}_j$即可。下面简记为$\mathbf{A}$。

解的存在唯一性

$\mathbf{A}{\mathbf{u}}_j={\mathbf{b}}_j$的解存在唯一等价于$\mathbf{A}{\mathbf{u}}_j=0$有且仅有零解。

这对应到原方程，相当于已知$f=0$和边界条件$a=0$时，能否推出解$u=0$。

1. 先证两端点为0。由测试函数的任意性，取$v{|}_{I_1}=u_h{|}_{I_1}\in {\mathbf{P}}^k(I_1)$。代入式（1）即为：
   $$-{\int }_{I_j}{u}_h(u_h{)}_x\mathrm{d}x+(u_h{)}_{x_{\frac{3}{2}}}^{-}(u_h{)}_{x_{\frac{3}{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$，显然有$(u_h{)}_{x_{\frac{3}{2}}}^{-}=(u_h{)}_{x_{\frac{1}{2}}}^{+}=a=0$。

2. 再证在整个单元内为0。由上一步知道有${\int }_{I_1}{u}_h{v}_x\mathrm{d}x=0$。而由条件知道近似解$u_h$多项式有$(x-x_{\frac{1}{2}})$的因子，可写为$u_h(x)=(x-x_{\frac{1}{2}})\hat{u}(x)$，其中$\hat{u}(x)\in {\mathbf{P}}^{k-1}(I_1)$。由测试函数的任意性，取$v(x)={\int }_{x_{\frac{1}{2}}}^x\hat{u}(\xi )\mathrm{d}\xi$。则因为$(x-x_{\frac{1}{2}})$恒大于0，积分${\int }_{I_1}(x-x_{\frac{1}{2}})\hat{u}(x)v_x\mathrm{d}x=0$当且仅当$\hat{u}(x)=0$，故整个单元$I_1$均为0.

3. 以此类推逐个单元。

误差量度

$L_2$范数误差：
$$L_2:=||u-u_h|{|}_{L_2}={\left({\int }_0^1{\left(u(x)-u_h(x)\right)}^2\mathrm{d}x\right)}^{\frac{1}{2}}$$

如何说明一个格式的阶数$r$呢？假定在$h$宽度下计算得到$e_h\approx C{h}^r$，加密一倍得到$e_{2h}\approx C(2h{)}^r$，则

$$\frac{e_{2h}}{e_h}=2^r\Rightarrow r={\log }_2\left(\frac{e_{2h}}{e_h}\right)$$