计算方法数值分析C语言源程序.docxVIP

第1章 线性方程组的直接算法 求解线性方程组的直接算法是基于矩阵分解的算法。常见的矩阵分解有两种： １．矩阵的三角分解 矩阵的三角就是把一个矩阵分解成两个三角形矩阵的乘积。比如： 简单的三角分解：，这里，是单位下三角矩阵，是上三角矩阵。 列主元三角分解：，这里，是初等置换阵，是单位下三角矩阵且各元素的模不超过1，是上三角矩阵。 全主元三角分解：，这里，，是初等置换阵，是单位下三角矩阵且元素的模不超过，是上三角矩阵。 ２．矩阵的正交三角分解 正交化三角化就是把一个矩阵分解成一个正交矩阵和一个上三角矩阵的乘积．即 ，这里 是正交矩阵，是上三角矩阵． 1.1 矩阵的三角分解 1.1.1 功能 把实矩阵分解成单位下三角形矩阵和上三角形矩阵的乘积．即 ． 该算法适用于各阶顺序主子式不等于的矩阵． 1.1.2 算法概述 所谓三角分解就是把阶方阵作如下分解 其中是单位下三角矩阵，上三角矩阵。 当时，构造Gauss变换 则 以此类推，只要对角线上的元素，便可以一直这样做下去。直到将其化为上三角矩阵为止。即有 ． 其中 ，且 从而有 而且 ． 综上所述，我们可以将两个矩阵因子继续存储在原矩阵的存储空间上．的主对角线上的１不予存储． 算法5.3（计算三角分解：Gauss消去法） 1.1.3 算法程序 1.1.3.1.1 参数说明 **a n阶矩阵； n 矩阵的阶； 1.1.3.1.2 Ｃ程序 bool GaussLU(double **a,int n)//n阶矩阵的LU分解 { for(int k=0;kn-1;k++) { if(a[k][k]==0) return false; for(int i=k+1;in;i++) { a[i][k]/=a[k][k]; for(int j=k+1;jn;j++) a[i][j]-=a[i][k]*a[k][j]; } } return true; } 1.1.4 例题 求矩阵 1 2 3 4 1 4 9 16 1 8 27 64 的三角分解，结果如下： 1 2 3 4 1 2 6 12 1 3 6 24 1 7 6 24 　 1.2 列主元三角分解 1.2.1 功能 用矩阵的列主元三角分解，分解矩阵： ， 这里，是初等置换阵，是单位下三角矩阵且各元素的模不超过1，是上三角矩阵。 　 1.2.2 算法概述 列主元三角分解法和普通三角分解法基本上类似，所不同的是在构造Gauss变换前，先在对应列中选择绝对值最大的元素（称为列主元），然后实施初等行交换将该元素调整到矩阵对角线上。 例如第步变换叙述如下： 选主元：确定使； 行交换：将矩阵的第行和第行上的元素互换位置。即． 实施Gauss变换：通过初行变换，将列主对角线以下的元素消为零．即 算法5.2（计算三角分解：列主元Gauss消去法） 1.2.3 程序设计 1.2.3.1 参数说明 **a：n阶矩阵； *p：记忆分解过程中进行的行交换； n：矩阵的阶． 1.2.3.2 Ｃ程序 int RowGaussLU(double **a,int *p,int n)//n阶矩阵的列主元LU分解 { int i,j, k; double max; for(i=0;in;i++) p[i]=i; for(i=0;in-1;i++) { max=a[i][i]; p[i]=i; for(j=i+1;jn;j++) if(fabs(a[j][i])fabs(max)) { max=a[j][i]; p[i]=j; } if(max==0) return i-1; else if(p[i]!=i) { for(j=0;jn;j++) swap(a[i][j],a[p[i]][j]); } for(j=i+1;jn;j++) a[j][i]/=a[i][i]; for(j=i+1;jn;j++) for(k=i+1;kn;k++) a[j][k]-=a[i][k]*a[j][i]; } if(a[n-1][n-1]==0) return n-1; else return n; } 1.2.4 例题 用列主元高期分解矩阵 结果， P=(0,3,2,3) 1.3 全主元三角分解 1.3.1 功能 将矩阵进行如下三角分解 ， 其中：是初等置换阵，单位下三角阵（各元素绝对值不超过１），上三角阵． 1.3.2 算法概述 全主元三角分解法和普通三角分解法基本上类似，所不同的是在构造Gauss变换前，先在对应列中选择绝对值最大的元素（称为列主元），然后实施初等行交换将该元素调整到矩阵对角线上。 例如第步变换叙述如下： 选主元：确定，使； 行列交换：将矩阵的第行和第行上的元素互换位置，再将第列和第列上的元素互换。即；． 实施Gauss变换：通过初行变换，将列主对角线以下的元素消