20190730个人笔记(川大)

官网课程视频地址

笔记主页

DG方法一维

以一维为例子(1973)

ux=f,x[0,1];u(0)=a.u_x = f,\quad x\in [0,1];\quad u(0) = a.

用差分格式

目标是求解U0,U1,,UNU_0, U_1, \cdots, U_N

迎风格式Uxx=xj=UjUj1hU_x|_{x=x_j} = \frac{U_{j} - U_{j-1}}{h}

UjUj1h=f(xj),U0=a, \frac{U_{j} - U_{{j-1}}}{h} = f(x_j),\quad U_0 = a,

U1=U0+hf(x1),U2=U1+hf(x2),U_1 = U_0 +hf(x_1), U_2 = U_1 + hf(x_2),\cdots

缺点是只有一阶精度,

二阶差分格式:

Uj+1Uj12h=f(xj),U0=a, \frac{U_{j+1} - U_{{j-1}}}{2h} = f(x_j),\quad U_0 = a,

但是U1U_{-1}未知,并且比较依赖等距网格

DG方法

逼近空间

Ij=[xj12,xj+12]I_j = [x_{j-\frac{1}{2}}, x_{j+\frac{1}{2}}]

uhVh,Vh={vvIjPk(Ij)}u_h \in V_h, V_h = \{ v| v|_{I_j} \in \mathcal P^k(I_j)\}

dim(Vh)=(k+1)N\dim(V_h) = (k+1)N

乘以测试函数vv分部积分,得到弱格式

(u,vx)j+Uj+12Vj+12Uj12Vj12=(f,v)j(*) -(u, v_x)_j + U_{j+\frac{1}{2}}V_{j+\frac{1}{2}} -U_{j-\frac{1}{2}}V_{j-\frac{1}{2}} = (f, v)_j \tag{*}

离散格式:求uhVhu_h\in V_h

