一、描述问题：
题目：求解线性方程组Ax=b，写成函数。其中，A为n×n的N阶矩阵，x为需要求解的n元未知数组成的未知矩阵，b为n个常数组成的常数矩阵。即
运行程序时的具体实例为：
转化为矩阵形式（为检验程序的可靠性，特意选取初对角线元素为0的矩阵方程组）即为：
二、分析问题并找出解决问题的步骤：
由高等代数知识可知，解高阶线性方程组有逆矩阵求解法、增广矩阵求解法等，而在计算机C语言中，有高斯列主消元法、LU分解法、雅克比迭代法等解法。
为了与所学的高等代数知识相一致，选择使用“高斯简单迭代消元法”，与高等代数中的“增广矩阵求解法”相一致。以下简述高斯消元法的原理：
算法基本原理：
首先，为了能够求解N阶线性方程组（N由用户输入），所以需要定义一个大于N维的数组a[dim+1][dim+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.#include 2.#include 3.#include 4.#define dim 10 //定义最大的维数10，为防止计算值溢出 5.double a[dim+1][dim+1],b[dim+1],x[dim+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(n<10):"); 14.scanf("%d",&n); 15. /*判断矩阵阶数是否超过界定值*/ 16. if(n>dim) 17. { 18. printf("错误：元数超过初设定的值%d，请重启程序重新输入\\n",dim); 19. exit(0); 20. } 21. 22. /*输入系数矩阵和常数矩阵（即增广矩阵）的元素*/ 23. getarray(n); 24. 25. /*使对角线上的主元素不为0*/ 26. for(j=1;j<=n-1;j++) 27. { 28. if(a[j][j]==0) 29. for(i=j+1;i<=n;i++) 30. { 31. if(a[i][j]!=0) 32. { 33. /*交换增广矩阵的第i行与第j行的所有元素*/ 34. for(k=1;k<=n;k++) 35. { 36. a[i][k]+=a[j][k]; 37. a[j][k]=a[i][k]-a[j][k]; 38. a[i][k]-=a[j][k]; 39. } 40. b[i]+=b[j]; 41. b[j]=b[i]-b[j]; 42. b[i]-=b[j]; 43. } 44. continue; //找到第j列第一个不为0的元素即跳回第一层循环 45. } 46. } 47. /*开始用高斯简单迭代消元法进行求解计算*/ 48. for(j=1;j<=n-1;j++) 49. { 50. /*使系数矩阵转化为上三角矩阵，常数矩阵相应进行变换*/ 51. for(i=j+1;i<=n;i++) 52. { 53. temp=a[i][j]/a[j][j]; 54. b[i]=b[i]-temp*b[j]; 55. for(k=1;k<=n;k++) 56. a[i][k]=a[i][k]-temp*a[j][k]; 57. printf("\\n通过初等行变换增广矩阵矩阵C化为：\\n"); 58. /*输出进行初等行变换的过程*/ 59. printf("C="); 60. for(p=1;p<=n;p++) 61. { 62. for(q=1;q<=n;q++) 63. printf("\%.3f",a[p][q]); . printf("\%.3f\\n",b[p]); 65. } 66. printf("\\n"); 67. } 68. } 69. 70. /*输出最终的增广矩阵C*/ 71. showarray(n); 72. 73. /* 开始按顺序反代求解x[i](i=n,n-1,n-2,…,2,1)*/ 74. x[n]=b[n]/a[n][n]; 75. for(j=n-1;j>=1;j--) 76. { 77. x[j]=b[j]; 78. for(k=n;k>=j+1;k--) 79. x[j]=x[j]-x[k]*a[j][k]; 80. x[j]=x[j]/a[j][j]; 81. } 82. printf("\\n原方程组的唯一一组实数解为：\\n"); 83. for(j=1;j<=n;j++) 84. printf("x[%d]= %.3f\\n",j,x[j]); 85.} 86. 87./*定义矩阵输入函数getarray(n)并打印以作检查*/ 88.double getarray(int n) .{ 90.printf("\\n请输入该矩阵各行的实数（以空格隔开）\\n"); 91.for(i=1;i<=n;i++) 92. { 93.printf("\\n第%d行:\",i); 94.for(j=1;j<=n;j++) 95. { 96. scanf("%lf",&a[i][j]); 97. printf("a[%d][%d]= %.3f",i,j,a[i][j]); 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("\%.3f",a[i][j]); 106.printf("\\n"); 107. } 108. printf("\\n"); 109. /*输入常数矩阵的各个数*/ 110. for(i=1;i<=n;++i) 111. { 112. printf("请输入常数b[%d] = ",i); 113. scanf("%lf",&b[i]); 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("\%.3f",a[i][j]); 126. printf("\%.3f",b[i]); 127. printf("\\n"); 128. } 129. 130. temp=1; 131. for(i=1;i<=n;i++) 132. temp*=a[i][i]; 133. printf("\\n矩阵的行列式|A|=%f\\n",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，其精确解为 所以各个解均在误差范围内。下载本文
