用平方根法和追赶法求解线性方程组
摘要：在自然科学和工程技术中，很多问题的解决常常归结为求解线性方程组。例如，电学中的网络问题，解非线性方程组问题，应用有限元法解结构力学问题，船体数学放样中建立三次样条函数问题，解常微分方程边值问题等，都导致求解线性方程组。应用有限元法解结构力学问题时，其系数矩阵大多具有对称正定性质。平方根法就是利用对称正定矩阵的三角分解而得到的求解对称正定方程组的一种有效方法，平发根法数值稳定计算量也只有直接分解法的一半，目前在计算机上广泛应用平方根法解此类方程组。而解常微分方程边值问题，解船体数学放样中建立三次样条函数中，都会要求解系数矩阵为对角占优的三角方程组，追赶法就是利用矩阵的对角占优来求解三对角方程组的一种计算量较小的方法，且追赶法计算不会出现中间结果数量级的巨大增长和摄入误差的严重累积。本文通过编写程序用平方根法和追赶法求解线性方程组，还应用不同的接线性方程组的方法进行求解，并将结果进行了分析和比较。最后对课程设计进行了评价。
关键词：线性方程组；对称正定矩阵；对角占优的三对角矩阵；平方根法；追赶法；
解线性方程组的迭代法 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的下三角阵
求解线性方程组的直接解法 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的元素,不需逐步 计算存储.
在求解线性方程组时,直接解法有顺序高斯消元法、列主元高斯消元法、全主元高斯消元法、高斯约当消元法、消元形式的追赶法、LU分解法、矩阵形式的追赶法，当我们遇到对称正定线性方程组时，我们就要用到平方根法(对称LLT 分解法)来求解，为了熟悉和熟练运用平方根法求解线性方程组，下面对运用平方根法求解线性方程组进行解析。 一、运用平方根法求解线性方程组涉及到的定理及定义 我们在运用平方根法求解线性方程组时，要判定线性方程组Ax=b的系数矩阵A是否是对称正定矩阵，那么我们就要了解正定矩阵的性质和如下定理及定义： 1、由线性代数知，正定矩阵具有如下性质： 1) 正定矩阵A是非奇异的 2) 正定矩阵A的任一主子矩阵也必为正定矩阵 3) 正定矩阵A的主对角元素均为正数 4) 正定矩阵 A的特征值均大于零 5) 正定矩阵A的行列式必为正数 定义一线性方程组Ax=b的系数矩阵A是对称正定矩阵，那么Ax=b是对称正定线性方程组。 定义二如果方阵A满足A=AT，那么A是对称阵。 2.1.4 平方根法和改进的平方根法 如果A是n阶对称矩阵，由定理2还可得如下分解定理： 定理2 若A为n阶对称矩阵，且A的各阶顺序主子式都不为零，则A可惟一分解为：A＝LDLT，其中L为单位下三角阵，D为对角阵。 证明因为A的各阶顺序主子式都不为零，所以A可惟一分解为：A＝LU 因为，所以可将 U分解为： 其中 D为对角矩阵，U1为单位上三角阵．于是：A＝LDU1＝L(DU1) 因为A为对称矩阵，所以，A＝AT＝U1TDTLT＝U1T(DLT)，由 A的 LU分解的惟一性即得：L＝U1T，即U1＝LT，故A＝LDLT。 工程技术中的许多实际问题所归结出的线性方程组，其系数矩阵常有对称正定性，对于具有此类特殊性质的系数矩阵，利用矩阵的三角分解法求解是一种较好的有效方法，这就是对称正定矩阵方程组的平方根法及改进的平方根法，这种方法目前在计算机上已被广泛应用。 定理3 对称矩阵A为正定的充分必要条件是A的各阶顺序主子式大于零。 2 对称正定矩阵的三角分解 定理 (Cholesky分解)设A为n阶对称正定矩阵，则存在惟一的主对角线元素都是正数的下三角阵L，使得：A＝LLT。 分解式A＝LLT称为正定矩阵的Cholesky分解，利用Cholesky分解来求解系数矩阵为对称正定矩阵的方程组AX＝b的方法称为平方根法。 设A为4阶对称正定矩阵，则由定理 4知，A＝LLT，即： 将右端矩阵相乘，并令两端矩阵的元素相等，于是不难算得矩阵L的元素的计算公式为：
§5 平方根法 一、教学设计 1．教学内容：对称正定矩阵的Cholesky 分解法、三对角线矩阵分解的追赶法。 2．重点难点：Cholesky 分解法、追赶法。 3．教学目标：掌握对称正定矩阵的Cholesky 分解的计算过程，掌握三对角线矩阵分解的追赶法。 4．教学方法：讲授与讨论。 二、教学过程 §5 平方根法 在工程计算中，常遇到求解解对称再正定线性方程组问题，如应用有限元法解结构力学问题，应用差分方法解椭圆型偏微分方程等，最后都归结为求解系数矩阵为对称正定阵的线性方程组。根据系数矩阵的特殊性，是否有更好的解决方案（在存贮空间上的好处是显而易见的），算法上是否有所简化？ 5－0对称正定矩阵及性质复习 定义：设n n R A ?∈，如果A 满足条件 （1）A A T =；（2）对任意非零向量n R ∈x ，有0>x x A T ，则称A 为对称正定矩阵。 定理1 （对称正定矩阵的性质）如果n n R A ?∈为对称正定矩阵，则 （1）A 为非奇异阵，且1-A 亦是对称正定阵； （2）记k A 为A 的顺序主子阵，则k A 亦是对称正定阵),,2,1(n k =； （3）A 的特征值),,2,1(0)(n i A i =>λ； （4）A 的顺序主子式都大于零，即),,2,1(0)det(n k A k =>。 定理2 设n n R A ?∈为对称矩阵（判据）
（1）若A 的特征值),,2,1(0)(n i A i =>λ，则A 为对称正定矩阵； （2）若A 的顺序主子式都大于零，即),,2,1(0)det(n k A k =>，则A 为对称正定阵。 5－1 对称正定矩阵的三角分解 由前述定理 3.1知，若n 阶方阵A 的顺序主子式)1,,2,1 ()d e t (-=n k A k 均不为零，则A 有唯一的三角分解LU A =，其中L 为单位下三角阵，U 为上三角阵。n 阶对称正定阵A 的顺序主子式都大于零，当然有LU 分解，进一步地，此时U L ,之间有什么关系？这对解方程组有用处。由LU A L U A T T T ===及分解的唯一性，想到若U 的主对角元素皆为1，就有可能获得一些结果。为此，再将U 分解 DR u u u u u u u u u u u u u u u U n n nn nn n n ≡??? ?????? ???????? ?????? ??? ? ? ? ?=????????? ?? ?=111222********* 11222 11211 易知),,2,1(0n i u ii => （用k k k U L A ,,分别记矩阵U L A ,,的k 阶 顺序主子阵，容易验证k k k U L A =于是 ii k i i ii k i k k k k k k u a U U L U L A ∏∏ =======1 )(1det det det )det(det ） 于是LDR LU A ==，所以 A DR L LU DL R LDR A T T T T =====)()()(， 即 )()(DR L DL R A T T == 由分解的唯一性知：T R L =，R L T =，于是T LDL A = 自然地，若记
作业六：分别编写用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
直接法解线性方程组 实习题目： 仿照三对角方程组的追赶法解五对角方程组，其中系数矩阵为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);
笔算开平方方法 一． 拿出一个数，以小数点为分界，两位为一节，从最高位开始开平方。 我们就拿256吧 两位一节，先看最高的是2,那最大开平方就是1，写下1,剩余1。 第二步就是重点了！ 再取两个下来，也就是56。前面还有1,组合成156。 将第一次的开平方数1,先扩大20倍，得到２０,加上可以取的最大值，这个最大值是什么最大呢？也就是x*(20+x)<=156的最大x,可以取６,也正好是６,所以开平方的结果是１６。 再拿个比较大的数：15625 这个数，我们还是两位一节，看最高位１，那就写１，没剩余。 第二步：再取两个下来，也就是５６,我们先将１扩大２０倍，再用刚才的方法，取最大的x，可以取２，那就写２,剩余56-2*(20+2)=56-44=12 第三步：再取两个下来，也就是２５,和刚才剩余的12组成1225,那我们再对刚才的开平方数12，再扩大20倍，得到240,再求最大的开平方数，正好是５,没有剩余。 所以结果是125 如果有剩余，那小数点后也是两位两位地加，也就是一次加两个0,方法和前面一样，对前面已开出来的先扩大20倍，再取最大开方数，一直到你所要的准确度。 二． 1．将被开方数的整数部分从个位起向左每隔两位划为一段，用撇号分开（竖式中的11’56），分成几段，表示所求平方根是几位数； 2．根据左边第一段里的数，求得平方根的最高位上的数（竖式中的3）； 3．从第一段的数减去最高位上数的平方，在它们的差的右边写上第二段数组成第一个余数（竖式中的256）； 4．把求得的最高位数乘以20去试除第一个余数，所得的最大整数作为试商（20×3除256，所得的最大整数是4，即试商是4）； 5．用所求的平方根的最高位数的20倍加上这个试商再乘以试商．如果所得的积小于或等于余数，试商就是平方根的第二位数；如果所得的积大于余数，就把试商减小再试（竖式中（20×3+4）×4=256，说明试商4就是平方根的第二位数）； 6．用同样的方法，继续求平方根的其他各位上的数． 如遇开不尽的情况，可根据所要求的精确度求出它的近似值． 例如求的近似值（精确到0.01），可列出上面右边的竖式，并根据这个竖式得到 笔算开平方运算较繁，在实际中直接应用较少，但用这个方法可求出一个数的平方根的具有任意精确度的近似值． 实例 例如，A=5：5介于2的平方至3的平方;之间。我们取初始值2.1，2.2，2.3，2.4，2.5，2.6，2.7，2.8，2.9都可以，我们最好取中间值2.5。 第一步：2.5+（5/2.5-2.5)1/2=2.2；即5/2.5=2，2-2.5=-0.5，-0.5×1/2=-0.25，2.5+(-0.25)=2.25，取2位数2.2。 第二步：2.2+(5/2.2-2.2)1/2=2.23；即5/2.2=2.27272，2.27272-2.2=-0.07272，-0.07272×1/2=-0.03636，2.2+0.03636=2.23。取3位数2.23。 第三步： 2.23+(5/2.23-2.23)1/2=2.236。即5/2.23=2.2421525，,2.2421525-2.23=0.0121525，,0.0121525×1/2=0.00607，,2.23+0.006=2.236.，取4位数。每一步多取一位数。这个方法又叫反馈开方，即使你输入一个错误的数值，
求解线性方程组——超松弛迭代法 #include
cout<<"输入松弛因子w (1 追赶法求解三对角线性方程组 一 实验目的 利用编程方法实现追赶法求解三对角线性方程组。 二 实验内容 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 ?????????????????????????=??????---?????????????????? …………… 笔算开平方法的计算步骤如下： 1．将被开方数的整数部分从个位起向左每隔两位划为一段，用撇号分开，分成几段，表示所求平方根是几位数；小数部分从最高位向后两位一段隔开，段数以需要的精度+1为准。 2．根据左边第一段里的数，求得平方根的最高位上的数。 3．从第一段的数减去最高位上数的平方，在它们的差的右边写上第二段数组成第一个余数。 4．把求得的最高位数乘以20去试除第一个余数，所得的最大整数作为试商 5．用商的最高位数的20倍加上这个试商再乘以试商．如果所得的积小于或等于余数，试商就是平方根的第二位数；如果所得的积大于余数，就把试商减小再试，得到的第一个小于余数的试商作为平方根的第二个数。 6．用同样的方法，继续求平方根的其他各位上的数。 如遇开不尽的情况，可根据所要求的精确度求出它的近似值． 笔算开平方运算较繁，在实际中直接应用较少，但用这个方法可求出一个数的平方根的具有任意精确度的近似值． 我国古代数学的成就灿烂辉煌，早在公元前一世纪问世的我国经典数学著作《九章算术》里，就在世界数学史上第一次介绍了上述笔算开平方法．据史料记载，国外直到公元五世纪才有对于开平方法的介绍． 手工开根号法,只适用于任何一个整数或者有限小数开二次方. 因为网上写不出样式复杂的计算式,所以只能尽量书写,然后通过口述来解释: 假设一个整数1456456,开根号首先要从个位开始,每两位数做个标记,这里用'表示,那么标记后变成1'45'64'56.然后根据你要开的小数位数在小数点后补0,这里的举例开到整,则补2个0,(原因等明白该做法后自会理解),解法如下: 解法中需要说明的几个问题: 1,算式中的....没有意义,是因为网上无法排版,为了能把版式排得整齐点而加上的 2,为了区别小数点，所以小数点用。表示，而所有的.都是为了排版需要 3、除了1'45'64'56中的'有特殊意义，在解题中有用处外，其他的'都是为了排版和对起位置，说明数字来源而加的，取消没有任何影响 ...........1..2..0..6。8 .........----------------------- .....1../..1'45'64'56.00.. (1) (1) ............-------- .......22..|.45.. (2) (44) ..............-------- ........240.|.1'64.. (3) 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 ，，，；，，， 三． 程序框图 《数值分析》课程设计追赶法解三对角方程组 院（系）名称信息工程学院 专业班级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 日 习题 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 手动开平方方法（最新方案） 虽然现在开方可以直接用符号表示，但考试中如果出一道开方让你写数值的题目怎么办呢？在最新的数学研究中，有一种最新的开平方法。 如有下题： 1522756=（） 开方步骤如下： （一）分位 把一个平方数分为几段。 1.从最低位（个位）开始。 2.每两个数为一位。 3.最高位可以是一位数。 1522756分为：1|52|27|56 分位后，1522756被分为了4段，开方结果为四位数（这里是完全平方数，没有小数）（二）开方 开方运算和除法类似，每运算1次都有一个递减过程。运算时也是从高位至低位。 如1|52|27|56先算1，再算52…… 格式如下： 平方根 52 | |1 56 | 27 运算过程 和除法类似，平方根写在横线上面，运算过程写在下面。 平方定义，12=1 所以如下： 1 52 | |1 56 | 27 1 ——————— 5 2 这第一步与除法佷像，但是是一次落2位，也就是1段。 下面的运算就与除法有些差别了，这是计算中非常麻烦的部分。 这一步骤叫：造数 首先，将已开出的平方根部分×2，得到1×2=2 然后，我们须要假设下一个我们要开出的平方根是A，A的范围是0~9中任何一个自然数。下面就需要我们去试一试了，我们要在0~9中找出一个数作为A的值，前提是：要使前面一步算出的2与A合为一个新数，就是以A为个位，2为十位，合成2A（注意：这里不指2和A相乘，如果A=6，那么这个数为26），并且2A×A最接近而不超过前面落下的52。下一步就是试数，经试验A=2合适，也就得到22×2=44。 这一步的44就是结果了，下一位平方根为A，也就是2，得到： 解：由图可知a<0，b>0，a-b<0 ∴ () 2a b a b a b a b a =----=---+=- 其实平方根与立方根是可以笔算算出来的，当你身边没有计算机的时候，掌握此类的算法十分有用。 至于怎样算，可以归纳为如下两条公式：平方根，20m+n ；立方根， 300m^2+30mn+n^2。 怎样去理解呢，很简单。模板是按除法的模式。以开平方为例，譬如要求72162的平方根，先要从个位开始将它分块，每两位一块，即7，21，62这样分。然后开始试商，从最高为试起，先来7，什么数的平方小于7的呢？明显是2。然后用7减去2的平方，得出的数字3为余数，将要在下一步与后两位数字合起来用来进行下一步运算。第二步，此时被除的变成了321，此时公式开始派上用场，上一步试出来的商2即为m ，至于n 呢，当然是第二步要试的商啦，而除数就是公式20m+n ，切记商与除数的积不要大过被除数。具体到刚才的数字，除数是321，而被除数则是20×2+n，即40几，要n×（20×2＋n ）小于等于321，最合适的就是n=6，即46×6=276，再用321减去276得出结果45用于第三步的试商。第三步，也像第二步一样试商，只不过此时的被除数变成4562，除数m=20×26+n，n 是第三步要试的商。由n×（20×26+n）小于等于4562得出第三步的试商n=8，第四步开始棘手了，因为个位之前的已经试完了，此时，应从小数点之后的十分位开始，如一开始一样，每两位分成一块，这之后，就可以按前面的方法一直试下去了。 至于立方根，也是与平方根一样的思路，只不过比平方根复杂一点。与平方根的区别主要有三点，一、分块变为每三位一块，如刚才的72162，要分为72，162；二、除数变成300m^2+30mn+n^2；三、余数的区别，平方根的余数肯定要比除数小的，不然说明试的商不合适，例如上面的题目，第二步余数45小于除数46，第三步余数338小于除数528；而立方根就有点不同，它在第二步开始试商的时候，得出来的余数是有可能比除数大的，而且经实践得出，这可能性不低，至于到了第三步，余数又开始回归正常了，即必定小于除数，否则试商有误。 数值分析实验报告 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追赶法求解三对角线性方程组
笔算开平方法的计算步骤
Gauss-Seidel迭代法求解线性方程组
追赶法解三对角方程组
第三章 解线性方程组的直接方法
迭代法解线性方程组
手动开平方方法(最新方案)
算平方根的简便方法
实验解线性方程组的基本迭代法实验
相关主题
文本预览