线性方程组的数值算法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语言实现(附代码)
线性方程组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 存在一个三角分解。
C语言LU分解法实现解线性方程组
C语言LU分解法实现解线性方程组C语言LU分解法实现解线性方程组#include#include//LU分解法实现解线性方程组double sumU(double L[5][5] ,double U[5][5], int i, int j ){ double sU = 0.0;for (int k = 1; k <= i-1 ; k++){sU += L[i-1][k-1] * U[k-1][j-1];}return sU;}//计算求和1double sumL(double L[5][5] ,double U[5][5], int i, int j ){ double sL = 0.0;for (int k = 0; k <= j-1; k++){sL += L[i-1][k-1] * U[k-1][j-1];}return sL;}//计算求和2double sumY(double L[5][5] ,double y[5],int i){double sY=0.0;for (int k = 1; k <= i - 1; k++){sY += L[i-1][k-1] * y[k-1];}return sY;}//计算求和3double sumX(double U[5][5] ,double x[5],int i ,int m){double sX = 0.0;for (int k = i+1; k <= m; k++){sX += U[i-1][k-1] * x[k-1];}return sX;}//计算求和4int main(){double a[5][5] = {4,5.3,-5.6,-3,-3.4,5,-2.1,3.2,4,-8,2,-4,-7.2,-5,-2.4,5,-3,-8,2.3,3,4.2,-3,0,0,-2};//将系数存入二维数组double L[5][5] = {0};double U[5][5] = {0};//初始化部分double b[5] = {100.16,-75.72,98.2,57.1,3.72};int n = 5;//n阶//输出[Ab]printf("[A]:\n");for (int i = 1; i <= n; i++){for (int j = 1; j <= n; j++){printf("%f\t", a[i-1][j-1]);}printf("\n");}//计算L,Ufor (int i = 1; i <= n; i++){L[i-1][i-1] = 1;//对角线元素为1for (int j = i; j <= n; j++)//由于数组下标从0开始所以i-1,j-1U[i-1][j-1] = a[i-1][j-1] - sumU(L,U,i,j);if(j+1 <= n) L[j][i-1] = (a[j][i-1] - sumL(L,U,j+1,i))/U[i-1][i-1];//i变j+1，j变i }}//输出Uprintf("U:\n");for (int i = 1; i <= n; i++){for (int j = 1; j <= n; j++){printf("%f\t",U[i-1][j-1]);}printf("\n");}//输出Lprintf("L:\n");for (int i = 1; i <= n; i++){for (int j = 1; j <= n; j++){printf("%f\t",L[i-1][j-1]);}printf("\n");}//由Ly=b 求ydouble y[5] = {0.0};y[0] = b[0];//y(1) = b(1);for (int i = 2; i <= n; i++)y[i-1] = b[i-1] - sumY(L,y,i);}//由Ux=y 求xdouble x[5] = {0.0};for (int i = n; i >= 1; i--){x[i-1] =( y[i-1] - sumX(U,x,i,n))/ U[i-1][i-1]; } //输出yprintf("y:\n");for (int i = 0; i < n; i++){printf("%f\n",y[i]);}printf("\n");//输出xprintf("x:\n");for (int i = 0; i < n; i++){printf("%f\n",x[i]);}printf("\n");system("pause");return 0;}。
高斯方法解线性方程组c程序
高斯消去法和高斯主元消去法解线性方程组：高斯消元法：#include<stdio.h>#include<math.h>main(){int gauss(int n,double a[],double b[]); int i;double a[3][3]={{3,-1,4},{-1,2,-2},{2,-3,-2}}; double b[3]={7,-1,0};if(gauss(3,&a[0][0],b)!=0)for(i=0;i<=2;i++)printf("\nx[%d]=%f\n",i,b[i]);}int gauss(int n,double a[],double b[]) {int i,k,j,p,q;double d,t;for(k=0;k<=n-2;k++){d=a[k*n+k];if(d==0)return(0);for(j=k+1;j<=n-1;j++){p=k*n+j;a[p]=a[p]/d;}b[k]=b[k]/d;for(i=k+1;i<=n-1;i++){for(j=k+1;j<=n-1;j++){p=i*n+j;a[p]=a[p]-a[i*n+k]*a[k*n+j];}b[i]=b[i]-a[i*n+k]*b[k];}}d=a[(n-1)*n+n-1];if(fabs(d)+1.0==1.0){printf("fail\n");return(0);}b[n-1]=b[n-1]/d;for(k=n-2;k>=0;k--){t=0.0;for(j=k+1;j<=n-1;j++)t=t+a[k*n+j]*b[j];b[k]=b[k]-t;}return (1);}⎪⎩⎪⎨⎧=---=-+-=+-0232122743321321321x x x x x x x x x结果：x1=2，x2=1，x3=0.5高斯全选主元法：#include<stdio.h>#include<math.h>#include<stdlib.h>main(){int gauss(int n,double a[],double b[]);int i;double a[3][3]={{3,-1,4},{-1,2,-2},{2,-3,-2}}; double b[3]={7,-1,0};if(gauss(3,&a[0][0],b)!=0)for(i=0;i<=2;i++)printf("\nx[%d]=%f\n",i,b[i]);}int gauss(int n,double a[],double b[]){int *js,i,j,L,k,is,p,q;double d,t;js=malloc(n*sizeof(int));L=1;for(k=0;k<=n-2;k++){d=0.0;for(i=k;i<=n-1;i++)for(j=k;j<=n-1;j++){t=fabs(a[i*n+j]);if(t>d){d=t;is=i;js[k]=j;}}if(d+1.0==1.0)L=0;else{if(js[k]!=k)for(i=0;i<=n-1;i++){p=i*n+k;q=i*n+js[k];t=a[p];a[p]=a[q];a[q]=t;}if(is!=k){for(j=k;j<=n-1;j++){p=k*n+j;q=is*n+j;t=a[p];a[p]=a[q];a[q]=t;}t=b[k];b[k]=b[is];b[is]=t;}}if(L==0){free(js);printf("fail\n");return(0);}d=a[k*n+k];for(j=k+1;j<=n-1;j++){p=k*n+j;a[p]=a[p]/d;}b[k]=b[k]/d;for(i=k+1;i<=n-1;i++){for(j=k+1;j<=n-1;j++){p=i*n+j;a[p]=a[p]-a[i*n+k]*a[k*n+j];}b[i]=b[i]-a[i*n+k]*b[k];}}d=a[(n-1)*n+n-1];if(fabs(d)+1.0==1.0){free(js);printf("fail\n");return(0);}b[n-1]=b[n-1]/d;for(i=n-2;i>=0;i--){t=0.0;for(j=i+1;j<=n-1;j++)t=t+a[i*n+j]*b[j];b[i]=b[i]-t;}js[n-1]=n-1;for(k=n-1;k>=0;k--)if(js[k]!=k){t=b[k];b[k]=b[js[k]];b[js[k]]=t;} free(js);return(1);}结果：x1=2，x2=1，x3=0.5。
数值计算方法课程设计(C语言)
数值计算方法课程设计(C语言)数值计算方法课程设计姓名学号成绩课程实际报告实验一：秦九韶算法题目用选列主元高斯消去法解线性方程组⎪⎪⎩⎪⎪⎨⎧=+-=-+-=-+-=--02 02 0 21 34343232121x x x x x x x x x x算法语言：利用c 语言的知识编写该算法程序算法步骤叙述：秦九昭算法的基思路是v[0]=a[0]*x+a[1] v[i]=v[i-1]*x+a[i+1];利用秦九昭算法计算多项式函数。
程序清单：#include <stdio .h >void main(){float a[5],x,sum;int i;printf("presase input the value of x=");scanf("%f",&x);for (i =5;i >=0;i --){printf("please input the value of a%d=",i);scanf("%f",&a[i]);}sum=a[5];for(i=5;i>=1;i--){sum=sum*x+a[i-1];}printf("f(x)=%f/n",sum); }输出结果计算：实验总结：通过运用C 语言，解决了秦九韶算法手写的复杂。
为以后的雪地打下基础。
实验二：用选列主元高斯消去法解线性方程组题目用选列主元高斯消去法解线性方程组⎪⎪⎩⎪⎪⎨⎧=+-=-+-=-+-=--02 02 0 21 34343232121x x x x x x x x x x算法步骤叙述第一步消元——在增广矩阵（A,b ）第一列中找到绝对值最大的元素，将其所在行与第一行交换，再对（A,b ）做初等行变换使原方程组的第一列元素除了第一行的全变为0；第二步消元——在增广矩阵（A,b ）中第二列中（从第二行开始）找到绝对值最大的元素，将其所在行与第二行交换，再对（A,b ）做初等行变换使原方程组的第二列元素除了第一和第二行的全变为0；第三步消元——在增广矩阵（A,b ）中第三列中（从第三行开始）找到绝对值最大的元素，将其所在行与第三行交换，再对（A,b ）做初等行变换使原方程组的第三列第四行元素为0；第四，按x4-x3-x2-x1的顺序回代求解出方程组的解,x[n]=b[n]/a[n][n],x[i]=(b[i]-Σa[i][j]x[j])/a[i][i],i=n-1,…,2,1 程序清单：#include<iostream>#include<math>#define N 4static double A[N][N] = {-3，-1,0,0，-1,2，-1,0，0，-1,2，-1，0,0，-1,2};static double B[N]={1,0,0,0};static double X[N];int i,j,k;void main(){for(k = 0; k < N-1 ;k++){int index = k;for(i = k; i< N ;i++){if(fabs(A[index][k]) < fabs(A[i][k])){index = i;}}double temp;for( i = k ; i < N ;i++ ){temp = A[index][i];A[index][i] = A[k][i];A[k][i] = temp;}temp = B[index];B[index] = B[k];B[k] = temp;for(i = k+1; i<N; i++){double T = A[i][k]/A[k][k];B[i] = B[i] - T * B[k];for ( j = k+1 ; j < N ; j++ ){A[i][j] = A[i][j] - T * A[k][j];}}}X[N-1] = B[N-1]/A[N-1][N-1];for (i = N-2; i >=0 ; i--){double Temp = 0;for (int j = i+1; j<N ;j++)Temp = Temp + A[i][j] * X[j];X[i] = (B[i] - Temp) /A[i][i];}cout << "线性方程组的解（X1，X2，X3......Xn ）为："<<endl;for( i = 0; i < N ;i++){cout << X[i] <<" ";}} 实验总结：通过c++语言的编写过程掌握高斯消元法及选列主元元素的技术，掌握了简单的c++程序编辑语言编写算法程序实验五：二分法与牛顿法题目用二分法和Newton 迭代法求下列方程的正根：要求结果的误差限为6105.0-⨯,05.01)1ln(22=---+-x x x x x1.二分法算法语言：C 语言算法思路：算法思路先给定区间[a,b]，要求f(a)与f(b)是异号，保证区间内与x 轴有交点，求x=(a+b)/2，求f(x)，检查f(x)与f(a)是否同号，如果是同号，把x 当成新的a ，否则把x 当成新的b ，得到新的区间，重复求a 和b 的中点的值，判断与f(a)是否同号，不断循环下去，直到达到精度为止。
线性方程组求解.cpp
#include <stdio.h>#include <math.h>#include <stdlib.h>#define MAX 20 //定义最大的方程组中方程个数double a[MAX+1][MAX+2] ; //方程组系数存放的增广矩阵存放数组double b[MAX+1][MAX+2] ;int n ; //方程的实际个数int xxx = 0; //是否初始化void init()//输入线性方程组{xxx= 1;printf("\n\t请输入方程组的方程个数: ");scanf ("%d",&n);if( n > MAX )//输入的数超过了定义的最大方程个数，退出{printf("sorry,can not deal !\n\n");printf("最大能处理的方程组为%d阶方程组\n\n,MAX");exit(0);}printf("\n请依次输入方程组的方程的系数(顺序从上到下，从左至右)：\n"); for(int i = 1 ; i < n+1 ; i++ )//n行n+1列的增广矩阵for(int j = 1 ; j < n+2 ; j++ )scanf ("%lf",&a[i][j]);for(int i = 1 ; i < n+1 ; i++ ) //备份for(int j = 1 ; j < n+2 ; j++ )b[i][j]=a[i][j];}void reset()//将数组a还原成初始状态{for(int i = 1 ; i < n+1 ; i++ )for(int j = 1 ; j < n+2 ; j++ )a[i][j]=b[i][j];}void display()//输出增光矩阵{for(int i = 1 ; i < n+1 ; i++ ){for(int j = 1 ; j < n+2 ; j++ )printf("%.3f ",a[i][j]);printf("\n");}}void print()//输出线性方程组{int count=0;//记录某一行中系数为0的个数for(int i = 1 ; i < n+1 ; i++ )for(int j = 1 ; j < n+2 ; j++ ){if( a[i][j] == 0 && j != n+1)count ++ ;else{if( j == 1 || count == j-1)//系数为第一项或者该项之前的所有项的系数都是0printf ("%.3f X%d ",a[i][j],j);//小数后面保留三位elseif( j == n+1 ){printf(" = %.3f\n",a[i][n+1]);count =0;}elsea[i][j]>0 ? printf (" + %.3f X%d ",fabs(a[i][j]),j): printf (" - %.3f X%d ",fabs(a[i][j]),j);}}}void selectMainEle(int k)//选取第k列的主元{double mainEle = fabs(a[k][k]) ,temp=0;int j = k ;//用变量j记住主元初始位置for(int i = k+1 ; i < n+1 ; i++ )//选出绝对值最大的数if( mainEle < fabs(a[i][k])){mainEle = fabs(a[i][k]) ;j = i; //j记录主元（列中的最大数）的行数}if( k != j ) //初始值不是最大值，k,j换行for(int i = 1 ; i < n+2 ; i++ ){temp = a[k][i] ;a[k][i] = a[j][i] ;a[j][i] = temp ;}}void printit()//输出方程组的根{printf("\n\n求解的方程组的根为：\n");for(int i=1;i<n+1;i++){if(i%4==0) printf("\n");//每一行输出三个解 printf("X%d=%.3f ",i,a[i][n+1]);}printf("\n");}void gause()//高斯列主消元法核心{printf("线性方程组为：\n");print();//输出输入的方程组for( int k=1 ; k < n ; k ++)//消元过程{selectMainEle(k);for(int i = k+1 ; i < n+1 ; i++ )a[i][k] /= a[k][k];for(int i = k+1 ; i < n+1 ; i++ )for(int j = k+1 ; j < n+2 ; j++ )a[i][j] -= a[i][k] * a[k][j] ;printf("\n消元过程如下：\n");for(int i=1;i<n+1;i++){for(int j=1;j<n+1;j++)if(j<k+1&&i>j)printf("%.4f ",0.0);elseprintf("%.4f ",a[i][j]);printf("\n");}}a[n][n+1] /= a[n][n] ;//回代过程for(int i=n-1 ; i>0 ; i--){double sum=0.0;for(int j=i+1 ;j<n+1 ;j++)sum += a[i][j] * a[j][n+1];a[i][n+1] = (a[i][n+1] - sum)/a[i][i];}printit();}void gauseyuedang()//高斯-约当列主元消去法，高斯列主元消去法的改进{printf("线性方程组为：\n");print();//输出输入的方程组double temp;for( int k=1 ; k < n+1 ; k++)//消元，只保留主元{selectMainEle(k);temp = a[k][k];//用temp记住主元的值for(int i=k ; i < n+2 ; i++)//将主元变为1a[k][i] /= temp;//开始消元，出去主元第k列元素全部消去变为零for(int i = 1 ; i < n+1 ; i++ )//行控制if( i != k ) //行变换不包括当前主元所在行for(int j = k+1 ; j < n+2 ; j++ )//列控制a[i][j]=a[i][j]-a[k][j]*a[i][k];printf("\n消元过程如下：\n");for(int i=1;i<n+1;i++){for(int j=1;j<n+1;j++)if(j<k+1&&i!=j)printf("%.4f ",0.0);elseprintf("%.4f ",a[i][j]);printf("\n");}}printit();//输出方程组的解}void LU()//LU分解法{double sum ;printf("线性方程组为：\n");print();//输出输入的方程组for(int i=2;i<n+1;i++)a[i][1] /= a[1][1];for(int r=2 ; r<n+1 ; r++){for(int i=r ; i<n+1 ; i++)//U变换{sum = 0;for(int k=1;k<r;k++)sum += a[r][k]*a[k][i];a[r][i] -= sum;}for(int i=r+1 ; i<n+1 ; i++)//L变换{sum = 0;for(int k=1 ; k<r ; k++)sum += a[i][k]*a[k][r];a[i][r] = (a[i][r] - sum)/a[r][r];}}printf("\n\nLU分解中的L为：\n");for(int i=1 ; i < n+1 ; i++){for(int j=1 ;j < n+1 ; j++)if(i < j) //上部为0printf("%.3f ",0.0);elsei==j ? printf("%.3f ",1.0): printf("%.3f ",a[i][j]); //主对角线为1 printf("\n");}printf("\n\nLU分解中的U为：\n");for(int i=1 ; i < n+1 ; i++){for(int j=1 ;j < n+1 ; j++)if(i > j) //上部为0printf("%.3f ",0.0);elseprintf("%.3f ",a[i][j]);printf("\n");}for(int i=2;i<n+1;i++)//第一次回带求解{sum=0;for(int k=1;k<i;k++)sum+=a[i][k]*a[k][n+1];a[i][n+1] -= sum;}printf("\n\n第一次求解的方程组(LY=b)的根为：\n");for(int i=1;i<n+1;i++){if(i%4==0) printf("\n");//每一行输出三个解printf("Y%d=%.3f ",i,a[i][n+1]);}a[n][n+1]/=a[n][n];//Xn的值for(int i=n-1;i>0;i--){sum=0;for(int k=i+1;k<n+1;k++)sum+=a[i][k]*a[k][n+1];a[i][n+1] = (a[i][n+1] - sum)/a[i][i];}printf("\n\n第二次求解的方程组(UX=y)的根，即原方程组的解为：\n"); for(int i=1;i<n+1;i++){if(i%4==0) printf("\n");//每一行输出三个解printf("X%d=%.3f ",i,a[i][n+1]);}printf("\n");}void Jacobi()//雅可比迭代法{printf("线性方程组为：\n");print();//输出输入的方程组int count=0;///迭代次数double latest[n+1],old[n+1],temp[n+1],sum,precision,max,temp2;for(int i=1;i<n+1;i++)if(a[i][i] == 0){printf("sorry，主对角线中有数据为0，不能用Jacobi迭代法处理！\n请返回重新选择");return ;}printf("\n输入求解精度值: ");scanf("%lf",&precision);for(int i=1 ; i < n+1 ; i++) //初始状态解为0{latest[i] = 0.0;//新值，下一次迭代结果old[i] = 0.0;//旧值，上一次的迭代结果}printf("Jacobi迭代法迭代过程如下:\n");do{ max =0;printf("\n\n第%d次迭代的结果: ",count++);for(int i=1 ; i<n+1 ; i++){if(i%4==0) printf("\n");//每一行输出三个解printf("X%d=%.4f ",i,old[i]);}for(int i = 1 ; i < n+1 ; i++)//迭代过程,完成新的一组解{sum =0.0;for(int j = 1 ; j < n+1 ; j++)if(j != i) sum = sum + a[i][j] * old[j];latest[i]=(a[i][n+1] - sum)/a[i][i];temp[i] = latest[i];//Y用来记录本次迭代的值temp2 = fabs(latest[i] - old[i]); //新值减旧值if(max < temp2) max = temp2 ;}printf("最大的偏差为：%.4f\n",max);for(int i=1 ; i<n+1 ; i++)old[i] = temp[i]; //将新值保留到old数组} while( max > precision && count <= 100);//前后两次解的偏差满足精度要求或者迭代次数超过100次printf("\n\n第%d次(最后的)求解的方程组的根为：\n",count);for(int i=1;i<n+1;i++){if(i%4==0) printf("\n");//每一行输出三个解printf("X%d=%.4f ",i,old[i]);}printf("\n");}void GS()//高斯-赛德尔迭代法{printf("线性方程组为：\n");print();//输出输入的方程组double precision,sum,X[n+1],max,temp,old;int count=0;///迭代次数for(int i=1;i<n+1;i++)if(a[i][i] == 0){printf("sorry，主对角线中有数据为0，不能用高斯赛德尔迭代法处理！\n请返回重新选择");return ;}printf("\n输入求解精度值: ");scanf("%lf",&precision);for(int i = 1 ; i < n+1 ; i++)// 解初始化为0X[i] = 0.0 ;printf("高斯赛德尔迭代法迭代过程如下:\n");do{printf("\n\n第%d次迭代的结果: ",count++);for(int i=1 ; i<n+1 ; i++){if(i%4==0) printf("\n");//每一行输出三个解printf("X%d=%.5f ",i,X[i]);}for(int i = 1 ; i < n+1 ; i++)//迭代过程,完成新的一组解{sum = 0.0;max = 0.0;old =X[i];for(int j = 1 ; j < n+1 ; j++)if(j != i) sum += a[i][j] * X[j] ;X[i] = (a[i][n+1] - sum) /a[i][i];temp = fabs(X[i]-old);if(max < temp) max = temp;}printf("最大的偏差为：%.5f\n",max);} while(max > precision && count <= 100);//前后两次解的偏差满足精度要求或者迭代次数超过100次printf("\n\n第%d次(最后的)求解的方程组的根为：\n",count);for(int i=1;i<n+1;i++){if(i%4==0) printf("\n");//每一行输出三个解printf("X%d=%.5f ",i,X[i]);}printf("\n");}int main(){int xz;while(1){printf("\t\t\t*****************************************\n");printf("\t\t\t**** \t 0.更改线性方程组 ****\n");printf("\t\t\t**** \t 1.高斯列主消元法 ****\n");printf("\t\t\t**** \t 2.高斯-约当列主消去法 ****\n");printf("\t\t\t**** \t 3.LU分解法 ****\n");printf("\t\t\t**** \t 4.雅各比迭代法 ****\n");printf("\t\t\t**** \t 5.高斯-赛德尔迭代法 ****\n");printf("\t\t\t**** \t 6.退出 ****\n");printf("\t\t\t*****************************************\n");printf("\t请根据相应的代号，从上述方法中选择其一： ");scanf("%d",&xz);switch(xz){case 0: init(); break;case 1: if(xxx){reset();gause();}else{printf("sorry,还没有输入线性方程组!\n");init();}break;//高斯列主消元法case 2: if(xxx){reset();gauseyuedang();}else{printf("sorry,还没有输入线性方程组!\n");init();}break;//高斯-约当列主元消去法case 3: if(xxx){reset();LU();}else{printf("sorry,还没有输入线性方程组!\n");init();}break;case 4:if(xxx){reset();Jacobi();}else{printf("sorry,还没有输入线性方程组!\n");init();}break;case 5: if(xxx){reset();GS();}else{printf("sorry,还没有输入线性方程组!\n");init();}break;case 6: exit(0);default:printf("\n\t错误输入o(︶︿︶)o 唉！\n\n");break;}system("pause");system("cls");}}。
线性方程组C语言编程
}
}
void main()
{
int code[100]; /* 输入符号标记 */
float b[100]; /* 方程右值 */
Input(b,code); /* 初始化 */
Sstart(b,code); /* 化标准型 */
}
return temp;
}
void JustArtificial()
{
int i;
for(i=m+indexe+indexl;i<s;i++)
if(fabs(x[i])>=0.000001){
printf("No Answern");
return;
}
}
int Check(int in)
for(i=0;i<m;i++)
cin>>matrix[n][i];
if(type==1)
for(i=0;i<m;i++)
matrix[n][i]=-matrix[n][i];
}
void Xartificial()
{
int i,j,k;
if(indexg!=0){
for(i=m+indexe+indexl;i<s;i++){
int m,n,s,type; /* 方程变量,约束数,求最大最小值的类型,0:最小 1:最大 */
int indexe,indexl,indexg; /* 剩余变量,松弛变量,人工变量 */
用C语言解决线性规划问题(用单纯形法解)
用C语言解决线性规划问题（用单纯形法解）#include<stdio.h>#include<math.h>#include<iostream.h>#define BORDER -0.00001#define M 100int main(){int k; //初始变量的个数int m; //约束条件的个数;cout<<"输入初始变量的个数"<<endl;cin>>k;cout <<"输入约束条件的个数"<<endl;cin>>m;cout<<"输入约束条件的种类"<<endl;cout<<"0 means <="<<endl;cout<<"1 means >="<<endl;cout<<"2 means ="<<endl;int NGET=0;int NLET=0;int NET =0;int I =0;//人工变量；int *code=new int[m]; //the array of the >= <= =;for(int i=0;i<m;i++){cin>>code[i];if(code[i]==0)NGET++;if(code[i]==1)NLET++;I++ ;if(code[i]==2)NET++;I++ ;}int n; //变量总和；n=k+2*NLET+NGET+NET;float *Index=new float[n+1]; //目标函数的系数;float *c=new float[n+1];int NTYPE;for(i=0;i<n+1;i++)Index[i]=0.0;//为语法要求定初值；cout<<"输入目标函数的系数"<<endl;for(i=0;i<k;i++)cin>>Index[i];for(i=k+NGET+NLET;i<n;i++)//人工变量所在的列；Index[i]=-M; //the initionalization of indexes of manul variable;cout<<"输入所求函数的类型是求最大还是最小"<<endl;cout<<"0 means min"<<endl;cout<<"1 means max"<<endl;cin>>NTYPE;if(NTYPE) //the initionalization of the c[i];{for(i=0;i<n+1;i++)c[i]=-Index[i];//一般情况下只要将目标函数系数的相反数输入；}if(!NTYPE){for(i=0;i<n+1;i++)if(i<k+NGET+NLET)c[i]=Index[i];else c[i]=-Index[i];}delete []Index;float **a=new float*[m+1]; //the array of all the variable to compute;for(i=0;i<m+1;i++)a[i]=new float[n+1];int INDEXG=k;int INDEXL=k+NGET;int INDEXE=k+NGET+NLET;int *ARTV=new int[I]; //保存人工变量；for(i=0;i<m;i++) //the analization of the code[](<= >= = );{if(code[i]==0){a[i][INDEXL]=1.0;INDEXL++;}if(code[i]==1){a[i][INDEXE]=1.0;INDEXE++;a[i][INDEXG]=-1.0;INDEXG++;ARTV[I]=i;I++;}if(code[i]==2){a[i][INDEXE]=1.0;INDEXE++;ARTV[I]=i;I++;}}if( (INDEXG!=k+NLET) || (INDEXL!=k+NGET+NLET) || (INDEXE!=n) )//excption {return -1;}cout<<"输入约束表达式左边的系数"<<endl;for(i=0;i<m;i++)for(int j=0;j<k;j++)cin>>a[i][j];float *b=new float[m];cout<<"输入约束表达式右边的值" <<endl;for(i=0;i<m;i++)cin>>b[i];for(i=0;i<m;i++)a[i][n]=b[i];float *temp=new float[n+1];if(I){for(i=0;i<I;i++)for(int j=0;j<n+1;j++){temp[j]=-a[ARTV[i]][j];c[j]+=M*temp[j];}}for(i=0;i<n+1;i++)a[m][i]=c[i];for(i=0;i<m+1;i++){for(int j=0;j<n+1;j++)cout<<a[i][j]<<" ";cout<<endl;}int flag=0;float temp1;float temp2;int K,J;int index;for(i=0;i<n;i++)if(a[m][i]<0)flag=1; //检验系数;while(flag) //Using a[][] to compute the result; {temp1=0;for(i=0;i<n;i++){if(temp1>a[m][i]){temp1=a[m][i];K=i;}}temp2=M;for(i=0;i<m;i++){if(a[i][K]>0&&(a[i][n]/a[i][K])<temp2){temp2=a[i][n]/a[i][K];J=i;}}if(temp2==M){cout<<"无解!"<<endl;return -1;}float temp3=a[J][K];for(i=0;i<n+1;i++){a[J][i]=a[J][i]/temp3;}for(i=0;i<m+1;i++){if(i!=J){float temp4=a[i][K];for(int j=0;j<n+1;j++){a[i][j]=a[i][j]- a[J][j]*temp4;}}cout<<endl;}flag=0;for(i=0;i<n+1;i++){if(a[m][i]<BORDER)flag=1;}cout<<"**************************************"<<endl;for(i=0;i<m+1;i++){for(int j=0;j<n+1;j++)cout<<a[i][j]<<" ";cout<<endl;}cout<<"***************************************"<<endl;getchar();}if(NTYPE)cout<<"The answer is :"<<a[m][n];if(!NTYPE)cout<<"The answer is :"<<-a[m][n];getchar();return 0;}。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
线性方程组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 存在一个三角分解。
而且，L 的对角线元素为1，U 的对角线元素非零。
得到L 和U 后，可通过以下步骤得到X ：(1)、利用前向替换法对方程组LY B =求解Y 。
(2)、利用回代法对方程组UX Y =求解X 。
4、雅可比迭代：考察一般形式的线性方程组：,1,2,3,N ij j i ja xb i N ==∑L L (2)设从(2)中分离变出变量i X 将它改写成：111,2,3,Ni i ij jj iij i x b a x i N a =≠⎛⎫ ⎪=-= ⎪ ⎪⎝⎭∑L L (3)据此建立雅可比迭代公式：(k 1)(k)111,2,3,N i i ij jj ii j i x b a x i N a +=≠⎛⎫ ⎪=-= ⎪ ⎪⎝⎭∑L L (4)5、高斯-赛德尔迭代：通过对雅可比迭代的观察，由于通常1k x +被认为是比k x 更好的x 的近似值，所以在计算1k y +时用1k x +来替换k x 是合理的。
所以对雅可比迭代进行了改进，得到了高斯-赛德尔迭代，其收敛速度更快。
它的迭代公式为： (k 1)(k 1)111,2,3,Ni i ij jj ii j i x b a x i N a ++=≠⎛⎫ ⎪=-= ⎪ ⎪⎝⎭∑L L (5)三、 实验内容：1、许多科学应用包含的矩阵带有很多零。
在实际情况中很重要的三角形线性方程组有如下形式: 111211122232223334322111111N N N N N N N N N N NNd x c x b a x d x c x b a x d x c x b a x d x c x b a x d x b --------+=++=++=++=+=OO O M (6)构造一个程序求解三角形线性方程组。
可假定不需要行变换，而且可用第k 行消去第k+1行的k x 。
2、使用C 语言编写一个程序求解线性方程组AX B =，其中：1357213500252631A -⎡⎤⎢⎥-⎢⎥=⎢⎥⎢⎥---⎣⎦ 1234B ⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎣⎦3、使用2中的程序求解线性方程组AX B =，其中=1j ij N N A a i -⨯⎡⎤==⎣⎦ ；而且1ij N B b ⨯⎡⎤=⎣⎦，11b N =当2i ≥时，1(i 1)/(i 1)N i b =--。
对3,7,11N =的情况分别求解。
精确解为[]1111X '=L。
对得到的结果与精确解的差异进行解释。
4、修改2中的程序，使得它可以通过重复求解N 个线性方程组J J AC E = 其中1,2,,J N =L L (7)来得到1A - ，则1212[C C C ][E E E ]N N A =L L (8)而且112[C C C ]N A -=L (9)5、设有如下三角线性方程组，而且系数矩阵具有严格的对角优势：111211122232223334322111111N N N N N N N N N N NNd x c x b a x d x c x b a x d x c x b a x d x c x b a x d x b --------+=++=++=++=+=OO O M (10)（1）、根据解线性方程组的迭代法原理，设计一个算法来求解上述方程组，算法必须有效利用系数矩阵的稀疏性。
(2)、用（1）中的程序解下列的三角线性方程组：(a)12123234345484950495043+43+43+434343m m m m m m m m m m m m m m m m +=+=+=+=++=+=M M M M (b)12123234345484950495041+42+41+424142m m m m m m m m m m m m m m m m +=+=+=+=++=+=M M M M6、利用高斯-赛德尔迭代法求解下列带状方程：1231234123452345646474849504748495048495012252122521225212521225212252125x x x x x x x x x x x x x x x x x x x x x x x x x x x x x -+=-+-+=-+-+=-+-+=-+-+=-+-=-+=M MMM M M四、计算结果及分析：1、此题由于没有提供具体的实例计算，所以,以《数值方法(MATLAB 版)》(第四版)第107页的11、12题为例，进行代码的计算验证。
此题以高斯消去法作为算法的核心，并依此来进行程序的编写。
实例：第11题:求解下列线性方程组。
121232343427239423102412x x x x x x x x x x +=+-=++=-=第12题: 求解下列线性方程组。
1212323434525934219262x x x x x x x x x x +=-+=--+=+=第11题的正确解为:[]1322x '=- 第12题的正确解为:[]2321x '=-精确的结果。
2、此题以带选主元的三角分解法作为算法的核心，以此进行程序的编写。
验算:将x 代入AX B =中，算出[]1.0000 2.0001 3.0000 4.0000B '= 与实际的B 进行比较，此算法在保留四位小数的时候，得出的解比较的正确。
3经过程序计算可看出，在3,7N =时得出的结果比较的精确，与所给的答案一致。
但在11N =时，得出的结果与实际的结果相差很大。
可能由于计算时迭代次数较多，机器误差累计造成。
4、此题由于没有提供具体的实例计算，所以,以《数值方法(MATLAB 版)》 (第四版)第119页的第4(b)题为例，进行代码的计算验证。
J 1-27A =42125-2⎡⎤⎢⎥⎢⎥⎢⎥⎣⎦ 1E 11J ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦经验算，可推出[]1012A C C C -=。
5、此题用高斯-赛德尔迭代法作为算法的核心，以此编写程序。
通过对给出的两个实例的计算：6、用5中的程序，通过高斯-赛德尔迭代求解此带状方程解为：五、实验结论：经过C程序计算的得出的数据，经验算后，都比较好的符合了计算的要求，这同时也证明了算法的正确性。
通过本次实验，了解科学计算中的解线性方程组的一些基本方法，通过亲自设计算法，对这些科学的方法有了更加深刻的认识。
对以后的学习有很大帮助。
附件：一、(高斯消去法)：#include<stdio.h>#include<stdlib.h>//头文件声明//float main(){float *a,*b,*c,*d,*x,**A;FILE *fp;//定义指向文件的指针//float temp,sum;int i,j,N;fopen_s(&fp,"D:\\vs程序\\3.4.6-1.12.txt","w" );//确定文件指针的指向// printf("请输入此方程组的元的个数：\n");scanf_s("%d",&N);a=(float *)malloc((N-1)*sizeof(float));//动态定义一维数组，下同//b=(float *)malloc(N*sizeof(float));c=(float *)malloc((N-1)*sizeof(float));d=(float *)malloc(N*sizeof(float));x=(float *)malloc(N*sizeof(float));A=(float **)malloc(N*sizeof(float *));//动态定义二维数组//for (i=0;i<N;i++)A[i]=(float *)malloc(N*sizeof(float));for(i=0;i<N;i++){for(j=0;j<N;j++)A[i][j]=0;//先使系数矩阵为一个0矩阵//}printf("请输入对角系数d[i]:\n");for(i=0;i<N;i++){scanf_s("%f",&d[i]);//输入对角元素的值//}printf("请输入系数a[i]:\n");for(i=0;i<N-1;i++){scanf_s("%f",&a[i]);//输入对角线下的元素系数//}printf("请输入系数c[i]:\n");for(i=0;i<N-1;i++){scanf_s("%f",&c[i]);//输入对角线下的元素系数//}for(i=0;i<N;i++){A[i][i]=d[i];//给系数矩阵的对角元素重新赋值//}for(i=0;i<N-1;i++){A[i+1][i]=a[i];A[i][i+1]=c[i];}//给系数矩阵重新赋值，使之成为三角形线性方程组//printf("请输入系数b[i]:\n");for(i=0;i<N;i++){scanf_s("%f",&b[i]);//输入题中的b矩阵//}for(i=0;i<N-1;i++){temp=(*(A[i+1]+i))/(*(A[i]+i));//对角元素所在列的下一行元素比上此对角元素，求出其比例关系//b[i+1]=b[i+1]-b[i]*temp;//根据上条指令求出的比例关系，进行行变换//for(j=0;j<N;j++)A[i+1][j]=A[i+1][j]-A[i][j]*temp;//进行行变换，已达到消元的目的// }x[N-1]=b[N-1]/A[N-1][N-1];//经过上面的行变换，求出只有一个未知数的那一行的未知数大小//for(i=2;i<=N;i++){sum=0;sum+=A[N-i][N+1-i]*x[N+1-i];x[N-i]=(b[N-i]-sum)/A[N-i][N-i];}//此循环根据书上的回代公式和上面求出的未知数的值，求出输入的方程组的解//if (fp!=NULL){fprintf(fp,"此方程组的解为：");for(i=0;i<N;i++)fprintf(fp,"\nx[%d]=%9.4f",i,x[i]);//将计算结果，输出到文件指针所指向的文件中去//fclose(fp);//关闭文件//}free(a);//释放动态定义指针所指向的空间，便于下次使用。