JACOBI迭代法 python

当我们计算线性方程组的解时，

可以把方程组变换为上面的下面这张图所示，即左右两边都有x，改图是写成矩阵的形式，如果我们写成方程组的形式，变为

我们手写计算通常都是待定系数法，如计算下面：

a+b=10;

a+2b=16;

待定系数即可算出a=4，b=6；但是计算机计算这些数量巨大的方程式使用迭代法，尤其当矩阵是大型稀疏矩阵（矩阵中有大部分元素都为0）时。

Jacobi迭代法也是依据上图方程式变形得到，

即：

这里的每一个x的下标i表示第i个x值，上标k+1表示迭代到了k+1次。

迭代求解的思想：

简单迭代法的求根过程分成两步，第一步先提供根的某个猜测值，即所谓迭代初值，然后将迭代初值逐步加工成满足精度要求的根。迭代法的设计思想是：f (x) = 0等价变换成 然后由迭代公式 逐步球的满足精度的解。实际迭代中不同迭代函数的求解可能影响求的精确解的运算量，甚至可能因为函数发散而无法求解。解题时可通过对导函数的判断而判断函数是否发散，而编写代码时可以通过判断循环次数——即循环过多次而不能从循环中出来时就判断为死循环，无法求得正解。

//c++代码
/*雅可比算法的代码实现*/
#include<iostream>
#include<iomanip>
#include<string>
#include<vector>
using namespace std;      

//函数求数组中的最大值
double MaxOfList(vector<double>x){
    double max=x[0];
    int n=x.size();
    for(int i=0;i<n;i++)
        if(x[i]>max) max=x[i];
    return max;
}

//雅可比迭代公式
void Jacobi(vector<vector<double> > A,vector<double> B,int n){
    vector<double> X(n,0);
    vector<double> Y(n,0);
    vector<double> D(n,0);
    int k=0; //记录循环次数
    do{ 
        X=Y;//注意到这一行代码与高斯赛德尔迭代法的区别
        for(int i=0;i<n;i++){
            double tem=0;
            for(int j=0;j<n;j++){
                if(i!=j) tem += A[i][j]*X[j];
            }
            Y[i]=(B[i]-tem)/A[i][i];
            cout<<left<<setw(8)<<Y[i]<<" ";
        }
        cout<<endl;
        k++;
        if(k>100){
            cout<<"迭代失败！（可能是函数不收敛）"<<endl;
            return ;
        }

        for(int a=0;a<n;a++){
            D[a]=X[a]-Y[a];
        }
    }while( MaxOfList(D)>0.00001 || MaxOfList(D)<-0.00001);

    return ;
}

int main(){

    int n;
    cout<<"请输入方程组未知数的个数n：";
    cin>>n;
    cout<<endl;

    vector<vector<double> >A(n,vector<double>(n,0));
    vector<double>B(n,0);

    cout<<"请输入方程组的系数矩阵："<<endl;
    for(int i=0;i<n;i++){
        for(int j=0;j<n;j++){
            cin>>A[i][j];
        }
    }
    cout<<endl;

    cout<<"请输入方程组的值向量："<<endl;
    for(int k=0;k<n;k++){
        cin>>B[k];
    }
    cout<<endl;

    cout<<"您输入的方程组为："<<endl;
    for(int a=0;a<n;a++){
        for(int b=0;b<n;b++){
            cout<<A[a][b]<<" ";
        }
        cout<<"    "<<B[a]<<endl;
    }
    cout<<endl;
    cout<<"由雅可比迭代公式求的方程组的解为："<<endl;
    Jacobi(A,B,n);
    return 0;
}