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

合集下载

线性方程组AX=B的数值计算方法实验

线性方程组AX=B的数值计算方法实验
printf("请输入矩阵各项:\n");
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语言程序 线性方程求解 有解释

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语言实现(附代码)

线性方程组的数值算法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语言程序

线性方程组求解 高质量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语言解线性方程的四种方法

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的求解方法

矩阵方程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-),男,汉族,辽宁沈阳人,硕士,讲师,主要从事应用数学的教学和研究工作。

线性方程组AX=B的数值计算方法实验

线性方程组AX=B的数值计算方法实验

线性方程组AX=B的数值计算方法实验学号:姓名:梁哲豪一、实验描述在自然科学和工程技术中很多问题的解决常常归结为解线性代数方程组。

例如电学中的网络问题,船体数学放样中建立三次样条函数问题,用最小二乘法求实验数据的曲线拟合问题,解非线性方程组问题,用差分法或者有限元法解常微分方程,偏微分方程边值问题等都导致求解线性方程组,而且后面几种情况常常归结为求解大型线性方程组。

线性代数方面的计算方法就是研究求解线性方程组的一些数值解法与研究计算矩阵的特征值及特征向量的数值方法。

关于线性方程组的数值解法一般有两类:直接法:若在计算过程中没有舍入误差,经过有限步算术运算,可求得方程组的精确解的方法。

迭代法:用某种极限过程去逐步逼近线性方程组精确解的方法。

迭代法具有占存储单元少,程序设计简单,原始系数矩阵在迭代过程中不变等优点,但存在收敛性及收敛速度等问题。

上三角线性方程组的求解:基本算法:高斯消元法:将原方程组化为三角形方阵的方程组:(k=1,2,…,n-1; i=k+1,k+2, …,n ;j=k+1,k+2, …,n+1)由回代过程求得原方程组的解:LU分解法:将系数矩阵A转化为A=L*U,L为单位下三角矩阵,U为普通上三角矩阵,然后通过解方程组l*y=b,u*x=y,来求解x。

二、实验内容1、许多科学应用包含的矩阵带有很多零。

在实际情况中很重要的三角形线性方程组有如下形式:……构造一个程序求解三角形线性方程组。

可假定不需要变换。

而且可用第k 行消去第k+1行的x。

k核心代码:#include<iostream.h>#include<math.h>#include<iomanip.h>#define N 4//矩阵阶数void ColPivot(double c[N][N+1],double[]);//函数声明void main(){int i,j;double x[N];double c[N][N+1]={1,3,5,7,1,2,-1,3,5,2,0,0,2,5,3,-2,-6,-3,1,4};cout<<"----------------------------------------"<<endl;cout<<"系数矩阵为: \n";for(i=0;i<N;i++){for(j=0;j<N;j++)cout<<setw(10)<<c[i][j];cout<<endl;}cout<<"右侧矩阵 y 为: \n";for(i=0;i<N;i++)cout<<setw(10)<<c[i][N];cout<<endl;cout<<"----------------------------------------"<<endl;ColPivot(c,x);//调用函数,进行高斯消去法变换cout<<"变换后得到的三角矩阵: \n";for(i=0;i<N;i++){for(j=0;j<N;j++)cout<<setw(10)<<c[i][j];cout<<endl;}cout<<"变换后的右侧矩阵 y 为: \n";for(i=0;i<N;i++)cout<<setw(10)<<c[i][N];cout<<endl;cout<<"----------------------------------------"<<endl; cout<<"方程的解为: \n";for(i=0;i<N;i++)cout<<" x["<<i<<"]= "<<x[i]<<endl;cout<<"----------------------------------------"<<endl; }void ColPivot(double c[N][N+1],double x[]){int i,j,k;double p,max;double t[N];for(i=0;i<=N-2;i++){max=0;k=i;for(j=i+1;j<N;j++)if(fabs(c[j][i])>max){k=j;max=fabs(c[j][i]);//选主元}if(k!=i)for(j=i;j<=N;j++){p=c[i][j];c[i][j]=c[k][j];//选出主元后进行交换c[k][j]=p;}for(j=i+1;j<N;j++){p=c[j][i]/c[i][i];for(k=i;k<=N;k++)c[j][k]-=p*c[i][k];//高斯消去,进行计算}}for(i=0;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=核心代码:#include <stdio.h>#include <math.h>#define L 30double a[L][L],b[L],l[L][L],u[L][L],x[L],y[L];int main(){int n,i,j,k,r;printf("请输入矩阵元次:\n");scanf("%d",&n);printf("请输入矩阵各项:\n");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){for(j=1;j<=n;++j){l[i][j]=0;u[i][j]=0.0;}}for(k=1;k<=n;++k){for(j=k;j<=n;++j){u[k][j]=a[k][j];for(r=1;r<k;++r){u[k][j]-=l[k][r]*u[r][j];}}for(i=k+1;i<=n;++i){l[i][k]=a[i][k];for(r=1;r<k;++r){l[i][k]-=l[i][r]*u[r][k];}l[i][k]/= u[k][k];}l[k][k]=1.0;}for(i=1;i<=n;++i){y[i]= b[i];for(j=1;j<i;++j){y[i]-=l[i][j]*y[j];}}for(i=n;i>0;--i){x[i]= y[i];for(j=i+1;j<=n;++j){x[i]-=u[i][j]*x[j];}x[i]/= u[i][i];}for(i=1;i<=n;++i){printf("%0.2lf\n",x[i]);}return0;}运行结果:3、使用程序3.3求解线性方程组AX=B,其中,A= [a ij] N×N= i j-1,而且B=[b ij] N×1, b11=N,当i≥2时,b i1=(i N-1)/(i-1),对N=3,7,11的情况分别求解。

