DG方法2

1. 个人笔记(川大)

笔记主页

1.2. 复习回顾

$$\begin{cases} u_t + f(u)_x = 0\\[1em] u(x,0) = u^0(x) \end{cases}$$

半离散DG方法

$V_h = \{ v:V|_{I_j} \in P^k(I_j)\}$

离散格式：求$u_h \in V_h$使得$\forall v\in V_h$
$$((u_h)_t, v)_j - (f(u_h), v_x)_j + U^-_{h,j+\frac{1}{2}}V^-_{j+\frac{1}{2}}- U^-_{h,j-\frac{1}{2}}V^+_{j-\frac{1}{2}}=0$$

---

$\widehat{f}_{j+\frac{1}{2}}$数值通量

$$\widehat{f}_{j+\frac{1}{2}} = \widehat f(U^-_{h,j+\frac{1}{2}}, U^+_{h,j+\frac{1}{2}}) = U^+_{h,j+\frac{1}{2}}$$

于是有

$$((u_h)_t, v)_j - (f(u_h), v_x)_j + \widehat{f}_{j+\frac{1}{2}} V^-_{j+\frac{1}{2}}- \widehat{f}_{j-\frac{1}{2}} V^+_{j-\frac{1}{2}}=0$$

向前欧拉$(u_h)_t = \frac{u_h^{n+1} - u_h^n}{\Delta t}$

1.3. 理论分析

1.3.1. 模型问题微分方程的特性

• 存在间断解，即使初值$u^0$是光滑的

• 不可能定义强解，因此定义弱解，弱解存在，但是没有唯一性

• 加一个熵条件(Entripy Solution)可保证弱解(熵解)唯一性

  • 熵条件有很多种定义办法，例如：

  • $U(u)_t + F(u)_x \leq 0$, $U''\geq 0$

  • 分析

  • $U'(u)(u_t + f(u)_x) = 0$

  • $U(u)_t + U'(u)f'(u)u_x=0$

  • $U(u)_t +F(u)_x = 0$

1.3.2. 局部熵条件与$L^2$稳定性

$U(u) = \frac{1}{2}u^2$

取$v= u_h \in V_h$

要注意
$$\frac{d}{dt}\int_{I_j} (u_h)_t u_h dx = \frac{d}{dt}\int_{I_j} \frac{1}{2}(u_h)^2dx$$
有
$$\frac{d}{dt}\int_{I_j} U(u_h) dx - (f(u_h), {u_h}_x)_j + \widehat{f}_{j+\frac{1}{2}} {u_h}^-_{j+\frac{1}{2}}- \widehat{f}_{j-\frac{1}{2}} {u_h}^+_{j-\frac{1}{2}}=0$$
定义
$$g(u_h)_x = f(u_h){u_h}_x$$
第二项分部积分，再定义
$$\hat{F}_{j+\frac{1}{2}} =- g(u_h)^-_{j+\frac{1}{2}}+ \widehat{f}_{j+\frac{1}{2}} {u_h}^-_{j+\frac{1}{2}}$$
得到
$$\frac{d}{dt}\int_{I_j} \! U(u_h) dx + \hat{F} _{j+\frac{1}{2}} - \hat{F} _{j-\frac{1}{2}} + \hat{F} _{j-\frac{1}{2}} + g(u_h)^+_{j-\frac{1}{2}}+ \widehat{f}_{j-\frac{1}{2}} {u_h}^+_{j-\frac{1}{2}}\!=\! 0$$

