用 C 语言求解 N 阶线性矩阵方程 Ax=b 的简单解法一、描述问题：题目：求解线性方程组 Ax=b，写成函数。其中，A 为 nn 的 N 阶矩阵，x 为需要求解的 n元未知数组成的未知矩阵，b 为 n 个常数组成的常数矩阵。即运行程序时的具体实例为：转化为矩阵形式（为检验程序的可靠性，特意选取初对角线元素为 0 的矩阵方程组）即为：二、分析问题并找出解决问题的步骤：由高等代数知识可知，解高阶线性方程组有逆矩阵求解法、增广矩阵求解法等，而在计算机 C 语言中，有高斯列主消元法、LU 分解法、雅克比迭代法等解法。为了与所学的高等代数知识相一致，选择使用“高斯简单迭代消元法”，与高等代数中的“增广矩阵求解法”相一致。以下简述高斯消元法的原理：算法基本原理：首先，为了能够求解 N 阶线性方程组（N 由用户输入），所以需要定义一个大于 N 维的数组 adim+1dim+1（dim 为设定的最大维数，防止计算量溢出），当用户输入的阶数 N 超过设定值时提示重启程序重新输入。进而，要判断方程组是否有解，无解提示重启程序重新输入，有解的话要判断是有无数不定解还是只有唯一一组解，在计算中，只有当原方程组有且只有一组解时算法才有意义，而运用高等代数的知识，只有当系数矩阵对应的行列式 |A|0 时，原方程组才有唯一解，所以输入系数矩阵后要计算该系数矩阵的行列式 |A|（定义了 getresult(n)函数计算），当行列式 |A|=0 时同样应提示重启程序重新输入， |A|0 时原方程组必然有且仅有唯一一组解。判断出方程组有且仅有唯一一组解后，开始将系数矩阵和常数矩阵（合并即为增广矩阵）进行初等行变换（以 a11 为基元开始，将第 j 列上 j 行以下的所有元素化为 0），使系数矩阵转化为上三角矩阵。这里要考虑到一种特殊情况，即交换到第 j-1 列后，第 j 行第 j列元素 ajj=0 ，那此时不能再以 ajj 为基元。当变换到第 j 列时，从 j 行 j 列的元素 ajj 以下的各元素中选取第一个不为 0 的元素，通过第三类初等行变换即交换两行将其交换到 ajj 的位置上，然后再进行消元过程。交换系数矩阵中的两行，相当于两个方程的位置交换了。再由高斯消元法，将第 j 列元素除 ajj 外第 j 行以下的其他元素通过第二种初等行变换化为 0，这样，就能使系数矩阵通过这样的行变换化为一个上三角矩阵，即，当系数矩阵 A 进行初等行变换时，常数矩阵也要进行对应的初等行变换，即此时 那么有接下来，进行“反代”，由 可求出 ，再往上代入 即可求出 以此类推，即可从 xn推到 xn-1 ，再推到 xn-2 直至 x1 。至此，未知矩阵 x 的所有元素就全部求出，即求出了原方程组有且仅有的唯一一组解。基本原理示意图：三、编写程序1. #include2. #include3. #include4. #define dim 10 /定义最大的维数 10，为防止计算值溢出5. double adim+1dim+1,bdim+1,xdim+1; /定义双精度数组6. double temp;7. double getarray(int n); /定义输入矩阵元素的函数8. double showarray(int n); /定义输出化简系数矩阵过程的函数9. int n,i,j,k,p,q;10. double main()11. 12. 13. printf(请输入系数矩阵的阶数 n(ndim)17. 18. printf(错误：元数超过初设定的值 %d，请重启程序重新输入 n,dim);19. exit(0);20. 21.22. /*输入系数矩阵和常数矩阵（即增广矩阵）的元素*/23. getarray(n);24. 25. /*使对角线上的主元素不为 0*/26. for(j=1;j=1;j-)76. 77. xj=bj;78. for(k=n;k=j+1;k-)79. xj=xj-xk*ajk;80. xj=xj/ajj;81. 82. printf(n 原方程组的唯一一组实数解为：n);83. for(j=1;j=n;j+)84. printf(x%d= %.3fn,j,xj);85. 86.87. /*定义矩阵输入函数 getarray(n)并打印以作检查*/88. double getarray(int n)89. 90. printf(n 请输入该矩阵各行的实数（以空格隔开）n);91. for(i=1;i=n;i+)92. 93. printf(n 第%d 行:t,i);94. for(j=1;j=n;j+)95. 96. scanf(%lf,&aij);97. printf(a%d%d= %.3f,i,j,aij);98. printf(n);99. 100. 101. printf(nA=);102.for(i=1;i=n;i+)103. 104.for(j=1;j=n;j+)105.printf(t%.3f,aij);106.printf(n);107. 108. printf(n);109. /*输入常数矩阵的各个数*/110. for(i=1;i=n;+i)111. 112. printf(请输入常数 b%d = ,i);113. scanf(%lf,&bi);114. 115.116.117./*定义增广矩阵 C 输出函数 showarray(n)*/118.double showarray(int n)119.120.printf(n 通过初等行变换最终增广矩阵矩阵 C 化为：n);121. printf(C=);122. for(i=1;i=n;i+)123. 124. for(j=1;j=n;j+)125. printf(t%.3f,aij);126. printf(t%.3f,bi);127. printf(n);128. 129.130. temp=1;131. for(i=1;i=n;i+)132. temp*=aii;133. printf(n 矩阵的行列式|A|=%fn,temp);134. /*判断原线性方程组是否有唯一解*/135. if(temp=0)136. 137. printf(n 该方程组无唯一解，请重新启动程序输入n);138. exit(0);139. 140.复制代码程序执行结果：四、误差分析由程序执行结果图可知，该 C 语言程序所求得的该 N 阶矩阵方程即 N 维线性方程组的解为即由于程序中所有变量除了增广矩阵的角标以外都定义为 double 型，而 double 型变量的精确度是 16 位，所以程序运行过程中变量的有效数字至多有 15 位，而为了程序执行时界面的清爽，将每个变量的有效数字只取了小数点后 3 位，就运行的具体程序来说，这小数点的后三维数字均为有效数字，所以本程序的误差至多为 0.001 即小数点后三位。而在该具体的 5 维线性方程组中，用克拉默法则计算出系数行列式 |A|=665，其精确解为所以各个解均在误差范围内。