《 数值线性代数 》 课件

1. 《数值线性代数》课件 河南师大数学与信息科学学院

2. 本章讨论解线性方程组 高斯(Gauss)消去法。 第三章 线性代数方程组解法 • 解下三角形方程组的前代法 • 解上三角形方程组的回代法 • Gauss变换 • 矩阵的三角分解法 • Gauss消去法 • 解三对角方程组的追赶法 • 解正定方程组的平方根法 • 矩阵的列主元三角分解法

3. 下三角形方程组形如 • 前代法的解析表达式 解下三角形方程组的前代法

4. 前代法的程序实现: void Lower(double**L,double *b,int n){ for(int i=0;i<n;i++){ for(j=1;j<i;j++) b[i]-=L[i][j]*b[j]; b[i]/=L[i][i]; } } 注意：程序中的下标起始值是从0开始的，这是C语言规则所致，组织数组L和b应当注意。

5. 前代法的消元改进 • 假定已将方程组化为下列同解形式（用增广矩阵表示） 运用消元方法，计算： 则方程组又同解于

6. 前代法程序2： void Lower(double**L,double*b,int n) { for(int i=0;i<n;i++){ b[i]/=L[i][i]; for(int j=i+1;j<n;j++) b[j]-=b[i]*L[j][i]; } }

7. 上三角方程组的回代法 上三角形方程组形如 • 前代法的解析表达式

8. 回代法的程序实现: void Uper(double**U,double *b,int n){ for(int i=n;i>0n;i--){ for(j=i+1;j<n;j++) b[i]-=U[i][j]*b[j]; b[i]/=U[i][i]; } } 注意：程序中的下标起始值是从0开始的，这是C语言规则所致，组织数组U和b应当注意。

9. 回代法的消元改进 假定已将方程组化为下列同解形式（用增广矩阵表示）

10. 回代法程序2： void Uper(double**U,double*b,int n) { for(int i=n-1;i>=0;i--){ b[i]/=U[i][i]; for(int j=0;j<i;j++) b[j]-=b[i]*U[j][i]; } }

11. Gauss变换 • 定义：形如 的单位下三角矩阵称为Gauss变换。它的向量表示形式为

12. Gauss变换的性质 • （1）变换性质

13. (2)Gauss变换的逆变换

14. 矩阵的LU分解

16. 这种计算三角分解的方法称作Gauss消去法.

17. Gauss消去法算法设计 • 根据Ｌ和Ｕ的形状以及它们的形成顺序，完全可以把计算结果保存在原矩阵Ａ的存储区域中，即 void Gauss(double**A,int n){ for (int i=0;i<n-1;i++) for (int j=i+1;j<n;j++){ A[j][i]/=A[i][i]; for (int k=i+1;k<n;k++) A[j][k]-=A[j][i]*A[i][k]; } }

18. 算法的理论分析 该定理从理论上给出了矩阵的三角分解能否继续进行的一个预先判定方法． 该定理矩阵的三角分解存在的充分必要条件． 注意：对一个矩阵进行三角分解，虽然可以预先判定三角分解是否存在．但在实际计算中，通常并不对其预先判定，而是进行每一步消元前，判定主元是否为零，若为零便退出算法．算法程序可改为

19. bool Gauss(double**A,int n) { for (int i=0;i<n-1;i++) { if(A[i][i]==0) return false; for (int j=i+1;j<n;j++) { A[j][i]/=A[i][i]; for (int k=i+1;k<n;k++) A[j][k]-=A[j][i]*A[i][k]; } } return true; }

20. 求解方程组的Gauss算法

21. 改进Gauss消去法算法摸拟

22. 解三对角方程组——追赶法 ( Crout分解 )

24. , 故有 (3.1)

25. 解 (3.2) 解 （3.3） (3.1) —(3.3) 叫追赶法，工作量小，非常有效。

26. 平方根法

28. 算法程序如下: • void Cholesky4(double **A,int n){ • for(int k=0;k<n;k++){ • for(int i=0;i<k;i++) • A[k][k]-=A[k][i]*A[k][i]; • A[k][k]=sqrt(A[k][k]); • for(i=k+1;i<n;i++){ • for(int j=0;j<k;j++) • A[i][k]-=A[k][j]*A[i][j]; • A[i][k]/=A[k][k]; • } • } • }

29. 用Gauss算法实现Cholesky分解 以四阶正定矩阵为例说明之.

32. 算法程序如下: void Cholesky(double **A,int n){ for(int k=0;k<n;k++){ A[k][k]=sqrt(A[k][k]); for(int i=k+1;i<n;i++) A[i][k]/=A[k][k]; for(i=k+1;i<n;i++) for(int j=i;j<n;j++) A[j][i]-=A[j][k]*A[i][k]; } }

33. 改进平方根算法(即LDLT) 算法程序如下页:

34. void Cholesky4(double **A,int n){ double *v=new double[n]; for(int j=0;j<n;j++){ for(int i=0;i<j;i++) v[i]=A[j][i]*A[i][i]; for(i=0;i<j;i++) A[j][j]-=A[j][i]*v[i]; for(i=j+1;i<n;i++){ for(int k=0;k<j;k++) A[i][j]-=A[i][k]*v[k]; A[i][j]/=A[j][j]; } } }

35. 列主元素三角分解 列主元消去法是对Gauss消去法的改进。这种改进一则可以保证非奇异矩阵的消去法可以顺利进去到底，二则在求解方程组时能提高数值稳定性。 列主元消去法就是在Gauss消去法的每步消元之前,先对包含对应未知量的剩余方程中选择未知量系数绝对值最大者为主元,然后将该方程调换到前面的位置上.这要用到初等置换阵 Pst用作行变换时,交换第s和第t行,列变换时,交换第s和第t列.

36. 列主元消去法的第K步消元的剖析

39. 求解线性方程组的列主元消去法 综合上述讨论,得到求解线性方程组Ax=b的列主元消去法下列算法: • 算法的实现 • 考虑算法实现时,需要考虑以下两点:

40. 选主元，并考虑置换阵的保存 • m=a[i][i];p[i]=I; • for(int j=i+1;j<n;j++) • if(fabs(m)<fabs(a[j][i])){ • m=a[j][i]; p[i]=j; } • 交换主元行 • if(p[i]!=i){ • for(int j=0;j<n;j++){ • m=a[i][j]; a[i][j]=a[p[i]][j]; a[p[i]][j]=m; } • } • 计算Pb • for(int i=0;i<n-1;i++){ • if(p[i]!=i){ • t=b[i];b[i]=b[p[i]];b[p[i]]=t; }

41. 算法程序 void Gauss(double**A,int *p,int n){ for (int i=0;i<n-1;i++){ //选主元 double m=a[i][i]; p[i]=i; for(int j=i+1;i<n;i++) if(fabs(m)<fabs(a[j][i])){ m=a[j][i]; p[i]=j; } if(A[i][i]==0) return false; if(p[i]!=i) //考虑行交换 for(j=0;j<n;j++){ m=a[i][j]; a[j][j]=a[p[i]][j]; a[p[i]][j]=m; } //消元过程 for (int j=i+1;j<n;j++){ A[j][i]/=A[i][i]; for (int k=i+1;k<n;k++) A[j][k]-=A[j][i]*A[i][k]; } } }