解线性方程组（二）——Jacobi迭代法求解（C++）

news2024/11/27 7:37:41

迭代法

相比于直接法求解，迭代法使用多次迭代来逐渐逼近解，其精度比不上直接法，但是其速度会比直接法快很多，计算精度可控，特别适用于求解系数矩阵为大型稀疏矩阵的方程组。

Jacobi迭代法

假设有方程组如下：
$\begin{cases} a_{11}x_1+a_{12}x_2+\cdots+a_{1n}x_n=b_1\\ a_{21}x_1+a_{22}x_2+\cdots+a_{2n}x_n=b_2\\ \cdots \qquad \qquad\cdots \qquad \qquad \cdots \\ a_{n1}x_1+a_{n2}x_2+\cdots+a_{nn}x_n=b_n\\ \end{cases}$
将其转换为矩阵形式
$A\vec{x}=\vec{b}$
$\begin{bmatrix} {a_{11}}&{a_{12}}&{\cdots}&{a_{1n}}\\ {a_{21}}&{a_{22}}&{\cdots}&{a_{2n}}\\ {\vdots}&{\vdots}&{\ddots}&{\vdots}\\ {a_{m1}}&{a_{m2}}&{\cdots}&{a_{mn}}\\ \end{bmatrix} \begin{bmatrix} {x_{1}}\\ {x_{2}}\\ {\vdots}\\ {x_{n}}\\ \end{bmatrix}= \begin{bmatrix} {b_{1}}\\ {b_{2}}\\ {\vdots}\\ {b_n} \end{bmatrix}$
对于是否可以使用Jacobi迭代法，需要满足以下条件之一：

A为行对角优阵，即 $|a_{ii}|>\sum_{j \neq i}|a_{ij}|(i=1,2,\cdots,n)$
A为行列角优阵，即 $|a_{jj}|>\sum_{j \neq i}|a_{ij}|(j=1,2,\cdots,n)$
A的元素满足 $\sum_{i \neq j}\frac{|a_{ij}|}{|aii|}<1(j,1,2,\cdots,n)$
若矩阵A满足上述条件之一，则可以使用Jacobi迭代法求解方程组。
首先将上述的方程组转为如下形式：
$\begin{cases} x_1=\frac{1}{a_{11}}(-a_{12}x_2-\cdots -a_{1n}x_n+b_1)\\ x_2=\frac{1}{a_{22}}(-a_{21}x_1-\cdots -a_{2n}x_n+b_2)\\ \cdots \qquad \qquad\cdots \qquad \qquad \cdots \\ x_n=\frac{1}{a_{nn}}(-a_{n1}x_1-\cdots -a_{nn-1}x_{n-1}+b_n)\\ \end{cases}$
写成矩阵形式可以得到Jacobi迭代式：
$(D+L+u)\vec{x}=\vec{b}\\ D\vec{x}=-(L+U)\vec{x}+\vec{b}\\ \vec{x}^{(k+1)}=-D^{-1}(L+U)\vec{x}^{(k)}+D^{-1}\vec{b}$
其中 $D$ 为对角矩阵， $L$ 为下三角矩阵- $D$ ， $U$ 为上三角矩阵- $U$ ， $D + L + U$ 为矩阵A。

代码实现

由于这个过程涉及大量的矩阵操作，整个算法分为两个源文件：Matrix.cpp实现矩阵操作，main.cpp实现Jacobi迭代法。
首先是Matrix.cpp的代码，其中矩阵求逆的原理参考：