矩阵方程AXB的解法

矩阵方程AXB的解法

2x2 -3x3 =3,x4 +2x5 -3x6 = -1,x1 +x2 -x3 =2,x4 +x5 -x6 = 0,解 得 x2 = 3-2x1,x3 = 1-x1x5 = 1- 2x4,x6 = 1-x4,于 是 所 求 矩 阵 的 解 为
X = éëêê31--x21xx11
x4 1-2x4 1-x4
方面都有非常广泛的应用.矩阵方程作为代数方程的一类,
在上述应用过程中,涉及核心问题必然是 矩 阵 方 程 解 的 判 定
和求解.文献已对矩 阵 方 AX =B 解 的 判 定 做 了 探 讨,本 文 来探讨矩阵方程 AX=B 的 求 解 方 法. 对 于 AX =B 的 求 解 方法,在我们熟知的线性代数教材中都局限于系数矩阵A 可 逆的情形,系数矩阵A 不可逆的情形基本上没有涉及.本文
例2 设 A =
éëêê341
-1 -3 3
230ùûúú ,B =
éëêê317
9 11 5
777ùûúú ,且 AX
= B,求 X .
解 :对 (A ,B)施 行 初 等 行 变 换 ,得 到 (A ,B)=
éëêê3 4 1
-1 -3 3
2 3 0
3 1 7
9 11 5
777ùûúú ~
éëêê100
给出了系数矩阵 A 不 可 逆 时 的 AX =B 的 两 种 普 适 性 求 解 方法,并通过例子进行了说明.文中出 现 的 数 学 符 号 都 是 标
准 的 ,可 参 见 文 献 .
二 、求 解 方 法 1——— 待 定 法 所谓的待定法,首先设出 未 知 矩 阵 X 中 的 待 定 元 素,然
后转化为线性方 程 组,用 解 方 程 组 的 方 法 求 出 矩 阵 X 的 待 定元素.

矩阵方程ax=b的解法

矩阵方程ax=b的解法

矩阵方程ax=b的解法
矩阵方程是数学学科中最常用的问题之一,它需要用矩阵构建相关的方程,来求解所求问题。

矩阵方程ax=b,就是矩阵a乘以一个矩阵x(其中x的元素均可以是实数)等于矩阵b的方程,其中的矩阵a、x和b都是矩阵,解矩阵方程就是确定这个方程中x的值。

解矩阵方程的一般方法有两种:一种是用逆矩阵的方法求解,另一种则是将矩阵方程拆解为一组线性方程组,再利用一些方法求解线性方程组求解。

首先来看用逆矩阵法解矩阵方程,逆矩阵是可以得到的,它乘以原矩阵能得到单位矩阵。

将a -1 乘以a、乘以b可以得到矩阵x,从而解出矩阵方程ax=b。

另一种方法就是将矩阵方程ax=b化为一组线性方程组,那么解线性方程组便可以得到每个元素。

解线性方程组的常用方法有高斯消元、高斯-约旦消元法、矩阵变换求解形式等。

在求解数学问题中,矩阵方程是比较常用的一种方法,能够有效地求解n元线性方程组问题,并且在科学计算中应用广泛。

矩阵方程的解法有许多,以上介绍的是两种较为常见的解法,希望对大家能够有所帮助。

