用C语言求解N阶线性矩阵方程Ax=b的简单解法

用C语言求解N阶线性矩阵方程Ax=b的简单解法一、描述问题：题目：求解线性方程组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.#includestdio.h2.#includestdlib.h3.#includemath.h4.#definedim10//定义最大的维数10，为防止计算值溢出5.doublea[dim+1][dim+1],b[dim+1],x[dim+1];//定义双精度数组6.doubletemp;7.doublegetarray(intn);//定义输入矩阵元素的函数8.doubleshowarray(intn);//定义输出化简系数矩阵过程的函数9.intn,i,j,k,p,q;10.doublemain()11.{12.13.printf(请输入系数矩阵的阶数n(n10):);14.scanf(%d,&n);15./*判断矩阵阶数是否超过界定值*/16.if(ndim)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(\t%.3f,a[p][q]);64.printf(\t%.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.doublegetarray(intn)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,&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(\t%.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.doubleshowarray(intn)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,a[i][j]);126.printf(\t%.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，其精确解为所以各个解均在误差范围内。