用C语言求解N阶线性矩阵方程Ax=b的简单解法
线性方程组AX=B的数值计算方法实验
for(i=1;i<=n;++i)
{
for(j=1;j<=n;++j)
{
scanf("%lf",&a[i][j]);
}
}
printf("请输入方程组的常数项：\n");
for(i=1;i<=n;++i)
{
scanf("%lf",&b[i]);
}
for(i=1;i<=n;++i)
t[i]=c[i][N];
for(i=N-1;i>=0;i--)//利用回代法求最终解
{
for(j=N-1;j>i;j--)
t[i]-=c[i][j]*x[j];
x[i]=t[i]/c[i][i];
}
}
运行结果：（以具体方程组为例）
2、（PA=LU：带选主元的分解法）求解线性方程组AX=B,其中：
A= B=
#define N 4//矩阵阶数
voidColPivot(doublec[N][N+1],double[]);//函数声明
voidmain(){
inti,j;
doublex[N];
doublec[N][N+1]={1,3,5,7,1,
2,-1,3,5,2,
0,0,2,5,3,
-2,-6,-3,1,4};
（3）[max1,j]=max(abs(A(p:N,p)));C=A(p,:);A(p,:)=A(j+p-1,:);A(j+p-1,:)=C;d=t(p);t(p)=t(j+p-1);t(j+p-1)=d;
c语言程序 线性方程求解 有解释
#include <stdio.h> //标准输入输出头文件#include <stdlib.h> //不是C标准库中的头文件包含了的C语言标准库函数的定义定义了五种类型、一些宏和通用工具函数#include <malloc.h> //包含内存分配函数的函数声明定义常量和表现由堆例程使用的类型#include <math.h> //数学函数库//下面代码为定义为双精度浮点数整形int GS(int,double**,double *,double); //求解函数声明double **TwoArrayAlloc(int,int); //创造二维数组函数声明void TwoArrayFree(double **); //析构二维数组函数声明//注：析构函数当对象脱离其作用域时（对象所在的函数已调用完毕）系统自动执行析构函数。
往往用来做“清理善后”的工作void main() //主函数{int i,j,n; //定义性变量double ep,**a,*b; //定义浮点型变量二维指针一维指针ep = 1e-4; //定义精确度之后用来判断printf("你要解几元线性方程组：\n"); //输出解方程的初步要求scanf("%d",&n); //输出函数输入整形数字//注：求解n元线性方程（即有n个未知数），需要有n个方程组，所以系数总数为n*n,每个方程各有一个常数，一共有n个常数//下面代码分别为系数和常数分配总内存a = TwoArrayAlloc(n,n); //为方程组分配系数所占内存:n*n*sizeof(double)b = (double *)calloc(n,sizeof(double)); //为方程组各个方程分配常数内存:n*sizeof(double)if(b == NULL) //判断内存分配情况{printf("内存分配失败\n"); //输出判断结果exit(1); //内存分配失败，退出程序}//下面代码将用户输入的系数和常数分别一一对应的写入上一步所分配的内存中for(i=0;i<n;i++) //输入n元方程组的系数{printf("请输入第%d行相应的系数：\n",i+1); //说明第几行的方程系数for(j=0;j<n;j++) //循环输入对应行的系数{printf("a[%d][%d]: ",i,j); //输入数在数组中的相对位置scanf("%lf",a[i]+j); //输入系数fflush(stdin); //清空缓存}printf("请输入第%d行相应的常数：\n",i+1); //说明第几行方程组的常数printf("b[%d]: ",i); //循环输入对应行的系数scanf("%lf",b+i); //满足上一个条件后数在数组中的对应位置fflush(stdin); //清空缓存}//注:向下循环如果一直满足j<n则一直向下运行下面步骤一直运行j+1的情况直到不能满足j<n//下面代码根据输入的系数和常数，在屏幕上打印出方程组的具体表达式printf("方程组：\n"); //显示方程组提示语for(i=0;i<n;i++) //打印方程组{for(j=0;j<n;j++) //输入对应行的系数{if(a[i][j]>0) //判断方程组的系数为正{if(j>0)printf(" + "); //从第二行开始如果系数为正输出+号if(a[i][j]!=1) //如果系数不等于1printf("%lfX%d",a[i][j],j+1); //输出系数与Xn的乘积else //否则printf("X%d",j+1); //输出Xn}if(a[i][j]<0) //判断系数小于0{if(j>0)printf(" - "); //如果系数为负输出-号if(a[i][j]!=-1) //如果系数不等于-1printf("%lfX%d",fabs(a[i][j]),j+1); //输出系数的绝对值与Xn的乘积else //否则printf("X%d",j+1); //输出Xn}//注:向下循环如果一直满足i<n 则一直向下运行下面步骤一直运行i+1的情况直到不能满足i<n注：x【n】意思是x包含了n个数据的数组}printf("= %lf\n",b[i]); //打印等号以及常数项}//下面代码判断是否可以用高斯消去法求解，具体逻辑请参考高斯消去法求解线性方程if(!GS(n,a,b,ep)) //判断是否可以用高斯消去法求解，{printf("不可以用高斯消去法求解\n"); //如果不能输出提示语exit(0); //不能用高斯消去法求解，退出程序}printf("该方程组的解为：\n"); //输出该方程组的解打印出结果for(i=0;i<n;i++) //输出结果printf("x%d = %.2f\n",i+1,b[i]); //输出结果TwoArrayFree(a); //释放系数所占内存free (b); //释放常数所占内存}//下面代码为高斯消去法的判断依据函数int GS(int n,double **a,double *b,double ep) //求解函数主体为双精度{int i,j,k,l; //定义整型变量double t; //定义浮点型变量for(k=1;k<=n;k++) //判断每个系数的大小{for(l=k;l<=n;l++) //判断每个系数的大小if(fabs(a[l-1][k-1])>ep) //如果每列系数大于epbreak; //跳出else if(l==n) //如果l=nreturn(0); //返回0if(l!=k) //如果1不等于k{for(j=k;j<=n;j++) //交换系数{t = a[k-1][j-1]; //把第k-1行的每个系数赋给ta[k-1][j-1]=a[l-1][j-1]; //把第l-1行的每个系数赋给k-1行a[l-1][j-1]=t; //把t赋给第1-l行的每个系数}t=b[k-1]; //把每1行的常数项赋给tb[k-1]=b[l-1]; //把l-1行的常数项赋给k-1行b[l-1]=t; //把t的值赋给l-1行的常数}t=1/a[k-1][k-1]; //将t置为Xn分之一用来消元for(j=k+1;j<=n;j++) //消元a[k-1][j-1]=t*a[k-1][j-1]; //将每行的系数乘以该行确定的t 保证Xn的系数为1b[k-1]*=t; //把t的值赋给k-1行的常数for(i=k+1;i<=n;i++) //变换系数{for(j=k+1;j<=n;j++) //变换系数a[i-1][j-1]-=a[i-1][k-1]*a[k-1][j-1]; //把k-1行的每个系数分别赋给j-1和i+1行b[i-1]-=a[i-1][k-1]*b[k-1]; //把k-1行的常数项分别赋给i-1行和b }}for(i=n-1;i>=1;i--) //进行消元for(j=i+1;j<=n;j++) //进行消元b[i-1]-=a[i-1][j-1]*b[j-1]; //把j-1行的常数项分别赋给i-1行和breturn(1); //返回1的步骤}//注:向下循环如果一直满足i>=1则一直向下运行下面步骤一直运行i-1的情况直到不能满足i>=1//下面代码为分配内存的函数double **TwoArrayAlloc(int r,int c) //二维函数创建函数主体{double *x,**y; //定义浮点型变量二维指针一维指针int n; //定义变量x=(double *)calloc(r*c,sizeof(double)); //用x指向r*c个大小为double的内存y=(double **)calloc(r,sizeof(double*)); //用y指向r个double*内存，if(!x||!y) //如果不满足x和y{printf("内存分配失败\n"); //输出判断结果exit(1); //内存分配失败，退出程序}for(n=0;n<=r-1;++n) //输入对应取值y[n]=&x[c*n]; //定义为x在c*n处的地址赋给yreturn (y); //返回y}//注:向下循环如果一直满足n<=r 则一直向下运行下面步骤一直运行n+1的情况直到不能满足n<=r//下面代码为释放内存的函数void TwoArrayFree(double **x) //析构函数{free(x[0]); //释放内存free(x); //释放内存。
线性方程组的数值算法C语言实现(附代码)
线性方程组AX=B 的数值计算方法实验一、 实验描述：随着科学技术的发展，线性代数作为高等数学的一个重要组成部分，在科学实践中得到广泛的应用。
本实验的通过C 语言的算法设计以及编程，来实现高斯消元法、三角分解法和解线性方程组的迭代法（雅可比迭代法和高斯-赛德尔迭代法），对指定方程组进行求解。
二、 实验原理：1、高斯消去法：运用高斯消去法解方程组，通常会用到初等变换，以此来得到与原系数矩阵等价的系数矩阵，达到消元的目的。
初等变换有三种：(a)、(交换变换)对调方程组两行；(b)、用非零常数乘以方程组的某一行；(c)、将方程组的某一行乘以一个非零常数，再加到另一行。
通常利用(c)，即用一个方程乘以一个常数，再减去另一个方程来置换另一个方程。
在方程组的增广矩阵中用类似的变换，可以化简系数矩阵，求出其中一个解，然后利用回代法，就可以解出所有的解。
2、选主元：若在解方程组过程中，系数矩阵上的对角元素为零的话，会导致解出的结果不正确。
所以在解方程组过程中要避免此种情况的出现，这就需要选择行的判定条件。
经过行变换，使矩阵对角元素均不为零。
这个过程称为选主元。
选主元分平凡选主元和偏序选主元两种。
平凡选主元：如果()0p pp a ≠，不交换行；如果()0p pp a =，寻找第p 行下满足()0p pp a ≠的第一行，设行数为k ，然后交换第k 行和第p 行。
这样新主元就是非零主元。
偏序选主元：为了减小误差的传播，偏序选主元策略首先检查位于主对角线或主对角线下方第p 列的所有元素，确定行k ，它的元素绝对值最大。
然后如果k p >，则交换第k 行和第p 行。
通常用偏序选主元，可以减小计算误差。
3、三角分解法：由于求解上三角或下三角线性方程组很容易所以在解线性方程组时，可将系数矩阵分解为下三角矩阵和上三角矩阵。
其中下三角矩阵的主对角线为1，上三角矩阵的对角线元素非零。
有如下定理：如果非奇异矩阵A 可表示为下三角矩阵L 和上三角矩阵U 的乘积： A LU = (1) 则A 存在一个三角分解。
三种典型矩阵方程的简单解法
即对矩阵 … 施行初等列变换 ,当把 A 变成 E 时 ,B
B 就变成 X 。(f ) 式提供了一个具体解矩阵方程 XA = B 的
简单方法 。 例2 解下列矩阵方程 。
2005 年 6 月第 3 期 三种典型矩阵方程的简单解法 3 (i) X - 1 3 1 2 4 0 1 2 3 - 1 3 3 A - 1 = 3 1 2 4 2 1 0 0 1 2 = (0 2 3) ; 1 2 4 0 1 2 3 0 0 5 0 1 2 1 2 0 1 2 3 - 1 3 ,B = ( 0 2 3) 。 1 4 7 0 1 2 0 0 - 8 1 0 - 1 0 1 0 0 0 1
… ,
X
可得解矩阵方程 AXB = C 的简单解法 例3 解下列矩阵方程 。
1 (i) 2 3 2 2 4 3 1 X 3 2 5 1 3 1 1 = 2 3 0 ;
…
- 4
…
- 3 1 (ii) X 4 7
…
3 2 5 8 3
于是有 X = ( - 4 - 3 3) 。
1 4 2 5 8 2 5 8 3 1 3 6 1 2 1 5 3 2 。 3 1 4 1 4 7 2 1 5 3 2 。 3 0 - 3 - 6 0 - 6 - 20 6 = 0 1 1 7 1 4 A 7
- 1 - 1 X = PL PL - 1 … P1- 1B 。 证毕
,再左乘 B 即得 X。 ,再右乘 B 即得 X。
- 1 - 1
若 XA = B ,则有 XAA - 1 = BA - 1 ,即 X = BA - 1 。于是
- 1
(1) (2)
又若 AXB = C ,则有 A AXBB
CB
A B
线性代数方程组求解
线性代数方程组求解一、实验要求编程求解方程组：方程组1：方程组2：方程组3：要求：用C/C++语言实现如下函数：1.bool lu(double* a, int* pivot, int n);实现矩阵的LU分解。
pivot为输出参数，pivot[0,n) 中存放主元的位置排列。
函数成功时返回false，否则返回true。
2.bool guass(double const* lu, int const* p, double* b, int n);求线代数方程组的解设矩阵Lunxn 为某个矩阵anxn 的LU 分解，在内存中按行优先次序存放。
p[0,n)为LU 分解的主元排列。
b 为方程组Ax=b 的右端向量。
此函数计算方程组Ax=b 的解，并将结果存放在数组b[0,n)中。
函数成功时返回false ，否则返回true 。
3. void qr(double* a, double* d, int n);矩阵的QR 分解假设数组anxn 在内存中按行优先次序存放。
此函数使用HouseHolder 变换将其就地进行QR 分解。
d 为输出参数，d [0,n) 中存放QR 分解的上三角对角线元素。
4. bool hshld(double const*qr, double const*d, double*b, int n); 求线代数方程组的解设矩阵qrnxn 为某个矩阵anxn 的QR 分解，在内存中按行优先次序存放。
d [0,n) 为QR 分解的上三角对角线元素。
b 为方程组Ax=b 的右端向量。
函数计算方程组Ax=b 的解，并将结果存放在数组b[0,n)中。
函数成功时返回false ，否则返回true 。
二、问题分析求解线性方程组Ax=b ，其实质就是把它的系数矩阵A 通过各种变换成一个下三角或上三角矩阵，从而简化方程组的求解。
因此，在求解线性方程组的过程中，把系数矩阵A 变换成上三角或下三角矩阵显得尤为重要，然而矩阵A 的变换通常有两种分解方法：LU 分解法和QR 分解法。
线性方程组求解 高质量C语言程序
课题：线性方程组求解目录课题描述及要求 (2)项目分析 (2)算法流程 (2)方法说明 (3)源代码 (3)程序说明 (6)运行结果 (7)总结 (7)参考文献 (8)线性方程组求解05111114 陈龙一．课题描述和功能要求1.描述：求解线性方程组Ax=b，写成函数。
其中，A为n乘n阶矩阵，x为n元未知向量，b为n个常数组成的矩阵。
2.要求：采用高斯先列主元消元法（也可采用其他方法）求解线性方程组AX=b。
二．项目分析数学上，高斯消去法或称高斯－约当消去法，由高斯和约当得名（很多人将高斯消去作为完整的高斯－约当消去的前半部分），它是线性代数中的一个算法，用于决定线性方程组的解，决定矩阵的秩，以及决定可逆方矩阵的逆。
当用于一个矩阵时，高斯消去产生“行消去梯形形式”。
例如：一个二元一次方程组，设法对每个等式进行变形，使两个等式中的同一个未知数的系数相等，这两个等式相减，得到一个新的等式，在这个新的等式中，细数相等的未知数就被除去了（系数为0）。
同样的也适合多元多次方程组。
我们知道m*n矩阵（用大写字母表示）是一个m行n列的数阵，n维向量（用加粗的小写字母表示）是n个数的数组，也就是一个n*1矩阵（列向量。
我们不考虑行向量）。
另外，大家也都知道矩阵乘法。
因此一个m*n线性方程组可以表示为Ax=b，其中A是由系数aij组成的m*n矩阵即系数矩阵，x是n维的未知数向量，b是m维的结果向量。
如果把向量b写到A的右边得到m*(n+1)的矩阵，得到的新矩阵称为这个方程组的增广矩阵。
每一个方程组均对应于一个增广矩阵。
三．算法流程图四．方法说明（1）第1步消元——在增广矩阵（A，b）第一列中找到绝对值最大的元素，将其所在行与第一行交换（2）第2步消元——在增广矩阵（A，b）中的第二列中（从第二行开始）找到绝对值最大的元素，将其所在行与第二行交换（3）第3步消元——在增广矩阵（A，b）中的第三列中（从第三行开始）找到绝对值最大的元素，将其所在行与第二行交换（4）按x4- x3- x2- x1的顺序回代求解出方程组的解。
C语言解线性方程的四种方法
C语言解线性方程的四种方法发了好几天编了个解线性方程组的小程序，可第一次实战就大败而归。
经过半天的调试，仍找不出纠正的方法。
因为并不是算法的问题，而是因为自己对编译器处理浮点函数的方法不是很理解。
明明D=0的方阵解出来不等于0了，跟踪调试发现，计算过程程序对数据进行了舍去处理，导致最终结果不对。
不过如果没有浮点型的话，这个程序应该算不错了。
复制代码代码如下:#include<stdio.h>#include<math.h>#include<mem.h>#define NUM 100void print(void) /* 使用说明*/{ clrscr();printf("\n\n\n\n\n\t\t\t\t Introduction \n");printf("\t*--------------------------------------------------------------*\n");printf("\t* This program was design for compute linear equations. *\n");printf("\t* The way of use it is very simple. *\n");printf("\t* First : Input the number of the equation;(Input 0 to exit) *\n");printf("\t* Second: Input the coefficient of every eqution; *\n");printf("\t* Third : Input the constant of every eqution; *\n");printf("\t* Last : Chose the way you want use to solve the equtions; *\n");printf("\t* That's all, input any key to run it . . . *\n");printf("\t*-------------------------By__TJX------------------------------*\n");getch(); }void chose(void) /*选择计算方法*/{ clrscr();fflush(stdin);printf("\n\n\n\n\n\t\t**********Introduction********** \n");printf("\t\t* Chose the way,please. * \n");printf("\t\t* a : Gauss eliminant. * \n");printf("\t\t* b : Gauss_yd eliminant. * \n");printf("\t\t* c : Iterative way. * \n");printf("\t\t* d : Cramer way. * \n");printf("\t\t* e : exit. * \n");printf("\t\t*************By__TJX************ \n");printf("\t\tPlease choose number :\n");}void input(double **a1,double b1[],int num) /*数据输入*/{ int i,j,t;double *p;char de1,de2;do{printf("Please input array a[%d][%d]: \n",num,num);printf("Warn: The first number of the array mustn't contain zero! \n");for(i=1;i<=num;i++){printf("Please input array a[%d][]: \n",i);for(j=1;j<=num;j++){t=0;if(i==1&&j==1){ do{if(t==0) { scanf("%lf",&a1[i][j]); t++;}else {printf("The input is invalid,please input again:\n"); scanf("%f",&a1[i][j]);}}while(a1[i][j]==0);}else scanf("%lf",&a1[i][j]);}}printf(" \nPlease check the value of array a[%d][%d],press Y to input again.\n",num,num);do{de1=getch();}while(de1!='y'&&de1!='Y'&&de1!='n'&&de1!='N');}while(de1=='y'||de1=='Y');do{printf("Please input array b[%d]: \n",num);p=b1+1;for(i=1;i<=num;i++)scanf("%lf",p++);printf(" \nPlease check the value of array b[%d],press Y to input \n",num);do{de2=getch();}while(de2!='y'&&de2!='Y'&&de2!='n'&&de2!='N');}while(de2=='y'||de2=='Y');}int max(double *t1, double x1[],int n) /*迭代子函数*/{ int i,temp=0;for(i=1;i<=n;i++)if(fabs(x1[i]-t1[i])>1e-2) {temp=1;break;}/* printf(" %d ",temp); */return temp;}int ddcompute(double **a1,double b1[],double x1[],int n) /*迭代法计算*/{double *t;int i,j,k=0;double sum1=0.0,sum2=0.0;t=(double*)malloc(n*sizeof(double));printf("\nPlease Input The Initial Value of x:\n");for(i=1;i<=n;i++)scanf("%lf",&x1[i]);do{ k++;for(i=1;i<=n;i++)t[i]=x1[i];for(i=1;i<=n;i++){ sum1=0.0;sum2=0.0;for(j=1;j<=i-1;j++) sum1=sum1+a1[i][j]*x1[j]; /*printf(" sum1= %0.4f ",sum1);*/for(j=i+1;j<=n;j++) sum2=sum2+a1[i][j]*t[j]; /* printf(" sum2= %0.4f ",sum2);}*/if(a1[i][i]==0||fabs(sum1)>1e+12||fabs(sum2)>1e+12){printf(" \nWarning: These equtions can't be solve by this way!\n Press any Key to continue...");getch();free(t);return 0;}x1[i]=(b1[i]-sum1-sum2)/a1[i][i];}}while(max(t,x1,n));/* for(i=1;i<=n;i++){if(i%3==0) printf("\n");printf(" %.4f ",x1[i]);}*/free(t);return 1; }int gscompute(double **a1,double b1[],double x1[],int n) /*高斯消元法计算*/ {int i,j,k;double m,sum;for(k=1;k<=n-1;k++)for(i=k+1;i<=n;i++){ if(a1[k][k]==0) {printf(" \nThese equtions can't be solve is this \n Press any Key to continue...");getch();return 0; }if((m=0-a1[i][k]/a1[k][k])==0) {i++; continue;}else {for(j=k+1;j<=n;j++)a1[i][j]=a1[i][j]+a1[k][j]*m;b1[i]=b1[i]+b1[k]*m;}}/* yi xia ji suan x zhi */x1[n]=b1[n]/a1[n][n];for(i=n-1;i>=1;i--){sum=0.0;for(j=n;j>=i+1;j--)sum=sum+a1[i][j]*x1[j];x1[i]=(b1[i]-sum)/a1[i][i];}return 1; }int gs_ydcompute(double **a1,double b1[],double x1[],int n) /*高斯_约当法计算*/ {int i,j,k;double m,sum;for(k=1;k<=n;k++){i=1;while(i<=n){ if(a1[k][k]==0) {printf(" \nThese equtions can't be solve is this way.\n Press any Key to continue...");getch(); return 0;}if(i!=k){ if((m=0-a1[i][k]/a1[k][k])==0) {i++; continue;}else {for(j=k+1;j<=n;j++)a1[i][j]=a1[i][j]+a1[k][j]*m;b1[i]=b1[i]+b1[k]*m;}i++;}else i++; }}/* yi xia ji suan x zhi */for(i=n;i>=1;i--)x1[i]=b1[i]/a1[i][i];return 1;}double computed(double **a,int h,int l, int *c1,int n) /*计算系数行列式D值*/{ int i, j,p=1;double sum=0.0;if(h==n)sum=1.0;else {i=++h;c1[l]=0;for(j=1;j<=n;j++)if(c1[j])if(a[i][j]==0) p++;else {sum=sum+a[i][j]*computed(a,i,j,c1,n)*pow(-1,1+p); p++; } c1[l]=1; }return sum; }void ncompute(double **a,double b[],double x[],int n,int *c,double h) /*克莱姆法计算*/ {int i,j;double t[NUM];for(j=1;j<=n;j++){ for(i=1;i<=n;i++){t[i]=a[i][j];a[i][j]=b[i];}x[j]=computed(a,0,0,c,n)/h;for(i=1;i<=n;i++)a[i][j]=t[i]; }}main(){double x[NUM];double b[NUM];int i,j=2,n=0;int *c;double he;char m,decision;double **a;a=(double**)malloc(NUM*sizeof(double*));for (i=0; i<NUM; i++)a[i]=(double*)malloc(NUM*sizeof(double));print();do{clrscr();do{if(n>=NUM) printf("n is too large,please input again:\n");elseprintf("Please input the total number of the equations n(n<NUM): \n");scanf(" %d",&n);}while(n>NUM);if(n==0) {for(i=1; i<NUM; i++) free(a[i]);free(a); exit(1);}input(a,b,n);c=(int *)malloc((n+1)*sizeof(int));memset(c,1,(n+1)*sizeof(int));he=computed(a,0,0,c,n);if(fabs(he)>1e-4){Other: chose();do{m=getche();}while(m!='a'&&m!='b'&&m!='A'&&m!='B'&&m!='c'&&m!='C'&&m!='d'&&m!='D' &&m!='e'&&m!='E');switch(m){ case 'a': ;case 'A': j=gscompute(a,b,x,n); break;case 'b': ;case 'B': j=gs_ydcompute(a,b,x,n); break;case 'c': ;case 'C': j=ddcompute(a,b,x,n); break;case 'd': ;case 'D': j=1; ncompute(a,b,x,n,c,he); break;case 'e': ;case 'E': j=2; break;default: j=2; break;}if(j==1){ clrscr();printf("\n\n\n\n");printf(" D=%.4f \n",he);for(i=1;i<=n;i++){if(i%5==0) printf("\n");printf(" %.4f ",x[i]);}}else if(j==0){printf(" \nThese equtions can't be solve is this way.\nPlease chose the other way."); goto Other;}else {for(i=1; i<NUM; i++) free(a[i]);free(a);free(c);exit(1);}}else printf(" \n\n\tD=%.4f\n This linear equations hasn't accurate answer!",he);printf(" \n Do you want to continue?(Y/N) \n");do{decision=getchar();}while(decision!='y'&&decision!='Y'&&decision!='n'&&decision!='N');}while(decision=='y'||decision=='Y');for(i=1; i<NUM; i++) free(a[i]);free(a);free(c);}如对本文有所疑问，请点击进入脚本之家知识社区提问。
矩阵方程AX=B的求解方法
矩阵方程AX=B的求解方法作者：姜鹏李扬来源：《课程教育研究》 2020年第17期姜鹏李扬（沈阳化工大学数理系辽宁沈阳 110142）【摘要】本文叙述了当A为可逆矩阵以及A为不可逆矩阵或者不是方阵时，矩阵方程AX=B的求解方法。
【关键词】矩阵方程线性代数【中图分类号】O151 【文献标识码】A 【文章编号】2095-3089（2020）17-0126-01矩阵方程是以矩阵为未知量的方程。
在矩阵方程AX=B中，A、B为已知矩阵，X为未知矩阵。
矩阵方程AX=B的求解问题，是线性代数中的一种典型问题，常用的求解方法主要分为如下的两种类型。
一、A为可逆矩阵当A为可逆矩阵时，用A的逆矩阵A-1分别左乘矩阵方程AX=B的左右两端，可得其唯一解为X=A-1B。
这种类型的矩阵方程，可细分为下列的两种解法。
二、A为不可逆矩阵或者不是方阵实际上，在计算矩阵方程AX=B时，并不知道矩阵A是否是可逆矩阵。
在具体操作时，当A为方阵时，可以按照上述的做法，先求出|A|或者对（A:B）施以初等行变换。
如果|A|=0或者A 化成的行最简形矩阵不是单位矩阵E，这时就说明A为不可逆矩阵。
当A为不可逆矩阵或者不是方阵时，就需要将矩阵X中的所有元素都设为未知数，并将原来的矩阵方程转化为关于上述未知数的线性方程组。
这时，矩阵方程AX=B就不一定有解。
参考文献：[1]谢彦红，吴茂全.线性代数及其MATLAB应用[M]. 第二版. 北京:化学工业出版社,2017年.[2]毛纲源.线性代数解题方法技巧归纳[M]. 第二版. 武汉:华中科技大学出版社, 2000年.作者简介:姜鹏（1976-），男，汉族，辽宁沈阳人，硕士，讲师，主要从事应用数学的教学和研究工作。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
用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列元素a jj=0 ，那此时不能再以a jj 为基元。
当变换到第j列时，从j行j列的元素a jj 以下的各元素中选取第一个不为0的元素，通过第三类初等行变换即交换两行将其交换到a jj 的位置上，然后再进行消元过程。
交换系数矩阵中的两行，相当于两个方程的位置交换了。
再由高斯消元法，将第j列元素除a jj外第j行以下的其他元素通过第二种初等行变换化为0，这样，就能使系数矩阵通过这样的行变换化为一个上三角矩阵，即，当系数矩阵A进行初等行变换时，常数矩阵也要进行对应的初等行变换，即此时那么有接下来，进行“反代”，由可求出，再往上代入即可求出以此类推，即可从x n推到x n-1，再推到x n-2直至x1。
至此，未知矩阵x的所有元素就全部求出，即求出了原方程组有且仅有的唯一一组解。
基本原理示意图：三、编写程序1.#include<stdio.h>2.#include<stdlib.h>3.#include<math.h>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("\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.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",&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.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",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，其精确解为所以各个解均在误差范围内。