定义上式最后三项$\theta = \hat{F} + g^+ - \hat f^+$，加一项减一项$g^-$泰勒展开可得
$$\theta = (g'(\xi) - \hat f)(u_h^+ - u_h^-)$$

$$(f'(\xi) - f)( u^+ - u^-) = \Big(\hat f(\xi, \xi) - \hat f(u^-, u^+) \Big)(u^+ - u^-)$$
如果$\hat f$满足三性质，则有$\theta >0$

于是有
$$\frac{d}{dt}\int_{I_j} \! U(u_h) dx + \hat{F} _{j+\frac{1}{2}} - \hat{F} _{j-\frac{1}{2}} \leq 0$$
假设是周期边界条件，上式对$j$求和得
$$\frac{1}{2}\frac{d}{dt}\int_{I} \! (u_h) ^2dx \leq 0$$
即
$$\int_{I} \! (u_h(x,t)) ^2dx \leq\int_{I} \! (u_h(x, 0)) ^2dx$$

1.4. 误差分析

目标：
$$||u-u_h||\leq Ch^{k+1}$$

这里只给思路

记$e_h = u -u_h$

$$\int_{I_j}u_t v dx - \int_{I_j}f(u)v_x dx + f(u^-_{j+\frac{1}{2}}) v^-_{j+\frac{1}{2}}- f(u^+_{j-\frac{1}{2}})v^+_{j-\frac{1}{2}}=0$$
只考虑迎风格式，根据前面证明稳定性，很显然，下面的式子不成立
$$|| e(,t) || \leq ||e(,0)||， \text{not hold}$$
不成立的根本原因是$e_h$不是多项式空间，因此需要找个多项式空间里$e_h$足够逼近的测试函数，于是做分解
$$e_h = (u-P_u) - (u_h-Pu):=\eta - \xi$$
其中$Pu \in V_h$ 是$u$在$V_h$上的投影

投影要满足
$$\begin{cases} (u-Pu)^-_{j+\frac{1}{2}}=0\\ (\eta, v)_j = 0, \quad\forall v\in P^{k-1} \end{cases}$$
根据1975 P.Ciarlet的分析
$$|| Pu - u || \leq Ch^{k+1}.$$

于是取测试函数$v = \xi \in V_h$

然后把全部是$\xi$的项放左边，与$\eta$相关的项放右边

1.5. 编程实现

$$u_h(x, t) = \sum_{l=0}^ka_j^l(t)\phi_j^l(x)$$
代入方程，然后分别去$v = \phi_j^m(x)$，后面过程从略

python代码

笔记主页

1.6. 龙格库塔方法

模型问题

$$u_t = L(u,t)$$

向前欧拉法
$$u^{n+1} = u^n + \Delta t L(u^n)$$
对$k>0$的情形不太好用

二阶龙格库塔方法,
$$u^{n+1} = \frac{1}{2}u^n + \frac{1}{2}(u^{(1)} + \Delta tL(u^{(1)}))$$
对$k=1$比较好，$k>1$不好用,

高阶

$$\begin{aligned} u^{(1)} &=u^{n}+\Delta t L\left(u^{n}, t^{n}\right) \\ u^{(2)} &=\frac{3}{4} u^{n}+\frac{1}{4} u^{(1)}+\frac{1}{4} \Delta t L\left(u^{(1)}, t^{n}+\Delta t\right) \\ u^{n+1} &=\frac{1}{3} u^{n}+\frac{2}{3} u^{(2)}+\frac{2}{3} \Delta t L\left(u^{(2)}, t^{n}+\frac{1}{2} \Delta t\right) \end{aligned}$$

1.7. 作业

1.7.1. 作业1

RKDG求解，$k=0,1,2$
$$u_t + u_x = 0， x\in[0,1], t>0$$
满足周期边界条件，以及
(1)
$$u(x,0)=\sin(2\pi x)$$
(2)
$$u(x,0)=\begin{cases} 1, \quad\frac{1}{4} \leq x \leq \frac{3}{4}\\[1em] 0, \quad \text{else} \end{cases}$$

下面给出$k=1$情形的python代码

import numpy as np
nElements = np.array([10,20,40,80,160,320])
nElement = nElements[0]
h = 1.0/nElement
dt = h*h
T = 2.3

k = 1
ip = np.array([-1, 1])

eledof = 2 # since deg = 1 
A = np.array([[0.5, 0.5], [-0.5, 0.5]])
M = 0.5*h*np.array([[2.0,1.0],[1.0, 2.0]])/3.0
C = np.array([[0, -1], [0, 0]]).T[-1]
MI = np.mat(M).I
MIA = MI.dot(-1.0*A)
MIC = MI.dot(-1.0*C)

## matrxi apply vector
def mathbbLdot(uhn):
    uhnp1 = np.zeros(nElement*eledof)
    for ic in range(nElement):
        tmp = MIA.dot(uhn[ic*eledof:(ic+1)*eledof])
        tmp += MIC*uhn[ic*eledof-1]
        uhnp1[ic*eledof:(ic+1)*eledof] = tmp
    return uhnp1

## initialize
ux0 = lambda x: np.sin(2*np.pi*np.array(x))
uhn = np.zeros(nElement*eledof)
for ic in range(nElement):
    ux0ele = ux0([ic*h, (ic+1)*h])
    uhn[ic*eledof:(ic+1)*eledof] = ux0ele

## iteration
def ite(uhn):
    u1 =  uhn + dt*mathbbLdot(uhn)
    uhnp1 = 0.5*uhn + 0.5*(u1 + dt*mathbbLdot(u1))
    return uhnp1
step = 0
maxstep = int(T/dt)
while(step < maxstep):
    uhn = ite(uhn)
    step += 1

得到的可交互误差图(鼠标悬停在点上可查看误差)如下

1.8. 作业2

证明DG数值求解如下方程的稳定性
$$u_t + (a(x)u)_x = 0$$
通量
$$\hat f_{j+\frac{1}{2}}=\begin{cases} a(x_{j+\frac{1}{2}}){u_h}_{j+\frac{1}{2}}^-, \quad a(x_{j+\frac{1}{2}})\geq 0\\[1em] a(x_{j+\frac{1}{2}}){u_h}_{j+\frac{1}{2}}^+, \quad \text{else} \end{cases}$$