一维热传导方程的前向、紧差分格式

中南林业科技大学
本科课程论文
学院：理学院
专业年级：09信息与计算科学一班
课程：偏微分方程数值解法
论文题目：一维热传导方程的前向Euler和紧差分格式指导教师：陈红斌
2012年7月
学生姓名：唐黎学号：20093936分工：程序编写，数值例子
学生姓名：何雄飞学号：20093925分工：格式建立，资料收集
学生姓名：汪霄学号：20093938分工：文档编辑，资料整理
学生姓名：毛博伟学号：20093931分工：公式编辑，查找资料
学生姓名：倪新东学号：20093932分工：数据分析，查找资料
学生姓名：何凯明学号：20093924分工：数据分析，查找资料
目录
1引言 (1)
2物理背景 (1)
3网格剖分 (2)
4．1.1向前Euler格式建立 (2)
4.1.2差分格式的求解 (4)
4.1.3收敛性与稳定性 (4)
4.1.4 数值例子 (7)
4.2.1紧差分格式建立 (10)
4.2.2差分格式求解 (12)
4.2.3数值例子 (13)
总结 (17)
参考文献 (18)
附录 (19)
1 引言
本文考虑的一维非齐次热传导方程的定解问题：
22(,),0,0,u u
a f x t x l t T t x
∂∂-=<<<≤∂∂
(,0)(),0,u x x x l φ=≤≤ (0,)(),
(1,)(),
0.u t t u t t t T αβ==<≤
其中a 为正常数，(,),(),(),()f x t x t t ϕαβ为已知函数，(0)(0),(1)(0).ϕαϕβ==
目前常用的求解热传导方程的差分格式有前向Euler 差分格式、向后Euler 差分格式、Crank-Nicolson 格式、Richardson 格式[1,2,3]．本文将给出前向Euler 格式和紧差分格式，并给出其截断误差和数值例子．
2 物理背景
热传导是由于物体内部温度分布不均匀,热量要从物体内温度较高的点流向温度较低的点处.以函数(),,,u x y z t 表示物体在t 时刻,(),M M x y =处的温度,并假设
(),,u x y z 关于,,x y z 具有二阶连续偏导数,关于t 具有一阶连续偏导数.(),,k k x y z =是物体在(),,M x y z 处的热传导系数,取正值.设物体的比热容为(),,c c x y z =,密度为
(),,x y z ρ.根据Fourier 热传导定律,热量守恒定律以及Gauss 公式得
,u u u u c kx k k t x x y y z z ρ
⎛⎫∂∂∂∂∂∂∂⎛⎫⎛⎫
=++ ⎪ ⎪ ⎪∂∂∂∂∂∂∂⎝⎭⎝⎭
⎝⎭ 如果物体是均匀的,此时,k c 以及ρ均为常数.令2
k
a c ρ
=
,上式方程化为 2222
2222,t u u u u a a u x y
z ⎛⎫
∂∂∂=++=∆ ⎪∂∂∂⎝⎭
若考虑物体内有热源,其热源密度函数为(),,F F x y z =,则有热源的热传导方程为
()2,,,,t u a u f x y z t =∆+
其中F f c ρ
=
. 3 网格剖分
取空间步长N l h /=和时间步长M T /=τ，其中M N ,都是正整数．用两族平行直线),1,0(N j jh x j ==和),,1,0(M k k t k ==τ将矩形域}0,0{T t l x G ≤≤≤≤=分割成矩形网格，网格节点为),(k j t x ．记),(k j k j t x u u =．以h G 表示网格内点集合，即位于开矩形G 的网点集合；h G 表示所有位于闭矩形的网点集合；h h h G G -=Γ是网格界点集合.
引进如下记号:
0max i i m
ω
ω∞
≤≤=，
ω=
1ω=，
1ω=
分别称为无穷范式（一直范式）2范数（2L 范数，平均范数），差商的2范数（差商的2L 范式）和1H 范式.
4.1.1向前Euler 格式建立
定义h τΩ上的网格函数
{0,0},
k i U U i m k n =≤≤≤≤
其中
(,),k i i k U u x t = 0,i m ≤≤ 0 1.k n ≤≤-
在结点处考虑微分方程(3.1-1),有
22(,)(,)(,),
11,0 1.
i k i k i k u u
x t a x t f x t i m k n t
x ∂∂-=≤≤-≤≤-∂∂ (3.2)
将
22411224
1(,)[(,)2(,)(,)](,)12i k i k i k i k ik k u h u
x t u x t u x t u x t t x h x ξ-+∂∂=-+-∂∂
242
11
4
(,),12k x
i
ik k i ik i h u U t x x x δξξ-+∂=-<<∂
和
212
1(,)[(,)(,)](,)2i k i k i k i ik u u x t u x t u x t x t t τητ+∂∂=--∂∂
21
2
(,),2k
t
i
i ik k ik k u
DU x t t t τηη+∂=-<<∂
代入(3.2),得到
224224
(,)(,)(,),212k
k t i
x
i
i k i ik ik k u
ah u
DU a U f x t x t t x τδηξ∂∂-=+-∂∂
11,0 1.i m k n ≤≤-≤≤- (3.3-1) 注意到初边值条件(3.1-2)和(3.1-3),有
0(,0)(),0,i i i U u x x i m ϕ==≤≤ (3.3-2)
0(),
k k U t α=
(),
1.
k
m k U t k n β=≤≤ (3.3-3)
在(3.3-1)~(3.3-3)中略去小量项
224(1)
24(,)(,),212ik
i ik ik k u
ah u
R
x t t x τηξ∂∂=-∂∂ (3.4)
并用k i u 代替k i U ，得到如下差分格式
2(,),k k
t i x i i k D u a u f x t δ-= 11,0 1.i m k n ≤≤-≤≤- (3.5-1)
0(),
i i u x ϕ= 0,i m ≤≤ (3.5-2)
0(),
k k u t α= 1.k n ≤≤ (3.5-3)
称(1)k R 为差分格式的局部截断误差。

记
241240101001max{max ,max },
212x x t T t T
u a u
c t x ≤≤≤≤≤≤≤≤∂∂=∂∂ (3.6-1)
则有
(1)
21(),
ik R c h τ≤+ 11,0 1.i m k n ≤≤-≤≤- (3.6-2)
4.1.2 差分格式的求解
记2=/r a h τ，称r 为步长比。

差分格式(3.5-1)可写为
111(12)()(,),
k k k k i i i i i k u r u r u u f x t τ+-+=-+++ 01,0 1.i m k n ≤≤-≤≤-
上式表明第k+1层上的值由第k 层上的值显示表示出来。

若已知第k 层的值
{|0}k i u i m ≤≤，则由上式就可直接得到第k+1层上的值+1{|0}k i u i m ≤≤。

有时也称
(3.5-1)为古典显格式。

可把古典显格式写成矩阵形式
1111012221222111112(,)12(,)12(,)12(,)k k k
k k k k k k m m m k k k k m m m k m r r u u f x t ru r
r r u u f x t r
r r u u f x t r r u u f x t ru ττττ+++---+---⎛⎫-⎛⎫⎛⎫+⎛⎫ ⎪ ⎪
⎪
⎪- ⎪ ⎪ ⎪ ⎪ ⎪
⎪ ⎪ ⎪=+ ⎪ ⎪
⎪ ⎪- ⎪ ⎪
⎪ ⎪ ⎪ ⎪ ⎪ ⎪-+⎝
⎭⎝⎭⎝⎭⎝
⎭
4.1.3 收敛性与稳定性
收敛性
设{(,)|0,0}i k u x t i m k n ≤≤≤≤为定解问题（3.1-1）~（3.1-3）的解，
{|0,0}k i u i m k n ≤≤≤≤为差分格式（）（）的解，则当1/2r ≤时，有
210max (,)(),
0k i k i i m
u x t u c T h k n τ≤≤-≤+≤≤，
其中1c 由（3.6-1）定义
证明
记
(,),
0,0.k k i i k i e u x t u i m k n =-≤≤≤≤
将（3.3-1）~（3.6-1）分别于（3.5-1）~（3.5-3）相减，得到误差方程
2(1)
,
11,01,k k t i x i ik D e a e R i m k n δ-=≤≤-≤≤-
00,
0,i e i m =≤≤
00,
0,
1.k k
m e e k n ==≤≤
当1/2r ≤时，有
1
(1)
2211110
max ()(),1.k k ik i m i e
R c k h c T h k n ττττ-∞
≤≤-=≤≤⋅+≤⋅+≤≤∑
证明完毕.
稳定性
如果在应用差分格式(3.5-1)~(3.5-3)时，计算右端函数(x ,t )i k f 有误差g k i ，计算初值(x )i ϕ有误差i ψ，则实际得到的是如下差分方程的解。

2(,),
k k k t i x i i k i D v a v f x t g δ-=+ 11,i m ≤≤- 01,k n ≤≤-
(),i i i v x ϕψ=+ 0,i m ≤≤ (3.8) 0(),k
k v t α= (),k m k v t β= 1.k n ≤≤
令
,k k k
i i i v u ε=- 0,i m ≤≤ 0,k n ≤≤ 将(3.5-1)~(3.5-3)与(3.8)相减，可得摄动方程组
2,k k k t i x i i D a g εδε-= 11,i m ≤≤- 01,k n ≤≤-
0,
i i εψ= 0,i m ≤≤ (3.9)
00,k ε= 0,k m ε= 1.k n ≤≤
当1/2r ≤时，有
1
,k k l
l g ε
ψ
τ-∞
∞
=∞
≤+∑ 1.k n ≤≤ (3.10)
上式说明当ψ
∞
和1
n i
i g τ-=∞
∑很小时，误差1max k
k n
ε∞
≤≤也很小。

摄动方程组(3.9)和差分方程(3.5-1)~(3.5-3)的形式完全一样。

上述结果可叙述如下。

当1/2r ≤时，差分格式(3.5-1)~(3.5-3)关于初值和右端项在下述意义下是稳定的：设{|0,0}k i u i m k n ≤≤≤≤为差分方程组 2,k k k t i x i i D u a u f δ-= 11,01,i m k n ≤≤-≤≤-
0,
i i u ψ= 0,i m ≤≤
00,0,
k k
m u u == 1k n ≤≤
的解，则有
1
00
,
k k l
l u
u
f τ-∞
∞
=∞
≤+∑ 1.k n ≤≤
下面来考虑>1/2r 的情况。

此时必存在0h 当0h h ≤时，
2(1)1
sin .22m h r π
-⎛⎫>
⎪⎝⎭
于是
2(1)14sin 1.2m h r π
-⎛⎫
-<-
⎪⎝⎭
设
0,k
i g = 11,01,i m k n ≤≤-≤≤- ()sin (1),
i i m x ψεπ=- 0,i m ≤≤
则（3.9）为
20,
k k t i x i D a εδε-= 11,01,i m k n ≤≤-≤≤-
()0sin (1),
i i m x εεπ=- 0,i m ≤≤
00,k ε=
0,k m ε= 1.k n ≤≤ 可以验证其解为
20
(1)14sin ,
2k
k i i m h r πεε⎡-⎤⎛⎫=-⋅ ⎪⎢⎥⎝⎭⎣⎦ 0,i m ≤≤ 0.k n ≤≤
由此易知
20
(1)14sin ,
2k
k
m h r πεε-⎛⎫=-⋅ ⎪⎝⎭ 20
(1)14sin .
2k
k
m h r πεε∞∞
-⎛⎫=-⋅ ⎪
⎝⎭
由于当k →∞时，
2(1)14sin ,
2k
m h r π-⎛⎫
-→∞ ⎪⎝⎭
所以不论初始误差多么小，均会使解有较大的误差。

我们得到如下结论。

定理3.1.4 当>1/2r 时，差分格式(3.5-1)~(3.5-3)是不稳定的。

由定理3.1.3和定理3.1.4知，当步长比1/2r ≤时，差分格式(3.5-1)~(3.5-3)是不稳定的；当步长比>1/2r 时，差分格式(3.5-1)~(3.5-3)是不稳定的。

我们把这种稳定性称为条件稳定性。

实际计算时选取步长比必须满足1/2r ≤，即2/1/2a h τ≤
4.1.4 数值例子
应用向前Euler 格式(3.5-1)-(3.5-3)计算定解问题
220,01,01u u x t t x ∂∂-=<<<≤∂∂ (,0),
01x u x e x =≤≤
1(0,),(1,),01t t u t e u t e t +==<≤
上述定解问题的精确解为(,)x t u x t e +=.
部分节点处数值解、精确解和误差的绝对值(1/10,1/200)h τ==
取不同不长时数值解的最大误差(1/2)r =
==的误差曲面图（1/2
1/10,1/200
hτ
r=）
==的误差曲面图（1/2 hτ
1/20,1/800
r=）
1/10,1/200
|(,1)(,1)|h h u x u x ττ==-的误差曲线图
1/20,1/800
|(,1)(,1)|h h u x u x ττ==-的误差曲线
4.2.1紧差分格式建立
令
22,
u v x ∂=∂
则有
1[(,)].
u
v f x t a t ∂=-∂
定义网格函数
(,),k i i k U u x t = (,),k i i k V v x t = 0,i m ≤≤ 0.k n ≤≤
由Taylor 展开，有
224462
246(,)(,)(,)
12360k x
i
i k i k ik k u h u h u
U x t x t t x x x δξ∂∂∂=++∂∂∂
224626(,)(,)(,)
12360i k i k ik k h u h u
v x t x t t x x ξ∂∂=++∂∂
22446246
(,)(,)1212360k
k i x i ik k ik k h h v
h u V V t t x x
δηξ⎡⎤∂∂=+-+⎢⎥∂∂⎣⎦
6641166111(10)(,)(,),12360144k k k
i i i ik k ik k u u V V V t t h x x ξη-+⎡⎤∂∂=+++-⎢⎥∂∂⎣⎦
其中11,(,)ik ik i i x x ζη-+∈将上式中上标为k 和k+1的两个等式相加并除以2， 可得
()22112k k x i x i U U δδ++111
222111(10)12k k k i i i V V V +++-+=++
66466
11(,)(,)360144ik k ik k u u
t t h x x ξη⎡⎤∂∂+-⎢⎥∂∂⎣⎦ 111112221(,)10(,)(,)12i i i k k k v x t v x t v x t -++++⎡⎤=
++⎢⎥⎣⎦
2222112221(,)10(,)(,)812i ik i ik i ik v
v v
x x x t t t τθθθ-+⎡⎤∂∂∂+⋅⋅++⎢⎥∂∂∂⎣⎦ 66466
11(,)(,).360144ik k ik k u u t t h x x ξη⎡⎤∂∂+-⎢⎥∂∂⎣⎦
利用(3.42),得到
111122
22211111012k k k k x
i t i t i t i U U U U a δδδδ+
+++-+⎛
⎫=⋅++ ⎪⎝⎭
1111122211(,)10(,)(,)12i i i k k k f x t f x t f x t a -++++⎡
⎤-⋅++⎢⎥
⎣⎦
23331133311(,)10(,)(,)2412i ik i ik i ik u u u
x x x a t t t τθθθ-+⎡⎤∂∂∂-⋅⋅++⎢⎥∂∂∂⎣⎦
2444112222221(,)10(,)(,)812i ik i ik i ik u
u u
x x x x t x t x t τθθθ-+⎡⎤∂∂∂+⋅⋅++⎢⎥∂∂∂∂∂∂⎣⎦
66466
11(,)(,),360144ik k ik k u u t t h x x ξη⎡⎤∂∂+-⎢⎥∂∂⎣⎦
其中111,(,),,,(,)ik ik ik i i k ik k k x x t t t ζηθθ-++∈∈ 则有
1111
2222
2111(10)12
k k k k t i t i t i t i U U U a U δδδδ++++-+++- 2411111222
1
[(,)10(,)(,)]()12i i i k k k f x t f x t f x t O h τ-++++=
++++
11,01i m k n ≤≤-≤≤-
注意到初值条件，忽略小项24()O h τ+，并用k i u 代替k i U ，得到如下差分格式
1111
2
2222111(10)12
k k k k t i t i t i t i u u u a u δδδδ++++-+++- 11111222
1
[(,)10(,)(,)]12i i i k k k f x t f x t f x t -++++=
++
11,01i m k n ≤≤-≤≤-
4.2.2 差分格式求解
差分格式可写成
1111111511()()()1226122k k k i i i r u r u r u +++-+-+++- 1111511
()()()1226122k k k i i i r u r u r u -+=++-++ 11
1112
2
2
[(,)10(,)(,)]12
i i i k k k f x t
f x t
f x t
τ
-++
+
+
+
++
11,01i m k n ≤≤-≤≤-
将其写出矩阵形式
11121
211511612211511
1226
12211511122
6
122115122
6k k k m k m r r u r r r u u r r r u r r +++-+-⎛⎫+- ⎪
⎪⎛⎫ ⎪ ⎪
-+- ⎪ ⎪ ⎪ ⎪ ⎪ ⎪
⎪ ⎪-+- ⎪
⎪ ⎪
⎪⎝⎭
⎪-+ ⎪⎝⎭ 12215116122115111226
12211511122
6
122115122
6k k k m k m r r u r r r u u r r r u r r --⎛⎫
-+ ⎪
⎪⎛⎫ ⎪ ⎪
+-+ ⎪ ⎪ ⎪ ⎪= ⎪ ⎪
⎪ ⎪
+-+ ⎪ ⎪ ⎪
⎪⎝⎭
⎪+- ⎪⎝
⎭
1
00112212
212
1
1121111()()(,)122122(,)
(,)
1111()()(,)122122k k k k m k k k m m m k r u r u Af x t Af x t Af x t r u r u Af x t ττττ++
+-++-+⎛⎫+--+ ⎪ ⎪
⎪ ⎪ ⎪+ ⎪
⎪ ⎪ ⎪ ⎪+--+ ⎪⎝
⎭ 01k n ≤≤-
其中
1111112
222
1
(,)[(,)10(,)(,)].12i i i i k k k k Af x t
f x t f x t f x t τ-++
+++=++
在每一时间层上只要解一个三对角线性方程组.
4.2.3 数值例子
应用紧差分格式(3.47-1)-(3.47-3)计算定解问题
220,01,01u u
x t t x ∂∂-=<<<≤∂∂ (,0),
01x u x e x =≤≤
1(0,),(1,),01t t u t e u t e t +==<≤
上述定解问题的精确解为(,)x t u x t e +=. 最大误差
111(,)max (,)k
i k i i m k n
E h u x t u τ∞≤≤-≤≤=-
从表中可以看出，当空间步长缩小到原来的1/2，时间步长缩小到原来的1/4时，最大误差约缩小到原来的1/16.
部分节点处数值解、精确解和误差的绝对值（1/10,1/100h τ==）
取不同步长的数值解的最大误差
==的误差曲面图hτ
1/10,1/100
1/20,1/400h τ==的误差曲面图
1/10,1/100
|(,1)(,1)|h h u x u x ττ==-的误差曲线图
1/20,1/400|(,1)(,1)|h h u x u x ττ==-的误差曲线图
1/40,1/1600|(,1)(,1)|h h u x u x ττ==-的误差曲线图
总结：
本文采用差分格式和紧差分格式来求解抛物型方程.差分格式采用二层共三个点,条件稳定显格式,当12r >
时误差随着r 无限增长其稳定条件为网比12r ≤,因此我们给出了当12
r =时的最大误差.而紧差分格式却采用二层共六个点,无条件隐格式,在每一时间层上均只需解三对角方程组.
参考文献
[1] 孙志忠．偏微分方程数值解法．北京：科学出版社，2005.
[2] 李荣华．偏微分方程数值解法．北京：高等教育出版社，2005.
[3]
精品
附录A Euler程序代码
流程图：
clear;
n=input('n=');%空间剖分数
m=input('m=');%时间剖分数
x=(0:1/n:1);
t=(0:1/m:1);
r=n*n/m;%网比
l=2;
A=diag(ones(1,n-1)*(1-2*r))+diag(ones(1,n-2)*r,-1)+diag(ones(1,n-2)*r,1); %---------------------精确解求解---------------------
for i=1:n+1
for j=1:m+1
u(i,j)=exp(x(i)+t(j));
end
end
%----------------------------------------------------
%--------------------初值条件求解---------------------
for i=1:n+1
u1(i,1)=exp(x(i));
end
for j=1:m+1
u1(1,j)=exp(t(j));
end
for j=1:m+1
u1(n+1,j)=exp(x(n+1)+t(j));
end
%---------------------数值解求解----------------------
for j=1:m
u1(2:n,j+1)=A*u1(2:n,j)+[r*u1(1,j),zeros(1,n-3),r*u1(n+1,j)]'; end
%----------------------------------------------------
%for i=1:10
% M(i)=(u(6,i*20+1));
% N(i)=(u1(6,i*20+1));
%end
%abs(M-N)'
%----------------------误差图-------------------------
u-u1;
[x,t]=meshgrid(x,t);
surf(x,t,u'-u1');
grid on;
xlabel('x');
ylabel('t');
zlabel('|u(x,t)-u1(x,t)|');
%---------------------误差曲线-------------------------
%E=abs(u(:,2)-u1(:,2));
%plot(x,E,'-o');
%u(:,m+1)'-u1(:,m+1)';
%plot(x,u(:,m+1)'-u1(:,m+1)','-o');
%xlabel('x');
%ylabel('|u(x,1)-u1(x,1)|');
%-------------------最大误差求解-----------------------
%for i=2:n
% for j=2:m+1
% K(i)=max(u(i,j)-u1(i,j));
% end
%end
%K=max(K)
%-----------------------------------------------------
附录B紧差分格式程序代码
clear;
n=input('n=');%空间剖分
m=input('m=');%时间剖分
H=1/n;%空间步长
T=1/m;%时间步长
x=(0:H:1);
t=(0:T:1);
r=n*n/m;%网比
l=2;
A=diag(ones(1,n-1)*(5/6+r))+diag(ones(1,n-2)*(1/12-r/2),-1)+diag(ones(1,n-2)*(1/12-r/2),1);%k+1层的系数矩阵
B=diag(ones(1,n-1)*(5/6-r))+diag(ones(1,n-2)*(1/12+r/2),-1)+diag(ones(1,n-2)*(1/12+r/2),1);%k层的系数矩阵
%---------------------精确解求解---------------------
for i=1:n+1
for j=1:m+1
u(i,j)=exp(x(i)+t(j));
end
end
%----------------------------------------------------
%--------------------初值条件求解---------------------
for i=1:n+1
u1(i,1)=exp(x(i));%条件u(x,0)=exp(x)
end
for j=1:m+1
u1(1,j)=exp(t(j));%条件u(0,t)=exp(t)
end
for j=1:m+1
u1(n+1,j)=exp(x(n+1)+t(j));%条件u(1,t)=exp(1+t)
end
%-------------------数值解的计算-------------
for j=1:m
u1(2:n,j+1)=A\(B*u1(2:n,j)+[(1/12+r/2)*u(1,j)-(1/12-r/2)*u(1,j+1),zeros(1,n-3),(1/12+r/2)*u(n+1,j)-(1/1 2-r/2)*u(n+1,j+1)]');
end
%----------------------------------------------------
%for i=1:n
% N(i)=(u(6,i*10+1));
% M(i)=(u1(6,i*10+1));
%end
%abs(N-M)'
%-------------------误差曲面----------------------
%u-u1;
%[t,x]=meshgrid(t,x);
%surf(t,x,u1-u);
%grid on;
%xlabel('x');
%ylabel('t');
%zlabel('|u(x,t)-u1(x,t)|');
%------------------误差曲线---------------------
%E=abs(u(:,m+1)'-u1(:,m+1)');
%plot(x,E,'-o');
%xlabel('x');
%ylabel('|u(x,1)-u1(x,1)|');
%-----------------------------------------------------
%for i=1:n
% M(i)=(u(6,i*10+1))
% N(i)=(u1(6,i*20+1));
%end
%abs(M-N)'
%----------------------------------------------------
%-------------------最大误差求解-----------------------%for i=2:n
% for j=2:m+1
% K(i)=max(abs(u(i,j)-u1(i,j)));
% end
%end
%K=max(K)
%-----------------------------------------------------。