苏佳园：二维热传导方程有限差分法的MATLAB实现

l所构成，其中xj?j?x?h,?x?h?J;tn?n?t?n?.

3.2插值函数的选择

选择不同的插值函数对偏微分方程进行估计，可得到不同的差分方程，进而稳定性和精度会有所不同。

用Taylor级数展开方法是最常用的方法，下面建立差分格式的同时引入一些基本概念及术语。

我们主要从对流方程的初值问题

?u??u?a?0,x?R,t?0,? (3.1) ?x??t??u(x,0)?g(x),x?R,和扩散方程的初值问题

?u?2u?a2,?t?xu(x,0)?g(x),(其中a?0)进行讨论。

x?R,t?0, (3.2) x?R.

假定偏微分方程初值问题的解u(x,y)是充分官话的，由Taylor级数展开有

u(xj,tn?1)?u(xj,tn)?un????(?t),?t?tj?u(xj,tn?1)?u(xj,tn?1)2??un???(?t),2?t?tj?u(xj?1,tn)?u(xj,tn)?(3.3) ?un???(h),h?xj?u(xj,tn)?u(xj?1,tn)??un???(h),?h?xj?u(xj?1,tn)?u(xj?1,tn)2?un???(h),?2h?xj?u(xj?1,tn)?2u(xj,tn)?u(xj?1,tn)n2?2u????(h).22h?xj?

????????????n其中???n或用，表示看括号内的函数在节点(xj,tn)处取的值。利用(3.3)表达式???jj中的第1式和第3式有

u(xj,tn?1)?u(xj,tn)?

?au(xj?1,tn)?u(xj,tn)h6

?u?un?[?a]j??(??h). ?t?x如果u(x,t)是满足偏微分方程(3.1)的光滑解，则

[?u?u?a]nj?0. ?t?x由此看一看出，偏微分方程(3.1)在u(xj,tn)处可以近似的用下面的方程来代替

?1nun?ujjnun?uj?1j?n?ah?0, (3.4)

j?0,?1,?2,?,n?0,1,2?.

其中uj为u(xj,tn)的近似值。(3.4)式称为逼近微分方程(3.1)的有限差分方程或简称差分方程。

差分方程再加上初始条件的离散型式就可以按时间逐层推进，算出各层的值。差分格式(3.4)和初始条件的离散形式结合在一起构成了一个差分格式。

3.3方程组的建立

将离散后的差分方程转化为方程组的形式，便于求解。

3.4方程组的求解

利用矩阵的解法求解方程组，再用MATLAB对矩阵求解方法进行程序化，以便对以后类似的方程进行求解。

隐式差分格式方程矩阵化后，得到的矩阵是严格的对角占优三对角矩阵，我们可以根据线性方程组的求解方法对其求解。其中这要应用的是追赶法，追赶法对于此类线性方程组的求解非常方便，用MATLAB对追赶法进行编程，就可以轻松实现矩阵的求解，进而解出差分方程的近似解。

7

苏佳园：二维热传导方程有限差分法的MATLAB实现

第4章 二维热传导方程

4.1 网格剖分

在区域D:?(x,y,t)0?x?X,0?y?Y,0?t?T?中，我们设二维热传导方程的初始值和边界条件如下：

2u?2u??u???t?a(2?2),0?x?X,0?y?Y,0?t?T,?x?y???u(x,y,0)?u0(x,y),0?x?X,0?y?Y??u(x0,y,t)?f(y,t),u(xX,y,t)?g(y,t),??u(x,y0,t)?h(x,t)，u(x,yY,t)?k(x,t)，其中a为正常数。

(4.1)

通过已知方程，建立一个关于时间和步长的函数，这样就把初始区域划分为一个网格图。 先将定义域

D:?(x,y,t)0?x?X,0?y?Y,0?t?T?

剖分为网格

Dh?{(xj,yh,tn)xj?j?x,j?0,1,?,J,J?x?X;yl?l?y,l?0,1?,J,J?y?Y;tn?n?t,n?0}.其中?t?

