线性代数方程组求解20985.pdf
线性代数方程组求解 一、实验要求 编程求解方程组：方程组 1：方程组 2：方程组 3：要求：用语言实现如下函数：实现矩阵的分解为输出参数中存放主元的位置排列函数成功时返回否则返回 求线代数方程组的解设矩阵为某个矩阵的分解 在内存中按行优先次序存放为分解的主元排列为方程组的右端向量 此函数计算方程组的解 并将结果存放在数组中 函数成功时返回否则返回矩阵的分解假设数组在内存中按行优先次序存放 此函数使用变换将其就地进行分解为输出参数中存放分解的上三角对角线元素求线代数方程组的解设矩阵为某个矩阵的分解 在内存中按行优先次序存放为分解的上三角对角线元素为方程组的右端向量函数计算方程组的解 并将结果存放在数组中函数成功时返回否则返回二、问题分析 求解线性方程组 Ax=b,其实质就是把它的系数矩阵 A 通过各种变换成一个下三角或上三角矩阵,从而简化方程组的求解;因此,在求解线性方程组的过程中,把系数矩阵 A 变换成上三角或下三角矩阵显得尤为重要,然而矩阵 A 的变换通常有两种分解方法：LU 分解法和 QR 分解法;1、LU 分解法：将 A 分解为一个下三角矩阵 L 和一个上三角矩阵 U,即：A=LU,其中 L=10010012121nnlll,U=nnnnuuuuuu0000022211211 2、QR 分解法：将 A 分解为一个正交矩阵 Q 和一个上三角矩阵 R,即：A=QR 三、实验原理 解 Ax=b 的问题就等价于要求解两个三角形方程组：Ly=b,求 y;Ux=y,求 x.设A为非奇异矩阵,且有分解式A=LU,L为单位下三角阵,U为上三角阵;L,U 的元素可以有 n 步直接计算定出;用直接三角分解法解 Ax=b 要求 A 的所有顺序主子式都不为零的计算公式：),2,1(niaulili,11/ualilil,i=2,3,n.计算 U 的第 r 行,L 的第 r 列元素 i=2,3,n：11rkkirkririulau ,i=r,r+1,n;rrrkkrikiriruulal/)(11,i=r+1,n,且 rn.求解 Ly=b,Ux=y 的计算公式；:,3,2,1111niylbybyikkikii .1,2,1,/)(,/1nniuxuyxuyxiinikkikiinnnn 四、实验步骤 1将矩阵 A 保存进计算机中,再定义 2 个空矩阵 L,U 以便保存求出的三角矩阵的值;利用公式,将矩阵 A 分解为 LU,L 为单位下三角阵,U 为上三角阵;2可知计算方法有三层循环;先通过公式计算出 U 矩阵的第一行元素liu 和 L 矩阵的第一列元素ill;再根据公式和,和上次的出的值,求出矩阵其余的元素,每次都要三次循环,求下一个元素需要上一个结果;3先由公式,Ly=b :,3,2,1111niylbybyikkikii 求出 y,因为 L 为下三角矩阵,所以由第一行开始求 y.4再由公式,Ux=y.1,2,1,/)(,/1nniuxuyxuyxiinikkikiinnnn 求出 x,因为 U 为上三角矩阵,所以由最后一行开始求 x.五、程序流程图 1、LU 分解法 开始输入系数矩阵A，常数项b及ndet1K=1,n-1,1调选列主元子程序i=k+1,n,1aikaik/akk(即求乘数mik)j=k+1,n,1aikaij-aik*akjjbibiaik*bkidetakkdetkbnbn/am(即求出xn)i=n-1,1,-1s0j=i+1,n,1ss+aij*bjjbi(bi-s)/aij输出b及det结束idetamndet回代过程消元过程 由主程序转来dakk,p ki=k+1,n,1|aik|d|daik,piid=0P=kj=k,n,1tapj,apj akj,akj tjtbp,bp bk,bk tdetdet返回主程序返回主程序输出失败标志，|A|=0det0结束选列主元=2、QR 分解法 开始输入数据 A，b矩阵 A的转置Householder 变换转置矩阵相乘输出上三角阵R=H3*H2*H1*A输入正交阵 Q=-（H3*H2*H1）T解上三角阵输出结果结束 六、实验结果 1、LU 分解法 方程组 1：方程组 2：方程组 3：2、QR 分解法 方程组 1：方程组 2：方程组 3：七、实验总结 为了求解线性方程组,我们通常需要一定的解法;其中一种解法就是通过矩阵的三角分解来实现的,属于求解线性方程组的直接法;在不考虑舍入误差下,直接法可以用有限的运算得到精确解,因此主要适用于求解中小型稠密的线性方程组;1、三角分解法 三角分解法是将 A 矩阵分解成一个上三角形矩阵 U 和一个下三角形矩阵 L,这样的分解法又称为 LU 分解法;它的用途主要在简化一个大矩阵的行列式值的计算过程,求反矩阵和求解联立方程组;不过要注意这种分解法所得到的上下三角形矩阵并非唯一,还可找到数个不同 的一对上下三角形矩阵,此两三角形矩阵相乘也会得到原矩阵;2、QR 分解法 QR分解法是将矩阵分解成一个正规正交矩阵Q与上三角形矩阵R,所以称为QR 分解法;在编写这两个程序过程中,起初遇到不少麻烦虽然课上老师反复重复着：“算法不难的,Its so easy”但是当自己实际操作时,感觉并不是那么容易;毕竟是要把实际的数学问题转化为计算机能够识别的编程算法,所以在编写程序之前我们仔细认真的把所求解的问题逐一进行详细的分析,最终转化为程序段;每当遇到问题时,大家或许有些郁闷,但最终还是静下心来反复仔细的琢磨,一一排除了错误,最终完成了本次实验;回头一想原来编个程序其实也没有想象的那么复杂,只要思路清晰,逐步分析,就可以慢慢搞定了;附 源代码:include include include include include using namespace std;bool ludouble a,int pivot,int n;6ft,Bi;coutnGuass 解法的误差为:n;for int i=0;in;i+ifi=3 coutendl;printf%.16ft,Bi-1;expct=expct+Bi-1;printfn 误差期望值为%.16ftn,expct/6;coutendl;else cout矩阵 QR 分解:n;qrA,D,n;6ft,Bi;coutnHouseholder 解法的误差为:n;for int i=0;in;i+ifi=3 coutendl;printf%.16ft,Bi-1;expct=expct+Bi-1;printfn 误差期望值为%.16ftn,expct/6;coutendl;return 0;bool ludouble a,int pivot,int n/矩阵 LU 分解 int i,j,k;double max,temp;max=0;temp=0;for i=0;in-1;i+/依次对第 i 列进行处理 /选出 i 列的主元,记录主元位置 max=fabsani+i;pivoti=i;forj=i+1;jmax max=fabsanj+i;pivoti=j;/对第 i 列进行行变换,使得主元在对角线上 ifpivoti=i forj=i;jn;j+/ij 与 pivotij 换 只用对上三角进行处理 temp=ani+j;ani+j=anpivoti+j;anpivoti+j=temp;forj=i+1;jn;j+/Pi 部分下三角 L anj+i=anj+i/ani+i;forj=i+1;jn;j+/计算上三角 U fork=i+1;kn;k+anj+k=anj+k-anj+iani+k;/计算下三角 L fori=0;in-2;i+/i 行 k 列 fork=i+1;kn-1;k+temp=anpivotk+i;anpivotk+i=akn+i;akn+i=temp;return false;bool guassdouble const lu,int const p,double b,int n/求线性代数方程组的解 int i,j,k;double temp;/按 qivot 对 b 行变换,与 LU 匹配 fori=0;in-1;i+/貌似错误在这里哦 temp=bpi;bpi=bi;bi=temp;/Ly=b,将 y 的内容放入 b fori=0;in;i+forj=0;j=0;i-forj=n-1;ji;j-bi=bi-luni+jbj;bi=bi/luni+i;return false;void qrdouble a,double d,int n/矩阵的 QR 分解 int i,j,l,k;double tem,m;double temp;temp=double mallocsizeofdoublen;for i=0;in-1;i+/依次对第 i 列进行处理 /获得 tem 值 m=0;forj=i;j0 m=-sqrtm;else m=sqrtm;/获得 temp 放入矩阵,并存主元 d tem=0;di=m;ani+i=ani+i-m;forj=i;j=n-1;j+tem=tem+anj+ianj+i;tem=sqrttem;forj=i;j=n-1;j+anj+i=anj+i/tem;/调整矩阵 fork=i+1;kn;k+forj=i;jn;j+tem=0;forl=i;ln;l+tem=tem+anj+ianl+ianl+k;tempj=ajn+k-2tem;forj=i;jn;j+ajn+k=tempj;dn-1=an-1n+n-1;bool householderdouble constqr,double constd,doubleb,int n/求线性代数方程组的解 int i,j,l;double rem;double temp;temp=double mallocsizeofdoublen;fori=0;in-1;i+forj=i;jn;j+rem=0;forl=i;ln;l+rem=rem+qrln+iqrjn+ibl;tempj=bj-2rem;forj=i;j-1;j-forl=n-1;l=j;-l bj=bj-blqrjn+l;bj =bj/dj;return false;