矩阵ax=b求x

矩阵ax=b求x

矩阵ax=b求x
可以用这两种方法解答:
1、初等变换法:有固定方法,设方程的系数矩阵为A,未知数矩阵为X,常数矩阵为B,即AX=B,要求X,则等式两端同时左乘A^(-1),有X=A^(-1)B。

又因为(A,E)~(E,A^(-1)),所以可用初等行变换求A^(-1),从而所有未知数都求出来了。

2、逆矩阵求解法:求解方法:容易算出已知矩阵的行列式等于-1。

然后计算伴随阵,具体方法是对于编号为mn 的元素,划去原阵的第m行和第n列,原阵退化为n-1阶矩阵,求出这个n-1阶阵的行列式,然后填入伴随阵的第n行第m列位置,最后乘以-1的m+n次幂。

矩阵第m行与第n列交叉位置的那个值,等于第一个矩阵第m行与第二个矩阵第n列,对应位置的每个值的乘积之和。

对于矩阵方程,当系数矩阵是方阵时,先判断是否可逆。

如果可逆,则可以利用左乘或右乘逆矩阵的方法求未知矩阵,如果方阵不可逆或是系数矩阵不是方阵,则需要用矩阵的广义逆来确定矩阵方程有解的条件,进而在有解的情形求出通解。

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

用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列元素a jj=0 ,那此时不能再以a jj 为基元。

当变换到第j列时,从j行j列的元素a jj 以下的各元素中选取第一个不为0的元素,通过第三类初等行变换即交换两行将其交换到a jj 的位置上,然后再进行消元过程。

ax=b的方程解法

ax=b的方程解法

ax=b的方程解法
解方程ax=b的方法有很多种,下面列举几种常见的解法。

1. 代入法:
将b代入方程ax中,得到ax=b,使用代入法可以将方程转化
为一个变量的一次方程,然后求解这个方程即可。

2. 直接除法:
如果a不为0,那么可以直接两边除以a,得到x=b/a,即为方
程的解。

3. 移项法:
将方程ax=b通过移项变形为ax-b=0,然后使用因式分解、配
方法等方法将方程转化为一次方程的形式,然后求解这个方程即可。

4. 矩阵方法:
将方程ax=b表示成矩阵的形式,即[A][x]=[B],其中[A]是一
个n行n列的矩阵,[x]是一个n行1列的列向量,[B]是一个
n行1列的列向量。

使用矩阵求逆的方法可以解得x=[A]^-1[B]。

以上是解方程ax=b的几种常见方法,具体选择哪一种方法取
决于方程的形式和具体的计算要求。

求解N阶线性矩阵方程

求解N阶线性矩阵方程

求解N阶线性矩阵方程一:编程语言:C语言二:问题描述求解线性方程组Ax=b。

手动输入A,其中A为n*n的N阶矩阵,x为列向量,即需要求解的未知的变量,b为列向量,是非齐次线性的常数。

三:问题分析及解题思路首先,用户输入矩阵A和列向量b,这些用C语言很容易实现,然后判断方程是否有解,及解的个数。

运用数学知识知道矩阵A的行列式不等于0才有唯一解。

然后进行初等变换,使系数转换为上三角矩阵,这些代码第一次作业做过了,很容易实现。

再然后根据数学知识代入即可求得原方程唯一解。