T为时间步长，?x?X,?y?Y分别为x轴和y轴的空间步长。 NMxMy4.2 稳定性分析

利用有限差分格式进行计算时是按时间层逐层推进的。那么计算第n?1层上的值时要用到第n层上计算出来的结果值，而计算第n层结果值时的舍入误差必然会影响到第n?1层的值。从而就要分析这种误差传播的情况。希望误差不至于越来越大，以至掩盖差分格式的解的面貌，这便是稳定性问题。

我们先考虑一维差分格式

?1nnnun?u?a?(u?ujjj?1j) (4.2)

的稳定性，其中??

?t为网格比，假设a?0差分格式从初层开始计算，当初始数据?x8

存在误差时考察这个误差在以后计算中的在传播情况。为方便起见，不考虑计算过程中的舍入误差。及确定初始数据误差绝对值为?，则差分格式在(xj,tn)处的误差为

nnmm?n(?1)n?m??Cm(1?a?)(?a?)m?0nn???Cm(1?a?)m(a?)m?n m?0?(1?a?)n?.于是，对于固定网格比?及a?0的情况，差分格式的解的误差随时间步长的步数

n的增加而增加。初始数据的误差将必定掩盖了差分格式的解的面貌，所以我们认为

差分格式(4.2)时不稳定的。

差分格式的稳定性不仅与差分格式本身有关，而且还与网格比的大小有关。差分格式的稳定性在差分方法的研究中具有特别的意义，我们再做进一步的叙述。

定义4.1[8] 为了度量误差及其他应用，引入范数

nhu00n???n2???(uj)h?, ?j????n12设uj有一个误差?j，则uj就有误差?j。如果存在一个正常数K,使得当

?t??t0,n?t?T时，一致的有

?n?K?0,

则称差分格式是稳定的。

差分格式一旦具有稳定性，就可以用差分格式计算出偏微分方程的近似解来。 一维热传导方程的各种类型的差分格式可以推广到二维热传导方程，利用向前差分格式

?1nun?ujjnnun?2u?uj?1jj?1?t对(4.1)式进行离散，引入记号

?h2

nnn?x2un?u?2u?ujlj?1,ljlj?1,l,nnn?y2un?u?2u?ujlj,l?1jlj,l?1. (4.3)

其中unjl为差分方程在节点(j,l,n)的计算值。

差分格式

9

苏佳园：二维热传导方程有限差分法的MATLAB实现

?1nun?ujljl

?t?a(?x2unjl?x2?2n?yujl?y2), (4.4)

利用Taylor级数展开易得差分格式(4.4)的截断误差为?(?t??x2??y2)。

t为方便稳定性的判断，设?x??y?h，令???2，为网格比。即改写为：

hn?1n2n2n ujl?ujl?a?(?xujl??yujl), (4.5)

用Fourier方法来分析(4.4)式的稳定性。令

nik1jhik2lhun?vee, jl把此式带入(4.5)式中有

vn?1n?1?2a?(cosk1h?1)?2a?(cosk2h?1)v,

??因此差分格式(4.5)的增长因子是

2k1h2k2hG(?t,k)?1?4a?(sin?sin),

22其中k?(k1,k2)。如果a??件是

14，则有G(?t,k)?1由此得出差分格式(4.4)的稳定性条

a??14.

容易看出在二维情况下采用这样的显式格式是不适合的，为此我们再转向考虑隐式格式。用一维向后差分格式

nn?1uj?uj?tnnnuj?1?2uj?uj?1?2h

的直接推广是

n?1unjl?ujl?t?a(?x2unjl?x2??y2unjl?y2),

(4.6)

此格式的截断误差仍为?(?t??x2??y2)，仍用Fourier方法来分析这个格式的稳定性，仿前可以得出其增长因子是：

G(?t,k)?1,khkh2221?4a?(sin1?sin)2210