DG方法2
笔记主页
⎩⎨⎧ut+f(u)x=0u(x,0)=u0(x)
半离散DG方法
Vh={v:V∣Ij∈Pk(Ij)}
离散格式:求uh∈Vh使得∀v∈Vh
((uh)t,v)j−(f(uh),vx)j+Uh,j+21−Vj+21−−Uh,j−21−Vj−21+=0
fj+21数值通量
fj+21=f(Uh,j+21−,Uh,j+21+)=Uh,j+21+
于是有
((uh)t,v)j−(f(uh),vx)j+fj+21Vj+21−−fj−21Vj−21+=0
向前欧拉(uh)t=Δtuhn+1−uhn
U(u)=21u2
取v=uh∈Vh
要注意
dtd∫Ij(uh)tuhdx=dtd∫Ij21(uh)2dx
有
dtd∫IjU(uh)dx−(f(uh),uhx)j+fj+21uhj+21−−fj−21uhj−21+=0
定义
g(uh)x=f(uh)uhx
第二项分部积分,再定义
F^j+21=−g(uh)j+21−+fj+21uhj+21−
得到
dtd∫IjU(uh)dx+F^j+21−F^j−21+F^j−21+g(uh)j−21++fj−21uhj−21+=0
定义上式最后三项θ=F^+g+−f^+,加一项减一项g−泰勒展开可得
θ=(g′(ξ)−f^)(uh+−uh−)
(f′(ξ)−f)(u+−u−)=(f^(ξ,ξ)−f^(u−,u+))(u+−u−)
如果f^满足三性质,则有θ>0
于是有
dtd∫IjU(uh)dx+F^j+21−F^j−21≤0
假设是周期边界条件,上式对j求和得
21dtd∫I(uh)2dx≤0
即
∫I(uh(x,t))2dx≤∫I(uh(x,0))2dx
目标:
∣∣u−uh∣∣≤Chk+1
这里只给思路
记eh=u−uh
∫Ijutvdx−∫Ijf(u)vxdx+f(uj+21−)vj+21−−f(uj−21+)vj−21+=0
只考虑迎风格式,根据前面证明稳定性,很显然,下面的式子不成立
∣∣e(,t)∣∣≤∣∣e(,0)∣∣,not hold
不成立的根本原因是eh不是多项式空间,因此需要找个多项式空间里eh足够逼近的测试函数,于是做分解
eh=(u−Pu)−(uh−Pu):=η−ξ
其中Pu∈Vh 是u在Vh上的投影
投影要满足
{(u−Pu)j+21−=0(η,v)j=0,∀v∈Pk−1
根据1975 P.Ciarlet的分析
∣∣Pu−u∣∣≤Chk+1.
于是取测试函数v=ξ∈Vh
然后把全部是ξ的项放左边,与η相关的项放右边
uh(x,t)=l=0∑kajl(t)ϕjl(x)
代入方程,然后分别去v=ϕjm(x),后面过程从略
python代码
笔记主页
模型问题
ut=L(u,t)
向前欧拉法
un+1=un+ΔtL(un)
对k>0的情形不太好用
二阶龙格库塔方法,
un+1=21un+21(u(1)+ΔtL(u(1)))
对k=1比较好,k>1不好用,
高阶
u(1)u(2)un+1=un+ΔtL(un,tn)=43un+41u(1)+41ΔtL(u(1),tn+Δt)=31un+32u(2)+32ΔtL(u(2),tn+21Δt)
RKDG求解,k=0,1,2
ut+ux=0,x∈[0,1],t>0
满足周期边界条件,以及
(1)
u(x,0)=sin(2πx)
(2)
u(x,0)=⎩⎨⎧1,41≤x≤430,else
下面给出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
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)
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
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
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
得到的可交互误差图(鼠标悬停在点上可查看误差)如下
证明DG数值求解如下方程的稳定性
ut+(a(x)u)x=0
通量
f^j+21=⎩⎨⎧a(xj+21)uhj+21−,a(xj+21)≥0a(xj+21)uhj+21+,else