偏微分方程现代数值方法

7
2 3 4 9 10 11
2 3 6 7 8 11 12
8
4 5 6 11 12 13
3 4 7 8 9 12 13
9
6 7 8 13 14 15
4 5 8 9 10 13 14
12
10 11 12 17 18 19
7 8 11 12 13 16 17
13
12 13 14 19 20 21
8 9 12 13 14 17 18
令求解区域 Ω (x, y)|0 x 1,0 y 1 ，显然满足边界条件 u(x,t) 0 。
初始条件设为 u0 (x) xy(1 x)(1 y)
这样，题目中的 f , ,u0 就确定了，然后就可以根据缺定的 f , ,u0 ，来求数值解。
三、求偏微分方程的变分形式
NF
∑
vjΨj+
dxdy
Ω i=1
j=1
i=1
j=1
i=1
j=1
i=1
j=1
NF
= ∬ f ∑ vjΨj dxdy
Ω j=1
考虑到vj的任意性，上方程可化为
∬
NF
[∑
∂ui ∂t
ΨiΨj
+
NF
∑
ui
∂Ψi ∂x
∂Ψj ∂x
+
NF
∑
ui
∂Ψi ∂y
∂Ψj ∂y
+
NF
(∑
uiΨi)2Ψj]
dxdy
=
∬
fΨjdxdy
表 1 节点总体编码和局部编码对照表 从边界条件可以确定函数值的节点称为Ⅰ类节点，个数记为 NE。未知函数值的节点称
5
为Ⅱ类节点，个数记为 NF，是有限元空间的维数。在示例中，NE=16, NF=9。 对Ⅱ类节点，确定其影响元素集和影响节点集，如表 2 所示。
nP（节点）
Sbe （影响元素集）
Sbn （影响节点集）
图 2 网格剖分示意图 表 1 节点总体编码和局部编码对照表
一二三四五六七八九十十十十十十十 一二三四五六
Ⅰ 1 7 2 8 3 9 4 10 6 12 7 13 8 14 9 15 Ⅱ 2 6 3 7 4 8 5 9 7 11 8 12 9 13 10 14 Ⅲ 6 2 7 3 8 4 9 5 11 7 12 8 13 9 14 10
,
v)
+
a(u,
v)
+
(u2,
v)
=
(f,
v)
故可令试探函数空间和检验函数空间为
W = V = *v|v ∈ H1(Ω), v|∂Ω = 0+ 从而偏微分原方程的变分形式为：
求u ∈ V满足
4
四、区域剖分
∂u ( ∂t
,
v)
+
a(u,
v)
+
(u2,
v)
=
(f,
v),
∀v
∈
V
u(0, x) = u0(x)
a7,2u2 + a7,3u3 + a7,6u6 + a7,7u7+a7,8u8 + a7,11u11 + a7,12u12 = f7
节点基函数由ΨⅠ (λⅠ, λⅡ) = λⅠ, ΨⅡ (λⅠ, λⅡ) = λⅡ, ΨⅢ (λⅠ, λⅡ) = λⅢ 给出。其中，
λⅠ, λⅡ, λⅢ是面积坐标，由下式确定
,
Ω i=1
i=1
i=1
i=1
Ω
下面对时间进行离散 采用向前 Euler 格式，
ui(0) = ui(0)
∂ui ∂t
≈
ui(n+1) − ∆t
ui(n)
j = 1, … , NF
∬
NF
[∑
ui(n+1) − ∆t
ui(n)
ΨiΨj
+
NF
∑
ui(n)
∇Ψi∇Ψj
+
NF
(∑
ui(n)Ψi)2Ψj]
dxdy
h2 a7,11 = a7,3 = 12∆t
h2 a7,12 = a7,2 = 12∆t
e二
e三
K̂
a7,3
=
1 ∆t [ ∬ λⅡ λⅢdxdy +
∬ λⅢλⅡ dxdy]
=
2 ∆t ∬
λⅢ λⅡ|J|dλⅠdλⅡ
=
h2 12∆t
e三
e四
K̂
1
2
h2
a7,6 = ∆t [ ∬ λⅡ λⅠdxdy + ∬ λⅠλⅡ dxdy] = ∆t ∬ λⅠλⅡ |J|dλⅠdλⅡ = 12∆t
e二
十十十二二二二二二二二二二三三三 七八九十一二三四五六七八九十一二 Ⅰ 11 17 12 18 13 19 14 20 16 22 17 23 18 24 19 25 Ⅱ 12 16 13 17 14 18 15 19 17 21 18 22 19 23 20 24 Ⅲ 16 12 17 13 18 14 19 15 21 17 22 18 23 19 24 20
由φi(Ψj) = δij得
Σ = *φi(v) = v(ai), i = 1,2,3+ φi(uk) = uk(Pi)
NF
NF
φi(uk) = φi (∑ αjΨj) = ∑ αjφi(Ψj) = αi
j=1
j=1
αi = uk(Pi)
NF
uk = ∑ uiΨi
i=1
六、建立有限元方程
将上面uk，vk的表达式代入离散方程，得
在原偏微分方程两边同乘以检验函数 v，并在求解区域Ω上积分，得
∫
∂u ∂t
vdx
−
∫
∆u
∙
vdx
+
∫
u2vdx
=
∫
fvdx
Ω
Ω
Ω
Ω
应用 Green 公式并注意到边界条件
引入记号
∫
∂u ∂t
vdx
+
∫
∇u∇vdx
+
∫
u2vdx
=
∫
fvdx
Ω
Ω
Ω
Ω
a(u, v) = ∫ ∇u∇vdx
Ω
则变分形式可写为
∂u ( ∂t
∬
∂uk ∂t
vkdxdy
+
∬(∂∂uxk
∂vk ∂x
+
∂uk ∂y
∂∂vyk)dxdy
+
∬
uk2vkdxdy
=
∬
fvkdxdy
Ω
Ω
Ω
Ω
6
∬
NF
*∑
∂uiΨi ∂t
∙
NF
∑
vjΨj
+
NF
∑
ui
∂Ψi ∂x
∙
NF
∑
vj
∂Ψj ∂x
+
NF
∑
ui
∂Ψi ∂y
∙
NF
∑
vj
∂Ψj ∂y
+
NF
(∑
uiΨi)2
偏微分方程现代数值方法
———实验报告
学院：能动学院 班级：硕 0014 姓名：任小龙 学号：3110013007
2
目录
目录........................................................................................................................................... 3 一、题目 ........................................................................................................................... 4 二、具体化求解问题 ....................................................................................................... 4 三、求偏微分方程的变分形式 ....................................................................................... 4 四、区域剖分 ................................................................................................................... 5 五、构造有限原子空间 ................................................................................................... 6 六、建立有限元方程 ....................................................................................................... 6 七、编写计算程序及结果分析 ..................................................................................... 10 程序.................................................................................................................................12
e九
K̂
来自百度文库
a7,7
=
1 ∆t
[
∬
λⅠ2
dxdy
+
∬ λⅢ2 dxdy +
∬ λⅡ 2 dxdy +
∬ λⅡ 2 dxdy +
∬ λⅢ 2 dxdy +
∬ λⅠ 2 dxdy]
e二
e三
e四
e九
e十
e 十一
=
2 ∆t
∬
(λⅠ2
+
λⅡ2
+
λⅢ 2 )
|J|dλⅠdλⅡ
=
h2 2∆t
K̂
h2 a7,8 = a7,6 = 12∆t
14
14 15 16 21 22 23
9 10 13 14 15 18 19
17
18 19 20 25 26 27
12 13 16 17 18 21 22
18
20 21 22 27 28 29
13 14 17 18 19 22 23
19
22 23 24 29 30 31
14 15 18 19 20 23 24
xy
λⅠ
=
∆1 ∆
=
1 2∆
|xⅡ xⅢ
yⅡ yⅢ
λⅡ
=
∆1 ∆
=
1 2∆
xⅠ
|x xⅢ
yⅠ y yⅢ
{λⅠ + λⅡ + λⅢ = 1
各标准元上的积分为
1 1| 1
=
*(yⅡ
−
yⅢ)
x
−
(xⅡ
−
xⅢ) 2∆
y
−
(xⅡyⅢ
−
xⅢyⅡ)+
1
1| 1
=
*(yⅢ
−
yⅠ)
x
−
(xⅢ
−
xⅠ) 2∆
y
−
(xⅢyⅠ
=
∬
fΨjdxdy
,
Ω i=1
i=1
i=1
Ω
j = 1, … , NF
其中，n 时刻为已知，n+1 时刻为未知。下面就对 NF 个节点所列出代数方程： 由于边界节点上值均为 0，故所有 NF 个节点均为内部节点，以 N=4 情形下 j=7 的节点 为例说明。
图 3 节点 7 的影响元素集和影响节点集示意图（N=4） 该节点函数值满足的代数方程可写为（未知量为 n+1 时刻）
3
偏微分方程现代数值方法
一、题目
求解下列方程
������������ ������������
−
∆������
+
������2
=
������
,
������������ Ω
u(x,t) 0
, on Ω
{ ������(0, ������) = ������0(������) 其中求解区域Ω自己给定，初始值以及������(������, ������, ������)也给定并推出 f,编写程序计算该问题的
K̂
∬
0∂∂λxⅡ
∂λⅢ ∂x
+
∂λⅡ ∂y
∂∂λyⅢ1
|J|d
λⅠdλⅡ
=
0
K̂
∬
λⅠm λⅡn λⅢk dxdy
=
∬
λⅠ m λⅡ n λⅢ k |J|dλⅠdλⅡ
=
|J|
(m
m! n! k! +n+k+
2)!
∆a1a2a3
K̂
则各系数为
1
2
h2
a7,2 = ∆t [ ∬ λⅢ λⅠdxdy + ∬ λⅠλⅢ dxdy] = ∆t ∬ λⅢ λⅠ|J|dλⅠdλⅡ = 12∆t
数值解，并与真解比较进行误差分析。
二、具体化求解问题
不妨令 u xy(x 1)( y 1)e-t ，代入偏微分方程中，就可以得到
f x2 y2 xy x2 y xy2 2y2 2y 2x2 2x et x2 y2 (1 x2 )(1 y2 )e2t
−
xⅠyⅢ)+
|J| = 2∆= h2, h = N−1
7
∬
[.∂∂λxⅠ
2
/
+
.∂∂λyⅠ/2]
|J|d
λⅠdλⅡ
=
1
K̂
∬
0(∂∂λxi)2
+
(∂∂λyi
2
)1
|J|d
λⅠdλⅡ
=
1 2
,
i = Ⅱ, Ⅲ
K̂
∬
0∂∂λxⅠ
∂λj ∂x
+
∂λⅠ ∂y
∂∂λyj1
|J|d
λⅠdλⅡ
=
−
1 2
,
i = Ⅱ, Ⅲ
表 2 Ⅱ类节点影响元素集、影响节点集对照表
五、构造有限原子空间
下面构造有限元子空间Vk ⊂ V，其中Vk = span*Ψ1, Ψ2, … , ΨNF+。基函数Ψk由分片多项 式构成，且满足φi(Ψj) = δij，则变分问题的离散形式为：
求uk ⊂ Vk，使得
(∂∂utk , vk) + a(uk, vk) + (uk2, vk) = (f, vk), ∀vk ∈ Vk 其中，uk = ∑Ni=F1 αiΨi。 选取的三角形单元自由度集合为
本题的求解区域如下图 1 所示。
图 1.求解区域Ω
采用三角形单元对求解区域进行剖分。对 x, y 方向分别 N 等分，记h = N1。节点P(xi, yj)的 总体编码为np = (N + 1)(i − 1) + j, i, j = 1,2, … , N + 1，如图 2 所示。然后在每个单元上对 节点进行局部编码，生成总体编码和局部编码对照表，如表 1 所示。