当前位置： 首页 > news >正文

计算方法实验四 解线性方程组的间接方法

news 来源：原创 2025/5/3 6:02:40

【实验性质】

综合性实验。

【实验目的】

掌握迭代法求解线性方程组。 

【实验内容】

应用雅可比迭代法和Gauss-Sediel迭代法求解下方程组：

【理论基础】

线性方程组的数值解法分直接算法和迭代算法。迭代法将方程组的求解转化为构造一个向量序列，如果该向量序列存在极限，其极限就是方程组的解。

迭代法程序简单，但有时工作量较大，在有限步内，得不到精确解，适宜系数矩阵阶数较高的问题。

构造关于解向量的迭代序列，的常见方法有Jacobi迭代和Gauss-Seidel加速迭代。

设：

统一的迭代公式为：

−1

对于雅可比迭代，矩阵；对于高斯－赛德尔迭代，G = −（𝐷＋L） 𝑈，f =(𝐷＋L)−1𝑏。

实际的计算与编程应尽量避免矩阵求逆、矩阵相乘等运算，而使用如下公式：

【实验过程】

取三个不同初值，分别用雅可比法和高斯－塞得尔法求解，用表格记载求解过程，并分析迭代初值、迭代公式对迭代的影响。

附表：

 

雅可比法

高斯－塞得尔法

 

初值

近似根

迭代次数

近似根

迭代次数

1

 

 

 

 

 

 

 

2

 

 

 

 

 

 

 

3

 

 

 

 

 

 

 

代码：

主函数：

//实验四

#include <iostream>

#include <windows.h>

#include "colvector.h"

#include "matrix.h"

#include <windows.h>

#include "linearequations.h"

using namespace std;

int main()

{

    double a1[]={8,7,0,0,6,12,5,0,0,4,9,3,0,0,1,2};

    double b1[]={0,-2,8,6};

    double x1[]={0,0,0,0};

    double e=0.0001;

    int flag=1,max=100;

    Matrix A(4,4);

    ColVector b(4),x(4),x0(4);

    A.initMatrix(a1);

    b.initMatrix(b1);

    x0.initMatrix(x1);

    LinearEquations obj;

    obj.JacobiIterativeMethod(A,b,e,x0,max,flag,x);

    cout<<x<<endl;

    return 0;

}

头文件：

#ifndef LINEAREQUATIONS_H
#define LINEAREQUATIONS_H
#include "matrix.h"
#include "colvector.h"
#include  <math.h>
class LinearEquations
{
public:
    void GaussElimin(Matrix A,ColVector b,int &flag,ColVector &x);
    void colPivotGaussElimin(Matrix A,ColVector b,int &flag,ColVector &x);
    void LUDecomposition(Matrix &mat,int &flag,Matrix &L,Matrix &U,int type=1);
    void TrigDecompositionMethod(Matrix &A,ColVector &b,int &flag,ColVector &x,int type=1);
    void Catchup(ColVector&a,ColVector &b,ColVector &c,ColVector &f,int &flag,ColVector &x,int type=1);
    LinearEquations();
    //雅可比迭代
    void JacobiIterativeMethod(Matrix &A,ColVector &b,double e,ColVector &x0,int MAX, int &flag,ColVector &x);
    //高斯-赛德尔迭代
    void Gauss_SeidelIterativeMethod(Matrix &A,ColVector &b,double e,ColVector &x0,int MAX, int &flag,ColVector &x);
 

};
 

#endif // LINEAREQUATIONS_H
 

 

源文件：
 

#include "linearequations.h"
 

LinearEquations::LinearEquations()
 

