矩阵迭代(Jacobi)

//本程序是用Jacobi迭代法求一个系数矩阵非奇异的方程组的近似解


#include<stdio.h>
#include<stdlib.h>
#include<math.h>

#define R 3 //宏定义矩阵行数
#define N 4 //宏定义矩阵列数
#define P 0.0000001 //宏定义计算的精度

float G[R][N-1]; //定义G[R][N-1]二维数组用来存放变形后的系数矩阵
float b[R]; //定义b[R]用来存放变形后的常数项
float Answer[R]; //定义外部数组用来存放Jacobi()传回的方程组的解

void Mat_Input(float a[R][N]); //自定义矩阵输入函数
void Jacobi(float a[R][N]); //自定义Jacobi迭代函数
void Mat_Change(float a[R][N]); //自定义矩阵变换函数
float Mat_Norm(float a[R][N]); //自定义求矩阵无穷范数函数
float Judge(float A[N]); //自定义找出相邻两个解向量的差的无穷范数
void Display(float G[R][N-1],float b[R] ); //自定义显示G和b的函数


int main (void) //主函数
{
int i;
float a[R][N]; //定义需要输入的矩阵
Mat_Input(a); //调用输入函数按要求输入矩阵
if(R!=(N-1)) //第一次判断，若系数矩阵的行数与列数不相等，则无法用此方法求解
{
printf("错误：该程序不能求解该方程组!"); //显示错误
}
else //否则执行下面的操作
{
Mat_Change(a); //调用矩阵变换函数，进行变形，x=G|x|+f,并把G传到 mat_0[R][N-1]，把f传到X_0[R]，以方便调用
Display(G,b); //调用Display函数显示G和b
if(Mat_Norm(a)>=1) //调用Mat_Norm()函数判断变形后的系数矩阵是否不小于1，若是则无法用此方法求解
{
printf("错误：该程序不能求解该方程组!"); //若无法求解，则显示错误提示
}
else //否则，就调用Jacobi函数进行求解
{
Jacobi(a);
}
printf("利用Jacobi迭代求得的方程组的解为:\n"); //显示最后计算结果
for(i=0;i<R;i++)
printf("x_%d=%f\n",(i+1),Answer[i]);
}
system("PAUSE");
return 0;
}


void Mat_Input(float a[R][N]) //实现矩阵的输入函数
{
int i,j;
printf("请输入矩阵，其输入格式如下:\n依次输入a11->回车->a12->回车->....->ann->回车。\n"); //提示按提示要求输入矩阵
for(i=0;i<R;i++) //设置循

环逐个输入元素
{
for(j=0;j<N;j++)
{
scanf("%f",&a[i][j]);
}
}

printf("
\n输入的增广矩阵为:\n"); //显示输入的矩阵
for(i=0;i<R;i++)
{
for(j=0;j<N;j++) //设置循环按格式显示出输入的矩阵
{
printf("%.6f\t ",a[i][j]);
}

printf("\n"); //用于矩阵的一行输入完后换行
}
}


void Mat_Change( float a[R][N]) //实现矩阵的变形函数
{
int i,j;
for(i=0;i<R;i++) //设置行的循环
{
G[i][i]=0; //变形后的矩阵G的对角线上的元素值都为零，即mat_x.mat_0[i][i]
for(j=0;j<i;j++) //计算出非对角线元素的值
{
G[i][j]=(-1*a[i][j]/a[i][i]); //j<i时mat_x.mat_0[i][j]=(-1*a[i][j]/a[i][i])
}
for(j=(i+1);j<(N-1);j++) //j>i时mat_x.mat_0[i][j]=(-1*a[i][j]/a[i][i])
{
G[i][j]=(-1*a[i][j]/a[i][i]);
}

}
for(i=0;i<R;i++) //算出变形后的常数项并传到mat_x.X_0[i]中
{
b[i]=(a[i][N-1]/a[i][i]);
}

}


float Mat_Norm(float a[R][N]) //实现求变形后的矩阵G的范数
{
float c[R],sum,max;
int i,j;
for(i=0;i<R;i++) //计算出每一行元素的绝对值和并存到数组c中
{
for(j=0,sum=0;j<N-1;j++)
{
sum+=fabs(G[i][j]);
}
c[i]=sum;
}
for(i=0;i<R;i++) //找出最大值，即求出了无穷范数
{
max=c[0];
if(max<fabs(c[i]))
{
max=fabs(c[i]);
}
}
return max; //返回最大值即返回无穷范数
}

float Judge(float A[R]) //实现找出相邻两个解向量的差的无穷范数
{
float q;
int i;
q=fabs(A[0]);
for(i=0;i<R;i++) //找出两个解向量的差的无穷范数
{
if(q<fabs(A[i]));
q=fabs(A[i]);
}
return q; //返回求得的值
}



void Jacobi(float a[R][N]) //实现Jacobi迭代
{
float X_1[R],X_2[R],A[R],m;
int i,j;
for(i=0;i<R;i++) //给出始解向量赋值为0
{
X_1[i]=0;
}
for(i=0;i<R;i++) //给第一次迭代的的解赋值
{
for(j=0,m=0;j<N-1;j++)
{
m+=X_1[i]*G[i][j]; //系数矩阵的第i行和解向量相乘的值
}
X_2[i]=m+b[i]; //m+改行常数项

}
for(i=0;i<R;i++)
{
A[i]=X_2[i]; //A[i]代表相邻解向量的差，由
于第一个解向量为0，那么第一个相邻解向量的差即第一次迭代出的解向量
}
while(Judge(A)>P) //设置循环开始迭代（当A[i]的无穷范数大于最大允许的误差，则继续迭代）
{
for(i=0;i<R;i++) //下一次迭代开始时把X_2[i]赋值给X_1[i],而X_2[i]用来装迭代出来的最新解向量
{
X_1[i]=X_2[i];
}
for(i=0;i<R;i++)
{
for(j=0,m=0;j<N-1;j++)
{
m+=X_1[j]*G[i][j]; //第i行与解向量相乘
}
X_2[i]=m+b[i]; //m+变形后的常数项
}
for(i=0;i<R;i++) //每次计算完成后计算相邻解向量的差
{
A[i]=(X_2[i]-X_1[i]);
}
}
for(i=0;i<N;i++) //直至迭代出符合要求的解，并把解向量传值到外部数组Answer中
{
Answer[i]=X_2[i];
}
}



void Display(float G[R][N-1],float b[R] )
{
int i,j;
printf("改写后的生成的矩阵G为:\n"); //显示矩阵G
for(i=0;i<R;i++)
{
for(j=0;j<N-1;j++)
{
printf("%f\t",G[i][j]);
}
printf("\n");
}
printf("改写后的生成的向量b为:\n");
for(i=0;i<R;i++)
{
printf("%f\t",b[i]);
}
printf("\n");
}