直接法求两点边值问题

课程设计（论文）任务书
数学与计算科学学院学院信息与计算科学专业班
课程名称科学仿真实验五
题目直接法求解两点边值问题（一）
任务起止日期：2014 年 6 月23 日～2014年7月 6 日
学生姓名学号
指导教师
教研室主任年月日审查
课程设计（论文）任务
注：1. 此任务书由指导教师填写。

如不够填写，可另加页。

2. 此任务书最迟必须在课程设计（论文）开始前下达给学生。

学生送交全部材料日期
学生（签名）指导教师验收（签名）
直接法求解两点边值问题（二）
摘要
线性方程组的数值解法可以分为直接法和迭代法两类。

所谓直接法，就是不考虑舍入误差，通过有限步骤四则运算即能求得线性方程组准确解的方法。

如克莱姆法则，但通过第一章的分析，我们知道用克莱姆法则来求解线性代数方程组并不实用，因而寻求线性方程组的快速而有效的解法是十分重要的。

本章讨论计算机上常用而有效的直接解法――高斯消去法和矩阵的三角分解等问题。

为方便计，设所讨论的线性方程组的系数行列式不等于零。

高斯消去法是解线性方程组最常用的方法之一，它的基本思想是通过逐步消元，把方程组化为系数矩阵为三角形矩阵的同解方程组，然后用回代法解此三角形方程组得原方程组的解。

关键词：线性方程组；直接解法；高斯消去法
DIRECT METHOD SOLVING TWO-POINT BOUNDARY VALUE PROBLEMS
(2)
ABSTRACT
Numerical algorithm of linear equations can be divided into two categories, direct method and iterative method. The so-called direct method, is not considered rounding error, through limited steps arithmetic which can obtain the accurate solution of linear equations method. Such as cramer's rule, but through the analysis of the first chapter, we know that cramer's rule is used to solve the linear algebraic equations is not practical, thus seeking quick and effective solutions of systems of linear equations solution is very important.
This chapter discuss computer commonly used and effective direct solution - gaussian elimination and triangle decomposition of matrices. For the convenience of meter, discussed the coefficient determinant of linear equations is not equal to zero.
Gauss elimination method is one of the most commonly used method of solving linear equations, the basic idea is to pass a gradual elimination, to coefficient matrix of the triangular matrix equations with solutions of the equations, then by back substitution method solving the triangle equations to the solution of the original equations.
Key words:linear equations; Direct method; Gaussian elimination
目录
1问题的提出 (1)
2 理论基础 (1)
2.1 高斯消去法 (2)
2.2 列主元消去法 (5)
2.3 矩阵的三角分解法 (6)
2.3.1 算法介绍 (6)
2.3.2 定理结论 (7)
2.3.3 计算公式 (9)
2.4解三对角方程组的追赶法 (10)
3 问题的求解 (12)
3.1顺序消去法 (12)
3.2 列主元消去法 (13)
3.3Doolittle分解法 (14)
3.4 追赶法 (15)
4 计算结果 (16)
参考文献 (20)
附录 (21)
1 问题的提出
考虑两点边值问题：
()()⎪⎩⎪⎨⎧==<<=+.
11,00,10,22y y a a dx dy dx y d ε 容易知道它的精确解为：
.1111ax e e a y x +⎪⎪⎭
⎫ ⎝⎛
---=
--εε
为了把微分方程离散，把[]1,0区间n 等分，令n
h 1
=，ih x i =，,1,,2,1-=n i 得到差分方程：
,212
11a h y y h
y y y i
i i i i =-++-++-ε
简化为：
()(),2211ah y y h y h i i i =++-+-+εεε
从而离散后得到的线性方程组的系数矩阵为：
()()()()⎥⎥
⎥
⎥⎥
⎥⎦
⎤
⎢⎢⎢⎢⎢⎢⎣⎡+-++-++-++-=h h h h h h h A εεεεεεεεεε2222
对1=ε，2
1
=
a ，100=n ，分别用顺序消去法、列主元消去法、Doolittle 分解法和追赶法求解线性方程组，然后比较与精确解的误差，对结果进行分析。

改变
n ，讨论同样问题。