{
 

 

}
 

 void LinearEquations::GaussElimin(Matrix A,ColVector b,int &flag,ColVector &x)
 

 {
 

      flag=1;
 

      int n=A.getRowSize();
 

      if(n!=A.getColSize() || n!=b.getRowSize())
 

      {
 

          flag =0;
 

          return;
 

      }
 

      cout<<"初始曾广矩阵"<<Matrix::merge(A,b)<<endl;
 

     for(int k=1;k<=n-1;k++)
 

     {
 

         for(int i=k+1;i<=n;i++)
 

         {
 

             if(A(k,k)==0)
 

             {
 

                 flag =0;
 

                 return;
 

             }
 

             double m=A(i,k)/A(k,k);
 

             for(int j=1;j<=n;j++)
 

             {
 

                 A(i,j)=A(i,j)-m*A(k,j);
 

             }
 

             b[i]=b[i]-m*b[k];
 

         }
 

         cout<<"第"<<k<<"次消元"<<Matrix::merge(A,b)<<endl;
 

     }
 

     cout<<"解 x"<<endl;
 

     x=ColVector(n);
 

     if(A(n,n)==0)
 

     {
 

         flag=0;
 

         return;
 

     }
 

     x[n]=b[n]/A(n,n);
 

     cout<<x<<endl;
 

     for(int i=n-1;i>=1;i--)
 

     {
 

         double sum=0;
 

         for(int j=i+1;j<=n;j++)
 

         {
 

             sum+=A(i,j)*x[j];
 

         }
 

         x[i]=(b[i]-sum)/A(i,i);
 

         cout<<"解 x"<<endl;
 

         cout<<x<<endl;
 

     }
 

 }
 

void LinearEquations::colPivotGaussElimin(Matrix A,ColVector b,int &flag,ColVector &x)
 

{
 

    flag=1;
 

    int n=A.getRowSize();
 

    if(n!=A.getColSize() || n!=b.getRowSize())
 

    {
 

        flag =0;
 

        return;
 

    }
 

    cout<<"初始曾广矩阵"<<Matrix::merge(A,b)<<endl;
 

   for(int k=1;k<=n-1;k++)
 

   {
 

       double max=fabs(A(k,k));
 

       int p=k;
 

       for(int w=k+1;w<=n;w++)
 

       {
 

           if(max< fabs(A(w,k))){
 

               max=fabs(A(w,k));
 

               p=w;
 

           }
 

       }
 

       if(max==0){
 

           flag=0;
 

           return;
 

       }
 

       if(p!=k){
 

           double temp=0;
 

           for(int s=1;s<=n;s++)
 

           {
 

               temp =A(k,s);
 

               A(k,s)=A(p,s);
 

               A(p,s)=temp;
 

           }
 

           temp=b[k];
 

           b[k]=b[p];
 

           b[p]=temp;
 

       }
 

       for(int i=k+1;i<=n;i++)
 

       {
 

           if(A(k,k)==0)
 

           {
 

               flag =0;
 

               return;
 

           }
 

           double m=A(i,k)/A(k,k);
 

           for(int j=1;j<=n;j++)
 

           {
 

               A(i,j)=A(i,j)-m*A(k,j);
 

           }
 

           b[i]=b[i]-m*b[k];
 

       }
 

       cout<<"第"<<k<<"次消元"<<Matrix::merge(A,b)<<endl;
 

   }
 

  cout<<"解 x"<<endl;
 

   x=ColVector(n);
 

   if(A(n,n)==0)
 

   {
 

       flag=0;
 

       return;
 

   }
 

   x[n]=b[n]/A(n,n);
 

   cout<<x<<endl;
 

   for(int i=n-1;i>=1;i--)
 

   {
 

       double sum=0;
 

       for(int j=i+1;j<=n;j++)
 

       {
 

           sum+=A(i,j)*x[j];
 

       }
 

       x[i]=(b[i]-sum)/A(i,i);
 

       cout<<"解 x"<<endl;
 

       cout<<x<<endl;
 

   }
 

}
 

void LinearEquations::LUDecomposition(Matrix &mat,int &flag,Matrix &L,Matrix &U,int type)
 

