DG方法初步
笔记主页
以一维为例子(1973)
u x = f , x ∈ [ 0 , 1 ] ; u ( 0 ) = a . u_x = f,\quad x\in [0,1];\quad u(0) = a. u x = f , x ∈ [ 0 , 1 ] ; u ( 0 ) = a .
目标是求解U 0 , U 1 , ⋯ , U N U_0, U_1, \cdots, U_N U 0 , U 1 , ⋯ , U N
迎风格式U x ∣ x = x j = U j − U j − 1 h U_x|_{x=x_j} = \frac{U_{j} - U_{j-1}}{h} U x ∣ x = x j = h U j − U j − 1
U j − U j − 1 h = f ( x j ) , U 0 = a , \frac{U_{j} - U_{{j-1}}}{h} = f(x_j),\quad U_0 = a, h U j − U j − 1 = f ( x j ) , U 0 = a ,
U 1 = U 0 + h f ( x 1 ) , U 2 = U 1 + h f ( x 2 ) , ⋯ U_1 = U_0 +hf(x_1), U_2 = U_1 + hf(x_2),\cdots U 1 = U 0 + h f ( x 1 ) , U 2 = U 1 + h f ( x 2 ) , ⋯
缺点是只有一阶精度,
二阶差分格式：
U j + 1 − U j − 1 2 h = f ( x j ) , U 0 = a , \frac{U_{j+1} - U_{{j-1}}}{2h} = f(x_j),\quad U_0 = a, 2 h U j + 1 − U j − 1 = f ( x j ) , U 0 = a ,
但是U − 1 U_{-1} U − 1
逼近空间
I j = [ x j − 1 2 , x j + 1 2 ] I_j = [x_{j-\frac{1}{2}}, x_{j+\frac{1}{2}}] I j = [ x j − 2 1 , x j + 2 1 ]
u h ∈ V h , V h = { v ∣ v ∣ I j ∈ P k ( I j ) } u_h \in V_h, V_h = \{ v| v|_{I_j} \in \mathcal P^k(I_j)\} u h ∈ V h , V h = { v ∣ v ∣ I j ∈ P k ( I j )}
dim ( V h ) = ( k + 1 ) N \dim(V_h) = (k+1)N dim ( V h ) = ( k + 1 ) N
乘以测试函数v v v
− ( u , v x ) j + U j + 1 2 V j + 1 2 − U j − 1 2 V j − 1 2 = ( 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{*} − ( u , v x ) j + U j + 2 1 V j + 2 1 − U j − 2 1 V j − 2 1 = ( f , v ) j ( * )
离散格式：求u h ∈ V h u_h\in V_h u h ∈ V h
− ( u h , v h ′ ) j + U h , j + 1 2 V h , j + 1 2 = ( f , v h ) j + U h , j − 1 2 V h , j − 1 2 -(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}} − ( u h , v h ′ ) j + U h , j + 2 1 V h , j + 2 1 = ( f , v h ) j + U h , j − 2 1 V h , j − 2 1
对∀ v h ∈ V h \forall v_h\in V_h ∀ v h ∈ V h j = 1 , 2 , ⋯ , N j=1,2,\cdots,N j = 1 , 2 , ⋯ , N V h ∣ I j V_h|_{I_j} V h ∣ I j V h ∣ I j = span { ϕ m j : m = 1 , 2 , ⋯ , N } V_h|_{I_j} = \text{span}\{\phi_m^{j} :m=1,2,\cdots, N\} V h ∣ I j = span { ϕ m j : m = 1 , 2 , ⋯ , N } ϕ m j ( x n j ) = { 1 , m = n , 0 , m ≠ n . \phi_m^j(x_n^j)=
\begin{cases}
1, \quad m=n,\\
0, \quad m\neq n.
\end{cases} ϕ m j ( x n j ) = { 1 , m = n , 0 , m = n .
于是在单元I j I_j I j u h = ∑ m = 1 k + 1 U m j ϕ m j u_h = \sum_{m=1}^{k+1}U^j_{m}\phi_m^j u h = m = 1 ∑ k + 1 U m j ϕ m j
在单元I j I_j I j ϕ n j , n = 1 , 2 , ⋯ , k + 1 \phi_n^j, n=1,2,\cdots,k+1 ϕ n j , n = 1 , 2 , ⋯ , k + 1 A j U j + B j U j = F j + S j , 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} A j U j + B j U j = F j + S j , j = 1 , 2 , ⋯ N . ( a )
A j = [ a m n ] ( k + 1 ) × ( k + 1 ) , a m n = ∫ I j ϕ n ϕ m ′ d x . \mathbb A^j =[a_{mn}]_{(k+1)\times(k+1)}, \quad
a_{mn} = \int_{I_j}\phi_n\phi'_m d x. A j = [ a mn ] ( k + 1 ) × ( k + 1 ) , a mn = ∫ I j ϕ n ϕ m ′ d x .
B j = [ b m n ] ( k + 1 ) × ( k + 1 ) , b m n = { 1 , m = n = k + 1 0 , 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} B j = [ b mn ] ( k + 1 ) × ( k + 1 ) , b mn = { 1 , 0 , m = n = k + 1 else
未知函数U j = [ U 1 j , U 2 j , ⋯ , U k + 1 j ] . \bm U^j = [U_1^j, U_2^j,\cdots,U_{k+1}^j]. U j = [ U 1 j , U 2 j , ⋯ , U k + 1 j ] .
很显然有U h , j − 1 2 = U 1 j U_{h,j-\frac{1}{2}}=U_{1}^j U h , j − 2 1 = U 1 j U h , j + 1 2 = U k + 1 j U_{h,j+\frac{1}{2}}=U_{k+1}^j U h , j + 2 1 = U k + 1 j
右端项
F j = [ F m j ] ( k + 1 ) × 1 , F m j = ∫ I j f ϕ m d x . \bm F^j = [F_m^j]_{(k+1)\times 1},\quad
F_m^j = \int_{I_j} f \phi_m d x. F j = [ F m j ] ( k + 1 ) × 1 , F m j = ∫ I j f ϕ m d x .
右端项
S j = [ s m j ] ( k + 1 ) × 1 \bm S^j = [s_m^j]_{(k+1)\times 1} S j = [ s m j ] ( k + 1 ) × 1 s 1 j = { U k + 1 j − 1 , j = 2 , 3 , ⋯ , N a , j = 1. s_1^j = \begin{cases}
U_{k+1}^{j-1},\quad &j=2,3,\cdots,N\\
a,\quad &j=1.\\
\end{cases} s 1 j = { U k + 1 j − 1 , a , j = 2 , 3 , ⋯ , N j = 1. s m j = 0 , m > 1. s_m^j = 0, m>1. s m j = 0 , m > 1.
注意: 线性系统(a)是逐个单元求解。
先分析k = 0 k=0 k = 0
此时u h u_h u h
取v = { 1 , x ∈ I j 0 , x ∉ I j v =\begin{cases}
1, \quad x\in I_j\\
0,\quad x \notin I_j
\end{cases} v = { 1 , x ∈ I j 0 , x ∈ / I j U j + 1 2 V j + 1 2 − U j − 1 2 V j − 1 2 = ( f , 1 ) j ≈ f ( x j ) h U_{j+\frac{1}{2}}V_{j+\frac{1}{2}}- U_{j-\frac{1}{2}}V_{j-\frac{1}{2}}=(f,1)_j \approx f(x_j)h U j + 2 1 V j + 2 1 − U j − 2 1 V j − 2 1 = ( f , 1 ) j ≈ f ( x j ) h U j − U j − 1 h = f ( x j ) \frac{U_j-U_{j-1}}{h} = f(x_{j}) h U j − U j − 1 = f ( x j ) U j + 1 2 − V j + 1 2 − − U j − 1 2 − V j − 1 2 + = ( f , 1 ) j ≈ f ( x j ) h U^-_{j+\frac{1}{2}}V^-_{j+\frac{1}{2}}- U^-_{j-\frac{1}{2}}V^+_{j-\frac{1}{2}}=(f,1)_j \approx f(x_j)h U j + 2 1 − V j + 2 1 − − U j − 2 1 − V j − 2 1 + = ( f , 1 ) j ≈ f ( x j ) h
设Ω = [ 0 , 1 ] 2 \Omega=[0,1]^2 Ω = [ 0 , 1 ] 2
a u x + b u y = f in Ω au_x +bu_y = f \quad \text{ in } \Omega a u x + b u y = f in Ω
满足边界条件u ( x , 0 ) = g ( x ) , u ( 0 , y ) = h ( y ) u(x, 0) = g(x), u(0, y)=h(y) u ( x , 0 ) = g ( x ) , u ( 0 , y ) = h ( y )
应用迎风原理，容易计算
缺点： 对于如下非线性及方程组例子不太好处理
( u 2 2 ) x = f (\frac{u^2}{2})_x = f ( 2 u 2 ) x = f u u x = f uu_x =f u u x = f
A u x = f \mathbb A \bm u_x = \bm f A u x = f A = ( 0 1 1 0 ) \mathbb A=\Big(
\begin{matrix}
0&1\\
1&0
\end{matrix}\Big) A = ( 0 1 1 0 ) ( u + v ) x = f 1 + f 2 (u+v)_x = f_1 +f_2 ( u + v ) x = f 1 + f 2
{ u t + f ( u ) x = 0 u ( x , 0 ) = u 0 ( x ) \begin{cases}
u_t + f(u)_x = 0\\[1em]
u(x,0) = u^0(x)
\end{cases} ⎩ ⎨ ⎧ u t + f ( u ) x = 0 u ( x , 0 ) = u 0 ( x )
半离散或龙格库塔方法
V h = { v : V ∣ I j ∈ P k ( I j ) } V_h = \{ v:V|_{I_j} \in P^k(I_j)\} V h = { v : V ∣ I j ∈ P k ( I j )}
离散格式：求u h ( , t ) ∈ V h u_h(,t) \in V_h u h ( , t ) ∈ V h ∀ v ∈ V h \forall v\in V_h ∀ v ∈ V h ( ( u h ) t , v ) j − ( f ( u h ) , v x ) j + U h , j + 1 2 − V j + 1 2 − − U h , j − 1 2 − V j − 1 2 + = 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 (( u h ) t , v ) j − ( f ( u h ) , v x ) j + U h , j + 2 1 − V j + 2 1 − − U h , j − 2 1 − V j − 2 1 + = 0
f ^ j + 1 2 \widehat{f}_{j+\frac{1}{2}} f j + 2 1
f ^ j + 1 2 = f ^ ( U h , j + 1 2 − , U h , j + 1 2 + ) = U h , j + 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}} f j + 2 1 = f ( U h , j + 2 1 − , U h , j + 2 1 + ) = U h , j + 2 1 +
于是有
( ( u h ) t , v ) j − ( f ( u h ) , v x ) j + f ^ j + 1 2 V j + 1 2 − − f ^ j − 1 2 V j − 1 2 + = 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 (( u h ) t , v ) j − ( f ( u h ) , v x ) j + f j + 2 1 V j + 2 1 − − f j − 2 1 V j − 2 1 + = 0
向前欧拉( u h ) t = u h n + 1 − u h n Δ t (u_h)_t = \frac{u_h^{n+1} - u_h^n}{\Delta t} ( u h ) t = Δ t u h n + 1 − u h n
( u h n + 1 − u h n Δ t , v ) j − ( f ( u h ) , v x ) j + f ^ j + 1 2 V j + 1 2 − − f ^ j − 1 2 V j − 1 2 + = 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 ( Δ t u h n + 1 − u h n , v ) j − ( f ( u h ) , v x ) j + f j + 2 1 V j + 2 1 − − f j − 2 1 V j − 2 1 + = 0
U h ( , t ) = U j ( t ) , x ∈ I j U_h(,t) = U_j(t), x\in I_j U h ( , t ) = U j ( t ) , x ∈ I j
取v = { 1 , x ∈ I j 0 , x ∉ I j v =\begin{cases}
1, \quad x\in I_j\\
0,\quad x \notin I_j
\end{cases} v = { 1 , x ∈ I j 0 , x ∈ / I j
h U j ′ ( t ) + f ^ j + 1 2 − f ^ j − 1 2 = 0. hU'_j(t) + \widehat{f}_{j+\frac{1}{2}} - \widehat{f}_{j-\frac{1}{2}} = 0. h U j ′ ( t ) + f j + 2 1 − f j − 2 1 = 0.
这种格式其实就是三点差分格式，也是一种单调格式，单调格式很简单，需要加三个条件
连续性：f ^ ( u , u ) = f ( u ) \widehat f(u, u) = f(u) f ( u , u ) = f ( u )
f ^ \widehat{f} f f ^ \widehat{f} f
DG方法求解
u x = f , x ∈ [ 0 , 1 ] ; u ( 0 ) = a . u_x = f,\quad x\in [0,1];\quad u(0) = a. u x = f , x ∈ [ 0 , 1 ] ; u ( 0 ) = a .
其中f ( x ) = cos x , a = 0 f(x)=\cos x, a=0 f ( x ) = cos x , a = 0 k = 0 , 1 , 2 k=0,1,2 k = 0 , 1 , 2 h = 1 10 , 1 20 , 1 40 , ⋯ , 1 320 h = \frac{1}{10},\frac{1}{20},\frac{1}{40},\cdots,\frac{1}{320} h = 10 1 , 20 1 , 40 1 , ⋯ , 320 1 L 1 , L ∞ L^1,L^\infty L 1 , L ∞
点我下载python代码(k k k
注: python做数值计算的优势：时代主打，免费简单，跨平台，安装包体积小且可优盘携带，第三方库多，文档多，开发快，容易混编...
上面的python代码支持任意k > 0 k>0 k > 0 k = 1 , 2 k=1,2 k = 1 , 2 f f f u h u_h u h
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) )
可以看到最终数值结果有限元误差和插值误差，在三种范数下L 1 , L 2 , L ∞ L^1,L^2,L^\infty L 1 , L 2 , L ∞ k + 1 k+1 k + 1
证明：
笔记主页