#include <Matrix.h>
#include <iostream>
#include <cmath>
//矩阵与向量相乘，输入矩阵A，向量b，运算结果result和维数n
void matrix_multiply_vector(double **A,double *b,double * result,int n)
{
    for(int i=0;i<n;i++)
    {
        result[i]=0.0;
        for(int j=0;j<n;j++)
        {
            result[i]+=A[i][j]*b[j];
        }
    }
}
//矩阵乘法
void matrix_multiply_matrix(double **A,double **B,double **result,int n)
{
    for(int i=0;i<n;i++)
    {
        for(int j=0;j<n;j++)
        {
            result[i][j]=0.0;
            for(int k=0;k<n;k++)
            {
                result[i][j]+=A[i][k]*B[k][j];
            }
        }
    }
}
//矩阵加减法
void matrix_add_matrix(double **A,double **B,double **result,int n,bool isAdd)
{
    for(int i=0;i<n;i++)
    {
        for(int j=0;j<n;j++)
        {
            if(isAdd)
            {
                result[i][j]=A[i][j]+B[i][j];
            }
            else{
                result[i][j]=A[i][j]-B[i][j];
            }
            
        }
    }
}
//向量的加减法
void vactor_add_vector(double *A,double *B,double *result,int n,bool isAdd)
{
    for(int i=0;i<n;i++)
    {
        if(isAdd)
        {
            result[i]=A[i]+B[i];
        }
        else{
            result[i]=A[i]-B[i];
        }
    }
}
//判断向量误差范围，只要符合精度即可
bool vector_equal(double *A,double *B,int n,double error)
{
    for(int i=0;i<n;i++)
    {
        if(fabs(A[i]-B[i])>error)
        {
            return false;
        }
    }
    return true;
}
//向量赋值
void vector_copy(double *A,double *B,int n)
{
    for(int i=0;i<n;i++)
    {
        B[i]=A[i];
    }
}
//矩阵初始化
void matrix_init(double **A,int n)
{
    for(int i=0;i<n;i++)
    {
        A[i]=new double [n];
        for(int j=0;j<n;j++)
        {
            A[i][j]=0.0;
        }
    }
}
//判断矩阵A是否有收敛性
bool astringency(double **A,int n)
{
    double abs_row_sum=0.0;
    double abs_col_sum=0.0;
    double the_third_condition=0.0;
    bool RowOptimalMatrix=true;
    bool ColOptimalMatrix=true;
    for(int i=0;i<n;i++)//判断是不是行对角优阵
    {
        abs_row_sum=0.0;
        for(int j=0;j<n;j++)
        {
            if(i!=j)
            {
                abs_row_sum+=fabs(A[i][j]);
            }
        }
        if(abs_row_sum>A[i][i])//证明不是行对角优阵
        {
            RowOptimalMatrix=false;
            break;
        }
    }
    for(int j=0;j<n;j++)//判断是不是列对角优阵
    {
        abs_col_sum=0.0;
        for(int i=0;i<n;i++)
        {
            if(i!=j)
            {
                abs_col_sum+=fabs(A[i][j]);
            }
        }
        if(abs_col_sum>A[j][j])
        {
            ColOptimalMatrix=false;
            break;
        }
    }
    return ColOptimalMatrix or RowOptimalMatrix;
}
//矩阵交换某两行
void matrix_swap_row(double **A,int i,int j,int n)
{
    double temp;
    for(int k=0;k<n;k++)
    {
        temp=A[i][k];
        A[i][k]=A[j][k];
        A[j][k]=temp;
    }
}
//矩阵第i行=矩阵第i行-矩阵第j行*a
void matrix_minus_inner(double **A,double a,int i,int j,int n)
{
    for(int k=0;k<n;k++)
    {
        A[i][k]-=a*A[j][k];
    }
}
//矩阵求逆
void matrix_inverse(double **A,double **A_inverse,int n)
{
    double **A_E=new double*[2*n];
    //构建增广矩阵
    for(int i=0;i<n;i++)
    {
        A_E[i]=new double [n*2];
        for(int j=0;j<n*2;j++)
        {
            if(j<n)
            {
                A_E[i][j]=A[i][j];
            }
            else if((j-n)==i){
                A_E[i][j]=1;
            }
            else{
                A_E[i][j]=0;
            }
        }
    }
    //首先将矩阵化为上三角矩阵
    for(int i=0;i<n;i++)
    {
        if(A_E[i][i]==0)
        {
            for(int k=i+1;k<n;k++)
            {
                if(A_E[k][i]!=0)
                {
                    matrix_swap_row(A_E,i,k,n*2);
                    break;
                }
            }
        }
        for(int j=i+1;j<n;j++)
        {
            matrix_minus_inner(A_E,A_E[j][i]/A_E[i][i],j,i,2*n);
        }
    }
    //判断矩阵是否可逆
    for(int i=0;i<n;i++)
    {
        if(A_E[i][i]==0)
        {
            std::cout<<"矩阵不可逆"<<std::endl;
            exit(0);
        }
    }
    //将上三角转换为对角矩阵
    for(int j=1;j<n;j++)
    {
        for(int i=0;i<j;i++)
        {
            matrix_minus_inner(A_E,A_E[i][j]/A_E[j][j],i,j,2*n);
        }
    }
    for(int i=0;i<n;i++)
    {
        for(int j=n;j<2*n;j++)
        {
            A_inverse[i][j-n]=A_E[i][j]/A_E[i][i];
        }
    }
}