2 理论基础
许多科学技术问题要归结为解含有多个未知量x 1, x 2, …, x n 的线性方程组：
⎪⎪⎩⎪⎪⎨
⎧=+++=+++=+++n
n nn n n n n n n b x a x a x a b x a x a x a b x a x a x a
22112
222212********* （2.1） 这里a ij (i , j = 1, 2, …, n )为方程组的系数，b i （i = 1, 2, …, n ）为方程组自由项。

方程组（2.1）的矩阵形式为：
AX = b
其中：
⎪⎪⎪⎪⎪⎭
⎫ ⎝⎛=⎪⎪⎪
⎪⎪⎭⎫
⎝⎛=⎪⎪⎪
⎪⎪⎭⎫
⎝⎛=n n nn n n n n b b b b x x x X a a a
a a a a a a A 2121212222111211
线性方程组的数值解法可以分为直接法和迭代法两类。

所谓直接法，就是不
考虑舍入误差，通过有限步骤四则运算即能求得线性方程组（2.1）准确解的方法。

如克莱姆法则，但通过第一章的分析，我们知道用克莱姆法则来求解线性代数方程组并不实用，因而寻求线性方程组的快速而有效的解法是十分重要的。

本章讨论计算机上常用而有效的直接解法――高斯消去法和矩阵的三角分解
等问题。

为方便计，设所讨论的线性方程组的系数行列式不等于零。

高斯（Gauss ）消去法是解线性方程组最常用的方法之一，它的基本思想是通
过逐步消元，把方程组化为系数矩阵为三角形矩阵的同解方程组，然后用回代法解此三角形方程组得原方程组的解。

2.1 高斯消去法
一般形式的线性方程组的解法中，为叙述问题方便，将b i 写成a i , n +1，i = 1,
2,…,n 。

