当前位置:首页 >> 理学 >>

偏微分方程数值解


偏微分方程数值解 (Numerical Solution of Partial Differential Equations)
主讲:王曰朋 eduwyp@yahoo.com.cn
1

参考数目 1. George J. Haltiner, Roger Terry Williams, Numerical Prediction and Dynamic Meteorology(2nd Edition), the United States of America, 1979. 2. Curtis F.Gerald and Patrick O., Applied Numerical Analysis, Person Education, Inc., 2004. 3. Eugenia Kalnay, Atmospheric Modeling, Data Assimilation and Predictability, the press Syndicate of the University of Cambridge,2003. 4. Arieh Iserles, A First Course in the Numerical Analysis of Differential Equations, Cambridge University Press,1996. 5. 李荣华,冯国忱. 微分方程数值解. 北京:人民教育出版社,1980. 6. 徐长发,李红. 实用偏微分方程数值解法. 华中科技大学出版社,2003. 7. 沈桐立,田永祥等. 数值天气预报. 北京:气象出版社,2007.

2

数值天气预报—PDE数值解
1. 挪威气象学家V.Bjerknes(1904)提出数值预 报的思想:通过求解一组方程的初值问题可以 预报将来某个时刻的天气—思想; 2. L.F.Richardson(1922):开创了利用数值积分 进行预报天气的先例,由于一些原因(如,计 算稳定性问题“Courant,1928”)并没有取得预 期的效果—尝试; 3. Charney, Fjortoft, and Von Neumann(1950), 借助于Princeton大学的的计算机(ENIAC),利 用一个简单的正压涡度方程 (C.G.Rossby,1940)对500mb的天气形式作 3 了24小时预报---成功;

The Electronic Numerical Integrator and Computer (ENIAC).

4

常微分方程的数值解
大气科学中 常微分方程和偏微分方程的关系 1. 大气行星边界层(近地面具有湍流运动特性的大 气薄层,1—1.5km), 埃克曼(V.W.Ekman)(瑞 典)螺线的导出; 2. 1963年,美国气象学家Lorenz在研究热对流的 不稳定问题时,使用高截断的谱方法,由 Boussinesq流体的闭合方程组得到了一个完全确 定的三阶常微分方程组,即著名的Lorenz系统。
5

Lorenz系统
dx / dt = a (y - x) dy / dt = x (b - z) - y dz / dt = xy - c z
其中,a=10,(Prandtl number); b=28(Rayleigh number); c=8/3; (x,y,z)_0=(0.01;0.01;1e-10)

6

50 40 30 20 10 0 -10 -20 -30

0

5

10

15

20

25

30

35

40

45

50

7

50 40 30 20 10 0 20

0 0 -10 -20 -30
8

-20

30

20

10

9

10

Franceshini 将Navier-Stokes方程截断为五维的

截谱模型如下:

ì x1?= ? ? ? ? x?= ? 2 ? ? ? ?x =? 3 í ? ? x?= ? 4 ? ? ? x?= ? 5 ? ?

2 x1 + 4 x2 x3 + 4 x4 x5 9 x2 + 3x1 x3 5 x3 - 7 x1 x2 + Re 5 x4 - x1 x6 x5 - 3x1 x4
11

欧拉法—折线法
? 常微分方程能直接进行积分的是少数,而多数是 借助于计算机来求常微分方程的近似解; ? 有限差分法是常微分方程中数值解法中通 常有效 的方法; ? 建立差分算法的两个基本的步骤: 1. 建立差分格式,包括:a. 对解的存在域剖分; b. 采用不同的算法可得到不同的逼近误差—截断 误差(相容性);c.数值解对真解的精度—整体 截断误差(收敛性);d.数值解收敛于真解的速 度;e. 差分算法—舍人误差(稳定性).
12

2.差分格式求解 将积分方程通过差分方程转化为代数方程求 解,一般常用递推算法。 在常微分方程差分法中最简单的方法是 Euler方法,尽管在计算中不会使用,但从 中可领悟到建立差分格式的技术路线,下 面将对其作详细介绍:

13

差分方法的基本思想“就是以差商 代替微商”
考虑如下两个Taylor公式:

u (t ? h) ? u (t ) ? u?(t )h ?

1 1 1 u??(t )h2 ? u???(t ) h3 ? ? ? u ( n ) (t ) hn ? O( hn ?1 ) (1) 2! 3! n!

u (t ? h) ? u (t ) ? u?(t )h ?

1 1 1 u??(t )h2 ? u???(t )h3 ? ? ? u ( n ) (t ) hn ? O( hn ?1 ) (2) 2! 3! n!

从(1)得到:

u?(ti ) ?

u (ti ?1 ) ? u (ti ) ? O ( h) h

14

从(2)得到:

u (ti ?1 ) ? u (ti ) u?(ti ) ? ? O ( h) h
从(1)-(2)得到:

u (ti ?1 ) ? u (ti ?1 ) u ?(ti ) ? ? O(h 2 ) 2h
从(1)+(2)得到:

u??(ti ) ?

u (ti ?1 ) ? 2u (ti ) ? u (ti ?1 ) ? O(h 2 ) h2
15

对经典的初值问题
? du ? ? f (t , u ) ? dt ?u (0) ? u0 ?
满足Lipschitz条件

t ? (0, T )

f (t, u1 ) ? f (t, u2 ) ? L u1 ? u2
保证了方程组的初值问题有唯一解。

16

一、算法构造:
1. 在求解域上等距离分割: u

t0
0

tn
ti

ti ?1 ? ti ? h

T

t

微分方 2. 在 [ti , ti ?1 ]有: 程的精 确解 u (ti ?1 ) ? u (ti ) ? f (t , u ) ? O(h) h 差分方 程的精 确解

ui ?1 ? ui ? f (ti , ui ) h

ui ?1 ? ui ? hf (ti , ui )
17

3. 应用时采用如下递推方式计算:

t0
u0

t1
u1

t2
u2
u3

4. 例题

对初值问题

? y? ? 2 x ? y ? ? y (0) ? 1

x ? (0,1)

h 用Euler法求解,用 N ? 5, 即, ? 0.2
f ( x0 , y0 ) ? 1.000, y1 ? 1.200; f ( x1 , y1 ) ? 1.600, y1 ? 1.520; f ( x2 , y2 ) ? 2.320, y3 ? 1.984; f ( x3 , y3 ) ? 3.184, y4 ? 2.621; f ( x4 , y4 ) ? 4.221, y5 ? 3.465
18

5. Euler法的几何意义

ui ?1
O(h2 )
u(ti ?1 )

u?(ti )
u(ti ) ? ui

t0
0

ti

ti ?1 ? ti ? h

t

在递推的每一步,设定 u(ti ) ? ui , 方程为: ui ?1 即:

过点 (ti , ui ) 作u(t )的切线,该切线的

? ui ? u?(ti )(ti ?1 ? ti )
19

ui ?1 ? ui ? hf (ti , ui )

二、误差分析 构造算法后,这一算法在实际中是否可行呢?也就是说是否使计算机仿真 而不失真,这还需要进一步分析。 1. 局部截断误差--相容性 为了分析分析数值方法的精确度,常常在 ui ? u (ti ) 成立的假定下, 估计误差 ei ?1 ? u(ti ?1 ) ? ui ?1 这种误差称为“局部截断误差”,如图。

h2 u (ti ?1 ) ? u (ti ) ? hu?(ti ) ? u??(? ) 2 h2 ? u (ti ) ? hf (ti , u (ti )) ? u ??(? ) 2
局部截断误差是以点 ti 的精确解 u(ti ) 为出发值,用数值方法推进到下一个点

h2 ei ?1 ? u (ti ?1 ) ? ui ?1 ? u??(? ) ? O(h 2 ) 2

ti ?1 而产生的误差。

20

2.整体截断误差—收敛性 整体截断误差是以点 t 0 的初始值 u0 为出发值,用数值方法推进i+1步到点

ti ?1 ,所得的近似值ui ?1 与精确值 (t )的偏差: i ?1 ? u i ?1
称为整体截断误差。

? u(ti ?1 ) ? ui ?1

ui ?1 ? ui ? hf (ti , ui )

u(ti ?1 ) ? u(ti ) ? hf (ti , u(ti )) ? ei ?1