{
 

    flag=0;
 

    int n=mat.getRowSize();
 

    if(n!=mat.getColSize()||type<1||type>2){
 

        flag=-1;
 

        return;
 

    }
 

    L=Matrix(n,n);
 

    U=Matrix(n,n);
 

    switch(type){
 

    case 1:
 

        for(int i=1;i<+n;i++){
 

        U(1,i)=mat(1,i);
 

        }
 

        for(int i=1;i<=n;i++){
 

        L(i,1)=mat(i,1)/U(1,1);
 

        }
 

        for(int k=2;k<=n;k++)
 

        {
 

            for(int j=k;j<=n;j++){
 

                double sum=0;
 

                for(int r=1;r<=k-1;r++){
 

                sum+=L(k,r)*U(r,j);
 

                }
 

                   U(k,j)=mat(k,j)-sum;
 

            }
 

            for(int i=k;i<=n;i++){
 

            double sum=0;
 

            for(int r=1;r<=k-1;r++){
 

            sum+=L(i,r)*U(r,k);
 

            }
 

            L(i,k)=(mat(i,k)-sum)/U(k,k);
 

            }
 

        }
 

        break;
 

    case 2:
 

        for(int k=1;k<=n;k++)
 

        {
 

            for(int j=k;j<=n;j++)
 

            {
 

                double sum=0;
 

                if(k>1){
 

                    for(int r=1;r<=k-1;r++){
 

                        sum+=L(j,r)*U(r,k);
 

                    }
 

                }
 

                L(j,k)=mat(j,k)-sum;
 

            }
 

            for(int j=k;j<=n;j++)
 

            {
 

                if(L(k,k)==0){
 

                    flag=0;
 

                    return ;
 

                }
 

                double sum=0;
 

                if(k>1)
 

                {
 

                    for(int r=1;r<=k-1;r++){
 

                        sum+=L(k,r)+U(r,j);
 

                    }
 

                }
 

                U(k,j)=(mat(k,j)-sum)/L(k,k);
 

            }
 

        }
 

     break;
 

 }
 

}
 

void LinearEquations::TrigDecompositionMethod(Matrix &A,ColVector &b,int &flag,ColVector &x,int type)
 

{
 

    flag =1;
 

    int n=A.getRowSize();
 

    if(n!=A.getColSize()||n!=b.getRowSize()||type<1||type>2)
 

    {
 

        flag=0;
 

        return;
 

    }
 

    Matrix L,U;
 

    ColVector y(n);
 

    x=ColVector(n);
 

    switch(type){
 

    case 1:
 

        LUDecomposition(A,flag,L,U);
 

        if(flag== -1)
 

        {
 

            return;
 

        }
 

        y[1]=b[1];
 

        for(int i=2;i<=n;i++)
 

        {
 

            double sum=0;
 

            for(int k=1;k<=i-1;k++)
 

            {
 

                sum+=L(i,k)*y[k];
 

            }
 

            y[i]=b[i]-sum;
 

        }
 

        cout<<"y="<<y<<endl;
 

        if(U(n,n)==0)
 

        {
 

            flag=0;
 

            return;
 

        }
 

        x[n]=y[n]/U(n,n);
 

        for(int i=n-1;i>=1;i--)
 

        {
 

            if(U(i,i)==0){
 

            flag=0;
 

            return;
 

    }
 

            double sum=0;
 

            for(int k=i+1;k<=n;k++)
 

        {
 

            sum+=U(i,k)*x[k];
 

       }
 

           x[i]=(y[i]-sum)/U(i,i);
 

    }
 

break;
 

    case 2:
 

        LUDecomposition(A,flag,L,U);
 

        if(flag<1)
 

        {
 

            return;
 

        }
 

        if(L(1,1)==0)
 

        {
 

            flag=0;
 

            return;
 

        }
 

        y[1]=b[1]/L(1,1);
 

        for(int i=2;i<=n;i++)
 

        {
 

            if(L(i,i)<=0){
 

                flag=0;
 

                return;
 

            }
 

            double sum=0;
 

            for(int k=1;k<=i-1;k++)
 

            {
 

                sum+=L(i,k)*y[k];
 

            }
 

            y[i]=(b[i]-sum)/L(i,i);
 

        }
 

        x[n]=y[n];
 

        for(int i=n-1;i>=1;i--)
 

        {
 

            double sum=0;
 

            for(int k=i+1;k<=n;k++)
 

        {
 

            sum+=U(i,k)*x[k];
 

       }
 

           x[i]=y[i]-sum;
 

    }
 

 break;
 

}
 

}
 

void LinearEquations::Catchup(ColVector&a,ColVector &b,ColVector &c,ColVector &f,int &flag,ColVector &x,int type)
 