四:代码实现#include"stdio.h"#include"stdlib.h"#include"math.h"#define dim 10double a[dim+1][dim+1],b[dim+1],x[dim+1];double temp;double getarray(int n);double showarray( int n);int n,j,i,k,p,q;void main(){printf("请输入系数矩阵的阶数n(n<10):");scanf("%d",&n);if(n>dim){printf("错误,请输入系数矩阵的阶数n(n<10):");scanf("%d",&n);}getarray(n);for(j=1;j<=n-1;j++){if(a[j][j]==0)for(i=j+1;i<=n;i++){if(a[i][j]!=0){for(k=1;k<=n;k++){a[i][k]+=a[j][k];a[j][k]=a[i][k]-a[j][k];a[i][k]-=a[j][k];}b[i]=b[i]+b[j];b[j]=b[i]-b[j];b[i]-=b[j];}continue;}}for(j=1;j<=n-1;j++){for(i=j+1;i<=n;i++){temp=a[i][j]/a[j][j];b[i]=b[i]-temp*b[j];for(k=1;k<=n;k++)a[i][k]=a[i][k]-temp*a[j][k];printf("\n通过初等行变换增广矩阵矩阵C化为:\n");printf("C=");for(p=1;p<=n;p++){for(q=1;q<=n;q++)printf("\t%.3f",a[p][q]);printf("\t%.3f\n",b[p]);}printf("\n");}}showarray(n);x[n]=b[n]/a[n][n];for(j=n-1;j>=1;j--){x[j]=b[j];for(k=n;k>=j+1;k--)x[j]=x[j]-x[k]*a[j][k];x[j]=x[j]/a[j][j];}printf("\n原方程组的唯一一组实数解为:\n");for(j=1;j<=n;j++)printf("x[%d]=%.3f\n",j,x[j]);}double getarray(int n){printf("\n请输入该矩阵各行的实数(空格隔开)\n");for(i=1;i<=n;i++){printf("\n第%d行:\t",i);for(j=1;j<=n;j++){scanf("%lf",&a[i][j]);printf("a[%d][%d]=%.3f",i,j,a[i][j]);printf("\n");}}printf("\nA=");for(i=1;i<=n;i++){for(j=1;j<=n;j++)printf("\t%.3f",a[i][j]);printf("\n");}printf("\n");for(i=1;i<=n;++i){printf("请输入常数b[%d]=",i);scanf("%lf",&b[i]);}return 0;}double showarray(int n){printf("\n通过初等变换最终增广矩阵C化为:\n");printf("C=");for(i=1;i<=n;i++){for(j=1;j<=n;j++)printf("\t%.3f",a[i][j]);printf("\t%.3f",b[i]);printf("\n");}temp=1;for(i=1;i<=n;i++)temp*=a[i][i];printf("\n矩阵的行列式|A|=%f\n",temp);if(temp==0){printf("\n该方程组无唯一解,请启动程序重新输入");exit(0);}return 0;}五:结果截图。

c++非方阵下的矩阵方程求解

c++非方阵下的矩阵方程求解

c++非方阵下的矩阵方程求解在C++编程中,解决矩阵方程是一个常见的任务。

通常情况下,我们思考的是一个方阵的矩阵方程。

然而,有时我们面临的是一个非方阵(m x n)的矩阵方程,其中m和n不一定相等。

我们可以使用C++中的线性代数库来解决这个问题。

一个著名的线性代数库是Eigen。

Eigen 提供了一个简单而强大的框架,可以方便地进行线性代数运算。

首先,我们需要在C++代码中包含Eigen库的头文件。

```c++#include <Eigen/Dense>```接下来,我们可以定义非方阵A和向量b,它们分别表示矩阵方程Ax = b中的系数矩阵和常数向量。

```c++Eigen::MatrixXd A(m, n); // 定义一个非方阵AEigen::VectorXd b(m); // 定义一个向量b```然后,我们可以使用Eigen的线性方程求解器来求解矩阵方程。

通过调用solve()函数,我们可以得到方程的解x。

```c++Eigen::VectorXd x = A.colPivHouseholderQr().solve(b);```注意,我们使用Eigen的colPivHouseholderQr()函数来执行QR分解并求解矩阵方程。

这是一种数值稳定的方法,在大多数情况下效果良好。

最后,我们可以输出解x。

```c++std::cout << "Solution: " << x << std::endl;```下面是一个完整的示例代码,展示了如何使用Eigen库求解非方阵的矩阵方程。

```c++#include <iostream>#include <Eigen/Dense>int main() {int m = 3; // 矩阵A的行数int n = 2; // 矩阵A的列数Eigen::MatrixXd A(m, n);Eigen::VectorXd b(m);// 为矩阵A和向量b赋值A << 1, 2, 3, 4, 5, 6;b << 1, 2, 3;// 求解矩阵方程Ax = bEigen::VectorXd x = A.colPivHouseholderQr().solve(b);// 输出解xstd::cout << "Solution: " << x << std::endl;return 0;}```这是一个简单的例子,展示了如何使用C++和Eigen库解决非方阵下的矩阵方程。

线性代数:Ax=b的解

线性代数:Ax=b的解

线性代数:Ax=b的解
n列的矩阵A,当且仅当向量b是列空间C(A)的⼀个向量时,Ax=b有解。

C(A)的零空间是N(A),N(A)正交补是A的⾏空间C(T(A)),
依据上⼀章的结论,任何Rn向量可以表⽰为r+n,其中n属于N(A),r属于C(T(A))。

因此,任何⼀个Ax=b的解可以表⽰为 x=r+n
A(r+n) = Ar+An = Ar = b,可见r也是Ax=b的解。

那么A的⾏空间⾥⾯是否有多个解。

