目录引言 (3)物理背景 (3)网格剖分 (4)向前Euler格式建立 (4)差分格式的求解 (6)收敛性与稳定性 (6)数值例子 (9)紧差分格式建立 (12)差分格式求解 (14)数值例子 (15)总结 (19)参考文献 (20)附录 (21)1 引言本文考虑的一维非齐次热传导方程的定解问题：22(,),0,0,u ua 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 以及ρ均为常数.令2ka c ρ=,上式方程化为 22222222,t u u u u a a u xy z ⎛⎫∂∂∂=++=∆ ⎪∂∂∂⎝⎭若考虑物体内有热源,其热源密度函数为(),,F F x y z =,则有热源的热传导方程为()2,,,,t u a u f x y z t =∆+其中Ff 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ωω∞≤≤=，21()mi j h ωω==∑，2111()mi i j h hωωω-=-=∑，2211ωωω=+，分别称为无穷范式（一直范式）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 ux t a x t f x t i m k n t x ∂∂-=≤≤-≤≤-∂∂ (3.2)将224112241(,)[(,)2(,)(,)](,)12i k i k i k i k ik k u h ux t u x t u x t u x t t x h x ξ-+∂∂=-+-∂∂242114(,),12k xiik k i ik i h uU t x x x δξξ-+∂=-<<∂和2121(,)[(,)(,)](,)2i k i k i k i ik u u x t u x t u x t x t t τητ+∂∂=--∂∂212(,),2ktii ik k ik k uDU x t t tτηη+∂=-<<∂代入(3.2),得到224224(,)(,)(,),212kk t ixii k i ik ik k uah uDU 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(,)(,),212iki ik ik k uah uRx t t x τηξ∂∂=-∂∂ (3.4)并用k i u 代替ki U ，得到如下差分格式2(,),k kt 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)ik R 为差分格式的局部截断误差。

记241240101001max{max ,max },212x x t T t Tu a uc t x ≤≤≤≤≤≤≤≤∂∂=∂∂ (3.6-1)则有(1)21(),ik R c h τ≤+ 11,0 1.i m k n ≤≤-≤≤- (3.6-2)4.1.2 差分格式的求解记2r 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 kk 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.13 收敛性与稳定性收敛性设{(,)|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 mu 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 km e e k n ==≤≤当1/2r ≤时，有1(1)2211110max ()(),1.k k ik i m i eR c k h c T h k n ττττ-∞≤≤-=≤≤⋅+≤⋅+≤≤∑证明完毕.稳定性如果在应用差分格式(3.5-1)~(3.5-3)时，计算右端函数(,)i k f x t 有误差k i g ，计算初值()i x ϕ有误差i ψ，则实际得到的是如下差分方程的解。

2(,),k k kt i x i i k i D v a v f x t g δ-=+ 11,i m ≤≤- 01,k n ≤≤- 0(),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 ki i i v u ε=- 0,i m ≤≤ 0,k n ≤≤ 将(3.5-1)~(3.5-3)与(3.8)相减，可得摄动方程组2,k k kt 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 ll g εψτ-∞∞=∞≤+∑ 1.k n ≤≤ (3.10)上式说明当ψ∞和1n ll g τ-=∞∑很小时，误差1max kk nε∞≤≤也很小。

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

上述结果可叙述如下。

当1/2r ≤时，差分格式(3.5-1)~(3.5-3)关于初值和右端项在下述意义下是稳定的：设{}0,0k iui 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 ≤≤ 的解，则有100,k k l l uuf τ-∞∞=∞≤+∑ 1.k n ≤≤下面来考虑1/2r >的情况。

此时必存在0,h 当h h ≤时，2(1)1sin .22m h r π-⎛⎫>⎪⎝⎭于是2(1)14sin 1.2m h r π-⎛⎫-<-⎪⎝⎭设0,ki 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,km ε= 1.k n ≤≤可以验证其解为20(1)14sin ,2kk i i m h r πεε⎡-⎤⎛⎫=-⋅⎪⎢⎥⎝⎭⎣⎦ 0,i m ≤≤ 0.k n ≤≤由此易知20(1)14sin ,2kkm h r πεε-⎛⎫=-⋅ ⎪⎝⎭2(1)14sin .2kk m h r πεε∞∞-⎛⎫=-⋅ ⎪⎝⎭由于当k →∞时，2(1)14sin ,2km h r π-⎛⎫-→∞ ⎪⎝⎭所以不论初始误差多么小，均会使解有较大的误差。