{
 

    flag=1;
 

    int n=b.getRowSize();
 

    if(n!=a.getRowSize()||n!=b.getRowSize()||n!=f.getRowSize())
 

    {
 

        flag=0;
 

        return;
 

    }
 

    if(fabs(b[1])<=fabs(c[1])||fabs(b[n])<=fabs(a[n])||b[1]*c[1]==0||a[n]*b[n]==0)
 

    {
 

        flag=0;
 

        return;
 

    }
 

    for(int i=2;i<=n-1;i++)
 

        if(fabs(b[i])<fabs(a[i])+fabs(c[i])||a[i]*c[i]==0)
 

        {
 

            flag=0;
 

            return;
 

        }
 

    ColVector d(n),u(n),y(n),l(n);
 

    x=ColVector (n);
 

    switch(type){
 

    case 1:
 

        d=c;
 

        u[1]=b[1];
 

        for(int i=2;i<=n;i++)
 

        {
 

            l[i]=a[i]/u[i-1];
 

            u[i]=b[i]-l[i]*c[i-1];
 

        }
 

        y[1]=f[1];
 

        for(int i=2;i<=n;i++)
 

        {
 

            y[i]=f[i]-l[i]*y[i-1];
 

        }
 

         x[n]=y[n]/u[n];
 

         for(int i=n-1;i>=1;i--)
 

         {
 

           x[i]=(y[i]-c[i]*x[i+1])/u[i];
 

         }
 

        break;
 

    case 2:
 

        l[1]=b[1];
 

        for(int i=1;i<=n-1;i++)
 

        {
 

            u[i]=c[i]/l[i];
 

            l[i+1]=b[i+1]-a[i+1]*u[i];
 

        }
 

        y[1]=f[1]/l[1];
 

        for(int i=2;i<=n;i++)
 

        {
 

            y[i]=(f[i]-a[i]*y[i-1])/l[i];
 

        }
 

        x[n]=y[n];
 

        for(int i=n-1;i>=1;i--)
 

        {
 

          x[i]=y[i]-u[i]*x[i+1];
 

        }
 

        break;
 

    }
 

}
 

 

void LinearEquations::JacobiIterativeMethod(Matrix &A,ColVector &b,
 

                                            double e,ColVector &x0,int MAX, int &flag,ColVector &x){
 

    flag =1;
 

    int n=A.getRowSize();
 

    if(n!=b.getRowSize()||n!=A.getColSize()||n!=x0.getRowSize())
 

    {
 

        flag=0;
 

        return;
 

    }
 

 

    for(int i=1;i<=n;i++){
 

        if(A(i,i)==0){
 

            flag=0;
 

            return ;
 

        }
 

    }
 

    double eps=0;
 

    int k=0;
 

    int sum=0;
 

    int i,j;
 

    ColVector x1(n);
 

    ColVector dist(n);
 

    while(k<MAX){
 

        for(i=1;i<=n;i++){
 

            double sum=0;
 

            for(j=1;j<=n;j++){
 

                if(j==1){
 

                    continue;
 

                }
 

                sum+=A(i,j)*x0[j];
 

            }
 

            x1[i]=(b[i]-sum)/A(i,i);
 

        }
 

        k++;
 

        dist=x1-x0;
 

        eps=dist.norm(3);
 

        cout<<endl;
 

        cout<<"k="<<k<<"\nx1="<<x1<<"x0="<<x0<<endl;
 

        cout<<"dist="<<dist<<endl;
 

        cout<<"eps="<<eps<<endl;
 

        if(eps<=e){
 

            break;
 

        }
 

        else{
 

            x0=x1;
 

        }
 

    }
 

    if(k==MAX){
 

        flag=0;
 

        return;
 

    }
 

    else{
 

        x=x1;
 

    }
 

}
 

 

void LinearEquations::Gauss_SeidelIterativeMethod(Matrix &A,ColVector &b,double e,
 

                                                  ColVector &x0,int MAX, int &flag,ColVector &x){
 

    flag =1;
 

    int n=A.getRowSize();
 

    if(n!=b.getRowSize()||n!=A.getColSize()||n!=x0.getRowSize())
 

    {
 

        flag=0;
 

        return;
 

    }
 

 

    for(int i=1;i<=n;i++){
 

        if(A(i,i)==0){
 

            flag=0;
 

            return ;
 

        }
 

    }
 

    double eps=0;
 

    int k=0;
 

    ColVector x1(n);
 

    ColVector dist(n);
 

    while(k<MAX){
 

        for(int i=1;i<=n;i++){
 

            double sum=0;
 

            for(int j=1;j<=n;j++){
 

                if(j<i){
 

                    sum+=A(i,j)*x1[j];
 

                }
 

                else if(j==i){
 

                    continue;
 

                }
 

                else{
 

                    sum+=A(i,j)*x0[j];
 

                }
 

            }
 

            x1[i]=(b[i]-sum)/A(i,i);
 

        }
 

        k++;
 

        dist=x1-x0;
 

        eps=dist.norm(3);
 

        cout<<"k="<<k<<"\nx1="<<x1<<"x0="<<x0<<endl;
 

        cout<<"dist="<<dist<<endl;
 

        cout<<"eps="<<eps<<endl;
 

        if(eps<=e){
 

            break;
 

        }
 

        else{
 

            x0=x1;
 

        }
 

    }
 

    if(k==MAX){
 

        flag=0;
 

        return;
 

    }
 

    else{
 

        x=x1;
 

    }
 

}
 

 
 