(uh,vh)j+Uh,j+12Vh,j+12=(f,vh)j+Uh,j12Vh,j12 -(u_h, v'_h)_j + U_{h,j+\frac{1}{2}}V_{h,j+\frac{1}{2}} = (f, v_h)_j +U_{h,j-\frac{1}{2}}V_{h,j-\frac{1}{2}}

vhVh\forall v_h\in V_h, j=1,2,,Nj=1,2,\cdots,N成立。
VhIjV_h|_{I_j}的基函数为拉格朗日插值多项式函数,即
VhIj=span{ϕmj:m=1,2,,N}V_h|_{I_j} = \text{span}\{\phi_m^{j} :m=1,2,\cdots, N\}满足
ϕmj(xnj)={1,m=n,0,mn. \phi_m^j(x_n^j)= \begin{cases} 1, \quad m=n,\\ 0, \quad m\neq n. \end{cases}

于是在单元IjI_j上有
uh=m=1k+1Umjϕmj u_h = \sum_{m=1}^{k+1}U^j_{m}\phi_m^j

在单元IjI_j上,分别取测试函数ϕnj,n=1,2,,k+1\phi_n^j, n=1,2,\cdots,k+1于是有如下线性系统
AjUj+BjUj=Fj+Sj,j=1,2,N.(a) \mathbb A^j \bm U^j + \mathbb B^j \bm U^j = \bm F^j + \bm S^j,\quad j=1,2,\cdots N. \tag{a}
其中

Aj=[amn](k+1)×(k+1),amn=Ijϕnϕmdx. \mathbb A^j =[a_{mn}]_{(k+1)\times(k+1)}, \quad a_{mn} = \int_{I_j}\phi_n\phi'_m d x.

Bj=[bmn](k+1)×(k+1),bmn={1,m=n=k+10,else \mathbb B^j =[b_{mn}]_{(k+1)\times(k+1)}, \quad b_{mn} = \begin{cases} 1,\quad &m=n=k+1\\ 0,\quad &\text{else} \end{cases}

未知函数
Uj=[U1j,U2j,,Uk+1j]. \bm U^j = [U_1^j, U_2^j,\cdots,U_{k+1}^j].

很显然有Uh,j12=U1jU_{h,j-\frac{1}{2}}=U_{1}^j以及Uh,j+12=Uk+1jU_{h,j+\frac{1}{2}}=U_{k+1}^j.

右端项

Fj=[Fmj](k+1)×1,Fmj=Ijfϕmdx. \bm F^j = [F_m^j]_{(k+1)\times 1},\quad F_m^j = \int_{I_j} f \phi_m d x.

右端项

Sj=[smj](k+1)×1 \bm S^j = [s_m^j]_{(k+1)\times 1}
其中
s1j={Uk+1j1,j=2,3,,Na,j=1. s_1^j = \begin{cases} U_{k+1}^{j-1},\quad &j=2,3,\cdots,N\\ a,\quad &j=1.\\ \end{cases}
smj=0,m>1.s_m^j = 0, m>1.

注意: 线性系统(a)是逐个单元求解。

先分析k=0k=0情形,

理论分析

二维推广

Ω=[0,1]2\Omega=[0,1]^2,方程

aux+buy=f in Ω au_x +bu_y = f \quad \text{ in } \Omega

满足边界条件u(x,0)=g(x),u(0,y)=h(y)u(x, 0) = g(x), u(0, y)=h(y).

应用迎风原理,容易计算

缺点: 对于如下非线性及方程组例子不太好处理

特殊例子1

(u22)x=f (\frac{u^2}{2})_x = f

uux=f uu_x =f

特殊例子2

Aux=f \mathbb A \bm u_x = \bm f
其中
A=(0110) \mathbb A=\Big( \begin{matrix} 0&1\\ 1&0 \end{matrix}\Big)
两方程相加即
(u+v)x=f1+f2 (u+v)_x = f_1 +f_2

换一种思路

{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}

半离散或龙格库塔方法

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

离散格式:求uh(,t)Vhu_h(,t) \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}

(uhn+1uhnΔt,v)j(f(uh),vx)j+f^j+12Vj+12f^j12Vj12+=0 \Big( \frac{u_h^{n+1} - u_h^n}{\Delta t} ,v\Big)_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

k=0k=0情形

Uh(,t)=Uj(t),xIj U_h(,t) = U_j(t), x\in I_j


v={1,xIj0,xIj v =\begin{cases} 1, \quad x\in I_j\\ 0,\quad x \notin I_j \end{cases}

hUj(t)+f^j+12f^j12=0. hU'_j(t) + \widehat{f}_{j+\frac{1}{2}} - \widehat{f}_{j-\frac{1}{2}} = 0.

这种格式其实就是三点差分格式,也是一种单调格式,单调格式很简单,需要加三个条件

练习

练习1

DG方法求解

ux=f,x[0,1];u(0)=a.u_x = f,\quad x\in [0,1];\quad u(0) = a.

其中f(x)=cosx,a=0f(x)=\cos x, a=0 取多项式次数k=0,1,2k=0,1,2,步长
h=110,120,140,,1320h = \frac{1}{10},\frac{1}{20},\frac{1}{40},\cdots,\frac{1}{320},计算L1,LL^1,L^\infty误差及阶。

点我下载python代码(kk任意)

注: python做数值计算的优势:时代主打,免费简单,跨平台,安装包体积小且可优盘携带,第三方库多,文档多,开发快,容易混编…

上面的python代码支持任意k>0k>0,但是代码太过繁琐,为方便理解,下面的python代码只针对k=1,2k=1,2情形,并且所有的矩阵都手动计算,右端ff也用插值,因此可以不涉及到积分,问题归结于矩阵处理,给出不到30行的求解uhu_h的python代码:

import numpy as np
nElements = np.array([10,20,40,80,160,320])
nElement = nElements[2]
h = 1.0/nElement
k = 2
eledof = k + 1
if k == 1:
    ip = np.array([-1, 1])
    A = np.array([[0.5, 0.5], [-0.5, 0.5]])
    M = np.array([[2.0,1.0],[1.0, 2.0]])/3.0
if k == 2:
    ip = np.array([-1, 0, 1])
    A = np.array([[3,4,-1],[-4,0,4],[1,-4,3]])/6.0
    M = np.array([[4,2,-1], [2,16,2], [-1,2,4]])/15.0
f = lambda x: np.cos(np.array(x))
u = lambda x: np.sin(np.array(x))
uh = np.zeros(eledof*nElement)
ui = np.zeros(eledof*nElement)
vectorU = np.zeros(eledof)
vectorU[-1] = a = 0
for ic in range(nElement):
    L, R = ic*h, (ic+1)*h
    xele = 0.5*(L+R) + 0.5*(R-L)*ip
    vectorF = M.dot(f(xele)) * (0.5*h)
    vectorF[0] += vectorU[-1]
    vectorU = np.linalg.solve(A, vectorF)
    uh[ic*eledof:(ic+1)*eledof] = vectorU
    ui[ic*eledof:(ic+1)*eledof] = u(xele)

在上面的代码后面再加个求误差的函数

def error(uh, u):
    from scipy.interpolate import interp1d
    pts, wts = np.polynomial.legendre.leggauss(20)
    L1sum = L2sum = Lmax = 0.0
    kind = ["slinear", "quadratic"]
    for ic in range(nElement):
        L, R = ic*h, (ic+1)*h
        xele = 0.5*(L+R) + 0.5*(R-L)*pts
        uhele = uh[ic*eledof: (ic+1)*eledof]
        iuhele = interp1d(ip, uhele, kind[k-1])
        tmp = u(xele) - iuhele(pts)
        tmp = abs(tmp)
        L1sum += (tmp*wts).sum()*h*0.5
        L2sum += (tmp*tmp*wts).sum()*h*0.5
        if(Lmax <= tmp.max()):
            Lmax = tmp.max()
    return L1sum, np.sqrt(L2sum), Lmax
    
print("h, k = ", h, k)
print("ui L1, L2, Loo error")
print(error(ui, u))
print("uh L1, L2, Loo error")
print(error(uh, u))

可以看到最终数值结果有限元误差和插值误差,在三种范数下L1,L2,LL^1,L^2,L^\infty的阶都是k+1k+1

练习2

证明:


笔记主页