实验一线性方程组的直接解法-武汉大学数学与统计学院.docVIP

实验一 线性方程组的直接解法 （武汉大学数学与统计学院 向华） 一、实习要点：掌握列选主元Gauss消去法，对称正定问题的Cholesky分解法。 二、基本内容回顾 1. Gauss消去法 直接法求解线性方程组，其中非奇异，。如果是上三角阵或下三角阵，则用回代容易求解。Gauss消去过程就是将问题化为易求解的上三角矩阵。定义如下矩阵（其中上标表示消去的步骤）： 。 定义乘子，并计算 。 第一列对角线以下元素化为0，系数矩阵变为 。 下一步消去第二列对角线以下元素；定义乘子，并计算 ； 变换后矩阵为 。 一般地，设前次消元后得到的系数矩阵为 ， 下一步需消去第列对角线以下元素。为此定义乘子，并计算 。 这样步消去后所得矩阵为 。 经过步消元后系数矩阵变为上三角阵，右端项也做相应的变换，最后原方程化为如下形式： ， 其中系数矩阵记为，下面亦以记之，表示其为上三角阵。 以上过程可以用矩阵变换简洁地表达。用乘子定义Gauss向量，以及Gauss变换矩阵。则Gauss消去过程可表示为： ， 亦即，其中为单位上三角阵。 利用Gauss变换矩阵的性质： 可得，且的具体形式为： 。 这一过程不中断，即主元不为零的充要条件是的阶顺序主子式非奇异。Gauss消去过程实际上给出了系数矩阵的LU分解。 由于数值稳定性，需要选主元策略。定义置换矩阵，左乘表示交换第与第行。列选主元的Gauss消去可以表示为： 。 利用置换矩阵的性质可得 。 令，，上式即为 ，所以，简记为。 2. 根平方法 对对称正定矩阵，存在唯一的对角元全为正的下三角矩阵，使，其中称作Cholesky因子。下面的公式基于Gaxpy计算第列元素()： 为了避免开方运算，可用改进的根平方法，即，其中是单位下三角矩阵，是对角元为正的对角阵。设，，对， 三、MATLAB语言描述 不选主元Gauss消去过程可用两行MATLAB代码描述： Cholesky分解过程用MATLAB描述为： 四、C代码示范 下面的C代码子程序gepp对系数矩阵作列选主元LU分解，并返回方程组的解。调用前参数a存储系数矩阵，参数b存储方程组右端项，参数N是矩阵维数。调用后，原系数矩阵对角线以下被单位下三角因子（不包括对角元素1）覆盖，其余存储上三角因子，右端项b被方程组的解所覆盖。 void gepp(double **a, double *b, int N) { // a-coefficient matrix, b-RHS, N-matrix size int i, j, jc, pv; double tm; for (j = 0; j N-1; j++) { // find the pivot pv = j; for (i = j; i N; i++) if (fabs(a[pv][j]) fabs(a[i][j])) pv = i; if (pv != j){ for (jc = j; jc N; jc++){ tm = a[j][jc]; a[j][jc] = a[pv][jc]; a[pv][jc] = tm; } tm = b[j]; b[j] = b[pv]; b[pv] = tm; } // multipliers, updating for (i = j+1; i N; i++) { a[i][j] /= a[j][j]; for (jc = j+1; jc N; jc++) a[i][jc] -= a[i][j]*a[j][jc]; b[i] -= a[i][j]*b[j]; } } // backward substitute b[N-1] = b[N-1]/a[N-1][N-1]; for (i = N-2; i = 0; i--) { for (j = i+1; j N; j++) b[i] -= a[i][j]*b[j]; b[i] /=a[i][i]; } } 下面的C代码子程序chol对对称正定矩阵作Cholesky分解，并返回方程组的解。调用前参数a存储系数矩阵，参数b存储方程组右端项，参数N是矩阵维数。调用后，原系数矩阵对角线及其以下存储Cholesky因子，右端项b被方程组的解覆盖。 void chol(double **a, double *b, int N) { // a-SPD matrix, b-RHS, N-matrix size int i, j, k; for (k = 0; k N; k++) { for (j = 0; j = k-1; j++) a[k][k] -= a[k]