雅可比迭代法:
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#define eps 1e-4//题目要求精度
#define MaxSize 100//最大迭代次数
#define n 6//矩阵的维数
double a[n][n] = {
{4,-1,0,-1,0,0},
{-1,4,-1,0,-1,0},
{0,-1,4,-1,0,-1},
{-1,0,-1,4,-1,0},
{0,-1,0,-1,4,-1},
{0,0,-1,0,-1,4} };
//如果是一列的向量可以用一维数组存储
double b[n]={0,5,-2,5,-2,6};//b
double X[n]={0,0,0,0,0,0};//解向量
double f[n]={0,0,0,0,0,0};//Invdiag*b
double Invdiag[n][n]={0};//D的逆
double Mat_Jacobi[n][n] = { 0 };//B=Invdiag*(L+U)
double temp[n]={0,0,0,0,0,0};//保存旧的X
double add1[n]={0,0,0,0,0,0};//保存B*X
double d;
//计算系数矩阵的Jacobi迭代矩阵B
void Jacobi_Mat()
{
int i,j;
for(j=0;j<n;j++)
{
if(a[j][j] != 0)
{
d=a[j][j];
Invdiag[j][j] = 1/a[j][j];//D的逆矩阵
for(i=0;i<n;i++)
{
if(i!=j)
{
Mat_Jacobi[i][j]=(-a[i][j])/d;//B矩阵
}
else
{
Mat_Jacobi[i][j] = 0;
}
}
}
else
{
printf("对角元素存在零元,无法求取Jacobi迭代矩阵!\n");
return;
}
}
}
//由于列向量的特殊性,二范数相当于各分量的绝对值的平方和的开方
double GetNorm()
{
int i;
double sum=0;
double e;
for (i=0;i<n;i++)
{
sum=sum+fabs(X[i]-temp[i])*fabs(X[i]-temp[i]);//计算X1-X0
e=sqrt(sum);//计算X1-X0的二范数
}
return e;
}
//雅可比迭代公式
void Jacobi()
{
//i表示行,j表示列,k用来计算迭代次数
int i,j,k;
double max,sum;
//判断系数矩阵是否严格对角占优,雅可比迭代法是否收敛
for(i=0;i<n;i++)
{
sum=0;
max=fabs(a[i][i]);
for(j=0;j<n;j++)
{
if(i!=j)
{
sum=sum+fabs(a[i][j]);
}
}
if(sum>=max)
{
printf("雅可比迭代法发散!\n");
}
}
//计算f=(D逆)*b
for(i=0;i<n;i++)
{
f[i]=Invdiag[i][i]*b[i];
}
for(k=0;k<MaxSize;k++)
{
//计算add1=B*x
for(i=0;i<n;i++)
{
sum=0;
for(j=0;j<n;j++)
{
sum=sum+Mat_Jacobi[i][j]*X[j];
}
add1[i]=sum;
}
//保存前一次的X解向量到temp
for(i=0;i<n;i++)
{
temp[i]=X[i];
}
//计算新的X=add1+f
for(i=0;i<n;i++)
{
X[i]=add1[i]+f[i];
}
//打印第k+1次迭代的解向量
if(G