追赶法解三对角方程组
《数值分析》课程设计追赶法解三对角方程组
院（系）名称信息工程学院
专业班级10普本信计
学号100111014
学生姓名刘银朋
指导教师张荣艳
2013 年05 月31日
数值分析课程设计评阅书
题目追赶法解三对角方程组
学生姓名刘银朋学号100111014 指导教师评语及成绩
指导教师签名：
年月日答辩评语及成绩
答辩教师签名：
年月日
教研室意见
总成绩：
教研室主任签名：
年月日
课程设计任务书
2012—2013学年第二学期
专业班级：10普本信息与计算科学学号：100111014 姓名：刘银朋
课程设计名称：数值分析Ⅰ、Ⅱ
设计题目：追赶法解三对角方程组
完成期限：自2013 年05月21 日至2013年05 月31日共10天
设计依据、要求及主要内容：
一、设计目的
理解追赶法,掌握追赶法的算法设计以及关于追赶法的分析和综合应用,能
够较熟练的应用Matlab软件编写求解追赶法的程序和应用Matlab软件数据库软
件.
二、设计内容
（1）认真挑选有代表性的三对角方程组.
（2）认真梳理解三对角方程组的解题思路.
（3）比较追赶法和高斯消去法的计算精度.
三、设计要求
1．先用Matlab数据库中的相应的函数对选定的方程,求出具有一定精度的解. 2．然后使用所用的方法编写Matlab程序求解.
3.对于使用多个方程解同意问题的,在界面上要设计成菜单的形式.
计划答辩时间：2013年06 月 5 日
工作任务鱼工作量要求：
查阅文献资料不少于3篇,课程设计报告1篇不少于3000字.
指导教师（签字）：教研室主任（签字）：
批准日期：2013 年05 月20 日
追赶法解三对角方程组
摘要
本文主要通过运用追赶法来求解三对角方程组的问题.追赶法是用来求解三对角方程组的专用方法,对于三对角方程组,追赶法比Gauss消去法的计算量要小的多,本文主要介绍了追赶法的原理,并用Matlab编写求解程序,以实现对三对角方程组的求解,进一步解决实际中的问题.并且根据所得出的结果分析追赶法算法和高斯消去的法的计算精度.
关键词：追赶法,三对角方程组,追赶法的Matlab程序
目录
1.前言 (1)
2.解题思想和方法 (1)
2.1 追赶法解题思想 (1)
2.2追赶法解题原理 (4)
3.对追赶法的MATLAB求解 (4)
3.1实验程序 (4)
3.2 应用举例 (5)
4.与高斯消去法的精度比较 (5)
课程设计总结 (8)
参考文献 (8)
1.前言
当今很多科学与工程计算问题大都可以化为线性方程组的形式,所以有效的求解线性方程组在科学和工程计算中是非常重要的.虽然线性代数方程的求解方法和数值计算软件包均很成熟,但随着并行计算机的发展,问题的求解速度和解题规模都大大提高,因而使数值计算方法和响应的数学软件包都产生了变化,相应的线性方程组的有效并行求解也引起了人们的普遍关注.
追赶法是用来求解三对角方程组的专用方法,生活中很多实际问题,都归结为求解线性方程组.例如解常微分方程边值问题,解热传导方程以及船体数学放样中建立三次样条函数等,都会要求就解系数矩阵成三对角线性的方程组.而解三对角方程组的最简单方法是追赶法,公式简单,计算量小,所占用的存储单元少,所以在小机器上也能求解.追赶大事用来求解三对角方程组的专用方法,对于三对角方程组,追赶法比Gauss消去法的计算量要小得多.应用追赶法求解三对角线性方程,追赶法仍然保持LU分解特性,它是一种特殊的LU分解.充分利用了系数矩阵的特点,而且使之分解更简单,得到对三对角线性方程组的快速解法. 本文讨论使用追赶法解线性方程组.介绍追赶法的理论,求解线性方程组的追赶法的实现以及追赶法的应用.
2.解题思想和方法
三对角矩阵是一种具有特殊意义的带状矩阵.用差分法求解二阶常微分方程边值问题时,最后常规解为求解具有三对角系数矩阵的线性方程组.对三对角矩阵实行Doolittle（或Crout）分解,便得到求解三对角方程组的最有效方法---追赶法.
设矩阵A非奇异,A有Crout分解A=LU,其中Ｌ为下三角矩阵,Ｕ为单位上三角矩阵,可先依次求出L,U中的元素后,令Ux=y,先求解下三角方程组Ly=f得出y,再求解上三角方程组Ux=y,事实上求解三对角方程组的2追赶法将矩阵三角分解的计算与求解两个三角方程组的计算放在一起,使算法更为紧凑.
本文主要介绍追赶法的原理,以及编写Matlab程序实现在计算机上的应用并分析他们的计算精度,比较解三对角方程组的最优解的问题.
2.1 追赶法解题思想
在实际问题中,经常遇到以下形式的方程组
?
??
???
?????=+=++=++=++=+-------+-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 （2.11）
这种方程组的系数矩阵A 为三对角矩阵,即
???
??????
?
?
?=---n n
n n n k k k b a c b a c b a c b a c b A 111
22211 以下针对这种方程组的特点提供一种简便有效的算法—追赶法.追赶法实际上是高斯消去法的一种简化形式,它同样分消元与回代两个过程.
先将（2.11）第一个方程中x 1的系数化为1
112111b d
x b c x =+
1
1
11
11b d y b c r =
=
（2.12）
有1211y x r x =+注意到剩下的方程中,实际上只有第二个方程中含有变量x 1,因此消元手续可以简化.利用（2.12）可将第二个方程化为2312y x r x =+,这样一步一步地顺序加工（2.11）的每个方程,设第k – 1个方程已经变成
111---=+k k k k y x r x
（2.13）
再利用（2.13）从第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
将上式与（2.11）中第11个方程联立,即可解出x n = y n ,这里n
n n n
n n n a r b a y d y 11----=
于是,通过消元过程,所给方程组（2.11）可归结为以下更为简单的形式：
?????
????==+=++n
n k k k k y x y x r x y x r x 11211 （2.14） 这种方程组称作二对角型方程组,其系数矩阵中的非零元素集中分步在主对角线和一条次主对角线上
???
???
?
?
??-1111112
1n k r r r r 对加工得到的方程组（2.14）自下而上逐步回代,即可依次求出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 （2.15）
上述算法就是追赶法,它的消元过程与回代过程分别称作“追”过程与“赶”过程.综合追与赶的过程,得如下计算公式：
????
?????=--=-===---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,2111111111
（2.16）
??
?--=-==+1
,,2,11
n n k x r y x y x k k k k n n （2.17）
2.2追赶法解题原理 (1) 计算{}i β的递推公式
111/c b β=,/(),2,3,
,1i i i i i c b a i n ββ=-=-.
(2) 解Ly=f
111/y f b =, 1()/(),2,3,
,;i i i i i i i y f a y b a i n β-=--=
(3) 解Ux=y
11,,1,2,
2,1n n i i i i x y x y x i n n β-+==-=--.
我们将计算系数121n βββ-→→→及12n y y y →→
→的过程称为追的过程.
将计算方程组的解11n n x x x -→→
→的过程称为追赶的过程.
3.对追赶法的MATLAB 求解
3.1实验程序
function x=chase(a,b,c,f) %定义函数chase n=length(b); if n-1==length(a) for i=n-1:-1:1 a(i+1)=a(i); end end
%将a 设置为n 维向量 c(1)=c(1)/b(1);f(1)=f(1)/b(1); for i=2:n-1
b(i)=b(i)-a(i)*f(i-1)/b(i); c(i)=c(i)/b(i);
f(i)=(f(i)-a(i)*f(i-1))/b(i);
end
f(n)=(f(n)-a(n)*f(n-1))/(b(n)-a(n)*c(n-1));
for i=n-1:-1:1
f(i)=f(i)-c(i)*f(i+1);
end
x=f;
3.2 应用举例
例1用追赶法解三对角方程组
设
270
329
045
A
??
?
= ?
?
??
5
6
7
f
??
?
= ?
?
??
求解线性方程组Ax=F
解：（1）在Matlab中编写一个名为chase.m的M文件,依次输入数据如下；
>> A=[2 7 0;3 2 9;0 4 5];
>> a=[3 4];
>> b=[2 2 5];
>>c=[7 9];
>> f=[5 6 7];
>> x=chase (a,b,c,f)
x= -3.0140 1.5754 0.1397
得到输出结果x= -3.0140 1.5754 0.1397即为原线性方程组的解.
结果验证：
>>A*x’
ans=5.0000
6.0000
7.0000
得到ans=f即结果正确
4. 与高斯消去法的精度比较
事实上,追赶法的求解过程就是将系数矩阵分解两个简单的二对角线矩阵,从而归结为求解两个简单三角形方程组的过程. 高斯消去法是求解线性方程组的最基本方法之一.Gauss 消去法是针对一般的线性方程组,与线性代数中的初等变换解线性方程组方法类似.追赶法只是针对系数矩阵为三对角阵的方程组,因此是一种特殊的方程组.此方法效率较高,不过不适用于一般的线性方程组.下面对追赶法与高斯消去法计算量进行比较.
例2 编写高斯消去法解线性方程组的Matlab 程序,并求解下面的线性方程组
1231231
235282821361
x x x x x x x x x ++=??
+-=??--=? 解：编程如下
f unction [RA,RB,n,X]=gaus(A,b) B=[A b];
n=length(b);RA=rank(A); RB=rank(B);zhica=RB-RA; if zhica>0,
disp('请注意：因为RA~=RB,所以此方程组无解.') return end if RA==RB if RA==n
disp('请注意：因为RA=RB=n,所以此方程组有唯一解.') X=zeros(n,1); C=zeros(1,n+1); for p= 1:n-1 for k=p+1:n
m= B(k,p)/ B(p,p); B(k,p:n+1)= B(k,p:n+1)-m* B(p,p:n+1); end end
b=B(1:n,n+1);A=B(1:n,1:n); X(n)=b(n)/A(n,n);
for q=n-1:-1:1
X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)))/A(q,q);
end
else disp('请注意：因为RA=RB end end 用Matlab调用程序：在Matlab中编写一个名为gauss.m的M文件,依次输入数据如下: A=[5,2,1;2,8,-3;1,-3,-6];b=[8;21;1]; [RA,RB,n,X]=gauss（A,b) 运行结果 RA =3 RB =3 n = 3 X =1 2 -1 ans = 1.0000 2.0000 -1.0000 综上所述,追赶法的原理和高斯消去法相同,但考虑到方程组的特点,计算时会把大量零元素撇开,从而大大节省计算量.追赶法实际上就是把高斯消去法用到求解三对角方程组上的结果.这时由于A特别简单,因此使得求解的计算公式非常简单,而且计算量仅为2(n-1)+(n-1)+[1+2(n-1)]=5n-4次乘除法,而另外解一个方程组Ax=f仅增加3n-2次乘除运算,易见追赶法的计算量是比较小的. 课程设计总结 首先,这次的课程设计让我把以前学习到的非线性方程的求解知识得到巩固和进一步的提高,对已有知识有了更进一步的理解和认识．同时也发现了自身有许多的不足之处.再者,在这次课程设计中碰到了很多的问题,自己通过上网,查阅了相关的书籍以及和我组其他的成员讨论的方式来解决的.另外编写程序是我的一大弱项,在以前做实验时总是按照老师给出的程序不加思考的敲进去,得出来结果就行,但在做这次课程设计时,需要自己真正的去理解和编程． 通过这次课程设计,我从中学会了很多,也发现自己真的还有很多不足以及很多东西需要去学习．所以在以后的生活学习中要不断的扩大自己的视野,多学习一些与专业有关的知识,不能只满足于课本上的知识．所以在完成本专业的基础上,要不断涉猎,完善自我,希望自己在以后的课程中会得到更好的锻练．总的来说这次课程设计还是有很多的收获的,并且特别感谢我们组的成员在做课程设计的过程中对我的帮助． 参考文献 [1] 冯国忱,黄明游．数值分析(上册)[M]．北京：高等教育出版社,2007． [2] 宋叶志,贾东志．MATLAB数值分析与应用[M]．北京：机械工业出版社,2009． [3] 张德丰．MATLAB数值分析与应用[M]．北京：国防工业出版社,2007． [4] 周煦,计算机数值方法及程序设计,机械工业出版社,2004.10.01. [5] 张静,杨文平,石傅强等,MATLAB程序设计与实例应用,中国铁道出版社,2003. 一、算法理论 在一些实际问题中，例如解常微分方程边值问题，解热传导方程以及船体数学放样中建立三次样条函数等，都会要求解系数矩阵为对角占优的三对角线方程组 ??????? ? ??=???????? ?????????? ? ?-----n n n n n n n n n f f f f x x x x b a c b a c b a c b 121121111 22211 ， 简记为f Ax =. 求解f Ax =等价于解两个三角形方程组： y f Ly 求,=；x y Ux 求,=.从而得到解三对角线方程组的追赶法公式： (1)计算{}i β的递推公式 ();1,,3,2,/,/111-=-==n i a b c b c i i i i i βββ (2)… (3) 解f Ly = ()();,,3,2,/,/11111n i a b y a f y b f y i i i i i i i =--==--β (4)解y Ux = .1,2,2,1,,1 --=-==+n n i x y x y x i i i i n n β 我们将计算系数 的过程称为追的过程，及n n y y y →→→→→→- 21121βββ 将计算方程组的解 的过程称为赶的过程。11x x x n n →→→- — 二、算法框图 ; \ 三、 算法程序 #include <> #include <> #include<> #define N 20 double a[N], b[N], c[N-1], f[N], r[N]; int n; （1） void LUDecompose(); ???????? ??2100012100012100012100012A －－－－－－－－＝??? ?? ?? ? ??=00001b 回车。 （2） 显示出 请输入下三角元素 输入4个a 值：-1 -1 -1 -1，回车。 （3） 显示出 请输入主对角线元素 输入5个b 值：2 2 2 2 2 ，回车。 （4） ！ （5） 显示出 请输入上三角元素 输入4个c 值：-1 -1 -1 -1，回车。 （6） 显示出 请输入5个方程组右端顶：1 0 0 0 0，回车。 其解为????? ????166667 .0333333.0500000.0666667.0833333 .0 例2．用该程序计算三对角线方程组 《数值分析》课程设计追赶法解三对角方程组 院（系）名称信息工程学院 专业班级10普本信计 学号100111014 学生姓名刘银朋 指导教师张荣艳 2013 年05 月31日 数值分析课程设计评阅书 课程设计任务书 2012—2013学年第二学期 专业班级：10普本信息与计算科学学号：100111014 姓名：刘银朋 课程设计名称：数值分析Ⅰ、Ⅱ 设计题目：追赶法解三对角方程组 完成期限：自2013 年05月21 日至2013年05 月31日共10天 设计依据、要求及主要内容： 一、设计目的 理解追赶法,掌握追赶法的算法设计以及关于追赶法的分析和综合应用,能 够较熟练的应用Matlab软件编写求解追赶法的程序和应用Matlab软件数据库软 件. 二、设计内容 （1）认真挑选有代表性的三对角方程组. （2）认真梳理解三对角方程组的解题思路. （3）比较追赶法和高斯消去法的计算精度. 三、设计要求 1．先用Matlab数据库中的相应的函数对选定的方程,求出具有一定精度的解. 2．然后使用所用的方法编写Matlab程序求解. 3.对于使用多个方程解同意问题的,在界面上要设计成菜单的形式. 计划答辩时间：2013年06 月 5 日 工作任务鱼工作量要求： 查阅文献资料不少于3篇,课程设计报告1篇不少于3000字. 指导教师（签字）：教研室主任（签字）： 批准日期：2013 年05 月20 日 追赶法解三对角方程组 摘要 本文主要通过运用追赶法来求解三对角方程组的问题.追赶法是用来求解三对角方程组的专用方法,对于三对角方程组,追赶法比Gauss消去法的计算量要小的多,本文主要介绍了追赶法的原理,并用Matlab编写求解程序,以实现对三对角方程组的求解,进一步解决实际中的问题.并且根据所得出的结果分析追赶法算法和高斯消去的法的计算精度. 关键词：追赶法,三对角方程组,追赶法的Matlab程序 追赶法求解三对角线性方程组 一 实验目的 利用编程方法实现追赶法求解三对角线性方程组。 二 实验内容 1、 学习和理解追赶法求解三对角线性方程组的原理及方法； 2、 利用MATLAB 编程实现追赶法； 3、 举例进行求解，并对结果进行分。 三 实验原理 设n 元线性方程组Ax=d 的系数矩阵A 为非奇异的三对角矩阵 11222=(1)(n 1)()()a c b a c A a n c b n a n ??????????--?????? ………… 这种方程组称为三对角线性方程组。显然，A 是上下半宽带都是1的带状矩阵。设A 的前n-1个顺序主子式都不为零，根据定理2.5的推论，A 有唯一的Crout 分解，并且是保留带宽的。 其中L 是下三角矩阵，U 是单位上三角矩阵。利用矩阵相乘法，可以1112212(1)1u(n 1)()()1l u m l u A LU l n m n l n ????????????????==?????--????????????…………… 得到： 由上列各式可以得到L 和U 。 引入中间量y ，令 y Ux =，则有： 已知 L 和d ，可求得y 。 则可得到y 的求解表达式： 11/1 2,3,,()(1)*y()=()[()(1)]/y d l i n m i y i li i di y i di m i y i li ==-+=--… 1111111/1(2)(1)(1)u (1)(11)/(1)(1)(1)l a l u c u c l mi bi i n a i m i i l i i n ci li ui ui ci li l i a i b i ui =*===≤≤+=+++≤≤-=?=+=+-+Ax LUx Ly d Ly d ====1112222(1)(n 1)(n 1)()()(n)(n)l y d m l y d l n y d m n l n y d ?????????????????????????=??????---?????????????????? …………… 一、实验题目 用追赶法解线性三对角方程组： ????? ???????--=????????????????????????322141141141124321x x x x 二、实验目的 1.熟悉掌握追赶法法的基本原理和基本方法。 2.学会用追赶法法解简单的方程组。 三、实验原理 （1）f Ly =，求y ；（2）y Ux =，求x ； 从而得到解三对角线方程组的追赶法公式。 1. 计算的递推公式： 111/b c =β )/(1--=i i i i i a b c ββ 2. 解：f Ly = 111/b f y = a a b y a f y i i i i i i i )/()(11----=β 3. 解：y Ux = n n y x = 1+-=i i i i x y x β 四、实验内容及结果 原始数据:a=[1,1,1]; b=[2,4,4,4]; c=[1,1,1]; d=[1,-2,2,-3]; 追赶法解三对角方程组： 程序源代码 %machase.m function x=machase(a,b,c,d) %用途：追赶法解三对角方程组Ax=d %格式：x= machase(a,b,c,d) a为次下对角线元素向量，b主对角元素 % 向量，c为次上对角线元素向量，d为右端向量，x返回解向量 n=length(a); for k=2:n b(k)=b(k)-a(k)/b(k-1)*c(k-1); d(k)=d(k)-a(k)/b(k-1)*d(k-1); end x(n)=d(n)/b(n); for k=n-1:-1:1 x(k)=(d(k)-c(k)*x(k+1))/b(k); end 输入：a=[1,1,1]; b=[2,4,4,4]; c=[1,1,1]; d=[1,-2,2,-3]; machase(a,b,c,d) 结果：ans = 0.9615 -0.9231 0.7308 五、实验结果分析 追赶法为一种特殊的LU分解法。追赶法是求解三对角矩阵的常用方法，但从整体编程角度分析，其程序编写较迭代法复杂，但通用性较好。追赶法求解三对角矩阵不但节省存储单元，而且可以减少计算量。 《数值分析》实验3 一．实验名称：追赶法 二、实验目的： （1）学会用追赶法解三对角方程组。 三、实验要求 （1） 按照题目要求完成实验内容 （2） 写出相应的C 语言程序 （3） 给出实验结果 （4） 写出相应的实验报告 四、实验题目 1. 追赶法解线性方程组 ?????? ? ??=??????? ????????? ??------500021121121124321x x x x 2．用追赶法解下列n 阶方程组 123141271411514115141151415n n x x x x x ---????????????--????????????--=??????????????????--??????--?????? ?? 取n=10和n=100。 #include { l[i][i-1]=a[i][i-1]/u[i-1][i-1]; u[i][i]=a[i][i]-l[i][i-1]*a[i-1][i]; } y[1]=b[1]; for(i=2;i<=n;i++) y[i]=b[i]-l[i][i-1]*y[i-1]; x[n]=y[n]/u[n][n]; for(i=n-1;i>=1;i--) x[i]=(y[i]-a[i][i+1]*x[i+1])/u[i][i]; for(i=1;i<=n;i++) printf(" x[%d]=%f\n",i,x[i]); } x[1]=8.705771 x[2]=7.823085 x[3]=7.586571 x[4]=7.523196 x[5]=7.506216 x[6]=7.501665 x[7]=7.500446 x[8]=7.500119 x[9]=7.500031 x[10]=7.500008 x[11]=7.500002 x[12]=7.500000 x[13]=7.500000 x[14]=7.500000 3）三对角形线性方程组 123456789104100000000141000000001410000000014100000000141000000001410000000014100000000141000000001410000000014x x x x x x x x x x -????????--????????--????--????????--????--????????--????--???????--??????-???? 7513261214455????????-?? ?? ??=??-?? ???? -?? ?????? ???-?? *(2,1,3,0,1,2,3,0,1,1)T x =--- 二、数学原理 设系数矩阵为三对角矩阵 1 122233111000000000 000000 n n n n n b c a b c a b A a b c a b ---?? ? ? ?= ? ? ? ? ?? ? 则方程组Ax=f 称为三对角方程组。 设矩阵A 非奇异，A 有Crout 分解A=LU ，其中L 为下三角矩阵，U 为单位上三角矩阵，记 1 122 233 1 10 00010 000 0001000 000100,00000000 00 0001n n n n b L U γαβγββγβ--???? ? ? ? ? ? ??== ? ? ? ? ? ? ? ? ? ??? ? ? ? 可先依次求出L ，U 中的元素后，令Ux=y ，先求解下三角方程组Ly=f 得出y ，再求解上三角方程组Ux=y 。 事实上，求解三对角方程组的2追赶法将矩阵三角分解的计算与求解两个三角方程组的计算放在一起，使算法更为紧凑。其计算公式为： 1111, 1111 ,111 ,2,3,,,1,2,,1i i i i i i i i i i i i i i n n i i i i c f b y i n c a b a f y y x y i n n x y x βγββαβγγβαβγ--+? ===?? =?? ?==-= ??? -?=?? =??=--?=-??对对(*) 三、程序设计 function x=chase(a,b,c,f) %求解线性方程组Ax=f,其中A 是三对角阵 %a 是矩阵A 的下对角线元素a(1)=0 %b 是矩阵A 的对角线元素 %c 是矩阵A 的上对角线元素c(n)=0 %f 是方程组的右端向量 n=length(f); x=zeros(1,n);y=zeros(1,n); d=zeros(1,n);u= zeros(1,n); %预处理 d(1)=b(1); for i=1:n-1 u(i)=c(i)/d(i); d(i+1)=b(i+1)-a(i+1)*u(i); end %追的过程 y(1)=f(1)/d(1); for i=2:n y(i)=(f(i)-a(i)*y(i-1))/d(i); end %赶的过程 x(n)=y(n); for i=n-1:-1:1 x(i)=y(i)-u(i)*x(i+1); end 实验追赶法解三对角方程组 一、实验目的 学会用追赶法解三对角方程组,并应用该算法于实际问题. 二、实验要求 给定三对角方程组，应用追赶法解得方程组的解。 三、实验内容 1、追赶法 2、以课本数值试验2为实例 3、如果有错，修改直至运行成功，查看运行结果； 四、实验环境 matlab 五、实验步骤和方法 1、程序设计 2、带入实例 3、撰写实验报告。 六、实验预习要求 得到实例的解 一、[源程序] function x = my_zgf2(A,d,flag) %MY_ZGF2 Summary of this function goes here [m,n]=size(A); %计算矩阵的大小 if nargin==2; %输入变量等于2的时候，A中储存所有元素的值for i=1:n a(i)=A(i+1,i); b(i)=A(i,i); c(i)=A(i,i+1); end a(1)=0; %补充不足的值 b(n)=A(n,n); c(n)=0; else c=[A(1,:) 0]; %flag==1时 b=A(2,:); a=[0 A(3,:)]; end u(1)=b(1); for i=2:n %第一次追赶，得到上、下三角矩阵 l(i)=a(i)/u(i-1); u(i)=b(i)-c(i-1)*l(i); end y(1)=d(1); %解Ly=d for i=2:n y(i)=d(i)-l(i)*y(i-1); end x(n)=y(n)/u(n); %解Ux=y for i=n-1:-1:1 x(i)=(y(i)-c(i)*x(i+1))/u(i); end 二、带入实例 A = -2.0000 -2.0000 -2.0000 -2.0000 -2.0000 -2.0000 -2.0000 0 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 -2.0000 -2.0000 -2.0000 -2.0000 -2.0000 -2.0000 -2.0000 0 d= 8.1400 0 0 0 0 0 0 0 >> d=A(4,:); my_zgf2(A,d,1) ans = 2.0350 1.0174 0.5086 0.2541 0.1267 0.0626 0.0298 0.0119 >> #include"iostream" #include"cmath" using namespace std; double a[100][100],l[100][100],u[100][100],b[100],y[100],x[100]; int n,s,i; void main(){ cout<<"请输入三对角阵的规模："< 实验课程名称计算机数值计算 实验项目名称追赶法 年级10级 专业信计 学生姓名成富 学号1007010167 理学院 实验时间：2012 年 3月 26 日 学生实验室守则 一、按教学安排准时到实验室上实验课，不得迟到、早退和旷课。 二、进入实验室必须遵守实验室的各项规章制度，保持室内安静、整洁，不准在室内打闹、喧哗、吸烟、吃食物、随地吐痰、乱扔杂物，不准做与实验内容无关的事，非实验用品一律不准带进实验室。 三、实验前必须做好预习（或按要求写好预习报告），未做预习者不准参加实验。 四、实验必须服从教师的安排和指导，认真按规程操作，未经教师允许不得擅自动用仪器设备，特别是与本实验无关的仪器设备和设施，如擅自动用或违反操作规程造成损坏，应按规定赔偿，严重者给予纪律处分。 五、实验中要节约水、电、气及其它消耗材料。 六、细心观察、如实记录实验现象和结果，不得抄袭或随意更改原始记录和数据，不得擅离操作岗位和干扰他人实验。 七、使用易燃、易爆、腐蚀性、有毒有害物品或接触带电设备进行实验，应特别注意规范操作，注意防护；若发生意外，要保持冷静，并及时向指导教师和管理人员报告，不得自行处理。仪器设备发生故障和损坏，应立即停止实验，并主动向指导教师报告，不得自行拆卸查看和拼装。 八、实验完毕，应清理好实验仪器设备并放回原位，清扫好实验现场，经指导教师检查认可并将实验记录交指导教师检查签字后方可离去。 九、无故不参加实验者，应写出检查，提出申请并缴纳相应的实验费及材料消耗费，经批准后，方可补做。 十、自选实验，应事先预约，拟订出实验方案，经实验室主任同意后，在指导教师或实验技术人员的指导下进行。 十一、实验室内一切物品未经允许严禁带出室外，确需带出，必须经过批准并办理手续。 追赶法求解三对角方程组 要求：对于给定的三对角系数矩阵和右端项，可以求 解线性代数方程组 一、 追赶法的数学理论 设系数矩阵为三对角矩阵 112223311100000000 000000 n n n n n b c a b c a b A a b c a b ---?? ? ? ?= ? ? ? ? ?? ? 则方程组Ax=f 称为三对角方程组。 设矩阵A 非奇异，A 有Crout 分解A=LU ，其中L 为下三角矩阵，U 为单位上三角矩阵，记 1 122 233 1 10 00010 000 0001000 000100,00000000 00 0001n n n n b L U γαβγββγβ--???? ? ? ? ? ? ??== ? ? ? ? ? ? ? ? ? ??? ? ? ? 可先依次求出L ，U 中的元素后，令Ux=y ，先求解下三角方程组Ly=f 得出y ，再求解上三角方程组Ux=y 。 事实上，求解三对角方程组的2追赶法将矩阵三角分解的计算与求解两个三角方程组的计算放在一起，使算法更为紧凑。其计算公式为： 1111, 1111 ,111 ,2,3,,,1,2,,1i i i i i i i i i i i i i i n n i i i i c f b y i n c a b a f y y x y i n n x y x βγββαβγγβαβγ--+?===?? =?? ?==-= ??? -?=?? =??=--?=-??对对(*) 二、 追赶法的算法和流程图 1.预处理 生成方程组的系数i u 及其除数i d ，事实上，按式(*)可交替生成i d 与i u ： 1d →1u →2d →…→1-n u →n d 其计算公式为 ??? ??-=-===+++1,...,2,1, ,/c u b 111 i i 11n i u a b d d d i i i i i 2.追的过程 顺序生成方程组右端： i y →2y →…→n y 据式(*)的计算公式为 n i d y a f y d f y i i i i i ,...,3,2, /)(/1111=?? ? -==- 作业五：编写用追赶法求解三对角线性方程组的标准程序，并求下列方程组的解。matlab 1233212417259x x x -???????????? =-?????????????????? function [ x,flag ] =zgf(A,b) % A 为方程组的系数矩阵 % b 为方程组的右端项 % x 为方程组的解 % flag 为指标向量，flag=‘failure ’表示失败，flag=‘OK ’表示成功 [n,m]=size(A);n=length(b); [L,U,flag]=LU_decom(A); for i=1:n if i==1 y(i)=b(i)/L(i,i); else y(i)=(b(i)-L(i,i-1)*y(i-1))/L(i,i); end end x(n)=y(n)/U(n,n); for k=n-1:-1:1 y(k)=y(k)-U(k,k+1)*x(k+1); x(k)=y(k)/U(k,k); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [L,U,flag]=LU_decom(A) [n,m]=size(A); if n~=m error('The rows and columns of matrix A must be equal!'); return; end L=eye(n); U=zeros(n); flag='OK'; for k=1:n for j=k:n z=0; for q=1:k-1 z=z+L(k,q)*U(q,j); end U(k,j)=A(k,j)-z; end if abs(U(k,k)) 2013-2014（1）专业课程实践论文题目：三对角方程组的追赶法 一、算法理论 在一些实际问题中，例如解常微分方程边值问题，解热传导方程以及船体数学放样中建立三次样条函数等，都会要求解系数矩阵为对角占优的三对角线方程组 11112222211111n n n n n n n n n b c x f a b c x f a b c x f a b x f -----????? ? ??? ? ??? ? ??? ?= ??? ? ??? ? ??? ?????? ? ， 简记为Ax f =。 求解Ax f =：等价于解两个三角形方程组 ,Ly f y =求；,Ux y x =求.从而得到解三对角线方程组的追赶法公式： (1)计算{}i β的递推公式 ()111/,/,2,3,,1;i i i i i c b c b a i n βββ==-=- (2) 解Ly f = ()()11111/,/,2,3,,;i i i i i i i y f b y f a y b a i n β--==--= (3)解Ux y = 1,,1,2,2,1.n n i i i i x y x y x i n n β+==-=-- 我们将计算系数 12112n n y y y βββ-→→→→→→ 及的过程称为追的过程， 将计算方程组的解 11n n x x x -→→→ 的过程称为赶的过程。 #include 计算方法与实习实验报告 4．编写用追赶法解三对角线性方程组的程序，并解下列方程组：（1） 2X1-X2=5 -X1+2X2-X3=-12 -X2+2X3-X4=11 -X3+2X4=-1 #include 计算方法实验报告 实验：追赶法解三对角矩阵方程 1.实验目的 ①熟悉matlab 软件的使用 ②掌握追赶法的基本方法 ③编程实现追赶法解三对角矩阵方程 2.实验内容 題目：用追赶法解方程组 ????????????????--------4114114114114????????????????54321x x x x x =??????? ?????????100200200200100 3.求解： 1.追赶法代码 function [x]=ZhuiGan(a,b,c,f) r=size(a); m=r(2); r=size(b); n=r(2); if size(a)~=size(b)|(m~=n-1)|size(b)~=size(f) error('变量不匹配,检查变量输入情况!'); end p=ones(1,m); Y=ones(1,n); x=Y; p(1)=a(1)/b(1); Y(1)=f(1)/b(1); t=0; for i=2:m t=b(i)-a(i-1)*p(i-1); p(i)=c(i)/t; Y(i)=(f(i)-a(i-1)*Y(i-1))/t; end Y(n)=(f(n)-a(n-1)*Y(n-1))/(b(n)-a(n-1)*p(n-1)); x(n)=Y(n);for i=n-1:-1:1 x(i)=Y(i)-p(i)*x(i+1); end 2.输入 >> a=-1*ones(1,4); >> c=a; >> b=4*ones(1,5); >> f=[100 200 200 200 100]; >> [x]=ZhuiGan(a,b,c,f) 4.结果 x =46.15384615384615 84.61538461538461 92.30769230769231 84.61538461538461 46.15384615384615 3总结 追赶法具有计算量少，方法简单，算法稳定等优点。 9.用追赶法解三对角方程组A x=b，其中 A= ( 2?10 ?1 2?1 0?12 ?1 0 0?1 2?1 0 0 0 ?1 2) ，b= ( 1 0) . 解：对于对角占优的三对角线性矩阵可用追赶法求解采用Fortran编程计算如下： program main implicit none integer::i,j,in,nb real::x real::a(1:100),b(1:100),c(1:100),f(1:100) real::d(1:100),y(1:100),v(1:100) open(unit=100,file='a.dat') open(unit=101,file='b.dat') open(unit=102,file='c.dat') open(unit=103,file='f.dat') open(unit=104,file='v.dat') do while(.not.eof(100)) read(unit=100,fmt=*)in,x a(in)=x end do nb=0 do while(.not.eof(101)) read(unit=101,fmt=*)in,x b(in)=x nb=nb+1 end do do while(.not.eof(102)) read(unit=102,fmt=*)in,x c(in)=x end do do while(.not.eof(103)) read(unit=103,fmt=*)in,x f(in)=x end do d(1)=c(1)/b(1) do i=2,nb-1 d(i)=c(i)/(b(i)-a(i)*d(i-1)) end do y(1)=f(1)/b(1) do i=2,nb y(i)=(f(i)-a(i)*y(i-1))/(b(i)-a(i)*d(i-1)) end do v(nb)=y(nb) do i=nb-1,1,-1 v(i)=y(i)-d(i)*v(i+1) end do write(unit=104,fmt=201)'y(i)',',','x(i)' do i=1,nb write(unit=104,fmt=200)y(i),',',v(i) end do 200 format(f7.4,a1,f7.4) 201 format(a7,a1,a7) end program 其中a.dat、b.dat、c.dat和f.dat的内容分别如图1、图2、图3和图4所示；计算结果保存在v.dat中，如图5所示： 图1 图2 图3 图4 图5 追赶法构造过程 追赶法仍然保持LU 分解特性,它是一种特殊的LU 分解。追赶法充分利用了系数矩阵的三对角特点，而且使之分解更简单,得到对三对角线性方程组的快速解法。 解出。 及可由时， 当，表示 ，则三对角方程的矩阵若记的计算公式，为： 和的及的元素于是得计算，，，有： 由矩阵乘法及相等定义y Ux d Ly LU A d Ax d d d d p b q n k c q a p b q q U p L n k c b p q a q p b q T n k k k k k k k k k i i i k k k k k k k k k =====-========+==--------),,,() ,,3,2() ,,3,2(211 111 1 111111 1 γγγγγ????????????????????????????????=????????????????----n n n n n n n n q q q p p p b a c b a c b a c b Doolittle 1221132111222 111111γγγ 分解形式 矩阵 2) 追赶法算法 1.输入变量个数n 、系数矩阵对应的三个向量a,b,c 、常数项b 2. For k=2,3,…,n 2.1 如果b k-1=0,则输出“追赶法失败”提示并终止 2.2 a k ? a k /b k-1 2.3 b k ?b k - a k *c k-1 2.4 d k ?d k - a k *d k-1 3.For k=n,n-1,…,1 3.1 x k ? (d k - c k *x k+1)/b k 时不能进行。 消元法的缺点，即当在消元法，因此也存 来源于空间，但是因为追赶法节省了计算时间和存贮程， 。追赶法的特殊求解过法次数仅有较简单，计算量、乘除追赶法。 组的方法亦称为追赶法用这组公式解线性方程，，， 分解的计算公式： 综合以上，求解出计算公式为： 0451 ,,1)(,,3,21 ,,2,1)(,,3,211111111111 1=-?????????-=-==-==-====--=-===-==+---+-k k k k k k n n n k k k k k k k k k k k k k k k k n n n k k k k q Gauss Gauss n n k q x c y x q y x y p d y n k c p b q q a p d y b q Doolittle n n k q x c y x q y x n k y p d y d y 追赶法求解三对角线性方程组 一 实验目的 利用编程方法实现追赶法求解三对角线性方程组。 二 实验容 1、 学习和理解追赶法求解三对角线性方程组的原理及方法； 2、 利用MATLAB 编程实现追赶法； 3、 举例进行求解，并对结果进行分。 三 实验原理 设n 元线性方程组Ax=d 的系数矩阵A 为非奇异的三对角矩阵 11222=(1)(n 1)()()a c b a c A a n c b n a n ??????????--?????? ………… 这种方程组称为三对角线性方程组。显然，A 是上下半宽带都是1的带状矩阵。设A 的前n-1个顺序主子式都不为零，根据定理2.5的推论，A 有唯一的Crout 分解，并且是保留带宽的。 其中L 是下三角矩阵，U 是单位上三角矩阵。利用矩阵相乘法，可以1112212(1)1u(n 1)()()1l u m l u A LU l n m n l n ????????????????==?????--????????????…………… 得到： 由上列各式可以得到L 和U 。 引入中间量y ，令y Ux =，则有： 已知 L 和d ，可求得y 。 则可得到y 的求解表达式： 11/1 2,3,,()(1)*y()=()[()(1)]/y d l i n m i y i li i di y i di m i y i li ==-+=--… 1111111/1(2)(1)(1)u (1)(11)/(1)(1)(1)l a l u c u c l mi bi i n a i m i i l i i n ci li ui ui ci li l i a i b i ui =*===≤≤+=+++≤≤-=?=+=+-+Ax LUx Ly d Ly d ====1112222(1)(n 1)(n 1)()()(n)(n)l y d m l y d l n y d m n l n y d ?????????????????????????=??????---?????????????????? ……………追赶法(经典计算)
追赶法解三对角程组
追赶法求解三对角线性方程组
追赶法
实验3 追赶法解三对角方程组 2
【良心出品】MATLAB 追赶法求解三对角方程组的算法原理例题与程序
追赶法解三对角矩阵
追赶法解三角矩阵
数值计算实验解三对角线性方程组的追赶法
(完整word版)追赶法求解三对角方程及其算例
追赶法求解三对角线性方程组
三对角方程组的追赶法
追赶法解线性方程组报告
计算方法追赶法解三对角矩阵方程
追赶法
追赶法构造过程
追赶法求解三对角线性方程组