线性代数方程组求解.doc
精品文档，仅供学习与交流，如有侵权请联系网站删除线性代数方程组求解一、实验要求编程求解方程组： 方程组1： 方程组2：方程组3：要求：用C/C+语言实现如下函数：1. bool lu(double* a, int* pivot, int n);实现矩阵的LU分解。pivot为输出参数，pivot0,n) 中存放主元的位置排列。函数成功时返回false，否则返回true。2. bool guass(double const* lu, int const* p, double* b, int n);求线代数方程组的解设矩阵Lunxn为某个矩阵anxn的LU分解，在内存中按行优先次序存放。p0,n)为LU分解的主元排列。b为方程组Ax=b的右端向量。此函数计算方程组Ax=b的解，并将结果存放在数组b0,n)中。函数成功时返回false，否则返回true。3. void qr(double* a, double* d, int n);矩阵的QR分解假设数组anxn在内存中按行优先次序存放。此函数使用HouseHolder变换将其就地进行QR分解。d为输出参数，d 0,n) 中存放QR分解的上三角对角线元素。4. bool hshld(double const*qr, double const*d, double*b, int n); 求线代数方程组的解设矩阵qrnxn为某个矩阵anxn的QR分解，在内存中按行优先次序存放。d 0,n) 为QR分解的上三角对角线元素。b为方程组Ax=b的右端向量。函数计算方程组Ax=b的解，并将结果存放在数组b0,n)中。函数成功时返回false，否则返回true。二、问题分析求解线性方程组Ax=b，其实质就是把它的系数矩阵A通过各种变换成一个下三角或上三角矩阵，从而简化方程组的求解。因此，在求解线性方程组的过程中，把系数矩阵A变换成上三角或下三角矩阵显得尤为重要，然而矩阵A的变换通常有两种分解方法：LU分解法和QR分解法。1、 LU分解法：将A分解为一个下三角矩阵L和一个上三角矩阵U，即：A=LU,其中 L=， U=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的所有顺序主子式都不为零）的计算公式： , ,i=2,3,，n.计算U的第r行，L的第r列元素(i=2,3,n)： ， i=r,r+1,n; , i=r+1,n,且rn.求解Ly=b，Ux=y的计算公式；四、实验步骤1>将矩阵A保存进计算机中，再定义2个空矩阵L，U以便保存求出的三角矩阵的值。利用公式，将矩阵A分解为LU，L为单位下三角阵,U为上三角阵。2>可知计算方法有三层循环。先通过公式计算出U矩阵的第一行元素 和L矩阵的第一列元素。再根据公式和，和上次的出的值，求出矩阵其余的元素，每次都要三次循环，求下一个元素需要上一个结果。3>先由公式 ，Ly=b求出y，因为L为下三角矩阵，所以由第一行开始求y.4>再由公式，Ux=y求出x, 因为U为上三角矩阵，所以由最后一行开始求x.五、程序流程图1、LU分解法2、QR分解法六、实验结果1、 LU分解法方程组1 ：方程组2：方程组3：2、 QR分解法方程组1：方程组2：方程组3：七、实验总结为了求解线性方程组，我们通常需要一定的解法。其中一种解法就是通过矩阵的三角分解来实现的，属于求解线性方程组的直接法。在不考虑舍入误差下，直接法可以用有限的运算得到精确解，因此主要适用于求解中小型稠密的线性方程组。1、三角分解法 三角分解法是将A矩阵分解成一个上三角形矩阵U和一个下三角形矩阵L，这样的分解法又称为LU分解法。它的用途主要在简化一个大矩阵的行列式值的计算过程，求反矩阵和求解联立方程组。不过要注意这种分解法所得到的上下三角形矩阵并非唯一，还可找到数个不同 的一对上下三角形矩阵，此两三角形矩阵相乘也会得到原矩阵。 2、 QR分解法QR分解法是将矩阵分解成一个正规正交矩阵Q与上三角形矩阵R,所以称为QR分解法。 在编写这两个程序过程中，起初遇到不少麻烦！虽然课上老师反复重复着：“算法不难的，It's so easy!”但是当自己实际操作时，感觉并不是那么容易。毕竟是要把实际的数学问题转化为计算机能够识别的编程算法，所以在编写程序之前我们仔细认真的把所求解的问题逐一进行详细的分析，最终转化为程序段。每当遇到问题时，大家或许有些郁闷，但最终还是静下心来反复仔细的琢磨，一一排除了错误，最终完成了本次实验。回头一想原来编个程序其实也没有想象的那么复杂，只要思路清晰，逐步分析，就可以慢慢搞定了。附 源代码:#include <iostream>#include <stdio.h>#include <math.h> #include <stdlib.h>#include <conio.h> using namespace std;bool lu(double* a, int* pivot, int n);/矩阵LU分解bool guass(double const* lu, int const* p, double* b, int n);/求线性代数方程组的解void qr(double* a, double* d, int n); /矩阵的QR分解bool householder(double const*qr, double const*d, double*b, int n); int main() int n=0; int temp=0; bool flag = false; double expct=0;/误差期望值 double devsq=0;/误差的方差 int * P= NULL; double * D= NULL; double A= 1, 1/2.0, 1/3.0, 1/4.0, 1/5.0, 1/6.0, 1/2.0, 1/3.0, 1/4.0, 1/5.0, 1/6.0, 1/7.0, 1/3.0, 1/4.0, 1/5.0, 1/6.0, 1/7.0, 1/8.0, 1/4.0, 1/5.0, 1/6.0, 1/7.0, 1/8.0, 1/9.0, 1/5.0, 1/6.0, 1/7.0, 1/8.0, 1/9.0, 1/10.0, 1/6.0, 1/7.0, 1/8.0, 1/9.0, 1/10.0,1/11.0 double B= 1+ 1/2.0+ 1/3.0+ 1/4.0+ 1/5.0+ 1/6.0, 1/2.0+ 1/3.0+ 1/4.0+ 1/5.0+ 1/6.0+ 1/7.0, 1/3.0+ 1/4.0+ 1/5.0+ 1/6.0+ 1/7.0+ 1/8.0, 1/4.0+ 1/5.0+ 1/6.0+ 1/7.0+ 1/8.0+ 1/9.0, 1/5.0+ 1/6.0+ 1/7.0+ 1/8.0+ 1/9.0+ 1/10.0, 1/6.0+ 1/7.0+ 1/8.0+ 1/9.0+ 1/10.0+1/11.0 n = 6 ; P = (int*)malloc(sizeof(int)*n); D = (double*)malloc(sizeof(double)*n); for (int i=0; i<n; i+) Pi=Di=0; cout<<"Please choose the way to solve (0: lu,1:qr):" cin>>flag; if(!flag) cout<<"矩阵LU分解: n" lu(A,P,n);/矩阵LU分解 for (int i=0; i<n; i+) for(int j=0; j<n; j+) printf("%ft",Ai); cout<<endl; cout<<"矩阵LU分解的主元位置: n" for (int i=0; i<n; i+) cout<<Pi<<"t" cout<<"nGuass求线性代数方程组求解:n" guass(A,P,B,n);/求线性代数方程组的解 for (int i=0; i<n; i+) if(i=3) cout<<endl; printf("%.16ft",Bi); cout<<"nGuass解法的误差为:n" for (int i=0; i<n; i+) if(i=3) cout<<endl; printf("%.16ft",Bi-1); expct=expct + Bi-1; printf("n误差期望值为%.16ftn",expct/6); cout<<endl; else cout<<"矩阵QR分解: n" qr(A,D,n);/矩阵qr分解 for (int i=0; i<n; i+) for(int j=0; j<n; j+) printf("%dt",Ai); cout<<endl; cout<<"QR分解的上三角矩阵对角线元素: n" for (int i=0; i<n; i+) printf("%ft",Di); cout<<"nHouseholder线性代数方程组求解:n" householder(A,D,B,n);/求线性代数方程组的解 for (int i=0; i<n; i+) if(i=3) cout<<endl; printf("%.16ft",Bi); cout<<"nHouseholder解法的误差为:n" for (int i=0; i<n; i+) if(i=3) cout<<endl; printf("%.16ft",Bi-1); expct=expct + Bi-1; printf("n误差期望值为%.16ftn",expct/6); cout<<endl; return 0 ; bool lu(double* a, int* pivot, int n)/矩阵LU分解 int i,j,k; double max,temp; max = 0; temp = 0; for (i=0; i<n-1; i+)/依次对第i列进行处理 / 选出i列的主元,记录主元位置 max = fabs(an*i + i); pivoti=i; for(j=i+1; j<n; j+)/j表示第j行 if( fabs(an*j + i)>max) max= fabs(an*j + i) ; pivoti=j; / 对第i列进行行变换，使得主元在对角线上 if(pivoti!=i) for(j=i; j<n; j+)/ij与pivotij换 只用对上三角进行处理 temp=an*i + j; an*i + j=an*pivoti+ j; an*pivoti+ j=temp; for(j=i+1; j<n; j+)/Pi 部分下三角L an*j + i=an*j+i/an*i+i; for(j=i+1; j<n; j+)/计算上三角U for(k=i+1; k<n; k+) an*j + k=an*j+k- an*j+i*an*i+k; /计算下三角 L for(i=0; i<n-2; i+)/i行k列 for(k=i+1; k<n-1;k+) temp=an*pivotk + i; an*pivotk + i=ak*n + i; ak*n + i=temp; return false ; bool guass(double const* lu, int const* p, double* b, int n)/求线性代数方程组的解 int i,j,k; double temp; /按qivot对b行变换，与LU匹配 for(i=0; i<n-1; i+) /貌似错误在这里哦 temp = bpi; bpi = bi; bi=temp; /Ly=b,将y的内容放入b for(i=0; i<n; i+) for(j=0; j<i; j+) bi=bi-lun*i+j*bj; /Uy=x,将x的内容放入b for(i=n-1; i>=0; i-) for(j=n-1; j>i; j-) bi=bi-lun*i+j*bj; bi=bi/lun*i+i; return false; void qr(double* a, double* d, int n) /矩阵的QR分解 int i,j,l,k; double tem,m; double *temp; temp = (double *)malloc(sizeof(double)*n); for (i=0; i<n-1; i+)/依次对第i列进行处理 /获得tem值 m = 0 ; for(j=i; j<n; j+)/j表示第j行 m = m +an*j + i*an*j + i ; if(an*i + i >0) m = -sqrt(m); else m = sqrt(m); /获得temp放入矩阵，并存主元d tem = 0 ; di = m ; an*i +i = an*i +i - m; for(j=i; j<=n-1; j+) tem=tem + an*j +i*an*j +i; tem= sqrt(tem); for(j=i; j<=n-1; j+) an*j +i = an*j +i/tem ; / 调整矩阵 for(k=i+1;k<n;k+) for(j=i; j<n; j+) tem = 0 ; for(l=i; l<n; l+) tem =tem + an*j + i*an*l + i*an*l + k; tempj = aj*n+k - 2*tem; for(j=i; j<n; j+) aj*n+k = tempj; dn-1 = a(n-1)*n+n-1; bool householder(double const*qr, double const*d, double*b, int n)/求线性代数方程组的解 int i,j,l; double rem; double *temp; temp = (double *)malloc(sizeof(double)*n); for(i=0; i<n-1; i+) for(j=i; j<n; j+) rem = 0; for(l=i;l<n; l+) rem = rem + qrl*n+i*qrj*n+i*bl; tempj = bj - 2*rem; for(j=i; j<n; j+) bj = tempj; for(j=n-1; j>-1; j-) for(l=n-1; l!=j;-l) bj =bj - bl*qrj*n+l; bj = bj /dj; return false;【精品文档】第 10 页
