追赶法求解三对角线性方程组
追赶法求解三对角线性方程组
一 实验目的
利用编程方法实现追赶法求解三对角线性方程组。
二 实验内容
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 ?????????????????????????=??????---??????????????????
……………
由y Ux =得：
111112221u(n 1)(n 1)(n 1)1(n)(n)u x y u x y x y x y ?????????????????????????=??????---??????????????????
………… 可得到X 的求解表达式：
()()
1,2,,1
()()u()(1)
x n y n i n n x i y i i x i ==--=-+…
从而得到Ax=d 的解x 。
四 Matlab 编程
根据以上的实验原理，在Matlab 中编程如下函数x=Trid （A,d ）
在Matlab中新建Trid.m,其中程序为如上虚线框内代码，放在工
作目录。在Command Window输入以下语句：
clear all;clc;
fprintf('输入非奇异三对角系数矩阵A\n');
A=input('A matrix='); % 输入系数矩阵
fprintf('系数矩阵');A
if det(A)==0 % 判断系数矩阵是否奇异
fprintf('系数矩阵A奇异！！！请重新输入！\n');
fprintf('重新输入非奇异三对角系数矩阵A\n');
A=input('A matrix=');
fprintf('系数矩阵');A
end
fprintf('输入矩阵d\n');
d=input('d matrix=');
fprintf('矩阵d');d
x=Trid(A,d) % 调用Trid.m中的Trid函数进行求解
fid = fopen('Ax=d.txt', 'wt'); % 生成Ax=d.txt文件
fprintf(fid,'%s\r\n','利用三对角线追赶法求解Ax=d');
fprintf(fid,'%s\r\n','==================================');
for i=1:size(A,1)
fprintf(fid, '%.1f\t', A(i,:)); % 输出Ax=d，以上=为分隔符fprintf(fid, '%s', 'x');
fprintf(fid, '%u\t', i);
fprintf(fid, '%.1f\n', d(i));
end
fprintf(fid,'%s\r\n','==================================');
五 举例计算及分析
以课本（《数值分析》第4版，颜庆津，北京航空航天大学出版社）27页例3为例进行计算，输入系数矩阵A 和d ：
41141=4114A ????????????????
1411，10.5=132d ????????-???????? 调用x=Trid(A,d)后，并生成Ax=d.txt 文件，Command Window
同时也会输出解为0.20.2=0.50.80.3x ????????-????????
，与课本答案一致; Ax=d.txt 文件内容
如下图：
fprintf(fid,'%s\r\n','求解得到结果如下：'); for i=1:size(A,1)
fprintf(fid,'%s','x'); % 输出解向量x （i ）
fprintf(fid, '%u', i);
fprintf(fid, '%s\t', '='); fprintf(fid,'%.5f\r\n',x(i)); end fclose(fid)
再次对程序进行验证，输入矩阵如下：
12213=120.531A ???????????
?，21=13d ????-????????
利用matlab 计算Ax 得到Ax=d ，验证所得到的x 即为方程组的解。
由以上两组计算可表明，该程序能满足追赶法解三对角线性方程
组，其中要求系数矩阵满足要求，即非奇异，三对角且前（n-1）个顺序主子式都不为零。
在又换了一组系数矩阵后，系数矩阵如下：
A=[1 1 0 0 0;2 2 1 0 0;0 2 3 1 0;0 0 2 4 1;0 0 0 2 5]
即11221=23124125A ????????????????
，11=111d ???????????????? 运行程序后出错，查找问题发现A 进行Crout 分解时，
l1=1
m2=2
u1=1
l2=a2-m2*u1=2-2×1=0
再计算u2=c2/l2时无法计算，分析发现A 的2阶顺序主子式为零，无法进行Crout 分解，所以该程序还无法判断输入的系数矩阵的（n-1）阶顺序主子式是否为零，这是不足的地方，也无法判断系数矩阵是否是三对角矩阵，还需要在程序上进行修改完善。
六 实验总结
1112212(1)1u(n 1)()()1l u m l u A LU l n m n l n ????????????????==?????--????????????……………
通过本次的课程实验设计，对追赶法求解三对角线性方程组有了深刻的认识和了解，用matlab编程实现了该解法；同时在实验过程中明白了系数矩阵需要满足一定的要求才能使用该种追赶法方法。
MATLAB代码 解线性方程组的迭代法
解线性方程组的迭代法 1.rs里查森迭代法求线性方程组Ax=b的解 function[x,n]=rs(A,b,x0,eps,M) if(nargin==3) eps=1.0e-6;%eps表示迭代精度 M=10000;%M表示迭代步数的限制值elseif(nargin==4) M=10000; end I=eye(size(A)); n=0; x=x0; tol=1; %迭代过程 while(tol>eps) x=(I-A)*x0+b; n=n+1;%n为最终求出解时的迭代步数tol=norm(x-x0); x0=x; if(n>=M) disp('Warning:迭代次数太多，可能不收敛！'); return; end end 2.crs里查森参数迭代法求线性方程组Ax=b的解 function[x,n]=crs(A,b,x0,w,eps,M) if(nargin==4) eps=1.0e-6;%eps表示迭代精度 M=10000;%M表示迭代步数的限制值 elseif(nargin==5) M=10000; end I=eye(size(A)); n=0; x=x0; tol=1; %迭代过程 while(tol>eps) x=(I-w*A)*x0+w*b; n=n+1;%n为最终求出解时的迭代步数tol=norm(x-x0); x0=x;
if(n>=M) disp('Warning:迭代次数太多，可能不收敛！'); return; end end 3.grs里查森迭代法求线性方程组Ax=b的解 function[x,n]=grs(A,b,x0,W,eps,M) if(nargin==4) eps=1.0e-6;%eps表示迭代精度 M=10000;%M表示迭代步数的限制值 elseif(nargin==5) M=10000; end I=eye(size(A)); n=0; x=x0; tol=1;%前后两次迭代结果误差 %迭代过程 while(tol>eps) x=(I-W*A)*x0+W*b;%迭代公式 n=n+1;%n为最终求出解时的迭代步数tol=norm(x-x0); x0=x; if(n>=M) disp('Warning:迭代次数太多，可能不收敛！'); return; end end 4.jacobi雅可比迭代法求线性方程组Ax=b的解 function[x,n]=jacobi(A,b,x0,eps,varargin) if nargin==3 eps=1.0e-6; M=200; elseif nargin<3 error return elseif nargin==5 M=varargin{1}; end D=diag(diag(A));%求A的对角矩阵 L=-tril(A,-1);%求A的下三角阵
追赶法(经典计算)
一、算法理论 在一些实际问题中，例如解常微分方程边值问题，解热传导方程以及船体数学放样中建立三次样条函数等，都会要求解系数矩阵为对角占优的三对角线方程组 ??????? ? ??=???????? ?????????? ? ?-----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．用该程序计算三对角线方程组
求解线性方程组的直接解法
求解线性方程组的直接解法 5.2LU分解 ① Gauss消去法实现了LU分解 顺序消元结束时的上三角矩阵U和所用的乘数,严格下三角矩阵。 将下三角矩阵的对角元改成1,记为L,则有A=LU, 这事实是一般的,我们不难从消去的第k个元素时的矩阵k行及k列元素的 历史得到这一点.因为从消元的历史有 u kj=a kj-m k1u1j- m k2u2j -…- m k,k-1u k-1,j, j=k,k+1,…,n m ik=(a ik-m i1u1k- m i2u2k -…-m i,k-1u k-1,k>/u kk i=k+1,k+2,…,n 于是a kj=m k1u1j+m k2u2j+…+m k,k-1u k-1,j+u kj, j=k,k+1,…,n a ik=m i1u1k+m i2u2k+…+m i,k-1u k-1,k+m ik u kk i=k+1,k+2,…,n 从前面两个式子我们可以直接计算L和U(见下段>.将矩阵分解为单位下 三角矩阵和上三角矩阵之积称为矩阵的LU分解.顺序消元实现了LU分 解,同时还求出了g, Lg=b的解. ②直接LU分解 上段我们得到(l ij=m ij> u kj=a kj-l k1u1j-l k2u2j -…- l k,k-1u k-1,j, j=k,k+1,…,n l ik=(a ik-l i1u1k-l i2u2k -…-l i,k-1u k-1,k>/u kk i=k+1,k+2,…,n 2 诸元素对应乘积,只不过算L的元素时还要除以同列对角元.这一规律很 容易记住.可写成算法(L和U可存放于A>: for k=1:n-1 for j=k:n u kj=a kj-l k1u1j-l k2u2j -…- l k,k-1u k-1,j end for i=k+1:n l ik=(a ik-l i1u1k-l i2u2k -…-l i,k-1u k-1,k>/u kk end end 这一算法也叫Gauss消去法的紧凑格式,可一次算得L,U的元素,不需逐步 计算存储.
追赶法求解三对角线性方程组
追赶法求解三对角线性方程组 一 实验目的 利用编程方法实现追赶法求解三对角线性方程组。 二 实验内容 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分解法。追赶法是求解三对角矩阵的常用方法，但从整体编程角度分析，其程序编写较迭代法复杂，但通用性较好。追赶法求解三对角矩阵不但节省存储单元，而且可以减少计算量。
数值分析5-用Jacobi迭代法和Gauss-Seidel迭代法求解线性方程组
作业六：分别编写用Jacobi迭代法和Gauss-Seidel迭代法求解线性方程组Ax=B的标准程序，并求下列方程组的解。 可取初始向量 X(0) =(0，0，0)’； 迭代终止条件||x(k+1)-x(k)||<=10e-6 （1） = （2） = Jacobi迭代法： 流程图 开 始 判断b中的最大值 有没有比误差大 给x赋初值 进行迭代 求出x，弱到100次还没到，警告不收 结束
程序 clear;clc; A=[8,-1,1;2,10,01;1,1,-5]; b=[1;4;3]; e=1e-6; x0=[0;0;0]'; n=length(A); x=zeros(n,1); k=0; r=max(abs(b)); while r>e for i=1:n d=A(i,i); if abs(d)
程序结果（1）
（2）
Gauss-Seidel迭代法： 程序 clear;clc; %A=[8,-1,1;2,10,01;1,1,-5]; %b=[1;4;3]; A=[5,2,1;-1,4,2;2,-3,10]; b=[-12;20;3]; m=size(A); if m(1)~=m(2) error('矩阵A不是方阵'); end n=length(b); %初始化 N=0;%迭代次数 L=zeros(n);%分解A=D+L+U,D是对角阵，L是下三角阵，U是上三角阵U=zeros(n); D=zeros(n); G=zeros(n);%G=-inv(D+L)*U d=zeros(n,1);%d=inv(D+L)*b x=zeros(n,1); for i=1:n%初始化L和U for j=1:n if i
追赶法解三对角方程组
《数值分析》课程设计追赶法解三对角方程组 院（系）名称信息工程学院 专业班级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 日
直接法解线性方程组
直接法解线性方程组 实习题目： 仿照三对角方程组的追赶法解五对角方程组，其中系数矩阵为A,右端向量为：r。将A分解为LU。其中L为下三角，U为单位上三角。A为7*7阶的矩阵，其中对角元为4 5 6 7 8 9 10。上下次三角对角线元素为1 2 3 4 5 6 ；上下第二条对角线元素为1 2 3 4 5;右端项为：1 2 3 4 5 6 7. 要求：输出系数矩阵A，右端向量r,下三角矩阵L，单位上三角矩阵U,下三角矩阵Ly=b 的解向量y,单位上三角方程组Ux=y的解（即最终的解向量。保留七位小数。 实现方法：通过MATLAB编程实现。建立MATLAB脚本文件。 首先通仿照三对角方程组的追赶法得到五对角矩阵的实现算法。 然后又MATLAB编程实现。 实验结果（MATLAB截图）：
结果分析： 通过提供的计算数据得到最终的解向量x及中间过程产生的下三角矩阵L，单位上三角矩阵U,下三角矩阵Ly=b 的解向量y。 同时为了确保算法的正确性，我还通过MATLAB的左除运算检验得使用此算法的计算结果正确。 这里由于是用MATLAB，最终结果为分数形式，考虑到精确解一般比近似解更好，因此未化成七位小数形式。 算法实现分析： 首先计算L和U的元素。由于已知L和U的特定形式（及除了对角线和上下次对角线和上下第二条对角线外，其余为0。故通过矩阵的乘法即可得到LU中元素的计算公式。（具体算法见MATLAB程序） 算法优劣点：
1.解此题时看上去要用较多的存储单元，但实际上只需存储系数矩阵A的不为0的元素。 2.A分解为LU计算完成后，后续计算x和y的“追赶过程”运算量一般来说计算量比较小。 3.此题也可用之前的LU算法求解。但此处算法与一般的LU分解的解线性方程组的算法，相比计算量小了不少。 4.对于此处特定的对称的系数矩阵A，算法还可以进一步优化。 5.由于我在此算法中A.L U的各对角值均用一个列向量表示，一个缺点在于输出A，L,U时要重新组成矩阵形式。不过优点在于减少了存储单元。 6.另一缺点是，未能将结果封装成一个文件。 后附MATLAB代码： c=[4,5,6,7,8,9,10];d=[1,2,3,4,5,6,0];b=[0,1,2,3,4,5,6];e=[1,2,3,4,5,0,0];a=[0,0,1,2,3,4,5]; r=[1 2 3 4 5 6 7]; w=zeros(7,1);x=zeros(7,1);y=zeros(7,1);m=zeros(7,1);n=zeros(7,1);h=zeros(7,1); w(1)=c(1);m(1)=d(1)/c(1);n(1)=e(1)/c(1); h(2)=b(2);w(2)=c(2)-h(2)*m(1);m(2)=(d(2)-b(2)*n(1))/w(2);n(2)=e(2)/w(2); for k=3:5 h(k)=b(k)-a(k)*m(k-2); w(k)=c(k)-a(k)*n(k-2)-h(k)*m(k-1); m(k)=(d(k)-h(k)*n(k-1))/w(k); n(k)=e(k)/w(k); end h(6)=b(6)-a(6)*m(4); w(6)=c(6)-a(6)*n(4)-h(6)*m(5); m(6)=(d(6)-h(6)*n(5))/w(6); h(7)=b(7)-a(7)*m(5); w(7)=c(7)-a(7)*n(5)-h(7)*m(6); y(1)=r(1)/w(1);y(2)=(r(2)-h(2)*y(1))/w(2); for k=3:7 y(k)=(r(k)-a(k)*y(k-2)-h(k)*y(k-1))/w(k); end x(7)=y(7); x(6)=y(6)-x(7)*m(6);
求解线性方程组——超松弛迭代法(c)
求解线性方程组——超松弛迭代法 #include
cout<<"输入松弛因子w (1 Gauss-Seidel迭代法求解线性方程组 一． 问题描述 用Gauss-Seidel 迭代法求解线性方程组 由Jacobi 迭代法中，每一次的迭代只用到前一次的迭代值。使用了两倍的存储空间，浪费了存储空间。若每一次迭代充分利用当前最新的迭代值，即在计算第i 个分量 ) 1(+k i x 时，用最新分量 ) 1(1 +k x ， ???+) 1(2 k x ) 1(1 -+k i x 代替旧分量 ) (1 k x ， ???) (2 k x ) (1 -k i x ，可以起 到节省存储空间的作用。这样就得到所谓解方程组的Gauss-Seidel 迭代法。 二． 算法设计 将A 分解成U D L A --=，则b x =A 等价于b x =--U)D (L 则Gauss-Seidel 迭代过程 ) ()1()1(k k k Ux Lx b Dx ++=++ 故 ) ()1()(k k Ux b x L D +=-+ 若设1 )(--L D 存在，则 b L D Ux L D x k k 1)(1)1()()(--+-+-= 令 b L D f U L D G 11)()(---=-=， 则Gauss-Seidel 迭代公式的矩阵形式为 f Gx x k k +=+) () 1( 其迭代格式为 T n x x x x ) ()0()0(2)0(1)0(，，，???= (初始向量)， ) (1 1 1 1 1 )()1()1(∑∑-=-+=++--=i j i i j k j ij k j ij i ii i i x a x a b a x )210i 210(n k ???=???=，，，；，，， 或者 ?? ???--=???=???==?+=∑∑-=-+=+++) (1)210i 210(111 1)()1()1()()1(i j i i j k j ij k j ij i ii i i i k i k i x a x a b a x n k k x x x ，，，；，，， 三． 程序框图 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 习题 3.1 1. 求下列方阵的秩： (1)??? ?? ??--340313021201;(2)????? ??----174034301320;(3)??????? ? ?---------12433023221453334 311 ;(4)??????? ??------34732038234202173132. 2. 求下列方阵的逆矩阵： (1) ?? ? ?? ? ?323513123； (2) ????? ?? ??-----1210232112201023. 3. 解下列矩阵方程 (1) 设 ???? ? ??--=????? ??--=1322 31,113122214B A ,求X 使B AX =; (2) 设 ??? ? ??-=? ???? ??---=132 321,433312120B A ,求X 使B XA =； （3） ?? ??? ??-=????? ??-=????? ??-=112510324, 123011113,1120111111C B A ,求X 使C AXB =. 4. 求下列行列式 （1）? ? ? ??? ??????71 1 0251020214214 ；（2）????????????-260523211213 141 2；（3）?? ? ???????---ef cf bf de cd bd ae ac ab ； （4） ????????????---d c b a 100110011001. 5. 判断下列线性方程组解的情况，如果有唯一解，则求出解. ???????=+++-=----=+-+=+++;01123,2532,242,5)1(432143214 3214321x x x x x x x x x x x x x x x x ? ? ???????=+=++=++=++=+;15,065,065,065,165)2(545434323212 1x x x x x x x x x x x x x (3) ? ?? ??=-++=-+-=-+-;3222, 2353, 132432143214321x x x x x x x x x x x x (4) ?????=---=--+=+++.034,0222,022432143214321x x x x x x x x x x x x 习题 3.2 1. 用回代法解上三角形线性方程组 （1）??? ????==+-=-+=++;63,3,6333,8484443432321x x x x x x x x x （2）?? ???? ?-=-=+--=+--=-+.63,1032,92,9244343242 1x x x x x x x x x 2. 用回代法解下三角形线性方程组 迭代法解线性方程组作业 沈欢00986096 北京大学工学院，北京100871 2011年10月12日 摘要 由所给矩阵生成系数矩阵A和右端项b，分析系数矩阵A，并用Jacobi迭代法、GS迭代法、SOR(逐步松弛迭代法)解方程组Ax=b 1生成系数矩阵A、右端项b,并分析矩阵A 由文件”gr900900c rg.mm”得到了以.mm格式描述的系数矩阵A。A矩阵是900?900的大型稀 疏对称矩阵。于是，在matlaB中，使用”A=zeros(900,900)”语句生成900?900的零矩阵。再 按照.mm文件中的描述，分别对第i行、第j列的元素赋对应的值，就生成了系数矩阵A，并 将A存为.mat文件以便之后应用。 由于右端项是全为1的列向量，所以由语句”b=ones(900,1)”生成。 得到了矩阵A后，求其行列式，使用函数”det(A)”,求得结果为”Inf”,证明行列式太大，matlaB无法显示。由此证明，矩阵A可逆，线性方程组 Ax=b 有唯一解。 接着，判断A矩阵是否是对称矩阵(其实，这步是没有必要的，因为A矩阵本身是对称矩阵，是.mm格式中的矩阵按对称阵生成的)。如果A是对称矩阵，那么 A?A T=0 。于是，令B=A?A T,并对B求∞范数。结果显示： B ∞=0,所以，B是零矩阵，也就是：A是对称矩阵。 然后，求A的三个条件数： Cond(A)= A ? A?1 所求结果是，对应于1范数的条件数为：377.2334;对应于2范数的条件数为：194.5739;对应 于3范数的条件数为：377.2334; 1 从以上结果我们看出，A是可逆矩阵，但是A的条件数很大，所以，Ax=b有唯一解并且矩阵A相对不稳定。所以，我们可以用迭代方法来求解该线性方程组，但是由于A的条件数太大迭代次数一般而言会比较多。 2Jacobi迭代法 Jacobi迭代方法的程序流程图如图所示： 图1:Jacobi迭代方法程序流程图 在上述流程中，取x0=[1,1,...,1]T将精度设为accuracy=10?3，需要误差满足： error= x k+1?x k x k+1 数值分析实验报告 0 a 12 K a 1,n 1 K a 2,n 1 U O M 则有： 第一步： Jacobi 迭代法 a 1n a 2n M , 则有： A D L U a n 1,n Ax b A A x D b L U (D L U)x b Dx (L U)x b x D (L U)x D b 令 J D (L U) 则称 J 为雅克比迭代矩阵 f D b 由此可得雅克比迭代的迭代格式如下： x (0) , 初始向量 x (k 1) Jx (k) f ,k 0,1,2,L 第二步 Gauss-Seidel 迭代法 Ax b (D L U )x b (D L)x Ux b x (D L) Ux (D L) b A D L U a 11 a 12 L a 1n a 11 A a 21 a 22 L a 2n a 22 M MM MO a n1 a n2 L a nn a 11 得到 D a 22 O a nn 由 a 21 0 M M O a n 1,1 a n 1,2 L 0 a nn a n1 a n2 L a n,n a 21 L M M O a n 1,1 a n 1,2 L a n1 a n2 L a n,n 1 a 12 K a 1,n 1 a 1n 0 K a 2,n 1 a 2n O M M a n 1,n 10 令 G (D L) U ，则称G 为Gauss-Seidel 迭代矩阵 f (D L) b 由此可得 Gauss-Seidel 迭代的迭代格式如下： x (0) , 初始向量 第三步 SOR 迭代法 w0 AD L U 1 ( D 1 wL ((1 w)D wU )) (D 1 wL) ((1 w)D wU ) w w w 令M w 1 (D wL), N 1 ((1 w)D wU )则有：A MN w w Ax b AM L W N M (M N )x b Mx Nx b x M Nx M b N M, 令W f Mb 带入 N 的值可有 L W ((1 w)D wU) (D wL) 1((1 w)D wU) (D wL) f 1 b w 1(D wL) 1b 1 (D wL) w 称 L W 为 SOR 迭代矩阵，由此可得 SOR 迭代的迭代格式如下： x (0) ,初始向量 二、算法程序 Jacobi 迭代法的 M 文件： function [y,n]=Jacobi(A,b,x0,eps) %************************************************* %函数名称 Jacobi 雅克比迭代函数 %参数解释 A 系数矩阵 % b 常数项 % x0 估计解向量 x (k 1) Gx (k) f ,k 0,1,2,L (k 1) f,k 0,1,2,L 实验追赶法解三对角方程组 一、实验目的 学会用追赶法解三对角方程组,并应用该算法于实际问题. 二、实验要求 给定三对角方程组，应用追赶法解得方程组的解。 三、实验内容 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<<"请输入三对角阵的规模："< 一． 问题描述 用Gauss-Seidel 迭代法求解线性方程组 由Jacobi 迭代法中，每一次的迭代只用到前一次的迭代值。使用了两倍的存储空间，浪 费了存储空间。若每一次迭代充分利用当前最新的迭代值，即在计算第i 个分量) 1(+k i x 时， 用最新分量) 1(1 +k x ，???+) 1(2 k x ) 1(1 -+k i x 代替旧分量)(1k x ，???) (2 k x ) (1-k i x ，可以起到节省存储 空间的作用。这样就得到所谓解方程组的Gauss-Seidel 迭代法。 二． 算法设计 将A 分解成U D L A --=，则b x =A 等价于b x =--U)D (L 则Gauss-Seidel 迭代过程 ) ()1()1(k k k Ux Lx b Dx ++=++ 故 )()1()(k k Ux b x L D +=-+ 若设1 )(--L D 存在，则 b L D Ux L D x k k 1)(1)1()()(--+-+-= 令 b L D f U L D G 11)()(---=-=， 则Gauss-Seidel 迭代公式的矩阵形式为 f Gx x k k +=+)()1( 其迭代格式为 T n x x x x )()0()0(2)0(1)0(，，，???= (初始向量)， )(1111 1 )() 1()1(∑∑-=-+=++--=i j i i j k j ij k j ij i ii i i x a x a b a x )210i 210(n k ???=???=，，，；，，， 或者 ?? ???--=???=???==?+=∑∑-=-+=+++) (1)210i 210(111 1)() 1()1()()1(i j i i j k j ij k j ij i ii i i i k i k i x a x a b a x n k k x x x ，，，；，，， 三． 程序框图 实验课程名称计算机数值计算 实验项目名称追赶法 年级10级 专业信计 学生姓名成富 学号1007010167 理学院 实验时间：2012 年 3月 26 日 学生实验室守则 一、按教学安排准时到实验室上实验课，不得迟到、早退和旷课。 二、进入实验室必须遵守实验室的各项规章制度，保持室内安静、整洁，不准在室内打闹、喧哗、吸烟、吃食物、随地吐痰、乱扔杂物，不准做与实验内容无关的事，非实验用品一律不准带进实验室。 三、实验前必须做好预习（或按要求写好预习报告），未做预习者不准参加实验。 四、实验必须服从教师的安排和指导，认真按规程操作，未经教师允许不得擅自动用仪器设备，特别是与本实验无关的仪器设备和设施，如擅自动用或违反操作规程造成损坏，应按规定赔偿，严重者给予纪律处分。 五、实验中要节约水、电、气及其它消耗材料。 六、细心观察、如实记录实验现象和结果，不得抄袭或随意更改原始记录和数据，不得擅离操作岗位和干扰他人实验。 七、使用易燃、易爆、腐蚀性、有毒有害物品或接触带电设备进行实验，应特别注意规范操作，注意防护；若发生意外，要保持冷静，并及时向指导教师和管理人员报告，不得自行处理。仪器设备发生故障和损坏，应立即停止实验，并主动向指导教师报告，不得自行拆卸查看和拼装。 八、实验完毕，应清理好实验仪器设备并放回原位，清扫好实验现场，经指导教师检查认可并将实验记录交指导教师检查签字后方可离去。 九、无故不参加实验者，应写出检查，提出申请并缴纳相应的实验费及材料消耗费，经批准后，方可补做。 十、自选实验，应事先预约，拟订出实验方案，经实验室主任同意后，在指导教师或实验技术人员的指导下进行。 十一、实验室内一切物品未经允许严禁带出室外，确需带出，必须经过批准并办理手续。Gauss-Seidel迭代法求解线性方程组
【良心出品】MATLAB 追赶法求解三对角方程组的算法原理例题与程序
第三章 解线性方程组的直接方法
迭代法解线性方程组
实验解线性方程组的基本迭代法实验
追赶法解三对角矩阵
追赶法解三角矩阵
Gauss-Seidel迭代法求解线性方程组
数值计算实验解三对角线性方程组的追赶法