main.cpp文件内容如下：

//Jacobi迭代法求解线性方程组
/*
5x1+2x2-2x3=1
x1+4x2+x3=2
x1-2x2+4x3=-1
*/
#include<iostream>
#include<cmath>
#include<Matrix.h>//自定义头文件
using namespace std;
int main()
{
    int n;
    cout<<"Enter the matrix dimension A: ";
    cin>>n;//输入数组维度
    double **A=new double *[n];
    cout<<"Enter the coefficient matrix:"<<endl;
    for(int i=0;i<n;i++)
    {
        A[i]=new double[n];
        for(int j=0;j<n;j++)
        {
            cin>>A[i][j];//每次输入一个数字都用空格隔开,输入样例
            //1 2 3\enter
            //4 5 6\enter
            //7 8 9\enter
        }
    }
    double *b=new double[n];
    cout<<"Input vectors b: ";
    for(int i=0;i<n;i++)
    {
        cin>>b[i];//输入方程组右边的向量，1 2 3\enter
    }
    bool isAstringency=astringency(A,n);//判断系数矩阵A是否具有收敛性
    if(isAstringency)
    {
        cout<<"矩阵A符合收敛性"<<endl;
    }
    else{
        exit(0);
        cout<<"矩阵A不符合收敛性"<<endl;
    }
    double *x=new double[n];//解向量X
    double *x_last=new double[n];//上一次的x
    for(int i=0;i<n;i++)
    {
        x[i]=0.0;//初始化x
    }
    double **A_L_U=new double*[n];//L+U
    double **A_D_inverse=new double*[n];//D的逆
    for(int i=0;i<n;i++)
    {
        A_D_inverse[i]=new double [n];
        A_L_U[i]=new double [n];
        for(int j=0;j<n;j++)
        {
            if(i==j)
            {
                A_L_U[i][j]=0.0;
                A_D_inverse[i][j]=1.0/A[i][j];//对角矩阵的逆为其倒数
            }
            else{
                A_L_U[i][j]=A[i][j];
                A_D_inverse[i][j]=0.0;
            }
        }
    }
    double **B=new double *[n];//公式前半段的矩阵
    matrix_init(B,n);
    matrix_multiply_matrix(A_D_inverse,A_L_U,B,n);//求D^(-1)(L+U)
    double *f=new double[n];
    matrix_multiply_vector(A_D_inverse,b,f,n);//求取D^-1 * b
    double *temp1=new double[n];
    do{
        vector_copy(x,x_last,n);
        matrix_multiply_vector(B,x_last,temp1,n);//计算公式前半段
        vactor_add_vector(f,temp1,x,n,false);
    }while(vector_equal(x,x_last,n,1e-6)==false);//判断向量在误差范围内相等
    cout<<"运行结果为："<<endl;
    for(int i=0;i<n;i++)
    {
        cout<<x[i]<<" ";
    }
    system("pause");
    return 0;
}

结果分析

代码运行结果如下：
在这里插入图片描述

当下一次的迭代结果与上一次的迭代结果的最大相差值小于1e-6时，认为迭代已经收敛，输出结果即可（当然也可以换成其它结束迭代方法，如：判断两个向量之差的二范数）。
与直接使用克拉默法则计算准确的解以及matlab计算结果比较，不难发现其 $x_1$ 和 $x_3$ 均不为0，只是是一个在我们设定的误差范围内接近0的数，符合迭代法的解的性质，只能在设定的误差范围内得到一个近似的解。