?i?1 ? ?i ? h[ f (ti , u(ti ) ? f (ti , ui )] ? Dh2
?i?1 ? ?i ? hL ?i ? Dh2
?i?1 ? (1? hL)i?1 ? 0 ? Dh2[1? (1? hL) ? (1? hL)2 ? ?? (1? hL)i ]
? (1 ? hL)
i ?1

Dh 2 ?0 ? [(1 ? hL)i ?1 ? 1] hL

21

? eTL ? 0 ?

Dh TL (e ? 1) L

特例,若不计初始误差,即 ? 0 则

?0

? i ?1


Dh TL ? (e ? 1) L

? i ?1 ? O(h)

3.舍入误差—稳定性 假设一个计算机仅表示4个数字(小数点后面), 那么
1 ? 0.3333, 3
x2 ? 1 9 1 x? 3

1 ? 0.1111 9

计算

x2 ?
?
x ? 0.3734

0.1112 ? 01111 ?1 0.3334 ? 0.3333

1 9 ? x?1 ? 0.6667 1 3 x ?0.3334 x? 3

22

我们的要求是:最初产生的小误差在以后的计算中虽然会传递下去,但不会无限制的

扩大,这就是稳定性所描述的问题。下面引进稳定性的概念:
设由初值 u0 得到精确解 若存在常数

un , 由初值 v0 得到精确解 vn , C 和充分小的步长h0 使得 un ? vn ? C u0 ? v0

则称数值方法是稳定的。

u

un
vn

u0
v0

0

tn

t
23

计算例题

2x ? dy ? ? y? y ? dx ? y(0) ? 1 ?
其解析解为:

x ? (0,1)

y ? 1 ? 2x
0.2000 x= 0.4000 0.6000 y= 1.3733 0.8000 1.0000

0

1.0000

1.2000

1.5315

1.6811

1.8269

24

1.9 1.8 1.7 1.6 1.5
y

1.4 1.3 1.2 1.1 1

0

0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

0.9

1

25

三、改进的Euler法 将微分方程 u?(t ) ?

f (t , u(t )) 在区间 [ti , ti ?1 ]
ti?1 ti

上积分,得到

u(ti ?1 ) ? u(ti ) ? ?
用梯形法计算积分的近似值,有

f (t , u(t ))dt

?
于是

ti ?1

ti

1 f (t , u (t ))dt ? [ f (ti , u (ti )) ? f (ti ?1 , u (ti ?1 ))] 2

1 u (ti ?1 ) ? u (ti ) ? h[ f (ti , u (ti )) ? f (ti ?1 , u (ti ?1 ))] 2
1 ui ?1 ? ui ? h[ f (ti , ui ) ? f (ti ?1 , ui ?1 )] 2

这是一个隐式格式,一般需要用迭代法来求ui ?1 ,而用显式的Euler法提供初值。
26

ui[0] ? ui ? hf (ti , ui ) ?1
k? k] ui[?1 1] ? ui ? h / 2[ f (ti , ui ) ? f (ti ?1, ui[?1 )]

为了简化计算的过程,在此基础上进一步变为如下算法:

ui ?1 ? ui ? hf (ti , ui )
ui?1 ? ui ? h / 2[ f (ti , ui ) ? f (ti ?1, ui ?1 )]
此式称为“改进的Euler法。 其局部截断误差为O(h3 ) 这个问题将在下节讨论。 接下来讨论其几何意义

预估 校正

27

u(t )
u

u(ti ?1 )

m1 ? f (ti ?1 , ui ?1 )

ui ?1 ui ?1

maverage ?

m1 ? m2 2

m0 ? f (ti , ui )

ui
0

ti

ti ?1

t

28

Euler法、改进的Euler法和解析解的比较

2 x x ? (0,1) ? dy ? ? y? y ? dx ? y(0) ? 1 y ? 1 ? 2x ?

1.9 1.8 1.7 1.6 1.5
y

1.4 1.3 1.2 1.1 1

29
0 0.1 0.2 0.3 0.4 0.5 x 0.6 0.7 0.8 0.9 1

四、(龙格-库塔)Runge-Kutta方法

简单的Euler法是建立在Taylor级数的一项展开;
改进的Euler法是以两项Taylor级数为基础建立的,如:

h2 u (ti ? h) ? u (ti ) ? hu?(ti ) ? u??(ti ) ? O(h3 ) 2 u(ti ? h) ? u(ti ) ? hu?(ti ) ? h2[(u?(ti ? h) ? u?(ti )) / h]/ 2 ? O(h3 )

u ?(ti ? h) ? u ?(ti ) ? u (ti ) ? h 2
f (ti ?1 , ui ?1 ) ? f (ti , ui ) ? ui ? h 2
如果我们截取Taylor级数的更多项会得到什么样的求解方法呢? 两个德国数学家(C.Runge & M.kutta)以这种思想为基础建立了求解微分 方程的龙格-库塔方法。它是常微分方程数值解法中使用最为广泛的方法之一。
30

一般地,一个K阶的Runge-Kutta方法可用下面的公式表示:
k ? ?ui ?1 ? ui ? ? w j K j j ?1 ? ? K1 ? hf (ti , ui ) ? ? ?? j ?1 ? ? K j ? hf (ti ? c j h, ui ? ? a jl Kl ) ? l ?1 ?

j ? 2,3,?, k

其中,

wj

是待定的加权系数, c2 , c3 ,?, ck , a21 ,?, ak ,k ?1 是待定的系数。

Euler法就是 k ? 1, w1 ? 1 的R-K法。 其系数的确定如下:将 ui ?1 展开成

h 的幂级数,并与微分方程的精确解 u(ti?1 )
p ? 1 项相同,这样确定的R-K法,

在点

ti 的Taylor展开式相比较,使两者的前

其局部截断误差为 O(h p?1 ) ,根据所得关于待定系数的方程组,求出它们的值后 代入公式,就成为一个 p 阶R-K方法。
31

例题 以二阶R-K法为例说明上述过程

?ui ?1 ? ui ? w1 K1 ? w2 K 2 ? ? K1 ? hf (ti , ui ) ? K ? hf (t ? c h, u ? a K ) i 2 i 21 1 ? 1

h2 ui ?1 ? ui ? hu?(ti ) ? u??(ti ) ?? 2
h2 ? ui ? hf (ti , ui ) ? f ?(ti , ui ) ?? 2
? ui ? hf (ti , ui ) ? h 2 [
把 K1 , K2 代入 ui ?1 中,有

1 1 f t (ti , ui ) ? f u (ti , ui ) f (ti , ui )]? 2 2

ui ?1 ? ui ? w1hf (ti , ui ) ? w2 hf [ti ? c2 h, ui ? a21hf (ti , ui )]

? ui ? w1hf (ti , ui ) ? w2h[ f (ti , ui ) ? c2 hft?(ti , ui ) ? a21hffu?(ti , ui )]
? ui ? h(w1 ? w2 ) f (ti , ui ) ? h2[w2c2 ft?(ti , ui ) ? w2a21 ffu?(ti , ui )]
32

经比较得到

? ? w1 ? w2 ? 1 ? 1 ? w2 c2 ? ? 2 ? 1 ? w2 a21 ? ? ? 2
取 c 2 为自由参数:

? w1 ? w2 ? 1 ? 1 ? w2 c2 ? ? 2 ? ?a21 ? c2 ?

1 2 c2 ? , ,1 2 3

1 3 1 1 ( w1 , w2 ) ? (0,1), ( , ), ( , ) 4 4 2 2
从而得到不同的但都是二阶的R-K方法,对应的有中点法、Heun(亨)法 以及改进的Euler法。

33

基于相同的过程,通过比较五次Taylor多项式,得到更加复杂的结果,给出了包含 13个未知数的11个方程。得到多组系数,其中常用的是以下四阶R-K法:

1 ? ui ?1 ? ui ? ( K1 ? K 2 ? K 3 ? K 4 ) ? 6 ? ? K1 ? hf (ti , ui ) ? 1 1 ? K 2 ? hf (ti ? h, ui ? K1 ) ? 2 2 ? 1 1 ? K 3 ? hf (ti ? h, ui ? K 2 ) ? 2 2 ? ? K 4 ? hf (ti ? h, ui ? K 3 ) ?
改进的Euler法、R-K法以及解析解的比较:

34

2x ? dy ? ? y? y ? dx ? y(0) ? 1 ?
1.8 1.7 1.6 1.5 1.4 1.3 1.2 1.1 1

x ? (0,1)
y ? 1 ? 2x

y

0

0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

0.9

1

35

五、线性多步(Linear Multistep Method)法 1. 预备知识:插值多项式 插值是离散函数逼近的重要方法,利用它可通过函数在有限个点处的取值状况, 估算出函数在其他点处的近似值。 从几何上理解:对一维而言,已知平面上n+1个不同点,要寻找一条n次多项式 曲线通过这些点。插值多项式一般常见的是拉格朗日插值多项式。 2. 气象应用 不均匀站点上的气象要素数据 3. 拉格朗日插值多项式 拉格朗日插值多项式逼近可能是求插值节点不均匀的插值多项式的最简单的方法。 实验观察结果或原始测量数据的分布通常是非均匀的。 例如,四个点可以确定一个三次多项式,其拉格朗日形式为: 插值

均匀网格点上的数据

P3 ( x) ?

( x ? x2 )( x ? x3 )( x ? x4 ) ( x ? x1 )( x ? x3 )( x ? x4 ) f1 ? f2 ( x1 ? x2 )( x1 ? x3 )( x1 ? x4 ) ( x2 ? x1 )( x2 ? x3 )( x2 ? x4 )
36

( x ? x1 )( x ? x2 )( x ? x3 ) ( x ? x1 )( x ? x2 )( x ? x4 ) ? f3 ? f4 ( x3 ? x1 )( x3 ? x2 )( x3 ? x4 ) ( x4 ? x1 )( x4 ? x2 )( x4 ? x3 )

4. Adams-Bashforth(阿达姆斯—贝雪福斯)公式 首先,用以下四个点对 f (t , u(t )) 进行三次Langrage插值:

(tn , f (tn , u(tn )),(tn?1, f (tn?1, u(tn?1 )), (tn?2 , f (tn?2 , u(tn?2 )),(tn?3 , f (tn?3 , u(tn?3 ))


f ( x, t ) ? ? 3 (t ) ? ?

(t ? tn ?1 )(t ? tn ?2 )(t ? tn ?3 ) f (tn , u (tn )) ? ? (tn ? tn ?1 )(tn ? tn ?2 )(tn ? tn ?3 )

(t ? tn )(t ? tn ?1 )(t ? tn ?2 ) f (tn ?3 , u (tn ?3 )) (tn ?3 ? tn )(tn ?3 ? tn ?1 )(tn ?3 ? tn ?2 )

于是,有 u (tn ?1 ) ? u (tn ) ?
tn?1 tn

?

tn?1

tn

f (t , u (t ))dt
*1 *2

? u (tn ) ? ? ? 3 (t )dt
容易算出,

?

xn?1

xn

? 3 ( x)dx ?

h [55F ( xn ) ? 59 F ( xn ?1 ) ? 37 F ( xn ? 2 ) ? 9 F ( xn ?3 )] 24
xn ? h

例如,我们可以算得

?

xn?1

xn

( x ? xn?1 )( x ? xn?2 )( x ? xn?3 ) 1 dx ? 3 ( xn ? xn?1 )( xn ? xn?2 )( xn ? xn?3 ) 6h

?

xn

( x ? xn ? h)( x ? xn ? 2h)( x ? xn ? 3h)dx
37

1 ? 3 6h

1 55 4 h ?0 (t ? h)(t ? 2h)(t ? 3h)dt ? 6h3 4 h ? 24 55
h

将(*2)代入(*1)得到Adams-Bashforth公式:

h un ?1 ? un ? [55 f ( xn , un ) ? 59 f ( xn ?1 , un ?1 ) ? 37 f ( xn ? 2 , un ? 2 ) ? 9 f ( xn ?3 , un ?3 )] 24
基于同样的计算过程,可以得到另外一个计算公式:

un ?1 ? un ?

h [9 f ( xn ?1 , un ?1 ) ? 19 f ( xn , un ) ? 5 f ( xn ?1 , un ?1 ) ? f ( xn ? 2 , un ? 2 )] 24

这称为Adams-Moulto公式。

un ?1 ? un ?

h [55 f ( xn , un ) ? 59 f ( xn ?1 , un ?1 ) ? 37 f ( xn ? 2 , un ? 2 ) ? 9 f ( xn ?3 , un ?3 )] 24
预估

un ?1 ? un ?

h [9 f ( xn ?1 , un ?1 ) ? 19 f ( xn , un ) ? 5 f ( xn ?1 , un ?1 ) ? f ( xn ? 2 , un ? 2 )] 24
校正
38

偏微分方程数值解
主讲:王曰朋

39

一、区域的离散 1.

2.

3.

40

则函数可表示为:

u ( x, y) ? u (ih, jk )
i ? 1, 2,3?

j ? 1, 2,3

二、1.(一维)一、二阶导数的有限差分近似表达式

41

2.(二维)一、二阶偏导数的有限差分近似

?u ?x

?
( xi , y j )

u ( xi ? h, y j ) ? u ( xi , y j ) h

? O ( h)

?u ?x
?u ?x

( xi , y j )

uij?1 ? uij ? h
? u ( xi ? h, y j ) ? u ( xi ? h, y j ) 2h ? O(h 2 )

( xi , y j )

?u ?x

( xi , y j )

uij?1 ? uij?1 ? 2h

? 2u ?x?y

( xi , y j )

1 uij??1 ? uij??1 ? uij??1 ? uij??1 1 1 1 ? ( 1 ) 2h 2k

O(h ) ? O(k )
2 2
42

3. 抛物型方程
(以热传导或磁扩散方程为例)
?u ? 2u ?? ?t ?x 2

初条:

u ( x ,0) ? f ( x )

(?? ? x ? ??, t ? 0)
??

初 值 问 题

精确解为 u (x , t ) ?

1 4??t

? (x ? ? ) 2 ? ? exp?? 4?t ? f (? )d? ? ? ??

不论初始分布如何集中,它 总在瞬间影响于无穷远,虽 该影响随距离按指数衰减, 然而它是以无限速度传播。 此乃抛物型方程解的特征。

43

三、热传导方程(抛物方程)
1. 热传导方程的介绍

? ?u ? 2u ? a2 2 ? ?t ?x ? ? ?u (0, t ) ? u( L, t ) ? 0 ?u ( x,0) ? f ( x) ? ? ?
j uN ? u(L, t ) ? 0

2. 离散化

u0j ? u(0, jk ) ? 0

ui0 ? u(ih,0) ? f (ih) ? fi
(1)向前差分格式:

ui ? ui ? 2ui ? u 2 u ?a k h2
j j i ?1 j

j ?1

j i ?1

44

计算:

ui

j ?1

? su ? (1 ? 2s)ui ? su
j i?1 j

j i ?1

这是一个显式格式(四点格式)

ka 2 s? 2 h
j ?1
j

i ?1

i

i ?1

每一层各个节点上的值是通过一个方程组求解得到的。这可以从

下面的计算过程看出来。
45

1 0 0 ?u1 ? su2 ? (1 ? 2s)u10 ? su0 ? 1 0 0 0 ?u2 ? su3 ? (1 ? 2s)u2 ? su1 ? 1 0 0 0 u3 ? su4 ? (1 ? 2s)u3 ? su2 ? ?? ? 0 0 0 ?u1 ?1 ? su N ? (1 ? 2s )u N ?1 ? su N ?2 ? N

系数矩阵为

s ?1 ? 2s ? ? ? s 1 ? 2s ? ? ? ? ? ? ? 1 ? 2s s ? ? ? s 1 ? 2s ? ? ?

46

计算实例:

? ?u ? u ? ? 2 ?t ?x ? ? ?u ( x,0) ? sin(? x) ?u (0, t ) ? u (1, t ) ? ?0 ? x ? 1;0 ? t ? 0.1 ?
2

47

3 2 1

u

0 -1 -2 -3 1 0.1 0.5 0.04 x 0 0.02 0 t 0.08 0.06

48

2. 向后差分格式

ui

j ?1

? ui 2 u ?a k
j
j

j ?1 i ?1

? 2ui ? u h2

j ?1

j ?1 i ?1

?suij??1 ? (1 ? 2s)uij?1 ? suij??1 ? uij 1 1
当知道第n层上的 ui 时,要确定第n+1层上各点值 ui 必须通过求解一个 线性代数方程组。
j ?1

?? su2j ?1 ? (1 ? 2 s )u1j ?1 ? su0j ?1 ? u1j ? j ?1 j ?1 j ?1 j ?? su3 ? (1 ? 2 s )u2 ? su1 ? u2 ? ? su4j ?1 ? (1 ? 2 s )u3j ?1 ? su2j ?1 ? u3j ? ?? ? j j j j ?? su N??11 ? (1 ? 2 s )u N?1 ? su N??11 ? u N ?

49

其矩阵表达式如下:

? u1j ?1 ? ? u1j ? ?s ?1 ? 2s ? ? ? ? u j ?1 ? ? u j ? ? ? s 1 ? 2s ?? 2 ? ? 2 ? ? ? ?? ? ? ?? ? ? ? ? ? j ?1 ? ? j ? ? s 1 ? 2s ? s ? ? u N ?1 ? ? u N ?1 ? ? ? ? ? ? ? ? s 1 ? 2 s ? ? u j ?1 ? ? u j ? ? ?? N ? ? N ?
j ?1
j

i ?1

i

i ?1

50

这是一个古典四点向后差分格式。计算实例

? ?u ? u ? ? 2 ?t ?x ? ? ?u ( x,0) ? sin(? x) ?u (0, t ) ? u (1, t ) ? ?0 ? x ? 1;0 ? t ? 0.1 ?
2
51

1 0.8 0.6

u
0.4 0.2 0 1 0.1 0.5 0.04 x 0 0.02 0 t 0.08 0.06

52

3. Crank-Nicolson格式,亦称六点对称格式

uij ?1 ? uij a 2 uij?1 ? 2uij ? uij?1 uij??1 ? 2uij ?1 ? uij??1 1 ? ( ? 1 ) 2 2 k 2 h h

?suij??1 ? (1 ? 2s)uij?1 ? suij??1 ? suij?1 ? (1 ? 2s)uij ? suij?1 1 1
? u1j ?1 ? ? 1 ? s ?s / 2 ? ? ? ? u j ?1 ? ? ?s / 2 1 ? s ?? 2 ? ? ? ?? ? ? ? ? ? ? j ?1 ? ? s / 2 1 ? s ? s / 2 ? ? u N ?1 ? ? ? ? ? j ?1 ? ? s 1 ? 2s ? ? u ? ? ? N ?
53

? u1j ? ?1 ? s s / 2 ? ? ??u j ? ? s / 2 1? s ?? 2 ? ? ? ?? ? ? ? ?? j ? s / 2 1 ? s s / 2 ? ? u N ?1 ? ? ? ?? j ? s / 2 1 ? 2s ? ? u ? ? ? N ?

j ?1
j

i ?1

i ?1

54

1 0.8 0.6

u
0.4 0.2 0 1 0.1 0.5 0.04 x 0 0.02 0 t 0.08 0.06

55

4. Richardson格式

ui

j ?1

? ui 2k

j ?1

u ? 2ui ? u ?a 2 h
2 j i ?1 j

j i ?1

uij?1 ? 2s(uij?1 ? 2uij ? uij?1) ? uij ?1
这是一个五点三层差分显式格式

j ?1
j

j ?1

i ?1

i

i ?1

56

讨论:

? uij 假若由于 ? 的作用,导致差分方程的近似解设为:
于是,我们可得到差分格式的误差方程如下:
j j ?ij?1 ? 2s(?i?1 ? 2?ij ? ?i?1) ? ?ij?1

t

x Richardson格式是不稳定的。
57

5. 稳定性判别 Von-Neumann 稳定性

在判断有限差分近似的稳定性方法中,以Von-Neumann方法使用较为
广泛,它仅适用于线性常系数的有限差分近似。其过程如下: 首先,要研究的差分方程可写为:

m?N1
如,

?a u
n j ?1

n?1 m j ?m

??b u
m?N0

n m j ?m

u
其次,对

n?1 j

? su
n j

? (1? 2s)u ? su
n j
n p

n j ?1

uij 进行变量分离:

2? p u ? ?V exp[i xj ] l p

58

最后将 u j

n

? V n exp[i? x j ]

代入所考察的有限差分方程。

V n?1 exp[i? x j ] ? sV n exp[i? x j ?1 ] ? (1 ? 2s)V exp[i? x j ] ? sV exp[i? x j ?1 ]
n n

V n?1 ? sV n exp[i? h] ?
定义为 放大系 数

(1 ? 2s)V n ? sV n exp[?i? h]
V n?1 ? s(cos ? h ? i sin ? h) ? (1 ? 2s) ? n V s(cos ? h ? i sin ? h) ? 1 ? 2s(1 ? cos ? h)

59

? 1 ? 4 s sin

2

?h
2
对所有的

V n ?1 下面说明,在什么条件下能使 ?1 n V

? 成立。

?1 ? 1 ? 4s sin
从上式,我们看出,

2

?h
2
2

?1

?1 ? 1 ? 4s sin

?h
2

s?

1 2sin
2

?h
2
60

1 s? 2
6. 交替显隐式格式 (1)显式预测隐式校正格式 在n+1/2层上用古典显式格式计算出过度值u 古典隐格式校正预测值,即:
n? j 1 2

,再在n+1层上用

? u n?1/ 2 ? u n u n?1 ? 2u n ? u n?1 j j j j ? j ? ? h2 ? ? 2 ? n?1 ?1 ?1 u j ? u n?1/ 2 u n?1 ? 2u n?1 ? u n?1 ? j j j ? j 2 ? ? h ? ? 2

61

(2). 跳点格式 首先将网格点(j, n)按j+n等于偶数或奇数分成两组,分别称为偶数 网点和奇数网点。从 t n 到 tn?1 的计算过程中,先在偶数网格点上

用古典显式格式计算,再在奇数网点上用古典隐格式计算,即:

? u n?1 ? u n u n?1 ? 2u n ? u n?1 j j j j ? j ? ? h2 ? ? n ? 1 ? j ? 偶数 ? n?1 ?1 ?1 u j ? un u n?1 ? 2u n?1 ? u n?1 ? j j j j ? ? ? h2 ? ? n ? 1 ? j ? 奇数 ?
62

三 双曲型方程
主要对象为 (a)一阶常系数线性双曲型方程 ?u ?u 其中a为常数 ?a ?0 ?t ?x (b)二阶常系数线性双曲型方程(波动方程) ? 2u ? 2u 其中a为常数 ? a2 2 ? 0 2 ?t ?x 这些方程的定解条件,可仅有初始条件, 也可以有初始条件和边界条件。
同椭圆型方程与抛物型方程相比,双曲 型方程差分格式的性质与定解问题解析解的 性质有更密切的关系。

63

1 一阶线性双曲型方程 (1) 初值问题
考虑
?u ?u ?a ? 0 ?? ? x ? ? ?t ?x u ( x,0) ? ? ( x) ? ? ? x ? ?

(3.1) (3.2)

由于u(x,t)沿x-t平面上方向为dx/dt=a的直线 x?at=C(C为常数)的变化率为0,即
du ?u ?u dx ?u dx ? ? ? ? ?a ?0 dt ?t ?x dt ?t dt

故沿x-t平面上任一条斜率为1/a的直线x?at=C, u(x,t)为常数。平行直线族x?at=C就是方程(3.1) 的特征线。
64

利用特征线,可以求出初值问题(3.1)、(3.2) 的解: u ( x, t ) ? ? ( x ? at )

由于u(x,t)在点 ( x0 , t0 ) 处的值依赖与?(x)在点 x0 ? at0 的值,故初始线t=0上的点 x0 ? at0 称为解 u(x,t)在点( x0 , t0 )的依赖区域。 与抛物型方程求解类似,对x-t平面进行矩 形网格剖分,x方向的步长为h,t方向的步长为?, 网点 ( x j 简记为(j,k)。 , tk )

65

(1) 偏心格式和中心差分格式 对方程(3.1),利用差商代替导数的方法,可得
u k ?1 ? u k j j u k ? u k?1 j j h u k?1 ? u k j j

?
u k ?1 ? u k j j

?a
?a

?0
?0

即 u k ?1 ? (1 ? ra)u k ? rau k?1 j j j
u k ?1 ? rau k?1 ? (1 ? ra)u k j j j
k ?1 j k j

?
u k ?1 ? u k j j

r u ? u ? a(u k?1 ? u k?1 ) ?a ?0 j j 2 ? 2h 其中 r ? ? / h j ? 0, ? 1, ? 2, ?, k ? 0, 1, 2, ? 前两个格式的局部截断误差阶为 O(? ? h) , 分别称为左、右偏心格式。

h u k?1 ? u k?1 j j

第三个格式的截断误差阶为 O(? ? h 2 ),它 称为中心差分格式。

66

从差分格式依赖区域和微分方程依赖区域的 关系,可以得到差分格式收敛的必要条件:
差分格式的依赖区域包含微分方程的依赖区 域(也称为CFL条件)。 对于左偏心格式, CFL条件为:a>0,且 ar ? 1 。 对于右偏心格式, CFL条件为:a<0,且 ar ? 1。 对于中心格式, CFL条件为: ar ? 1 。

可以证明:左、右偏心格式的稳定条件与其 CFL条件相同,但中心差分格式是恒不稳定的,
所以,当 a>0(或a<0)时,左(或右)偏心格 式才有实用值;中心差分格式无实用价值。
67

(2)Lax格式 对方程(3.1),可利用特征线构造格式。 P 设系数a>0, ? t k 上网 t t ? t k ?1 ? 点A(j?1,k),B(j,k),C(j+1,k) 处的解值已经算出,从点 P(j,k+1)作特征线,它与线 t ? tk ?? ? ? 段AB交于点 D。 C A D B h ? a? h ? a? 由u(p)=u(D),有 u ( D) ? u (C ) ? u ( A) 2h 2h 这样,得到Lax格式: 1 k ?1 u j ? [(1 ? ar )u k?1 ? (1 ? ar )u k?1 ] j j 2 j ? 0, ? 1, ? 2, ?, k ? 0, 1, 2, ? 当 ar ? 1,Lax格式稳定,截断误差阶为O(? ? h 2 ? h 2 ? )。
68

(3) Lax—Wendroff格式 对方程(3.1),利用特征线作二次插值,即可 得到Lax—Wendroff格式: ar k a 2r 2 2 k u k ?1 ? u k ? (u j ?1 ? u k?1 ) ? ?x uj j j j 2 2 当 ar ? 1 Lax—Wendroff格式稳定,它的截 , 断误差阶为 O(? 2 ? h 2 )。

69

(2) 初边值问题 应该注意:边值条件的给法与其它两类方程不同。 Y 如果 a>0,方程特 征线向右倾,只能在 x 变化区域的左边界上给 a>0 出边界条件: u(0, t ) ? ? 1 (t ) X 0 ? x< ? O 如果 a<0,方程特 Y 征线向左倾,只能在 x 变化区域的右边界上给 出边界条件,即使 x 的 a<0 变化区域为0 ? x ? d , 也只能给出边界条件: X 0 ? x ?d O u(d , t ) ? ? 2 (t ) 70

设常数a>0 ,考虑下面模型问题:
?u ?u ?a ?0 0?t ?T 0? x ?? ?t ?x u ( x,0) ? ? ( x) 0 ? x ? ? u (0, t ) ? ? (t ) 0 ? t ? T

前面建立的几个显格式,都适用于这个问题。 下面建立隐格式。

71

(1)最简隐格式
u
k ?1 j

?u

k j

h k ? 0, 1, 2, ?, j ? 1, 2, ?

?

?a

u

k ?1 j

?u

k ?1 j ?1

?0

( j ? 1, k ? 1) ( j, k ? 1)

?

? ?

( j, k )
k ?1 j

令r ? ? / h ,格式可改写为 u
k u 0 ? ? j , u0 ? ? k j

连同初始条件与边界条件:

ar k ?1 1 ? u j ?1 ? uk j 1 ? ar 1 ? ar

j ? 1, 2, ?, k ? 0, 1, 2, ?

该格式可在0 < x, t < ?内所有网点上显示地计算 解之近似值。

该格式的局部截断误差阶为 O(? ? h) 。
72

(2)Wendroff格式 1 1 在点 P( j ? , k ? ) 处,用 2 2 1 ? ? ?u ? ? ?u ? ? ? ?u ? ? ? ? ?t ? ? 2 ? ? ?t ? ? ? ?t ? ? ? ?P ? ? E ? ?G ? ?

D

?

F

? t ? t k ?1
C E

G

P

1 ? ? ?u ? ? ?u ? ? ?u ? ? ? t ? tk A? B H ? ?? ? ? ? ? ? ? ? ? ?x ? x ? jh ? ? P 2 ? ? ?x ? H ? ?x ? F ? x ? ( j ? 1)h 然后用中心差商逼近这些导数值,则可得到 Wendroff格式:
?1 ?1 u k ?1 ? u k u k?1 ? u k?1 u k ? u k?1 u k ?1 ? u k?1 1 j a j j j j j j j ( ? ) ( ? ? ) 0 ? 2 ? ? 2 h h k ? 0, 1, 2, ?, j ? 1, 2, ?
73

令r ? ? / h ,格式可改写为 1 ? ar k k ?1 k ?1 u j ? u j ?1 ? (u j ? u k?1 ) j 1 ? ar 连同初始条件与边界条件:
k u 0 ? ? j , u0 ? ? k j ? 1, 2, ?, k ? 0, 1, 2, ? j

该格式可在0 < x, t < ?内所有网点上显示地计 算解之近似值。 该格式的局部截断误差阶为 O(? 2 ? h2 ),且 无条件稳定。

74

2 二阶线性双曲型方程(波动方程) ? 2u ? 2u ? a 2 2 ? 0 ?? ? x ? ? a为正常数 (3.3) 考察 2 ?t ?x ?u (3.4) u ( x,0) ? ? ( x) ( x,0) ? ? (t ) ?? ? x ? ? ?t 对x-t平面进行矩形网格剖分,x方向的步长 为h,t方向的步长为?,网点( x j , t k ) 简记为(j,k)。 用二阶中心差商代替(3.3) k ?1 ? 中的二阶导数,则得到网点 ? ? ? k (j,k)处的差分方程: k ?1 ? 1 2 k a2 2 k ? uj ? 2 ?x uj j j ?1 j ?1 2 t ? h 或 u k ?1 ? a2r 2 (u k?1 ? u k?1 ) ? 2(1 ? a 2r 2 )u k ? u k ?1 j j j j j 其中 r ? ? / h; j ? 0, ? 1, ? 2, ?, k ? 0, 1, 2, ? 。
75

初始条件 u ( x,0) ? ? ( x) 离散: 取为 u 0 ? ? j j ?u ( x,0) ? ? (t ) 离散: 初始条件 ?t ?1 1 (u j ? u ?1 ) ? ? j j ? 2? ? u ?1 ,得 由 ? 消去 j 2 ? 1 ? 2u 0 ? a ? 2u 0 ?? 2 t j h 2 x j ? r2 2 u1j ? a [? j ?1 ? ? j ?1 ] ? (1 ? r 2 a 2 )? j ? ?? j 2 上述差分格式与初始条件的截断误差阶均 为 O(? 2 ? h 2 )。 该格式稳定的充分条件为 ar ? 1。
76

上述方法也可用于求解初边值问题:
? 2u ? 2u ? a 2 2 ? 0 0 ? x ? 1, 0 ? t ? T ?t 2 ?x ?u u ( x,0) ? ? ( x) ( x,0) ? ? (t ) 0 ? x ? 1 ?t u (0, t ) ? ? (t ) u(1, t ) ? ? (t ) 0 ? t ? T

3 交替方向隐格式

77

由于实际问题的具体特征、复杂性以 及算法自身的适用范围决定了应用中必须 选择、设计适合于自己特定问题的算法, 因而掌握数值方法的思想至关重要。
要在各种可能的求解方法中找到一种 统一地适用于计算材料学领域(或其它领 域)的理想方法,一般是不现实的。 任何模拟方法,都必须在最佳计算速 度和数值精度之间寻找平衡点。

科学计算(数值模拟)已经被公认为 与理论分析、实验分析并列的科学研究三 大基本手段之一,但三者之间相辅相成。
78

79

80

第三部分 偏微分方程的有限元方法
一 边值问题的变分原理
1 引论 (1)等周问题 在长度一定的所有平面封闭曲线中,求 所围面积为最大的曲线。
? dx ? ? dy ? 模型:在条件 ? ? ? ? ? ? ds ? l 下 s1 ? ds ? ? ds ? 1 s2 ? dy dx ? 求使得泛函 s( x, y ) ? ?s ? x ? y ?ds 2 1 ? ds ds ? 达到最大的函数 x( s), y ( s) 。
s2 2 2

81

定义:当求泛函在一个函数集合K中的极小 (或极大)问题,则该问题称为变分问题。 变分问题与微分方程的定解问题有一 定的联系。 (2)初等变分原理 ① 一元二次函数的变分原理 设 J ( x) ? Lx 2 ? 2 fx ( L ? 0, L, f为实常数) 考察J(x)的极值情况。 变分原理: 求 x0 ? R,使 J ( x0 ) ? min J ( x)
x?R

与求解方程 Lx = f 等价。
82

② 多元二次函数的变分原理

n ? 1 n n J ( x ) ? J ( x1 , x2 ,?, xn ) ? ?? aij xi x j ? ? bi xi 2 i ?1 j ?1 i ?1

求J(x)取极小值的驻点, 其中
? a11 a12 ? a1n ? ?a a22 ? a2 n ? ? 对称正定 A ? ? 21 ? ? ? ? ? ? ? ? ?an1 an 2 ? ann ? ? ? T 设 x ? ( x1 , x2 ,?, x n ) b ? (b1 , b2 ,?,b n )T ? ? ? 1 ? ? 则J(x)可表示为: J ( x ) ? ( Ax , x ) ? (b , x ) 2
83

变分原理:
设矩阵A对称正定,则下列两个命题等价: ? ? ? n (a) 求 x0 ? R ,使 J ( x0 ) ? minn J ( x ) x?R ? ? ? 1 ? ? 其中 J ( x ) ? ( Ax , x ) ? (b , x ) 2? ? ? (b) x0 是方程 Ax ? b 的解 上述两个例子表明: 求二次函数的极小值问题和求线性代数 方程(组)的解是等价的。

84

2 两点边值问题的变分原理 (1)弦平衡的平衡原理与极小位能原理 考察一根长为 l 的弦,两端固定在点 A(0,0) 和B(l,0)。当没有外力作用时,它的位置沿水平 方向与X轴重合。设有强度为f(x)的外荷载垂直 向下作用在弦上,于是弦发生形变。假定荷载 很小,因而发生的形变也很小。用 u(x) 表示在 荷载f(x)的作用下弦的平衡位置。 B(l ,0) A(0,0) ? X ?

u
85

平衡原理 求弦的平衡位置归结为求解两点边值问题: ?? Tu??( x) ? f ( x) 0 ? x ? l 其中T是弦的张力。 ? ?u (0) ? 0 u (l ) ? 0 极小位能原理: 弦的平衡位置 (记为 u? ? u? (x))将在满足边 值条件 u(0)=0,u(l)=0 的一切可能位置中,使 位能取极小值。 设弦处于某一位置u=u(x),可得到其总位能为 1 l J (u ) ? ? (Tu?2 ? 2uf )dx 2 0 弦的平衡位置 u? ? u? (x)是下列变分问题的解 J (u? ) ? min2 J (u )
u?C0
86

在数学上,要将某个微分方程的定解问题 转化为一个变分问题求解,必须针对已给的定 解问题构造一个相应的泛函,并证明定解问题 的解与泛函极值问题的解等价。 有限元方法正是利用这种等价性(边值问 题与变分问题的等价性),先将微分方程定解 问题转化为变分问题(或变分方程)的求解问 题,然后再设法近似求解变分问题(或变分方 程)。

87

(2)两点边值问题的变分原理 考察二阶常微分方程边值问题: ? d ? du ? ? Lu ? ? ? p ? ? qu ? f x ? (a, b) dx ? dx ? ? ?u (a) ? 0 u?(b) ? 0 ? ① 构造泛函 1 J (u ) ? ( Lu, u ) ? ( f , u ) 2 b b 1 b d ? du ? 2 ?? ? ? p ?udx ? ?a qu dx ? ?a fudx 2 a dx ? dx ? 1 b ? ? ( pu?2 ? qu 2 ? 2 fu)dx 2 a b du dv 引入泛函算子 ? (u, v) ? ? [ p ? quv]dx a dx dx 1 则 J (u ) ? ? (u, u ) ? ( f , u ) 2

88

② 变分问题

与前述二阶常微分方程边值问题相应的变分 问题是
1 求 u? ? H E ,使 J (u? ) ? min1 J (u )
u?H E

1 其中 H E [a, b] ? {u ( x) u ( x) ? H 1[a, b], u (a) ? 0} 1 J (u ) ? ? (u, u ) ? ( f , u ) 2 1 b ? ? ( pu?2 ? qu 2 ? 2 fu)dx 2 a

89

③ 变分原理(变分问题与边值问题的等价性) 设 f ? C 0 ( I ) , * ? C 2 是边值问题 u
? d ? du ? ? Lu ? ? ? p ? ? qu ? f dx ? dx ? ? ?u (a) ? 0, u?(b) ? 0 ? x ? ( a, b)

的解,则 u* 使 J(u) 达到极小值; 1 反之,若 u? ? C 2 ? H E使 J(u) 达到极小值, 则 u* 是边值问题的解。 1 其中 J (u ) ? ? (u, u ) ? ( f , u ) 2 u (a) ? 0是强制边界条件,u?(b) ? 0 是自然边 界条件,区别这两类边界条件在用有限元方法求 解边值问题时很重要。 90

(3)虚功原理
? d ? du ? ? Lu ? ? ? p ? ? qu ? f x ? (a, b) 对两点边值问题:? dx ? dx ? ?u (a) ? 0, u?(b) ? 0 ? 1 设 v ? H E ,以v乘方程两端,沿[a,b]积分, 并利用 v(a) ? 0, v?(b) ? 0 ,得变分方程 1 ? (u, v) ? ( f , v) ? 0 ?v ? H E b du dv 其中 ? (u, v) ? ? [ p ? quv]dx a dx dx 虚功原理 u* ? C 2 设 ,则 u*是边值问题解的充要条件是: 1 u* ? H E,且满足变分方程: 1 v ? HE 对任意 ? (u* , v) ? ( f , v) ? 0 在力学里, (u, v) ? ( f , v)表示虚功 91 ?

对于复杂的边界条件,边值问题的求解一 般是困难的。若将微分方程化为相应的变分问 题或变分方程,则只需处理强加边界条件,无 需处理自然边界条件(自然边界条件已包含于 变分问题中泛函的构造或已包含于给出的变分 方程之中)。这一特点对研究微分方程离散化 方法及其数值解带来了极大的方便。

92

3 二阶椭圆边值问题的变分原理 ? ? ? 2u ? 2u ? ?? ?u ? ?? 2 ? 2 ? ? f ( x, y ) ( x, y ) ? G ? ?x 模型方程 ? ?y ? ? ? ?u ? 0 ? |? 其中G是平面有界区域。
(1)极小位能原理 ① 构造泛函 1 J (u ) ? (??u, u ) ? ( f , u ) 2 ?u ?v ?u ?v ? )dxdy 引入泛函算子 ? (u, v) ? ?? ( ?x ?x ?y ?y G 1 则 J (u ) ? ? (u, u ) ? ( f , u ) 2

93

② 变分问题 与前述二阶椭圆边值问题相应的变分问题是
1 求 u? ? H 0 (G),使 J (u? ) ? min J (u ) 1
u?H 0 ( G )

1 其中 H 0 (G) ? {u( x, y) | u( x, y) ? H 1 (G), u( x, y) ? ? 0} 1 J (u ) ? (??u, u ) ? ( f , u ) 2 1 ? ? (u, u ) ? ( f , u ) 2

94

③ 变分原理(变分问题与边值问题的等价性) 设 f ? C 0 ( I ) ,* ? C 2 (G )是二阶椭圆边值问题 u 的解,则 u* 使 J(u) 达到极小值; 1 反之,若 u? ? C 2 (G ) ? H 0 (G) 使 J(u) 达到极 小值,则 u*是二阶椭圆边值问题的解。 1 其中 J (u ) ? ? (u, u ) ? ( f , u ) 2 ? 对第一边值问题,无论齐次或非齐次边界条 件,泛函是一样的,只是边界条件要作为强加 边值条件加在所取的函数类上。 ? 对第二、三类边值问题,无论齐次或非齐次 边界条件,二次泛函形式相对于第一边值问题 有所改变,但函数类的选取与边界条件无关。
95

(2)虚功原理
? ?? ?u ? f ( x, y ) ? ? ?u|?1 ? 0 ? ? ?u ? au|? ? 0 2 ? ?n ? ( x, y ) ? G

问题

a?0

1 设 v ? H E (G) ,以v乘方程两端后在G上积分, 并利用Green公式,得变分方程

? (u, v) ? ( f , v) ? 0

1 ?v ? H E (G)

? ?u ?v ?u ?v ? 其中 ? (u, v) ? ?? ? ? ?x ?x ? ?y ?y ?dxdy ? ??auvds ? 2 ? G ? 1 H E (G) ? {u ( x, y ) | u ( x, y ) ? H 1 (G), u ( x, y) ?1 ? 0}
96

虚功原理
u? ? C 2 (G )是边值问题的解,则对任意 设 1 u v ? H E (G) , ? 满足变分方程。

反之,若 u? ? C (G ) ? H (G),且对任意 1 v ? H E (G) 满足变分方程,则 u? 为边值问题的解。
2 1 E

在力学里, (u, v) ? ( f , v)表示虚功 ?

? 与极小位能原理类似,第一类边界条件为强 加边界条件,第二、三类边界条件为自然边界 条件。 ? 虚功原理比极小位能原理应用更广。
97

4 Ritz-Galerkin方法 目的:求解相应的变分问题或相应的变分方程。 Ritz方法是近似求解变分问题(即二次泛函极 小值)的算法。Galerkin方法是近似求解变分方程 的算法,这两种算法统称为Ritz-Galerkin方法。 Ritz-Galerkin方法的基本思想 用有限维空间的函数代替变分问题(或变分 方程)中无限维空间的函数,从而在有限维函数 空间中求变分问题(或变分方程)的近似解,并 要求当有限维空间的维数不断增加时,有限维 近似解逼近原变分问题(或变分方程)的解。 1 以下用V表示 H 0 , H E , H 1 等Sobolev空间, L表示微分算子,?(u,v)为由L及边值条件决定 的双线性泛函。

98

(1)Ritz方法
由极小位能原理得出的变分问题为: 求 u? ?V,使 J (u? ) ? min J (u ) u?V 1 其中 J (u ) ? ? (u, u ) ? ( f , u ) , 2 Ritz方法:求变分问题的近似解。 ? 设 Vn 是V 的n维子空间, 1 , ?2 ,?, ?n 是 Vn 的一 组基底(称为基函数) 。Vn 中任一元素 u n可表示为
un ? ? ci?i
i ?1 n

Ritz方法: 求 un ?Vn,使

J (un ) ? min J (v)
v?Vn
99

即选择适当的 c1 , c2 ,?, cn ,使 J (un ) 取极小值。

1 展开J (u u)n ) ? ? (un , un ) ? ( f , un ) J( n 2
n 1 n n ? ??? (?i , ? j )ci c j ? ? ( f , ? j )c j 2 i ?1 j ?1 j ?1



?J (un ) ?0 ?c j
n

j ? 1, 2, ?, n

则 c1 , c2 ,?, cn满足

?? (? ,? )c
i ?1 i j

i

? ( f ,? j )

j ? 1, 2, ?, n
n

解出 ci (i ? 1,?, n) 代入 u n,则得 un ? ? ci?i
i ?1
100

Ritz方法步骤为: ① 根据最小位能原理构造相应于微分方程或物 理问题的变分问题; ② 取 {?i }(i ? 1,?, n) 作为 Vn 的一组基底,即用 Vn ? span{?1 , ?2 ,?, ?n }近似代替无穷维空间V; ③ 根据二次函数取极值的必要条件,得到
un ? ? ci?i ? Vn 中 ci 所满足的方程组:
i ?1 n n

?? (? ,? )c
i ?1 i j

i

? ( f ,? j )

j ? 1, 2, ?, n

④ 求解关于 ci的线性代数方程组。
101

(2) Galerkin方法 由虚功原理得出的变分方程为:

? (u, v) ? ( f , v) ? 0
Galerkin方法:求变分方程的近似解。

? 设 Vn 是V 的n维子空间, 1 , ?2 ,?, ?n 是 Vn 的一 组基底(称为基函数) 。Vn 中任一元素 u n可表示为
un ? ? ci?i
i ?1 n

u Galerkin方法: 求 un ?Vn,使对 ?v ?Vn , n 满足

? (un , v) ? ( f , v) ? 0
102

将 un ? ? ci?i ? Vn 代入变分方程,则
i ?1
n

n

? (? ci?i , v) ? ( f , v)
i ?1

?v ?Vn

由 v ?Vn 的任意性,取?1 , ?2 ,?, ?n 作为v ,则得

?? (? ,? )c
i ?1 i j

n

i

? ( f ,? j )

j ? 1, 2, ?, n
n

解出 ci (i ? 1,?, n) 代入 u n,则得 un ? ? ci?i
i ?1

103

Galerkin步骤为: ① 根据虚功原理构造相应于微分方程或物理问 题的变分方程; ② 取 {?i }(i ? 1,?, n) 作为 Vn 的一组基底,即用 Vn ? span{?1 , ?2 ,?, ?n }近似代替无穷维空间V; ③ 取 ?1 , ?2 ,?, ?n作为v ,将 un ? ? ci?i ? Vn 代 入变分方程,得到 ci 满足的方程组:
i ?1 n

?? (? ,? )c
i ?1 i j

n

i

? ( f ,? j )

j ? 1, 2, ?, n

④ 求解关于 ci的线性代数方程组。
104

Ritz-Galerkin方法应用的困难 ① 基函数选取必须满足强加边界条件,因此 选取困难;
② 计算量、存储量巨大; ③ 方程组求解病态严重。 有限元法广泛应用的原因 ① 充分发挥了变分形式和Ritz-Galerkin方法的 优点; ② 摆脱了传统的基函数取法; ③ 各种问题的结构程序格式统一。
105

二 椭圆型方程的有限元方法
差分法解偏微分方程,解得的结果就是准确 解u在节点上的近似值;
Ritz-Galerkin方法得到近似的解析解,但对一 般区域,却往往难以实现。

有限元方法基于变分原理, 又具有差分方法 的一些特点,并且适于较复杂的区域和不同粗细 的网格。
有限元方法与传统Ritz-Galerkin方法的差别在 于有限维函数空间的构造方法。Ritz-Galerkin方法 选用的基函数在整个定解区域上整体光滑,有限 元则取分段或分片连续且局部非零的基函数。
106

1 一维问题的线性元 考虑两点边值问题: d du ? a? x?b ? Lu ? ? ( p ) ? qu ? f dx dx ? ?u (a) ? 0, u?(b) ? 0 ? 将区间[a,b]分割为n个子区间 [ xi ?1 , xi ](i ? 1,2,...n)。 第i个单元记为I i ? [ xi ?1 , xi ],其长度 hi ? xi ? xi ?1。 (1)试探函数与试探函数空间 设
1 Vh ? {uh ( x) uh ( x) ? H E ( I ),

且uh ( x)在I i 上为次数不超过m的多项式} u 则 Vh称为试探函数空间, h ?Vh 称为试探函数。
107

(2) 用单元形状函数表示试探函数

最简单的试探函数空间Vh 由分段线性函数组成。 设在节点上试探函数 u h 在节点上的一组值为 u0 ? 0, u1 , u2 ,..., un
在第i个单元 I i ? [ xi ?1 , xi ] 上的线性插值函数为 x ? xi x ? xi ?1 uh ( x) ? ui ?1 ? ui x ? I i xi ?1 ? xi xi ? xi ?1 即
xi ? x x ? xi ?1 uh ( x) ? ui ?1 ? ui hi hi x ? Ii

u 当 x ? I i 时, h (x) 的(线性)插值公式称 为(线性)单元形状函数。
108

把每个单元形状函数合并起来,就得到整 个区间[a,b] 上都有定义的函数 uh (x) :
x ? x0 ? x1 ? x x ? I1 ? h u0 ? h 1 ? 1 ??? ?? ?? ?? ? ? xi ? x ui ?1 ? x ? xi ?1 ui x ? Ii ? hi hi uh ( x) ? ? ? xi ?1 ? x u ? x ? xi u x ? I i ?1 i i ?1 ? hi ?1 hi ?1 ? ?? ?? ?? ?? ? ? xn ? x x ? xn ?1 un ?1 ? un x ? In ? hn ? hn ui ?1 ui u1 ui ?1

un

x0

x1

xi ?1

xi

xi ?1

xn

109

为使分段插值标准化,通常用仿射变换 x ? xi ?? hi 把 ? ? [0,1] 变到 x ?[ xi ?1 , xi ],令 则 变为
N0 (? ) ? 1 ? ? N1 (? ) ? ? xi ? x x ? xi ?1 uh ( x) ? ui ?1 ? ui hi hi uh ( x) ? N 0 (? )ui ?1 ? N1 (? )ui ~ u (? ) ? N (? )u ? N (? )u
h 0 i ?1 1

x ? Ii

x ? Ii


显然

i

? ?[0,1]

N0 (0) ? 1 N0 (1) ? 0 N1 (0) ? 0 N1 (1) ? 1
110

(3) 用节点基函数表示试探函数 定义基函数系 ?i (x)
? x ? xi ?1 ? ? hi ? xi ?1 ? x ? ?i ( x) ? ? ? hi ?1 ? ? 0 ? ? x ? [ xi ?1 , xi ] x ? [ xi , xi ?1 ] x ? [ xi ?1 , xi ?1 ] i ? 1, 2, ?, n ? 1

? x1 ? x ? h ? 1 ? 0 ( x) ? ? ? 0 ? ?

? x ? xn ?1 x ? [ x0 , x1 ] ? h ? n ? n ( x) ? ? ? 0 x ? [ x0 , x1 ] ? ?

x ? [ xn ?1 , xn ] x ? [ xn ?1 , xn ]

111

?0 ( x), ?1 ( x),..., ?n ( x) 线性无关,它们可组 成试探函数空间的基,常称为节点基函数。
几何形状如图 ? 0 ( x) ? i (x)
i ? 1,?, n ? 1

? n (x)

xi ?1 xi xi ?1 xn ?1 xn 任一试探函数 uh ?Vh 可表示为
uh ( x) ? ? ui?i ( x)
i ?0 n

a

x 0 x1

b

x ? [ a, b]

用这类插值型基函数,可以构造出适合各 种边界条件的试探函数。
112

节点基函数可用变量 i ?1 若借助前述放射变换 ? ??表示为
? ? N1 (? ) ? ? ?i ( x) ? ? N 0 (? ) ? ? 0 ? ?
x?x hi

x ? [ xi ?1 , xi ] x ? [ xi , xi ?1 ] 其它

??

x ? xi ?1 hi i ? 1,..., n ? 1 2,

x ? xi ?? hi ?1

? ? N 0 (? ) ? 0 ( x) ? ? ? 0 ? ? ? N1 (? ) ? n ( x) ? ? ? 0 ?

x ? [ x0 , x1 ] 其它 x ? [ xn ?1 , xn ] 别处

x ? x0 ?? h1 x ? xn ?1 ?? hn
113

(4)从Ritz方法出发形成有限元方程 ① 直接形成有限元方程 (a) 把表达式 uh ? ? ui?i ( x)代入泛函;
i ?1 n

(b) 将泛函表达式中积分区间[a,b]变到[0,1]; (c) 由 J (uh )达到极小值的条件 ?J (uh ) ?0 j ? 1,2,..., n ?u j 得到含 u j ?1 , u j , u j ?1 ( j ? 1, 2, ?, n)的有限元方程 a j ?1, j u j ?1 ? a jju j ? a j ?1, j u j ?1 ? b j ? 0 ( j ? 1, 2, ?, n) 这儿 u0 ? 0 an?1,n ? 0 (d) 解出有限元方程的数值解 ui (i ? 1,..., n),就 得到使二次泛函取极小的近似函数(有限元解)
uh ( x) ? ? ui?i ( x)
i ?1 n
114

有限元方程可用矩阵表示为 ? ? Ku ? b
其中
? a11 a21 ?a a22 a32 12 K ?? ? ? ? ? ? 称为总刚矩阵。 ? ? ? ? ? ann ?

? an ?1,n

115

② 用单元刚度分析形成有限元方程
工程中形成有限元方程时,通常先在每个 单元上形成单元矩阵(称为单元刚度矩阵), 然后由单元刚度矩阵形成总刚度矩阵(称为总 体合成)。 (a) 把 J (uh )按单元组织,则在第i个单元上,令 ( ? (i ) ? f i ?i1) ? ui ?1 ? ? (i ) ? (i ) (i ) ? ai ?1,i ? ai ,i ?1 u ? ? f ? ? (i ) ? ?u ? ?f ? ? i ? ? i ?
i) i) ? ai(?1,i ?1 ai(?1,i ? ? k (i ) ? ? (i ) (i ) ? ?a ? i ,i ?1 ai ,i ? k (i ) 称为单元刚度矩阵。各元素可计算 其中 得到。
116

再把 k (i )扩展成n?n矩阵,使其第i?1行、第i行 和第i?1列、第i列交叉位置的元素就是单元刚度矩 阵的四个元素,其余全为零( k (1) 只是第一行,第 一列元素非零)。即 ?? ?? ? ? i) i) ai(?1,i ?1 ai(?1,i ? K (i ) ? ? (i ) (i ) ? ? ai ,i ?1 ai ,i ? ? ?? ?? ? ? T 记 u ? (u1 , u2 ,..., un ) b ? (b1 , b2 ,..., bn )T ? ? 1 ?T ? ?T ? 1 ? ? 则 J (un ) ? u Ku ? u b ? ( Ku , u ) ? (b , u ) 2 2

其中 K ? ? K (i ) 称为总刚矩阵。
i ?1
117

n

(b) 由 J (uh ) 达到极小值的条件 ?J (uh ) ?0 j ? 1,2,..., n ?u j ? ? 得到有限元方程 Ku ? b 。
(c) 解出有限元方程的数值解 ui (i ? 1,..., n),就 得到使二次泛函取极小的近似函数(有限元解)
uh ( x) ? ? ui?i ( x)
i ?1 n

118

(5)从Galerkin方法出发形成有限元方程
把表达式 uh ? ? ui?i ( x)代入变分方程
i ?1 n

对前面的两点边值问题,变分方程变为

其中

?? (? ,? )u ? ? f? dx ? (u, v) ? ? ( pu?v? ? quv)dx
b i ?1 i j i a j

n

j ? 1,2,..., n

b

a

该方程即为Galerkin法形成的有限元方程。 ? 与Ritz方法相比, Galerkin方法形成的有限 元方程其系数矩阵就是总刚矩阵。 ? 由Galerkin方法推导有限元方程更加方便直 接,且适用面广。
119

2 一维问题的高次元

若希望在每个单元上提高逼近的精确度,则 可通过提高插值多项式次数来实现,
整个问题计算的全过程除分析单元插值外, 均与前面框架类似。

在单元 I i ? [ xi ?1 , xi ] 上可构造一、二、三及高 次插值多项式,其方法有两种:
① Lagrange型:在单元内部增加一些插值节点。

② Hermite型:在节点引进一阶、二阶乃至更高 阶导数。
120

① 线性元(Lagrange型) 要求:在每一个单元上是一次多项式,在单元 节点处连续。 插值条件:在单元的两个端点取指定值。 ② 二次元 (Lagrange型) 要求:在每一个单元上是二次多项式,在单元 节点处连续。 插值条件:在单元的两个端点及单元中点取指 定值。 ③ 三次元(Hermite型) 要求:在每一个单元上是三次多项式,在单元 节点处连续。 插值条件:在两个端点取指定的函数值和一阶 121 导数值。

采用高次元,有限元方程形成的方法和线性 元类似,但工作量增加。一是计算积分的复杂性 增加,二是矩阵的带宽增加。 高次元的主要优点是收敛阶高,且提高了函 数逼近的光滑性。

122

3 二维问题的矩形元

假定区域G可以分割成有限个矩形的和,且 每个小矩形(单元)的边和坐标轴平行。
采用矩形剖分后,任一个矩形 Rij ? {xi ? x ? xi ? ?x; y j ? y ? y j ? ?y} 通过仿射变换 ? ? ( x ? xi ) / ?x

? ? ( y ? yi ) / ?y

总可变成单位正方形 I ? I ? [0,1; 0,1] 如果在 Rij 上造出单元形状函数,就可得到试 探函数 uh (x)。而 Rij 上的形状函数可通过先在 I ? I 上造出形状函数,再通过仿射变化而得到。
123

在 I ? I 上构造形状函数,也采用Lagrange型 和Hermite型插值。
Lagrange型:根据若干插值节点处的函数 值决定插值函数。 Hermite型:根据若干插值节点处的函数值、 一阶偏导数乃至更高阶偏导数决定插值函数。

124

(1)Lagrange型公式 ① 双一次插值 ? ? 插值条件:给定 I ? I 顶点上的函数值 ? ? u00 , u01, u10 , u11 求:双线性函数 p11 (? ,? ) ? c0 ? c1? ? c2? ? c3?? 满足 p11 (0,0) ? u00 , p11 (0,1) ? u01, p11 (1,0) ? u10 , p11 (1,1) ? u11 设 p11(? ,? ) ? N00 (? ,? )u00 ? N01(? ,? )u11 ? N10 (? ,? )u10 ? N11(? ,? )u11 令
N 00 (0,0) ? 1, N 00 (0,1) ? N 00 (1,0) ? N 00 (1,1) ? 0 N 01 (0,1) ? 1, N 01 (0,0) ? N 01 (1,0) ? N 01 (1,1) ? 0 N10 (1,0) ? 1, N10 (0,0) ? N10 (0,1) ? N10 (1,1) ? 0 N11 (1,1) ? 1, N11 (0,0) ? N11 (0,1) ? N11 (1,0) ? 0

由 N ij (? ,? )为双线性函数,可求得 N ij (? ,? ) i, j ? 0,1

125

令 则

N 0 (? ) ? 1 ? ?

N1 (? ) ? ?

p11 (? ,? ) ? N 0 (? ) N 0 (? )u00 ? N 0 (? ) N1 (? )u11 ? N1 (? ) N 0 (? )u10 ? N1 (? ) N1 (? )u11

通过仿射变换消去?、? ,就得到 Rij上的 形状函数。把这些函数按单元叠加,即对所有 单元求和,就得到G上的试探函数。 实际计算时,并不消去中间变量?、? , 因为计算刚度矩阵元素(定积分)用?、?作 自变量更为方便。

126

② 双二次插值
插值条件:给定I?I上九个插 值节点(0,0)、(1/2,0)、(1,0)、 (0,1/2)、(1/2,1/2)、(1,1/2)、 (0,1)、(1/2,1)、(1,1)的函数值。 求:双二次函数
p22 (? ,? ) ? c0 ? c1? ? c2? ? c3? ? c4?? ? c5?
2 2

? ? ?

? ? ?

? ? ?

? c6? 2? ? c7?? 2 ? c8? 2? 2

满足

i j p22 ( , ) ? u i j , 2 2 2 2

i, j ? 0, 1, 2

127

设 p22 (? ,? ) ?

i , j ?0

?N

2

i 2

(? )N j (? )u i
2

j 2 2 ,

1 N 0 (0) ? 1, N 0 ( ) ? N 0 (1) ? 0 令 2 1 N 1 ( ) ? 1, N 1 (0) ? N 1 (1) ? 0 2 2 2 2 1 N1 (1) ? 1, N1 (0) ? N1 ( ) ? 0 2

由 N i (? ) 为二次函数,可求得 N 0 (? ) ? (2? ? 1)(? ? 1) N 1 (? ) ? 4? (1 ? ? ) N1 (? ) ? ? (2? ? 1) 故
p22 (? ,? ) ?
i , j ?0

?N

2

2
i 2

(? ) N j (? )u i
2

j 2 2 ,

通过仿射变换消去?、? ,就得到 Rij上的形状 128 函数。

③ 双三次插值

插值条件:给定I?I上十六个 插值节点(见图) 。
求:双三次函数

? ? ? ?

? ? ? ?

? ? ? ?

? ? ? ?

p33 (? ,? ) ? c0 ? c1? ? c2? ? c3? 2 ? c4?? ? c5? 2 ? c6? 3 ? c7? 2? ? c8?? 2 ? c9? 3 ? c10? 3? ? c11? 2? 2 ? c12??3 ? c13? 3? 2 ? c14? 2? 3 ? c15? 3? 3

满足 设

i j p33 ( , ) ? u i j , 2 2 3 3

i, j ? 0, 1, 2, 3

p33 (? ,? ) ?

i , j ?0

?N

3

i 3

(? ) N j (? )u i
3

j 3 3 ,

129



1 2 N 0 (0) ? 1, N 0 ( ) ? N 0 ( ) ? N 0 (1) ? 0 3 3 1 2 N 1 ( ) ? 1, N 1 (0) ? N 1 ( ) ? N 1 (1) ? 0 3 3 3 3 3 3 2 1 N 2 ( ) ? 1, N 2 (0) ? N 2 ( ) ? N 2 (1) ? 0 3 3 3 3 3 3 1 2 N1 (1) ? 1, N1 (0) ? N1 ( ) ? N1 ( ) ? 0 3 3

由 N i (? ) 为三次函数,可求得
9 1 2 27 2 N 0 (? ) ? ? (? ? 1)(? ? )(? ? ), N 1 (? ) ? ? (? ? 1)(? ? ), 2 3 3 2 3 3 27 1 9 1 2 N 2 (? ) ? ? ? (? ? 1)(? ? ), N1 (? ) ? ? (? ? )(? ? ), 2 3 2 3 3 3



p33 (? ,? ) ?

i , j ?0

?N

3

i 3

(? ) N j (? )u i
3

j 3 3 ,

130

(2) Hermite型公式 Lagrange型公式中不出现导数,这样的试探 函数只属于 C 0。为了得到属于 C 1 的试探函数,需 要Hermite型插值公式。 双三次多项式含有十六项: 1, x, y, x 2 , xy, y 2 , x 3 , x 2 y, xy2 , y 3 , x 3 y, x 2 y 2 , xy3 ,
x y ,x y ,x y
3 2 2 3 3 3

可以在四个顶点分别给定函数值、两个一 阶偏导数的值和二阶混合偏导数的值(共十六条 件),确定一个双三次多项式的十六个系数。

简单且常用的是不完全的双三次多项式插值。 它去掉双三次多项式中的 x 2 y 2 , x3 y 2 , x 2 y 3 , x 3 y 3 项。
131

插值条件:给定I?I上四个插值节点。

?

?

求:不完全双三次函数 ? 2 ? 2 2 2 ~( x, y) ? a ? a x ? a y ? a x ? a xy ? a y ? a x y ? a xy p 0 1 2 3 4 5 6 7
? a8 x 3 ? a9 y 3 ? a10 x 3 y ? a11xy3

满足
p 四个顶点处 ~ 的函数值等于 u 在该点的函数值; ?~ p 四个顶点处 的值等于 u x 在该点的值; ?x ?~ p 四个顶点处 的值等于u y 在该点的值。 ?y (y ? yj) ( x ? xi ) ?? 根据仿射变换 ? ? ?x ?y

则可将原插值问题转化为I?I上的插值问题。
132

插值条件:给定I?I上四个插值节点 (0,0)、(1,0)、 (0,1)、 (1,1) 。
求:不完全双三次函数 p(? ,? ) ? d 0 ? d1? ? d 2? ? d3? 2 ? d 4?? ? d5? 2 ? d 6? 2? ? d 7?? 2
? d8? 3 ? d9? 3 ? d10? 3? ? d11??3

满足

四个顶点处 p 的函数值等于 u 在该点的函数值; ?p 四个顶点处 ?? 的值等于 u x 在该点的值乘以?x; ?p 四个顶点处 的值等于u y 在该点的值乘以?y。 ?? 类似于Lagrange型公式的构造,可以求得 Rij 上的形状函数。
133

4 二维问题的三角形元 用有限元求解二维椭圆边值问题时,应用最 广的是三角形元。

在三角形元的有限元方法 中,先将定解区域G化分为若 干个小三角形(称作单元)。 然后在每个单元上构造插值型 函数,并用分片函数(但整体 连续的函数)代替变分问题或 变分方程中所需求解的函数。

134

(1)三角剖分 将定解区域化分成若干个小三角形单元 ? i 时应注意: 错误 ① 为了方便构造插值型函数, 要求每个单元的顶点是相邻 单元的顶点。 应避免 ② 为了保证有限元解有较好的 精度,每个单元中应尽量避免出现大的钝角。 ③ 为了保证有限元解的精确度和收敛性,并 避免其离散后代数方程组系数矩阵的病态性, 网格剖分中疏密的过渡不要太陡。
④ 单元顶点的编号顺序可以任意,但节点编 号顺序将影响有限元方程组系数矩阵的结构 (带宽)。

135

(2)面积坐标及有关公式 在三角形单元上构造插值型函数,并不简单 类同于矩形单元。 ① 面积坐标 考虑一个面积为S的三角形单元,其顶点 按反时针顺序记为i, j, k。在此单元内部任取一 点p(x,y),连接p和三个顶点 ,此单元则被分成 三个小三角形它们的面积记为 S i , S j 和 S k。
Sj Si Sk 记 Li ? , L j ? , Lk ? S S S 单元内任一点p(x,y)的位置与 三维数( Li , L j , Lk ) 一一对应, i 称 ( Li , L j , Lk ) 为面积坐标。
k

Sj

P(x,y)

Sk

Si
j
136

② 面积坐标与直角坐标的关系 面积坐标与坐标系无关,这是采用面积坐标的优点
面积坐标和直角坐标之间的关系:
2Si 1 ? ? Li ? 2 S ? 2 S [( x j yk ? y j xk ) ? ( y j ? yk ) x ? ( xk ? x j ) y ] ? 2S j 1 ? ? [( xk yi ? yk xi ) ? ( yk ? yi ) x ? ( xi ? xk ) y ] ?L j ? 2S 2S ? 2S k 1 ? ? Lk ? 2 S ? 2S [( xi y j ? yi x j ) ? ( yi ? y j ) x ? ( x j ? xi ) y ] ?

? x ? xi Li ? x j L j ? xk Lk ? ? ? y ? yi Li ? y j L j ? yk Lk ?

面积坐标有性质:Lm ( xn , yn ) ? ? mn

?1 ?? ?0

m?n m?n
137

③ 任意三角形到标准等腰直角三角形的变换 由于面积坐标满足 Li ? L j ? Lk ? 1 ,将其代入 ? x ? xi Li ? x j L j ? xk Lk ? Y ( xk , yk ) ? ? y ? yi Li ? y j L j ? yk Lk ? 得: (x j , y j ) ? x ? xi (1 ? L j ? Lk ) ? x j L j ? xk Lk ? ( xi , yi ) ? X o ? y ? yi (1 ? L j ? Lk ) ? y j L j ? yk Lk ? Lk 将 L j , Lk看作是某一平面 的坐标,则上式表明了一种 (0,1) k 变换。可以证明它把X-Y坐 标系的任意三角形单元映射 i j L 为 L j ? Lk 坐标系中的标准等 o (1,0) j 腰直角三角形单元。且 ( xi , yi ) ? (0,0) ( x j , y j ) ? (1,0) ( xk , yk ) ? (0,1) 138

这个变换的Jacobi行列式
?( x, y ) J? ? 2S ? 0 ?( L j , Lk ) 即该变换是仿射变换。 该变换除了能将三角单元仿射变换为标准 三角单元,还能将三角单元上的插值型函数变 换为标准三角单元上的同类型函数。因此,采 用面积坐标可使计算工作简单化、标准化。

另外,利用面积坐标表示的齐次多项式在 ?(i,j,k)上的积分也非常容易计算。即 p!q!r! p q r ??j ,kLi L j Lk dxdy ? 2S ( p ? q ? r ? 2)! ? (i , ) 其中p, q, r是任意非负整数。
139

(3)Lagrange型公式 两个变量x, y的高次多项 式可用Pascal三角形表示: 一般在三角形元 ?(1,2,3) 上,构造一个m次完全多项式 m
x x2 x3 x2 y

1 y xy xy2 y2 y3

p m ( x, y ) ?

i ? j ?0

?c x y
i ij

j

x 4 x3 y x 2 y 2 xy3 y 4 ??????????

逼近u(x,y)时,该多项式

pm ( x, y)

具有的项数为

1 1 ? 2 ? 3 ? ? ? (m ? 1) ? (m ? 1)( m ? 2) 2 1 故Lagrange型插值需要 (m ? 1)( m ? 2)个节点的 2 函数插值条件。
140

① 一次多项式 p1 ( x, y) 是一次多项式,插 值节点数是3。取?(1,2,3)的三个 顶点为插值节点,运用待定系 数法,易得 1 p1 ( x, y) ? L1u1 ? L2u2 ? L3u3 ② 二次多项式 p2 ( x, y)是二次多项式,插 值节点数是6。取?(1,2,3)的三个 顶点及三边中点为插值节点。设 运用待定系数法,可得

3

2

3 5?
? ?4

2

1

6

2 p2 ( x, y) ? c1L1 ? c2 L2 ? c3 L2 ? c4 L1L2 ? c5 L2 L3 ? c6 L3 L1 2 3

p2 ( x, y ) ? ? Li (2 Li ? 1)ui ? 4( L2 L3u4 ?L1 L3u5 ? L1 L2u6 )
i ?1

3

141

(4) Hermite 型公式 ① 二次多项式 在 ?(1,2,3)的顶点上取指定的 u1 , u2 , u3, ? ?u ? ? ?u ? ? ?u ? 在三边中点处的法向导数取指定的 ? ? , ? ? , ? ?

? ?n ? 4 ? ?n ?5 ? ?n ?6 2 设 p2 ( x, y) ? c1L1 ? c2 L2 ? c3 L2 ? c4 L1L2 ? c5 L1L3 ? c6 L2 L3 2 3

根据插值条件可得
3

3 5?
?4

? ? ( 2) ? ?u ? ? ( 2) ? p2 ( x, y ) ? ? ?? i ( x, y )ui ? ? i ( x, y )? ? ? 2 ? ?n ?3?i ? ? i ?1 ? 6 ? 1 ? ( 2) 1 1 1 ? 2 ?? i ( x, y ) ? Li ? (1 ? ? j ) Li L j ? (1 ? ? k ) Li Lk ? (2 ? ? j ? ? k ) L j Lk 2 2 2 ? 2s ? ( 2) ?? i ( x, y ) ? ? l Li ( L j ? Lk ) jk ? ?

? 142 其中 l jk 是 jk 边长度,i , ? j , ? k是顶点i, j, k的偏心率。

② 三次多项式 由三个顶点上的函数值和两个一阶偏导数, 加上在三角形形心上的函数值(共十个条件)确 3 定三次插值函数。
3 2 2 p3 ( x, y) ? c1L1 ? c2 L3 ? c3 L3 ? c4 L1 L2 ? c5 L1 L3 设 2 3 ? c6 L2 L1 ? c7 L2 L3 ? c8 L2 L1 ? c9 L3 L2 ? c10 L1L2 L3 2 2 3 3

?

2

可求得三次多项式为:
3 ( 3) 0

1

p3 ( x, y ) ? ? ( x, y )u0 ? ? ? i(3) ( x, y )ui ? ? i(3) ( x, y )(u x ) i ? ? i( 3) ( x, y )(u y ) i
i ?1 ( ? 03) ( x, y ) ? 27 L1 L2 L3

?

?

? i(3) ( x, y ) ? L3 ? 3L2 ( L j ? Lk ) ? 7 Li L j Lk i i ? i(3) ( x, y ) ? ( x j ? xi )( L2 L j ? Li L j Lk ) ? ( xk ? xi )( L2 Lk ? Li L j Lk ) i i ? i(3) ( x, y ) ? ( y j ? yi )( L2 L j ? Li L j Lk ) ? ( yk ? yi )( L2 Lk ? Li L j Lk ) i i
143

只利用插值节点处函数值得信息构造出来 的插值函数称作是Lagrange型插值函数。 除利用插值节点处的函数值信息,还利用 插值节点处的导数信息作为插值条件构造出来 的插值函数称作是Hermite型插值函数。 就整体而言,Hermite型插值函数的光滑程 度会有所提高。 一般情况下,高次元有限元解的精度比线性 有限元解的精度高。但在形成和求解有限元方程 组的过程中所花费的工作量也随之猛增。通常除 特殊需要,一般不必采用高次元,而是采用线性 元、网格细分、外推等措施。 除了三角形元和矩形元,还可考虑任意(甚至 是曲边)的四边形单元及曲边三角形单元。 144

5 用有限元方法解题主要步骤 ① 把边值问题化成等价的变分问题(或建立等价 的变分方程)。 ② 对求解域作网格。 ③ 构造基函数(或单元形状函数)。 ④ 形成有限元方程。

(a) 从变分原理及单元形状函数出发,先形成单元 刚度矩阵,再由单元刚度矩阵叠加成总刚矩阵。
(b) 从虚功原理及节点基函数出发形成有限元方程。

⑤ 求解线性代数方程组。
145

146


相关文章:
偏微分方程数值解第一次大作业
偏微分方程数值解第一次大作业_理学_高等教育_教育专区。清华大学偏微分方程数值解第一次大作业 偏微分方程数值解——第一次大作业 作业一:用有限差分法计算: ...
偏微分方程数值解法的MATLAB源码
[原创 偏微分方程数值解法的 MATLAB 源码【更新完毕】 原创]偏微分方程数值解法的 源码【更新完毕】 原创 说明:由于偏微分的程序都比较长,比其他的算法稍复杂一些...
偏微分方程数值解法
偏微分方程数值解法_数学_自然科学_专业资料。《偏微分方程数值解法》 课程设计 题姓学专班学 目: 六点对称差分格式解热传导方程的初边值 问题 名: 院: 业:...
第十章 偏微分方程数值解法
第十章 偏微分方程数值解法 偏微分方程问题,其求解十分困难。除少数特殊情况外,绝 大多数情况均难以求出精确解。因此,近似解法就显得更为重要。本章仅 介绍求解...
偏微分方程 数值解
中南林业科技大学 偏微分方程数值解法学生姓名: 学 号: 学 院: 专业年级: 设计题目:一维热传导方程几种差分格式 2012 年 07 月 一维热传导方程 ?u ? 2u ?...
偏微分方程数值解试题参考答案
偏微分方程数值解试题参考答案 - 偏微分方程数值解 1 ( Ax , x) ? (b, x) ( x ? R n ) ,证明下 2 一(10 分) 、设矩阵 A 对称正定,定义 J (...
偏微分方程数值解实验报告
偏微分方程数值解实验报告_数学_自然科学_专业资料。偏微分方程数值解实验报告 一、题目: 1 、用有限元方法求下列边值问题的数值解: ? ?2 ?x y = 2sin( ...
几种常见的偏微分方程数值求解问题
几种常见的偏微分方程数值求解问题_数学_自然科学_专业资料。一. 椭圆型问题 1.1 单位圆盘的泊松方程 泊松方程是最简单的椭圆型 PDE 问题。 该问题的公式为 该...
偏微分方程数值解总复习
偏微分方程数值解总复习 - 偏微分方程数值解总复习 一、考虑一维经典的初值问题: ? du ? = f (t , u ) ? dt ? u (0) = u 0 ? t ∈ (0, T)...
偏微分方程数值解
数学与计算科学学院 实验报告 实验项目名称 用 Eular 方法求解一阶常微分方程数值解 所属课程名称 实验类型实验日期 偏微分方程数值解 验证性 2015-3-26 班学姓...
更多相关标签: