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

偏微分方程数值解

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

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

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

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

232121000

3432310001121(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 ττττ----⎧=+-++⎪=+-++⎪⎪=+-++⎨

⎪⎪⎪=+-++⎩ (1.2.5) 若记

()

T

k

N k k k u u u 121,,,-= u ，()()()()T N x x x 1

21,,,-=ϕϕϕϕ ，()()()()T N x f x f x f 121,,,-=τττ f 则显格式(1.2.4)可写成向量形式

10

,0,1,,1

k k k M φ

+⎧=+=-⎨=⎩u Au f u (1.2.6)

其中

⎪⎪⎪

⎪⎪⎪

⎭

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

r r

21002100210021

A

而对于向前差分格式，当网比12r ≤

时稳定，当1

2

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 计算得到该方程的解析解，数值解以及相对误差如下：