求解线性方程组,用的是高斯列主元消元法
求解线性方程组#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define ROW 100
#define COL 100
int guass(double *[ROW], double [ROW],int, int);
void exchange(double *a[ROW], int flag, int row);
int main()
{
FILE *fp1,*fp2;
double a[ROW][COL],b[ROW];
int i=0,j=0,number=0,col=0,row=0;
char ch;
//打开文件
if((fp1=fopen("d:\\data1.txt","r"))==NULL)
{
printf("File open error!\n");
exit(0);
}
if((fp2=fopen("d:\\data2.txt","r"))==NULL)
{
printf("File open error!\n");
exit(0);
}
//将线性方程的数据读入到数组中
fseek(fp1,0,0);
for(i=0; i<100; i++)
{
for(j=0; j<100; j++)
{
fscanf(fp1,"%lf",&a[i][j]);
}
fscanf(fp1,"\n");
}
ch=fgetc(fp1);
while(!feof(fp1))
{
if (ch>='0'&&ch<='9')
number++;
ch=fgetc(fp1);
}
for(i=0; i<100; i++)
{
fscanf(fp2,"%lf",&b[i]);
fscanf(fp2,"\n");
}
ch=fgetc(fp2);//计算系数矩阵的行数
while(!feof(fp2))
{
if (ch>='0'&&ch<='9')
row++;
ch=fgetc(fp2);
}
col=number/row;//计算系数矩阵的列数
guass(a,b,row,col+1);
//进行回代求解
double x[row];
x[row-1]=b[row-1]/a[row-1][row-1]; //因矩阵化为了上三角矩阵,即可求出最后一行的解,接下来从最后一行往上一行进行回代。
int sum=0;
for (i=row-2; i>=0; i--)
{
for(j=i+1; j<col; j++)
sum+=a[i][j]*x[j]; //累加和。
x[i]=b[i]-sum;
x[i]=x[i]/a[i][i];
}
printf("solution is :\n");
for(i=0;i<row;i++)
printf("%lf\n",x[i]);
//关闭文件
if (fclose(fp1))
{
printf("Can not close the file!\n");
exit(0);
}
if (fclose(fp2))
{
printf("Can not close the file!\n");
exit(0);
}
return 0;
}
void exchange(double *a[ROW], int flag, int row)
{
int i;
double *temp;
for(i=flag+1; i<row; i++)
//同一列的数依次进行比较大小
if(fabs(*(a[flag]+flag))<fabs(*(a[i]+flag)))
{
//将每一列的最大值放至主对角线上
temp = a[flag];
a[flag] = a[i];
a[i] = temp;
}
return;
}
int guass(double *a[ROW],double b[ROW],int row,int col)
{
int temp=0,j;
for (int i=0; i<row-1; i++)
{
exchange(a,i,row);
for(j=i+1; j<row; j++)
temp=-*(a[j]+i)/ *(a[i]+i);
b[j]+=b[i]*temp;
for(int k=i; k<col; k++)
*(a[j]+k) += temp * *(a[i]+k);
}
return;
}
[此贴子已经被作者于2019-6-10 22:23编辑过]