20190731个人笔记(川大)

官网课程视频地址

笔记主页

复习回顾

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

半离散DG方法

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

离散格式:求uhVhu_h \in V_h使得vVh\forall v\in V_h
((uh)t,v)j(f(uh),vx)j+Uh,j+12Vj+12Uh,j12Vj12+=0 ((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


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

f^j+12=f^(Uh,j+12,Uh,j+12+)=Uh,j+12+ \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}}

于是有

((uh)t,v)j(f(uh),vx)j+f^j+12Vj+12f^j12Vj12+=0 ((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

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

理论分析

模型问题微分方程的特性

局部熵条件与L2L^2稳定性

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

v=uhVhv= u_h \in V_h

要注意
ddtIj(uh)tuhdx=ddtIj12(uh)2dx \frac{d}{dt}\int_{I_j} (u_h)_t u_h dx = \frac{d}{dt}\int_{I_j} \frac{1}{2}(u_h)^2dx

ddtIjU(uh)dx(f(uh),uhx)j+f^j+12uhj+12f^j12uhj12+=0 \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(uh)x=f(uh)uhx g(u_h)_x = f(u_h){u_h}_x
第二项分部积分,再定义
F^j+12=g(uh)j+12+f^j+12uhj+12 \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}}
得到
ddtIj ⁣U(uh)dx+F^j+12F^j12+F^j12+g(uh)j12++f^j12uhj12+ ⁣= ⁣0 \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

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

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

于是有
ddtIj ⁣U(uh)dx+F^j+12F^j120 \frac{d}{dt}\int_{I_j} \! U(u_h) dx + \hat{F} _{j+\frac{1}{2}} - \hat{F} _{j-\frac{1}{2}} \leq 0
假设是周期边界条件,上式对jj求和得
12ddtI ⁣(uh)2dx0 \frac{1}{2}\frac{d}{dt}\int_{I} \! (u_h) ^2dx \leq 0

I ⁣(uh(x,t))2dxI ⁣(uh(x,0))2dx \int_{I} \! (u_h(x,t)) ^2dx \leq\int_{I} \! (u_h(x, 0)) ^2dx

误差分析

目标:
uuhChk+1||u-u_h||\leq Ch^{k+1}

这里只给思路

eh=uuhe_h = u -u_h

IjutvdxIjf(u)vxdx+f(uj+12)vj+12f(uj12+)vj12+=0 \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)e(,0)not hold || e(,t) || \leq ||e(,0)||, \text{not hold}
不成立的根本原因是ehe_h不是多项式空间,因此需要找个多项式空间里ehe_h足够逼近的测试函数,于是做分解
eh=(uPu)(uhPu):=ηξ e_h = (u-P_u) - (u_h-Pu):=\eta - \xi
其中PuVhPu \in V_huuVhV_h上的投影

投影要满足
{(uPu)j+12=0(η,v)j=0,vPk1\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的分析
PuuChk+1. || Pu - u || \leq Ch^{k+1}.

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

然后把全部是ξ\xi的项放左边,与η\eta相关的项放右边

编程实现

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

python代码

笔记主页

龙格库塔方法

模型问题

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

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

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

高阶

u(1)=un+ΔtL(un,tn)u(2)=34un+14u(1)+14ΔtL(u(1),tn+Δt)un+1=13un+23u(2)+23ΔtL(u(2),tn+12Δt) \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

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

下面给出k=1k=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

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

作业2

证明DG数值求解如下方程的稳定性
ut+(a(x)u)x=0 u_t + (a(x)u)_x = 0
通量
f^j+12={a(xj+12)uhj+12,a(xj+12)0a(xj+12)uhj+12+,else \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}