⎪⎪⎩
⎪⎪
⎨
⎧=++++=++++=+++++++1,3322111
,223232221211,11313212111n n n nn n n n n n n n n n a x a x a x a x a a x a x a x a x a a x a x a x a x a
（2.1）
如果a 11 ≠ 0，将第一个方程中x 1的系数化为1，得：
)
1(1,1)1(12)1(121+=+++n n n a x a x a x
其中： )0(11
)
0()1(1a a a
ij
j
=
j = 1, …, n + 1
（记ij ij a a =)
0( i = 1, 2, …, n ; j = 1, 2, …, n + 1）
从其它n –1个方程中消x 1，使它变成如下形式
⎪⎪
⎩
⎪⎪⎨
⎧=++=++=++++++)1(1,)1(2)1(2)
1(1
,2)1(22)1(22)
1(1,1)1(12)1(121n n n nn n n n n n n n a x a x a a x a x a a x a x a x （2.2）
其中： n i a m a a ij
i ij ij ,,2)
1(1)1( =⋅-=
1,,3,211
)
1(1
1+==
n j a a m i i
由方程（2.2）到（2.3）的过程中，元素11a 起着重要的作用，特别地，把11a 称为主元素。

如果（2.3）中0)
1(22≠a ，则以)1(22a 为主元素，又可以把方程组（2.3）化为：
⎪⎪
⎪
⎩
⎪⎪⎪⎨⎧=++=++=+++=+++++++)2(1
,)2(3)2(3)
3(1,3)2(33)2(33)
2(1
,2)2(23)2(232)
1(1,1)1(12)1(121 n n n nn n n n n n n n n n n a x a x a a x a x a a x a x a x a x a x a x
（2.3）
针对（2.4） 继续消元，重复同样的手段，第k 步所要加工的方程组是：
⎪⎪
⎪⎪
⎪⎩
⎪⎪
⎪
⎪⎪⎨⎧=++=++=+++=+++=++++-+---+---+-----++)
1(1,)1()1()
1(1,)1()1()
1(1,1)1()1(11)
2(1
,2)2(23)2(232)
1(1,1)1(13)1(132)1(121 k n n n k nn k k nk k n k n k nn k k kk k n k n k kn k k k k n n n n n n a x a x a a x a x a a x a x a x a x a x a x a x a x a x a x
设0)
1(≠-k kk
a ，第k 步先使上述方程组中第k 个方程中x k 的系数化为1：
)(1,)()(1,k n k n k kn k k k k k a x a x a x ++=++
然后再从其它（n - k ）个方程中消x k ，消元公式为：
⎪⎪⎪
⎩⎪⎪⎪⎨⎧+=++=⋅-=++==----n
k i n k j a a a a n k k j a a a k kj
k ik k ij k ij k kk
k kj
k kj ,11,,11,,1,)
()1()1()()
1()
1()( （2.4） 按照上述步骤进行n 次后，将原方程组加工成下列形式：
⎪⎪
⎪⎩
⎪
⎪⎪⎨
⎧==+=+++=+++++-+---++)(1,)
1(1,1)1(1)
2(1
,2)2(23)2(232)1(1,1)1(13)1(132)1(121 n n n n n n n n n nn n n n n n n n a x a x a x a x a x a x a x a x a x a x 回代公式为：
⎪⎩
⎪⎨⎧-=-==∑+=++1
,,11
)
()
(1,)
(1
, n k x a a x a x n
k j j k kj
k n k k n n n
n （2.5）
综上所述，高斯消去法分为消元过程与回代过程，消元过程将所给方程组加
工成上三角形方程组，再经回代过程求解。

由于计算时不涉及x i , i = 1, 2, …, n ，所以在存贮时可将方程组AX = b ，
写成增广矩阵（A, b ）存贮。

下面，我们统计一下高斯消去法的工作量；在（2.5）第一个式子中，每执行
一次需要)(k n n --次除法，在（2.5）第二个式子中，每执行一次需要
)()]1([k n k n -⨯--次除法。

因此在消元过程中，共需要：
[]
)
12)(1(61
)1()1()()1(121++=+-=+-+-⨯+-∑∑==n n n k n k n k n k n n
k n
k 次乘作法。

此外，回代过程共有：)1(2
)(1-=
-∑=n n
k n n
k 次乘法。

汇总在一起，高斯消去法的计算量为：3
3)13(3232n
n n n n n -+=
-+次乘除法。

2.2 列主元消去法
在列主元消去法中，未知数仍然是顺序地消去的，但是把各方程中要消去的那个未知数的系数按绝对值最大值作为主元素，然后用顺序消去法的公式求解。

由于解方程组取决于它的系数，因此可用这些系数（包括右端项）所构成的“增广矩阵”作为方程组的一种简化形式。

列主元消去法计算步骤：
1、输入矩阵阶数 n ，增广矩阵 )1,(+n n A
2、对 n k ,,2,1 =
（1）按列选主元：选取l 使
0max ≠=≤≤ik n
i k lk a a
（2）如果 k l ≠，交换 )1,(+n n A 的 第k 行与底l 行元素 （3）消元计算 ：
n k i a a m kk
ik
ik ,,1 +=←
1,,1 ,,1, ++=+=-←n k j n k i a m a a kj ik ij ij
（4）回代计算：
1,,1, 1
1, -=-
←∑+=+n n i x a
a x n
i j j ij
n i i
（5）输出解向量：n x x x ,,,21 。

2.3 矩阵的三角分解法
下面我们进一步用矩阵理论来分析高斯消去法，从而建立矩阵的三角分解定
理，而这个定理在解方程组的直接解法中起着重要作用，我们将利用它来推导某些具有特殊的系数矩阵方程组的数值解法。

2.3.1 算法介绍
考虑线性方程组：
AX = b ，A ∈R n ⨯n
设解此方程组用高期消去法能够完成（不进行变换两行的初等变换），由于对A 施行初等变换相当于用初等矩阵左乘A ，于是，高斯消去法的求解过程用矩阵理论来叙述如下：
记：
⎪⎪⎪⎪⎪⎪⎪⎭
⎫
⎝⎛++=⎪⎪⎪⎪
⎪⎪⎪⎭
⎫
⎝⎛--=+-+1 1 111 1 11,,11
,,1k n k k k k n k k k l l L l l L
其中),,1()
()
(n i j a a l i ii
i ji ji +==
记 A A ≡)1(
于是：
⎪⎪⎪⎪⎪⎭
⎫ ⎝⎛⎪
⎪⎪⎪⎪⎪⎪⎪⎭
⎫
⎝⎛--=)1(1
)
1(22
)1(21)
1(1)1(12)
1(11)1(11
)1(1)
1(11
)1(21)
1(111n n n a A a a a a a a a a LA
)2()2(22)
2(12
)2(11)2()2(2)
2(2)2(22)
1(1)1(12)1(11)1
(1111)1(22
1)1(11
)1(2211
)
1(111)1(11
100
01A A A A a a a a a a a a r c A r a A c r a I a c nn n n n T T
T n =⎪⎪⎭
⎫
⎝⎛=⎪⎪⎪⎪⎪
⎭⎫
⎝
⎛=⎪⎪⎪⎭
⎫ ⎝⎛-=⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎭
⎫ ⎝⎛-=-
)3()3(22)
3(12
)3(11
)2(332
2)2(22
111
2)2(222
)2(20000001A A A A A c r a r a I a c A L T T n T T
=⎪⎪⎭
⎫ ⎝⎛=⎪⎪⎪⎪⎭
⎫
⎝
⎛⎪⎪⎪⎪⎪⎪⎭
⎫
⎝
⎛
-=- 如此继续下去，到n – 1步时有：
U a a a a a A L nn n n n n ≡⎪⎪
⎪⎪⎪
⎭
⎫ ⎝⎛=--)2()
2(2)2(22)
1(1)1(11)
1(1
也就是说：
U A L L L L n n =--1221
记：
⎪⎪⎪⎪⎪
⎪⎭
⎫ ⎝⎛==---1112132
312111
1
211n n n l l l l l L L L L 则有：A = LU 。

其中L 为单位下三角矩阵，U 为上三角矩阵。

2.3.2 定理结论
总结上述讨论得到重要性定理： 定理1：
（矩阵的三角分解）设A 为n ⨯ n 实矩阵，如果解AX = b 用高斯消去法能够
完成（限制不进行行的交换，即n k a k kk ,,2,1,0)
( =≠），则矩阵A 可分解为单位
下三角矩阵L 与上三角知阵U 的乘积。

A = LU
且这种分解是唯一的。

证明：由上述的讨论，存在性已得证，现在证唯一性。

设 A = L 1 U 1 = LU
其中L 1，L 1为单位下三角阵，U 1，U 为上三角阵，设11-U 存在，于是有
1111--=UU L L
上式右端为上三角矩阵，左边为单位下三角阵，故应为单位矩阵。

即
L 1 = L U 1 = U 。

定理2：
约化主元素0)
(≠i ii a （i = 1, 2, …, k ）的充要条件是矩阵A 的顺序主子式：
,0,022
21
12112111≠=≠=a a a a D a D
01111≠=kk
k k
k a a a a D
由上述讨论知，解 AX = b 的高斯消去法相当于实现了A 的三角分解，如果我
们能直接从矩阵A 的元素得到计算L ，U 的元素的公式，实现A 的三角分解，而不需要任何中间步骤，那么求解AX = b 的问题就等价于求解两个三角形矩阵方程组： （1）Ly = b 求y
（2）UX = y 求x
下面来说明L 、U 的元素可以由A 的元素直接计算确定。

显然，由矩阵乘法。

),,2,1(11n i u a i
i ==
得到U 的第一行元素；由1111u l a i i =得：
),,2,1(11
1
1n i u a l i i ==
即L 的第一列元素。

设已经求出U 的第1行~第r – 1行元素，L 的第1列~第r – 1列元素，由矩阵乘法可得：
),0(1
11k r l u u l u l a rk ri
ki r k rk ki n k rk ri <=+==∑∑-==
rr ir kr r k ik kr n
k ik ir u k u l u l a +==∑∑-==1
1
1
即可计算出U 的第r 行元素，L 的第r 列元素。

2.3.3 计算公式 综上所述，可得到用直接三角分解法解AX = b 的计算公式。

（1）n i a u i i ,,2,111 ==
（3.1）
n i u a l i i ,,2,111
11 ==
（3.2） 对于r = 2, 3, …, n 计算 （2）计算U 的第r 行元素
),,1,( 1
1n r r i u l a u r k ki rk ri ri +=-=∑-=
（3.3）
（3）计算L 的第r 列元素 (r ≠ n )
),,1()
(1
1
n r i u u l a l rr
r k kr ik ir ir +=-=
∑-=
（3.4）
（4） ⎪⎩
⎪⎨
⎧=-==∑-=)
,,3,2(111n i y l b y b y k
i k ik i i （3.5）
（5）⎪⎪⎩
⎪⎪⎨⎧
-=⎪⎭⎫ ⎝⎛-==∑+=)
1,2,,1(1 n i u x u y x u y x ii
k n
i k ik i i nn
n n （3.6）
（1）、（2）、（3）是矩阵A 的LU 分解公式,称为Doolittle 分解。

同理，可推出矩阵A = LU 分解的另一种计算公式，其中L 为下三角，U 为单位上三角，这种矩阵的分解公式称为矩阵的Crout 分解。

2.4 解三对角方程组的追赶法
在实际问题中，经常遇到以下形式的方程组：
⎪
⎪⎪
⎪⎩⎪
⎪⎪⎪⎨⎧=+=++=++=++=+-------+-n n n n n n n n n n n n k k k k k k k d x b x a d x c x b x a d x c x b x a d x c x b x a d x c x b
111112111232221212111 （4.1）
这种方程组的系数矩阵称为三对角矩阵。

⎪⎪⎪
⎪⎪⎪⎪⎪⎭
⎫
⎝
⎛=---n n
n n n k k k b a c b a c b a c b a c b A 111
22211 以下针对这种方程组的特点提供一种简便有效的算法—追赶法。

追赶法实际上是高斯消去法的一种简化形式，它同样分消元与回代两个过程。

先将（4.1）第一个方程中x 1的系数化为1：
1
12111b d x b c x =+
记：1
1
11
1
1b d y b c r =
=
（4.2）
有： 1211y x r x =+
注意到剩下的方程中，实际上只有第二个方程中含有变量x 1，因此消元手续可以简化。

利用（4.2）可将第二个方程化为：
2312y x r x =+
这样一步一步地顺序加工（4.1）的每个方程，设第k – 1个方程已经变成：
111---=+k k k k y x r x
（4.3）
再利用（4.3）从第k 个方程中消去x k -1，得：
k k k k k k k k k a y d x c x a r b 111)(-+--=+-
同除()k k k a r b 1--，得：
n k a r b a y d x a r b c x k
k k k
k k
k k k k k k ,,3,21111 =--=-+
--+-
记：
k
k k k
k k k k
k k k
k a r b a y d y a r b c r 111-----=
-=
则有：k k k k y x r x =++1
这样做n – 1步以后，便得到：
111---=+n n n n y x r x
将上式与（4.1）中第11个方程联立，即可解出：x n = y n 这里：
n
n n n
n n n a r b a y d y 11----=
于是，通过消元过程，所给方程组（4.1）可归结为以下更为简单的形式：
⎪⎪⎪⎩⎪
⎪⎪⎨⎧==+=++n
n k k k k y x y x r x y x r x 11211 （4.4）
这种方程组称作二对角型方程组，其系数矩阵中的非零元素集中分步在主对
角线和一条次主对角线上：
⎪⎪⎪
⎪⎪⎪⎭
⎫
⎝⎛-11111121n k r r r r 对加工得到的方程组（3.4）自下而上逐步回代，即可依次求出x n ，x n -1，…，
x 1，计算公式为：
⎩⎨
⎧--=-==+1
,,2,11
n n k x r y x y x k k k k n n （4.5）
上述算法就是追赶法，它的消元过程与回代过程分别称作“追”过程与“赶”
过程。

综合追与赶的过程，得如下计算公式：
⎪⎪⎪⎩
⎪⎪⎪⎨⎧=--=-===---n
k a r b a y d y a r b c r b d y b c r k k k k k k k
k
k k k k
,,3,211111
1111
（4.6）
⎩⎨
⎧--=-==+1
,,2,11
n n k x r y x y x k k k k n
n （4.7）
3 问题的求解
3.1 顺序消去法
function x = gaussElim(A,b) N = length(b); if A(1,1)==0
n_temp=find(A(:,1)~=0); c=A(n_temp(1),:); A(n_temp(1),:)=A(1,:); A(1,:)=c; end for j=2:N
m = A(j:N,j-1)/A(j-1,j-1)
A(j:N,:) = A(j:N,:) - m*A(j-1,:);
b(j:N) = b(j:N) - m*b(j-1);
if A(j,j)==0 &j~=N
n_temp=find(A(j:end,j)~=0);
c=A(j-1+n_temp(1),:);
A(j-1+n_temp(1),:)=A(j,:);
A(j,:)=c;
end
end
if A(N,N)==0
error('the coefficient matrix of the linear equations is sigular'); end
x = zeros(N,1);
x(N) = b(N)/A(N,N);
for j=N-1:-1:1,
x(j) = (b(j)-A(j,j+1:N)*x(j+1:N))/A(j,j);
end
3.2 列主元消去法
function x=f(A)
%列主元高斯消去法
[m,n]=size(A);
B=zeros(size(A(1,:)));
for i=1:(m-1)
for j=i:m
if A(j,i)>A(i,i)
B=A(j,:);A(j,:)=A(i,:);A(i,:)=B; %选取主元
end
end
for j=(i+1):m
A(j,:)=A(j,:)-A(i,:)*A(j,i)/A(i,i); % 消元法
end
end
%回代过程
x(m)=A(m,n)/A(m,m);
for i=(m-1):-1:1
x(i)=(A(i,n)-A(i,i+1:m)*x(i+1:m)')/A(i,i);
end
x
3.3 Doolittle分解法
function[x,l,u]=malu(A,b)
%用途：用LU分解法解方程组
%格式：[x,l,u]=malu(A,b),A为系数矩阵，b为右端向量，x返回解向量，l返回下三角矩阵，u返回上三角矩阵
format short %LU分解
n=length(b);
u=zeros(n,n);l=eye(n,n);
u(1,:)=A(1,:);l(2:n,1)/u(1,1);
for k=2:n
u(k,k:n)=A(k,k:n)-l(k,1:k-1)*u(1:k-1,k:n);
l(k+1:n,k)=(A(k+1:n,k)-l(k+1:n,k-1)*u(l:k-1,k))/u(k,k);
End %解下三角方程组Ly=b
y=zeros(n,1);
X(n)=y(n)/u(n,n);
For k=n-1:-1:1
x(k)=(y(k)-u(k,k+1:n)*x(k+1：n))/u(k,k); end
3.4追赶法
function [x,L,U]=thomas(e,p,n)
a=e.*ones(1,n-2);
b=-(2*e+1/n).*ones(1,n-1);
c=(e+1/n).*ones(1,n-2);
f=(p*(1/n)*(1/n)).*ones(1,n-1);
n=length(b);
%对A进行分解a,b,c,f
u(1)=b(1);
for i=2:n
if u(i-1)~=0
l(i-1)=a(i-1)/u(i-1);
u(i)=b(i)-l(i-1)*c(i-1);
else
break;
end
end
L=eye(n)+diag(l,-1);
U=diag(u)+diag(c,1);
x=zeros(n,1);
y=x;
%求解Ly=b
y(1)=f(1);
for i=2:n
y(i)=f(i)-l(i-1)*y(i-1);
end
%求解Ux=y if u(n)~=0
x(n)=y(n)/u(n); End
4 计算结果
（1）当e=1，a=0.5，n=100时，以上方法求解线性方程组解出的结果完全相同，但与真正的结果相差很大。

平均相对误差S1，S2更是高达1.1371.因为：
()(),2211ah y y h y h i i i =++-+-+εεε
5
210501.001.05.0-⨯=⨯⨯=ah 与01.101.01=+=+h ε，
01.201.02)2(-=--=+-h ε，1=ε相差很大。

（2）当e=1，a=0.5，n=10时，以上方法求解线性方程组解出的结果完全相同，与真正的结果相差无几。

平均相对误差S1，S2是0.0040.因为：
()(),2211ah y y h y h i i i =++-+-+εεε
3-21051.01.05.0⨯=⨯⨯=ah 与1.11.01=+=+h ε，1.21.02)2(-=--=+-h ε，
ε不在一个数量级。

=
1
（3）当e=1，a=0.5，n=200时，以上方法求解线性方程组解出的结果完全相同，与真正的结果相差无几。

平均相对误差S1，S2是 0.0021。

（4）当e=1，a=0.5，n=300时，以上方法求解线性方程组解出的结果完全相同，与真正的结果相差无几。

平均相对误差S1，S2是 0.0016。

参考文献
[1] 严蔚敏、吴伟民．数据结构(C语言版)[M].北京：清华大学出版社，1997.4.
[2] 李庆扬、王能超、易大义.数值分析[M].武汉.华中科技大学出版社，2006.7.
[3] 清华大学、北京大学计算方法编写组.计算方法[M].北京.科学出版社，1980.
附录
（1）命令文件
e=1;
a=0.5;
n=100;
x=(0+1/n):1/n:(1-1/n);
y=(1-a)*(1-exp(-x./e))/(1-exp(-1/e))+a.*x; %正确结果A=finput(e,a,n); %输入矩阵A
y1=f(A); % 列主元消去法结果
[y2,L,U]=thomas(e,a,n);%追赶法结果
plot(x,y,x,y1,'*',x,y2,'+')
s1=0;
s2=0;
for i=1:n-2
s1=s1+(y1(i)-y(i))/ y(i);
s2=s2+(y2(i)-y(i))/ y(i);
end
s1= s1/(n-2)
s2=s2/(n-2)
（2）输入矩阵A
function A=finput(e,a,n)
%输入矩阵
A=zeros(n-1,n);
h=1/n;
for i=1:n-1
for j=1:n
if j==i
A(i,j)=-(2*e+h);
end
if j==(i-1)
A(i,j)=e;
end
if j==(i+1)
A(i,j)=e+h;
end
if j==n
A(i,j)=a*h*h;
end
end
End
（3）顺序消去法
function x = gaussElim(A,b)
N = length(b);
if A(1,1)==0
n_temp=find(A(:,1)~=0);
c=A(n_temp(1),:);
A(n_temp(1),:)=A(1,:);
A(1,:)=c;
end
for j=2:N
m = A(j:N,j-1)/A(j-1,j-1)
A(j:N,:) = A(j:N,:) - m*A(j-1,:); b(j:N) = b(j:N) - m*b(j-1);
if A(j,j)==0 &j~=N
n_temp=find(A(j:end,j)~=0);
c=A(j-1+n_temp(1),:);
A(j-1+n_temp(1),:)=A(j,:);
A(j,:)=c;
end
end
if A(N,N)==0
error('the coefficient matrix of the linear equations is sigular'); end
x = zeros(N,1);
x(N) = b(N)/A(N,N);
for j=N-1:-1:1,
x(j) = (b(j)-A(j,j+1:N)*x(j+1:N))/A(j,j);
end
（4）列主元高斯消去法
function x=f(A)
%列主元高斯消去法
[m,n]=size(A);
B=zeros(size(A(1,:)));
for i=1:(m-1)
for j=i:m
if A(j,i)>A(i,i)
B=A(j,:);A(j,:)=A(i,:);A(i,:)=B;%选取主元
end
end
for j=(i+1):m
A(j,:)=A(j,:)-A(i,:)*A(j,i)/A(i,i); % 消元法
end
end
%回代过程
x(m)=A(m,n)/A(m,m);
for i=(m-1):-1:1
x(i)=(A(i,n)-A(i,i+1:m)*x(i+1:m)')/A(i,i);
end
X
（5）Doolittle分解法
function[x,l,u]=malu(A,b)
%用途：用LU分解法解方程组
%格式：[x,l,u]=malu(A,b),A为系数矩阵，b为右端向量，x返回解向量，l返回下三角矩阵，u返回上三角矩阵
format short %LU分解
n=length(b);
u=zeros(n,n);l=eye(n,n);
u(1,:)=A(1,:);l(2:n,1)/u(1,1);
for k=2:n
u(k,k:n)=A(k,k:n)-l(k,1:k-1)*u(1:k-1,k:n);
l(k+1:n,k)=(A(k+1:n,k)-l(k+1:n,k-1)*u(l:k-1,k))/u(k,k);
End %解下三角方程组Ly=b
y=zeros(n,1);
X(n)=y(n)/u(n,n);
For k=n-1:-1:1
x(k)=(y(k)-u(k,k+1:n)*x(k+1：n))/u(k,k);
End
（6）追赶法求解
function [x,L,U]=thomas(e,p,n)
a=e.*ones(1,n-2);
b=-(2*e+1/n).*ones(1,n-1);
c=(e+1/n).*ones(1,n-2);
f=(p*(1/n)*(1/n)).*ones(1,n-1); n=length(b);
%对A进行分解a,b,c,f
u(1)=b(1);
for i=2:n
if u(i-1)~=0
l(i-1)=a(i-1)/u(i-1);
u(i)=b(i)-l(i-1)*c(i-1); else
break;
end
end
L=eye(n)+diag(l,-1);
U=diag(u)+diag(c,1);
x=zeros(n,1);
y=x;
%求解Ly=b
y(1)=f(1);
for i=2:n
y(i)=f(i)-l(i-1)*y(i-1); end
%求解Ux=y
if u(n)~=0
x(n)=y(n)/u(n);
end。