Jacobi迭代法求线性方程组

#include<iostream>
using namespace std;


int num;//矩阵维度&未知数个数
double **A;//初始矩阵&迭代矩阵
double *b;//右端参数向量
double *Y;//初始向量
double e;//准许误差
int M;//最大迭代次数M


void print(int n, double**A,char *name) {
cout << "接下来打印矩阵" << name<<endl;
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
cout << A[i][j] << "   ";
}
cout << endl;
}
cout << "打印完毕" << endl;
}


void print(int n, double*s, char *name) {
cout << "接下来打印向量" << name<<endl;
for (int i = 0; i < n; i++) {
cout << s[i] << "   ";
}
cout << endl;
cout << "打印完毕" << endl;
}


void Init() {


printf("请输入未知数的个数：");
cin >> num;


A = new double*[num];
for (int i = 0; i < num; ++i)
A[i] = new double[num];


b = new double[num];
printf("请按行输入系数矩阵A：\n");
for (int i = 0; i < num; ++i)
for (int j = 0; j < num; ++j)
cin>>A[i][j];
print(num, A, "A");


printf("请输入n维结果向量：\n");
for (int i = 0; i < num; ++i)
cin>>b[i];
print(num, b, "b");


printf("请输入容许误差e：\n");
cin>>e;


printf("请输入最大容许迭代次数M：\n");
cin >>M;
}


//计算范数||X-Y||的值
double getNorm(double *x, double *y){
double sum = 0;
for (int i = 0; i < num; i++)
sum += abs(x[i] - y[i]);
return sum;
}


void Jacobi() {
//计算迭代矩阵B存在A中
for (int i = 0; i < num; i++) {
if (abs(A[i][i]) < e) {
cout << "求解失败" << endl;
exit(0);
}
double T = A[i][i];
for (int j = 0; j < num; j++)
A[i][j] = -(A[i][j]/T);
A[i][i] = 0;
b[i] = b[i] / T;
}


//初始化初始向量
Y = new double[num];
for (int i = 0; i < num; i++)
Y[i] = 0.0;


double *X = new double[num];
int k = 1;


//迭代
while (true) {
if (k > M) {
printf("迭代不收敛,已退出。");
return;
}
for (int i = 0; i < num; i++) {
double sum = 0;
for (int j = 0; j < num; j++)
sum += A[i][j] * Y[j];
X[i] = b[i] + sum;
}
if (getNorm(X, Y) < e) {
printf("迭代解得:");
for (int i = 0; i < num; ++i)
printf("x%d=%f ", i + 1, X[i]);
printf("\n");
return;
}
else {
k++;
for (int i = 0; i < num; i++)
Y[i] = X[i];
}
}
}


int main() {
Init();
Jacobi();
system("pause");
}