first DG方法：二维ODE的边值问题

考虑二维ODE的边值问题：
$$\{\begin{array}{c}a{u}_x+b{u}_y=f(x,y)\\ u(x,0)=g_1(x),0\le x\le 1\\ u(0,y)=g_2(y),0\le y\le 1\end{array}$$

格式导出

类似一维ODE的边值问题求法，做网格划分（mesh decomposition），其中$x$和$y$方向分别有$N_x$和$N_y$个单元：

$$\begin{gathered} \begin{array}{c}0=x_{\frac{1}{2}}<x_{\frac{3}{2}}<...<x_{N_x+\frac{1}{2}}=1\\ 0=y_{\frac{1}{2}}<y_{\frac{3}{2}}<...<y_{N_x+\frac{1}{2}}=1\end{array} \\ {I}_{ij}=\left[x_{j-\frac{1}{2}},x_{j+\frac{1}{2}}\right]\times \left[y_{j-\frac{1}{2}},y_{j+\frac{1}{2}}\right] \\ h:=\max \left\{\underset{i}{\max }\Delta x_i,\underset{j}{\max }\Delta y_j\right\} \end{gathered}$$

解空间定义为：
$${\mathbf{V}}_h^k=\left\{v:v{|}_{I_{ij}}\in {\mathbf{P}}^k(I_{ij}),1\le i\le N_x,1\le j\le N_y\right\}$$

测试函数则是${\mathbf{V}}_h^k$中的任意函数：$\forall v\in {\mathbf{V}}_h^k$.

记近似解为$u_h$，则原问题的弱形式为：

$${\int }_{I_{ij}}\left[a(u_h{)}_{x}v+b(u_h{)}_{y}v\right]\mathrm{d}x\mathrm{d}y={\int }_{I_{ij}}fv\mathrm{d}x\mathrm{d}y\text{\hspace{1em}(1)}$$

类似一维ODE边值问题时的处理，单独在每个方向上用分部积分：

$${\int }_{I_{ij}}(u_h{)}_{x}v\mathrm{d}x\mathrm{d}y={\int }_{y_{j-\frac{1}{2}}}^{y_{j+\frac{1}{2}}}\left[{\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}(u_h{)}_{x}v\mathrm{d}x\right]\mathrm{d}y$$

其中

$${\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}(u_h{)}_{x}v\mathrm{d}x=-{\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}{u}_h{v}_x\mathrm{d}x+(uv){|}_{x_{i+\frac{1}{2}}}-(uv){|}_{x_{i-\frac{1}{2}}}$$

代回原式有

$${\int }_{I_{ij}}(u_h{)}_{x}v\mathrm{d}x\mathrm{d}y=-{\int }_{y_{j-\frac{1}{2}}}^{y_{j+\frac{1}{2}}}{\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}{u}_h{v}_x\mathrm{d}x\mathrm{d}y+{\int }_{y_{j-\frac{1}{2}}}^{y_{j+\frac{1}{2}}}{u}_h(x_{i+\frac{1}{2}},y)v(x_{i+\frac{1}{2}},y)\mathrm{d}y-{\int }_{y_{j-\frac{1}{2}}}^{y_{j+\frac{1}{2}}}{u}_h(x_{i-\frac{1}{2}},y)v(x_{i-\frac{1}{2}},y)\mathrm{d}y$$

对含$(u_h{)}_y$项也如此处理，式（1）的左端项得到

$${\int }_{I_{ij}}\left[a(u_h{)}_{x}v+b(u_h{)}_{y}v\right]\mathrm{d}x\mathrm{d}y=-{\int }_{I_{ij}}(a{u}_h{v}_x+b{u}_h{v}_y)\mathrm{d}x\mathrm{d}y+{\int }_{y_{j-\frac{1}{2}}}^{y_{j+\frac{1}{2}}}{u}_h(x_{i+\frac{1}{2}},y)v(x_{i+\frac{1}{2}},y)\mathrm{d}y-{\int }_{y_{j-\frac{1}{2}}}^{y_{j+\frac{1}{2}}}{u}_h(x_{i-\frac{1}{2}},y)v(x_{i-\frac{1}{2}},y)\mathrm{d}y+{\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}{u}_h(x,y_{j+\frac{1}{2}})v(x,y_{j+\frac{1}{2}})\mathrm{d}y-{\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}{u}_h(x,y_{j-\frac{1}{2}})v(x,y_{j-\frac{1}{2}})\mathrm{d}y$$

可见单元$I_{ij}$的四个边界信息都会被用到。

对于Galerkin方法，测试函数$v\in {\mathbf{V}}_h^k$。

当$k=0$时，期望DG退化到一阶迎风FD方法，此时$u_h{|}_{I_{ij}}=u_{ij}$
$$a\frac{u_{ij}-u_{i-1,j}}{\Delta x}+b\frac{u_{ij}-u_{i,j-1}}{\Delta y}=f_{ij}$$
即
$$a(u_{ij}-u_{i-1,j})\Delta y+b(u_{ij}-u_{i,j-1})\Delta x=f_{ij}\Delta x\Delta y$$

由于测试函数是任意选取的${\mathbf{V}}_h^0$函数，所以将其取为$v=\{\begin{array}{c}1,I_{ij}\\ 0,\mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\end{array}$，则$v_x=v_y=0$。于是l类似于一维DG的推导，可以确定单元边界上的两侧值的取法：
$$-{\int }_{y_{j-\frac{1}{2}}}^{y_{j+\frac{1}{2}}}a{u}_h(x_{i-\frac{1}{2}}^{-},y)v(x_{i-\frac{1}{2}}^{+},y)\mathrm{d}y+{\int }_{y_{j-\frac{1}{2}}}^{y_{j+\frac{1}{2}}}a{u}_h(x_{i+\frac{1}{2}}^{-},y)v(x_{i+\frac{1}{2}}^{-},y)\mathrm{d}y-{\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}b{u}_h(x,y_{j-\frac{1}{2}}^{-})v(x,y_{j-\frac{1}{2}}^{+})\mathrm{d}x+{\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}b{u}_h(x,y_{j+\frac{1}{2}}^{-})v(x,y_{j+\frac{1}{2}}^{-})\mathrm{d}x$$

可以归结为，测试函数的取值永远朝向本单元的方向取，这样才能有值（=1）。而$u_h$的取向永远朝向迎风的方向取。

所以完整的DG格式为：

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

$$-{\int }_{I_{ij}}(a{u}_h{v}_x+b{u}_h{v}_y)\mathrm{d}x\mathrm{d}y+a{\int }_{y_{j-\frac{1}{2}}}^{y_{j+\frac{1}{2}}}\left[u_h(x_{i+\frac{1}{2}}^{-},y)v(x_{i+\frac{1}{2}}^{-},y)-u_h(x_{i-\frac{1}{2}}^{-},y)v(x_{i-\frac{1}{2}}^{+},y)\right]\mathrm{d}y+b{\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}\left[u_h(x,y_{j+\frac{1}{2}}^{-})v(x,y_{j+\frac{1}{2}}^{-})-u_h(x,y_{j-\frac{1}{2}}^{-})v(x,y_{j-\frac{1}{2}}^{+})\right]\mathrm{d}x={\int }_{I_{ij}}fv\mathrm{d}x\mathrm{d}y$$

求解方式

在单元$I_{ij}$上取基函数为${\mathbf{P}}^k(I_{ij})$中的函数做local bases，记为${\phi }_{ij}^l(x.y),l=1,...,K$，其中$K=\frac{(k+1)(k+2)}{2}$，共有$K$个基函数。例如，当$k=0$时$K=1$，取${\phi }_{ij}^1(x.y)=1$；当$k=1$时$K=3$，取${\phi }_{ij}^1(x.y)=1,{\phi }_{ij}^2(x.y)=x,{\phi }_{ij}^3(x.y)=y$。当然选取不是唯一的，也可以取Lagrange插值基函数。

假定$u_h{|}_{I_{ij}}={\sum }_{l=1}^K{u}_{ij}^l{\phi }_{ij}^l(x,y)$，需要求解${\mathbf{u}}_{ij}=[u_{ij}^1,u_{ij}^2,...,u_{ij}^K{]}^T$。测试函数是任意的，因此取$v={\phi }_{ij}^m(x,y),m=1,2,...,K$。所以完整的DG格式可以写成

$$-{\int }_{I_{ij}}\left[a\left(\sum _{l=1}^K{u}_{ij}^l{\phi }_{ij}^l(x,y)\right){\left({\phi }_{ij}^m\right)}_x+b\left(\sum _{l=1}^K{u}_{ij}^l{\phi }_{ij}^l(x,y)\right){\left({\phi }_{ij}^m\right)}_y\right]\mathrm{d}x\mathrm{d}y+a{\int }_{y_{j-\frac{1}{2}}}^{y_{j+\frac{1}{2}}}\left[\left(\sum _{l=1}^K{u}_{ij}^l{\phi }_{ij}^l(x_{i+\frac{1}{2}},y)\right){\phi }_{ij}^m(x_{i+\frac{1}{2}},y)-u_h(x_{i-\frac{1}{2}}^{-},y){\phi }_{ij}^m(x_{i-\frac{1}{2}},y)\right]\mathrm{d}y+b{\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}\left[\left(\sum _{l=1}^K{u}_{ij}^l{\phi }_{ij}^l(x,y_{j+\frac{1}{2}})\right){\phi }_{ij}^m(x,y_{j+\frac{1}{2}})-u_h(x,y_{j-\frac{1}{2}}^{-}){\phi }_{ij}^m(x,y_{j-\frac{1}{2}})\right]\mathrm{d}x={\int }_{I_{ij}}f{\phi }_{ij}^m(x,y)\mathrm{d}x\mathrm{d}y$$

整理出来未知的待求解的$u_{ij}^l$值有

$$\sum _{l=1}^K{u}_{ij}^l\left[-{\int }_{I_{ij}}a{\phi }_{ij}^l{\left({\phi }_{ij}^m\right)}_x-{\int }_{I_{ij}}b{\phi }_{ij}^l{\left({\phi }_{ij}^m\right)}_y+a{\int }_{y_{j-\frac{1}{2}}}^{y_{j+\frac{1}{2}}}{\phi }_{ij}^l(x_{i+\frac{1}{2}},y){\phi }_{ij}^m(x_{i+\frac{1}{2}},y)+b{\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}{\phi }_{ij}^l(x,y_{j+\frac{1}{2}}){\phi }_{ij}^m(x,y_{j+\frac{1}{2}})\right]={\int }_{I_{ij}}f{\phi }^m+a{\int }_{y_{j-\frac{1}{2}}}^{y_{j+\frac{1}{2}}}{u}_h(x_{i-\frac{1}{2}}^{-},y){\phi }_{ij}^m(x_{i-\frac{1}{2}},y)+b{\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}{u}_h(x,y_{j-\frac{1}{2}}^{-}){\phi }_{ij}^m(x,y_{j-\frac{1}{2}})$$

即${\sum }_{l=1}^K{a}_{ij}^{ml}{u}_{ij}^l=b_{ij}^m$，遍历$m=1,2,...,K$则组成一个线性方程组：${\mathbf{A}}_{ij}{\mathbf{u}}_{ij}={\mathbf{b}}_{ij}$，其中${\mathbf{A}}_{ij}=(a_{ij}^{ml}{)}_{K\times K}$，${\mathbf{b}}_{ij}={\left[b_{ij}^1,...,b_{ij}^K\right]}^T$。

所以只要知道了邻居单元$I_{i-1,j}$和$I_{i,j-1}$的值（它们贡献给右端向量），就能解出来$u_{ij}^l$，也就线性组合出$u_h{|}_{I_{ij}}={\sum }_{l=1}^K{u}_{ij}^l{\phi }_{ij}^l(x,y)$。计算的次序为：$I_{1,1}=>I_{1,2},I_{2,1}=>I_{1,3},I_{2,2},I_{3,1}=>...$。注意DG的stencil非常紧凑，无论单元内部多高阶，都只要邻居单元在边界上传递一个变量即可，而不需像有限差分那样需要传递冗长的stencil依赖。

DG的一些理论性质

优点：

1. 解的存在唯一性。证明类似一维。

2. $a{u}_x+b{u}_y=f(x,y)$对于任意的常量$a$或$b$均成立

3. 精度，有Lasaint & Ravart 在1974年证明stable accuracy。对于光滑的$u$，如果解空间选为${\mathbf{V}}_h^k=\left\{v:v{|}_{I_{ij}}\in {\mathbf{Q}}^k(I_{ij})\right\}$，有
   $$||u-u_h|{|}_{L_2}\le C{h}^{k+1}$$
   其中${\mathbf{Q}}^k(I_{ij})$内的函数为${\sum }_{m=0}^k{\sum }_{n=0}^k{a}_{mn}{x}^m{y}^n$，区别于${\mathbf{P}}^k(I_{ij})$内的函数为${\sum }_{m+n\le k}{a}_{mn}{x}^m{y}^n$。例如，当$k=1$时，${\mathbf{P}}^1$包含$1,x,y$，而${\mathbf{Q}}^1$包含$1,x,y,xy$。
   而如果采用${\mathbf{P}}^k$，则数值上可以观测到optimal order of accuracy，但没有证明。

4. 可推广到三角单元。此时对于光滑的$u$，如果解空间选为${\mathbf{V}}_h^k=\left\{v:v{|}_{I_{ij}}\in {\mathbf{P}}^k(I_{ij})\right\}$，有
   $$||u-u_h|{|}_{L_2}\le C{h}^{k+\frac{1}{2}}$$

5. $h$-自适应性。对任意的三角划分，可以局部加密，且允许hanging nodes（即一个三角形的一条边同时与另外多个三角形的边相挨）。

6. $p$-自适应性。不同单元可以使用不同阶数的多项式基函数。而这对于传统的连续FE是不允许的，因为FE要求单元边界时连续的，如${\mathbf{P}}^2$与${\mathbf{P}}^3$就无法在单元相挨边界处连续。

7. 方法可以推广到线性的初边值问题（非定常）。例如，将上述原问题的$y$空间项换成时间项$t$，整个PDE结构并没有改变。

局限性：

1. 难以处理非线性。例如$2u{u}_x+u_y=f$，意味着$u$的大小决定了风向，不知道$u$什么时候取什么符号。
2. 难以处理hyperbolic系统。例如$\{\begin{array}{c}{u}_x+v_y=0\\ v_x+u_y=0\end{array}$，变量替换$\xi =u+v$和$\eta =u-v$可以解耦得到$\{\begin{array}{c}{\xi }_x+{\xi }_y=0\\ {\eta }_x-{\eta }_y=0\end{array}$。于是两个scalar方程在$y$方向上的风向不一样，不知道如何初始化，因为有三个方向需要邻居的信息，所以不能one by one地去算。也就是说不能应对单元需要多于两个边界条件地情况。
3. time-dependent的初边问题。例如$\{\begin{array}{c}{u}_t+u_x=0\\ u(x,0)=g(x)\\ u(0,t)=u(1,t)\end{array}$的周期边界问题。

所以应对的办法是：RKDG（Runge-Kutta Discontinuous Galerkin）。