因此有G(?t,k)?1，即差分格式（4.6）是绝对稳定的。

为了提高精度，对微分方程（4.1）也可以用Crank-Nicolson型差分格式，这也是一维问题的直接推广。其格式可写为

?1nun?ujljl2n?1n2n?1n?ta?x(ujl?ujl)?y(ujl?ujl)?[?], （4.7） 222?x?y这也是二阶精度格式，其增长因子是

2k1h21?2a?sin?2a?sin2G(?t,k)?2k1h21?2a?sin?2a?sin2k2h2,k2h 2因此，对任何?都有G(?t,k)?1，所以(4.7)式也是绝对稳定的。

现在考虑一下隐式格式（4.6）式和（4.7）式的求解方法。我们知道，在一位格式形成的方程组是系数矩阵为三对角矩阵的线性代数方程组，因此用追赶法很容易求解。而对于（4.6）式和（4.7）式导出的系数矩阵不是三对角矩阵，因此求解就不容易了。

我们对于显式格式和隐式格式的分析知道，在实际使用上都受到限制，因此构造每层计算量不大的绝对稳定的格式就成为一个具有现实意义并很有兴趣的问题。在一维中，隐式格式是绝对稳定并可用追赶法很容易求解。由此产生了下面将使用的交替方程隐式格式。它具有绝对稳定、容易求解和有相当精度的特点。

4.3 建立方程组

22?u?u我们在构造微分方程（4.1）的隐式格式和显式格式中，对2和2做了同样的

?x?y2?u处理，即同时在第n层或第n?1层取值。为了构造一维形式的隐式格式，对二阶导数2?x2?u用u在第n?1层上用未知的二阶中心差分来代替，而2则用u在第n层上用已知的二

?y阶中心差分来代替，这样得到的方程组在仅x方向是隐式的。比较容易求解，用追赶法就可以了。同理，为对称起见，在下一时间层上重复上述步骤，即又仅在y方向是隐式的，对x方向是显式的。这样相邻的两个时间层合并起来构成一个差分格式。故用

11

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

多次追赶法就可以解出ujl了。

我们解向后差分格式的方程。 令rx?n?ta?ta,r?，则（4.6）式变为: y?x2?y2n?1un?u?a?t(jljl?x2unjl?x2?2n?yujl

?y2) （4.8）

2n?rx?x2un?r?jlyyujl,（4.8）在x方向是隐式时，代入（4.3），变形为

nnnunjl?ry(un?2u?u)?2ruj,l?1jlj,l?1xjl

（4.9）化为矩阵形式：

?unjl?rx(unj?1,l?unj?1,l),0?0 （4.9）

2rx?ry）?ry0?1?（?1?（2rx?ry）?ry??ry?0?ry1?（2rx?ry）??????000??n??uj1???n?0?0??uj2?????ry?0????

??????????n?u???ry1?（2rx?ry）????j,m?1??nnn??n?1u?r(u?u)?ru?j1xj?1,1j?1,1yj0?n?1nn??u?r(u?u)??j2xj?1,1j?1,1 ??. （4.10） ????????n?1nnn??uj,m?1?rx(uj?1,1?uj?1,1)?ryujm???（4.8）在y方向是隐式时，代入（4.3），可变形为

nnnnun?r(u?2u?u)?2rujlxj?1,ljlj?1,lyjl

（4.11）化为矩阵形式：

?u?ry(unjlnj,l?1?unj,l?1), （4.11）

12

2rx?ry）?rx0?1?（??rx1?（2rx?ry）?rx??0?rx1?（2rx?ry）??????000??n??u???1nl?0?0??u2l??????rx?0? ???????????n??u???rx1?（2rx?ry）???m?1,l?0?0?n?1nnrn?u?1l?ry(u1,l?1?u1,l?1)??

?un?1xu0l?2l?rnn?y(u1,l?1?u1,l?1)??????. ???????un?1rnnn?m?1,l?y(u1,l?1?u1,l?1)?ryuml??13

（4.12）

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

第5章 MATLAB编程

5.1 追赶法

对于前面求出的（4.10）和（4.12）矩阵，我们称之为三对角矩阵。在一些实际问题中，例如解常微分方程的边值问题，接热传导方程一级船体数学放样中建立三次样条函数等，都会要求解系数矩阵为对角占优的三对角方程组，解三对角矩阵我们有一种特殊的方法叫追赶法。

例如：求解线性方程组Ax?b，其中A为三对角矩阵。首先对矩阵A做如下三角分解：

A?LU

其中L为下三角矩阵，U为上三角矩阵，且主对角线元素全为1.

求解Ax?b等价于求解两个三角形方程组。先求解Ly?b的解y，再求解Ux?y的解x。

追赶法公式实际上就是把高斯消去法用到求解三对角方程组上去的结果。这是由于A特别简单，因此使得求解的计算公式非常简单，而且计算量很小。另外追赶法的计算公式中不会出现中间结果数量级的巨大增长和舍入误差的严重积累。最后在MATLAB中编程实现。

追赶法程序： 见附录。

例5.1 用追赶法求解线性方程组

?2x1?x2?6,??x1?3x2?2x3?1,???2x2?4x3?3x4??2,???3x3?5x4?1. ?解：在MATLAB命令窗口中输入求解程序：

>> A=[2 -1 0 0;-1 3 -2 0;0 -2 4 -3;0 0 -3 5]

A =

2 -1 0 0 -1 3 -2 0

0 -2 4 -3

14

0 0 -3 5

>>f=[6 1 -2 1]’; >>x=chase(A,f) x =

9 12 1 0.

输出计算结果为：

即追赶法求出了方程组的解为：

[x1,x2,x3,x4]?[9,12,1,0].

5.2 二维热传导方程有限差分解的MATLAB编程

程序见附录。

15