运行结果：
 

 

【实验心得】 
 

我在实验中使用了间接方法解线性方程组，获得了一些心得。首先，我发现使用间接方法可以简化计算过程，特别适用于方程组较大的情况。通过将方程组转化为矩阵形式，并利用矩阵的运算性质来求解，可以大大减少计算量。其次，我注意到间接方法的精度较高，能够得到比较准确的解。与直接方法相比，间接方法不容易出现舍入误差，因为矩阵运算过程中的舍入误差可以通过精确的结果进行修正。此外，我发现间接方法还可以应用于其他数学问题的求解中，例如线性最小二乘问题等。总体来说，间接方法是一种简单高效、精度较高的解线性方程组的方法，对于数学问题的求解也具有一定的推广价值。通过这次实验，我不仅学到了解线性方程组的间接方法，还深入了解了矩阵运算的原理和应用，对数学问题的求解能力也有了一定的提升。
 

 

得    分________
 

 
 

评阅日期________ 
 

教师签名________ 
 



                
                

                  

                    
                      
                  



                

                

                  
相关文章：

                    
                      
  • 
                        【Unity】使用XLua实现C#访问Lua文件
                      
    • 
                    
                      
  • 
                        JavaScript常规解密技术解析指南
                      
    • 
                    
                      
  • 
                        对称加密算法（AES、ChaCha20和SM4）Python实现——密码学基础(Python出现No module named “Crypto” 解决方案)
                      
    • 
                    
                      
  • 
                        Linux系统：进程程序替换以及相关exec接口
                      
    • 
                    
                      
  • 
                        大学生入学审核系统设计与实现【基于SpringBoot + Vue 前后端分离技术】
                      
    • 
                    
                      
  • 
                        记忆翻牌游戏：认知科学与状态机的交响曲
                      
    • 
                    
                      
  • 
                        关于在vscode终端不能执行npm
                      
    • 
                    
                      
  • 
                        React 组件prop添加类型
                      
    • 
                    
                      
  • 
                        【数据结构】String字符串的存储
                      
    • 
                    
                      
  • 
                        React memo
                      
    • 
                    
                      
  • 
                        C++八股--three day --设计模式之单例和工厂
                      
    • 
                    
                      
  • 
                        【数据结构】励志大厂版·初阶（复习+刷题）：栈与队列
                      
    • 
                    
                      
  • 
                        系统架构设计师：设计模式——结构型设计模式
                      
    • 
                    
                      
  • 
                        MCP智能体意图识别与工具路由：让AI自主决策调用链路
                      
    • 
                    
                      
  • 
                        Oracle无法正常OPEN(三）
                      
    • 
                    
                      
  • 
                        Webug4.0靶场通关笔记13- 第22关越权修改密码
                      
    • 
                    
                      
  • 
                        【Linux网络编程】http协议的状态码，常见请求方法以及cookie-session
                      
    • 
                    
                      
  • 
                        AE脚本 关键帧缓入缓出曲线调节工具 Flow v1.5.0 Win/Mac
                      
    • 
                    
                      
  • 
                        AI 生成内容的版权困境：法律、技术与伦理的三重挑战
                      
    • 
                    
                      
  • 
                        使用 Java 实现一个简单且高效的任务调度框架
                      
    • 
                    

  • 南京106亿元成交19宗涉宅地块：建邺区地块楼面单价重回4.5万元
    • 


  • 澎湃回声丨23岁小伙“被精神病8年”续：今日将被移出“重精”管理系统
    • 


  • 百年传承，再启新程，参天中国迎来2.0时代
    • 


  • 西班牙葡萄牙电力基本恢复
    • 


  • 对谈｜李钧鹏、周忆粟：安德鲁·阿伯特过程社会学的魅力
    • 


  • 老凤祥一季度净利减少两成，去年珠宝首饰营收下滑19%
    •