抛物形扩散方程的有限差分法及数值实例(完整资料).doc

此文档下载后即可编辑

偏微分方程数值解

所在学院：数学与统计学院

课题名称：抛物形扩散方程的有限差分法及数值实例

学生姓名：向聘

抛物形扩散方程的有限差分法及数值实例

1.1抛物型扩散方程

抛物型偏微分方程是一类重要的偏微分方程。考虑一维热传导方程：

22(),0u u

a f x t T t x

∂∂=+<≤∂∂

(1.1.1)

其中a 是常数，()f x 是给定的连续函数。按照初边值条件的不同给法，可将(1.1.1)的定解分为两类：

第一，初值问题（Cauchy 问题）：求足够光滑的函数()t x u ,，满足方程(1.1.1)和初始条件： ()()x x u ϕ=0,, ∞<<∞-x

(1.1.2)

第二，初边值问题（也称混合问题）：求足够光滑的函数()t x u ,，满足方程(1.1.1)和初始条件： ()()x x u ϕ=0,, 0x l <<

(1.1.3) 及边值条件

()()0,,0==t l u t u ， T t ≤≤0

(1.1.4)

假定()x f 和()x ϕ在相应的区域光滑，并且于()0,0，()0,l 两点满足相容条件,则上述问题有唯一的充分光滑的解。

1.2抛物线扩散方程的求解

下面考虑如下热传导方程

22()(0.)(,)0(,0)()u u

a f x t x u t u L t u x x ϕ⎧∂∂=+⎪∂∂⎪⎪

==⎨⎪=⎪⎪⎩

(1.2.1)

其中，0x l <<，T t ≤≤0,a (常数)是扩散系数。 取N

l

h =

为空间步长，M

T =

τ为时间步长，其中N ，M 是自然

数，用两族平行直线

jh

x x j ==, ()N j ,,1,0Λ=和

k t t k τ

==,

()M k ,,1,0Λ=将矩形域G {}T t l x ≤≤≤≤=0;0分割成矩形网格。其中

(),j

k

x t 表示网格节点；h

G

表示网格内点（位于开矩形G 中的网格

节点）的集合；h G 表示位于闭矩形G 中的网格节点的集合；h Γ表示h G -h G 网格边界点的集合。

k j u 表示定义在网点(),j k x t 处的待求近似解，

N j ≤≤0，M k ≤≤0。 现在对方程进行差分近似：

(一) 向前差分格式

=-+τ

k j

k j u u 111

2

2(())

k k k

j j j j j j u u u a

f f f x h +--++=

(1.2.2)

()

j j j x u ϕϕ==0，

k

u 0

=

k

N

u =0

(1.2.3) 计算后得：

111(12)k k k k j j j j j

u ru r u ru f τ++-=+-++

(1.2.4) 其中，2

a r h τ

=

，1,,1,0-=N j Λ，1,,1,0-=M k Λ。

显然，这是一个四点显示格式，每一层各个节点上的值是通过一个方程组求解到的。方程组如下：

1000

121011000

2321210003432310001121

(12)(12)(12)(12)N N N N N u ru r u ru f u ru r u ru f u ru r u ru f u ru r u ru f ττττ----⎧=+-++⎪=+-++⎪⎪=+-++⎨⎪⎪⎪=+-++⎩M

(1.2.5)

若记

()

T

k

N k k k u u u 1

21,,,-=Λu ，()()()()T N x x x 121,,,-=ϕϕϕϕΛ，()()()()T N x f x f x f 121,,,-=τττΛf

则显格式(1.2.4)可写成向量形式

10

,0,1,,1

k k k M φ

+⎧=+=-⎨=⎩L u Au f u

(1.2.6) 其中

⎪⎪⎪

⎪⎪⎪

⎭

⎫ ⎝⎛----=r r r r r r r r

r r

21002100210021ΛO M

O O O M O ΛA 而对于向前差分格式，当网比12r ≤时稳定，当12

r >时不稳定。这就意味着给定空间步长h 以后，时间步长τ必须足够小，才能保证稳定。

1.3抛物型热传导方程数值算例

对于(1.2.1)所描述的扩散方程，取1a =已知方程的精确解为

2

sin t u e x ππ-=

：

22(,0)sin()(0,)(1,)001,00.5

u u

t x u x x u t u t x t π⎧∂∂=⎪∂∂⎪⎪

=⎨⎪==⎪≤≤≤≤⎪⎩

(1.3.1)

设空间步长1/h M =，时间步长为0.5/N τ=，网格比2/r g h =。 向前格式：

111

2

2,1,...,1,1,...,k k k k k

j j

j j j u u u u u j M k N

h τ

++---+=

=-=

边值条件：

1100

101112

20,

,k k

k k k

k k u u u u u j u u h τ

++----+===， 1111112

2,

,k k k k k k k M M

M M M M M u u u u u j M u u h

τ

+++-+---+===. 初值条件：

(,0)sin ,0,1,...,j j u x x j M

π==

对时间和空间进行分割，令M=40，N=1600，通过Matlab 计算得到该方程的解析解，数值解以及相对误差如下：