假设存在r'使得Ar'=b, 那么有 A(r-r') = 0, r-r' 是N(A)的成员,由于r-r'⼜是C(T(A))的成员,有r-r'为零向量。

因此,A的⾏空间中有且有⼀个向量是Ax=b的解。

结合上⾯两个结论,Ax=b的解x=r+n,其中r是⼀个⾏空间中的⼀个向量,n是零空间的任意成员。

考虑解x的长度,||x||平⽅ = (r+n).(r+n) = r.r+n.n, Ax=b的解⾥⾯长度最短的向量就是r。

再具体求解⽅程的时候,得出的Ax=b的解的形式 x = v+ (t1.n1+t2.n2+..), 其中v是⼀个特解,n1,n2..是A零空间的⼀组基,t1,t2...是变参。

依据上⾯的结论 v = r+n,n=(v-r);r是⾏空间的向量,n是零空间的向量;求解r仍然是个困难。

C语言实现方阵的LU分解

C语言实现方阵的LU分解

对n 阶方阵A 的LU 分解实验目的掌握对n 阶方阵A 的LU 分解方法,并用C 语言编程实现;观察L 、U 矩阵的特征,并在此基础上得出求解方程方程组A X=b 的方法。

实验原理设矩阵*A ()ij n n a =的各阶顺序主子式均不为零,即:()0(1,2,3,)i ij a i k ≠= ,则由定理:设n 阶方阵A 的顺序主子式121,,,n D D D - 均不为零,则A 的LU 分解存在且唯一。

可得:A 有分解式 111212122212n n n n nn a a a a a a a a a ⎛⎫ ⎪ ⎪ ⎪ ⎪⎝⎭=2112111n n l l l ⎛⎫ ⎪ ⎪ ⎪ ⎪⎝⎭ 11121222n n nn u u u u u u ⎛⎫ ⎪ ⎪ ⎪ ⎪⎝⎭ ≡LU (1) 对比等式左边和右边乘积矩阵LU 的第r 行主对角元(含主对角元)的对应元素,得 1r ri rk ki k a lu ==∑ (,,;1,2,,i r n r n == (2) 再对比等式两边第r 列主对角元以下(不含主对角元)的对应元素,得1r ir ik kr k a l u ==∑ (1,,;1,2,,i r n r n =+=-(3) 当r=1时,由(2)和(3)式分别解出11i i u a = (1,2,,),i n = (4) 1111i i a l u = (2,3,,)i n = (5) 假设已求出U 的第1至r-1行,L 的第1至r-1列,由式(2)和式(3)分别得出计算U的第r 行、L 的第r 列元素的计算公式:11r ri ri rk ki k u a l u -==-∑ (,,;2,3,i r n r n== (6) 11r ir rk krk ir rr a l u l u -=-=∑ (1,,;2,3,i r n r n =+=- (7)称由(4)~(7)式所表示的矩阵分解为Doo l itt l e 分解,也称为LU 分解。

矩阵方程AX=0,AX=B有解的充要条件的一个简单证明

矩阵方程AX=0,AX=B有解的充要条件的一个简单证明

矩阵方程AX=0,AX=B有解的充要条件的一个简单证明
王跃恒;刘进波;赵人可
【期刊名称】《湖南科技学院学报》
【年(卷),期】2006(27)11
【摘要】矩阵方程 AX=0,AX=B,有解的充要条件通常用方程组理论进行证明,本文用矩阵或向量组的理论,给出一个简单的证明方法.
【总页数】2页(P56-57)
【作者】王跃恒;刘进波;赵人可
【作者单位】长沙理工大学,数学与计算机科学学院,湖南,长沙,410076;长沙理工大学,数学与计算机科学学院,湖南,长沙,410076;长沙理工大学,数学与计算机科学学院,湖南,长沙,410076
【正文语种】中文
【中图分类】O156.4
【相关文献】
1.解矩阵方程AX+YB=D与AX+XB=D的迭代方法 [J], 杨震;司步文
2.也谈四元矩阵方程AX-YB=C有解的条件及其解 [J], 明廷桥
3.矩阵方程AX=B有解的判定和解的结构 [J], 陈钦亚;姜德烁
4.矩阵方程AX-XTB=C可解的充要条件 [J], 王平心
5.矩阵方程AX=A+X有正定解和幂零解的充要条件 [J], 蒋永泉
因版权原因,仅展示原文概要,查看原文内容请购买。

  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 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,其精确解为所以各个解均在误差范围内。

相关文档
最新文档