1、LU分解法和QR分解法。1、LU分解法:将A分解为一个下三角矩阵L和一个上三角矩阵U,即:A=LU,其中 L=, U=2、QR分解法:将A分解为一个正交矩阵Q和一个上三角矩阵R,即:A=QR三、实验原理解A*=b 的问题就等价于要求解两个三角形方程组: Ly=b,求y; U*=y,求*.设A为非奇异矩阵,且有分解式A=LU, L为单位下三角阵,U为上三角阵。L,U的元素可以有n步直接计算定出。用直接三角分解法解A*=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
2、,U*=y的计算公式;四、实验步骤 1将矩阵A保存进计算机中,再定义2个空矩阵L,U以便保存求出的三角矩阵的值。利用公式,将矩阵A分解为LU,L为单位下三角阵,U为上三角阵。2可知计算方法有三层循环。先通过公式计算出U矩阵的第一行元素和L矩阵的第一列元素。再根据公式和,和上次的出的值,求出矩阵其余的元素,每次都要三次循环,求下一个元素需要上一个结果。3先由公式 ,Ly=b求出y,因为L为下三角矩阵,所以由第一行开场求y.4再由公式,U*=y求出*, 因为U为上三角矩阵,所以由最后一行开场求*.五、程序流程图1、LU分解法2、QR分解法六、实验结果方程组1 :方程组2:方程组1: 方程组3:七、
3、实验总结为了求解线性方程组,我们通常需要一定的解法。其中一种解法就是通过矩阵的三角分解来实现的,属于求解线性方程组的直接法。在不考虑舍入误差下,直接法可以用有限的运算得到准确解,因此主要适用于求解中小型稠密的线性方程组。1、三角分解法 三角分解法是将A矩阵分解成一个上三角形矩阵U和一个下三角形矩阵L,这样的分解法又称为LU分解法。它的用途主要在简化一个大矩阵的行列式值的计算过程,求反矩阵和求解联立方程组。不过要注意这种分解法所得到的上下三角形矩阵并非唯一,还可找到数个不同 的一对上下三角形矩阵,此两三角形矩阵相乘也会得到原矩阵。 2、 QR分解法QR分解法是将矩阵分解成一个正规正交矩阵Q与上三
4、角形矩阵R,所以称为QR分解法。在编写这两个程序过程中,起初遇到不少麻烦!虽然课上教师反复重复着:“算法不难的,Its so easy!但是当自己实际操作时,感觉并不是则容易。毕竟是要把实际的数学问题转化为计算机能够识别的编程算法,所以在编写程序之前我们仔细认真的把所求解的问题逐一进展详细的分析,最终转化为程序段。每当遇到问题时,大家或许有些郁闷,但最终还是静下心来反复仔细的琢磨,一一排除了错误,最终完成了本次实验。回头一想原来编个程序其实也没有想象的则复杂,只要思路清晰,逐步分析,就可以慢慢搞定了。附 源代码:*include stdio.hmath.hstdlib.hconio.husin
5、g 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; dou
6、ble e*pct=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、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
8、; P = (int*)malloc(sizeof(int)*n); D = (double*)malloc(sizeof(double)*n); for (int i=0; in; i+) Pi=Di=0; coutflag; if(!flag) 矩阵LU分解: n lu(A,P,n);/矩阵LU分解 for(int j=0; j j+)printf(%ft,Ai);endl;矩阵LU分解的主元位置:Pima*) ma*= fabs(an*j + i) ; pivoti=j; / 对第i列进展行变换,使得主元在对角线上 if(pivoti!=i) for(j=i; j+)/ij与pivoti
9、j换 只用对上三角进展处理 temp=an*i + j; an*i + j=an*pivoti+ j; an*pivoti+ j=temp; j+)/Pi 局部下三角L an*j + i=an*j+i/an*i+i; j+)/计算上三角U for(k=i+1; k=0; i-) for(j=n-1; j j-) 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(size
10、of(double)*n); /获得tem值 m = 0 ; 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;=n-1; tem=tem + an*j +i*an*j +i; tem= sqrt(tem); an*j +i = an*j +i/tem ; / 调整矩阵k for(l=i; l l+) tem =tem + an*j + i*an*l + i*an*l + k; tempj = aj*n+k - 2*tem; 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; rem = 0;l rem = rem + qrl*n+i*qrj*n+i*bl; tempj = bj - 2*rem; bj = tempj;-1; for(l=n-1; l!=j;-l) bj =bj - bl*qrj*n+l; bj = bj /dj;
copyright@ 2008-2022 冰豆网网站版权所有
经营许可证编号:鄂ICP备2022015515号-1