解线性方程组(二)——Jacobi迭代法求解(C++)

迭代法

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

Jacobi迭代法

假设有方程组如下:
{ a 11 x 1 + a 12 x 2 + ⋯ + a 1 n x n = b 1 a 21 x 1 + a 22 x 2 + ⋯ + a 2 n x n = b 2 ⋯ ⋯ ⋯ a n 1 x 1 + a n 2 x 2 + ⋯ + a n n x n = b n \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} a11x1+a12x2++a1nxn=b1a21x1+a22x2++a2nxn=b2an1x1+an2x2++annxn=bn
将其转换为矩阵形式
A x ⃗ = b ⃗ A\vec{x}=\vec{b} Ax =b
[ a 11 a 12 ⋯ a 1 n a 21 a 22 ⋯ a 2 n ⋮ ⋮ ⋱ ⋮ a m 1 a m 2 ⋯ a m n ] [ x 1 x 2 ⋮ x n ] = [ b 1 b 2 ⋮ b n ] \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} a11a21am1a12a22am2a1na2namn x1x2xn = b1b2bn
对于是否可以使用Jacobi迭代法,需要满足以下条件之一:

  1. A为行对角优阵,即 ∣ a i i ∣ > ∑ j ≠ i ∣ a i j ∣ ( i = 1 , 2 , ⋯ , n ) |a_{ii}|>\sum_{j \neq i}|a_{ij}|(i=1,2,\cdots,n) aii>j=iaij(i=1,2,,n)
  2. A为行列角优阵,即 ∣ a j j ∣ > ∑ j ≠ i ∣ a i j ∣ ( j = 1 , 2 , ⋯ , n ) |a_{jj}|>\sum_{j \neq i}|a_{ij}|(j=1,2,\cdots,n) ajj>j=iaij(j=1,2,,n)
  3. A的元素满足 ∑ i ≠ j ∣ a i j ∣ ∣ a i i ∣ < 1 ( j , 1 , 2 , ⋯ , n ) \sum_{i \neq j}\frac{|a_{ij}|}{|aii|}<1(j,1,2,\cdots,n) i=jaiiaij<1(j,1,2,,n)
    若矩阵A满足上述条件之一,则可以使用Jacobi迭代法求解方程组。
    首先将上述的方程组转为如下形式:
    { x 1 = 1 a 11 ( − a 12 x 2 − ⋯ − a 1 n x n + b 1 ) x 2 = 1 a 22 ( − a 21 x 1 − ⋯ − a 2 n x n + b 2 ) ⋯ ⋯ ⋯ x n = 1 a n n ( − a n 1 x 1 − ⋯ − a n n − 1 x n − 1 + b n ) \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} x1=a111(a12x2a1nxn+b1)x2=a221(a21x1a2nxn+b2)xn=ann1(an1x1ann1xn1+bn)
    写成矩阵形式可以得到Jacobi迭代式:
    ( D + L + u ) x ⃗ = b ⃗ D x ⃗ = − ( L + U ) x ⃗ + b ⃗ x ⃗ ( k + 1 ) = − D − 1 ( L + U ) x ⃗ ( k ) + D − 1 b ⃗ (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+u)x =b Dx =(L+U)x +b x (k+1)=D1(L+U)x (k)+D1b
    其中 D D D为对角矩阵, L L L为下三角矩阵- D D D U U U为上三角矩阵- U U U D + L + U D+L+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];//解向量Xdouble *x_last=new double[n];//上一次的xfor(int i=0;i<n;i++){x[i]=0.0;//初始化x}double **A_L_U=new double*[n];//L+Udouble **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 * bdouble *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_1 x1 x 3 x_3 x3均不为0,只是是一个在我们设定的误差范围内接近0的数,符合迭代法的解的性质,只能在设定的误差范围内得到一个近似的解。

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.mzph.cn/news/685918.shtml

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

QGIS004:【08图层工具箱】-导出到电子表格、提取图层范围

摘要&#xff1a;QGIS图层工具箱常用工具有导出到电子表格、提取图层范围等选项&#xff0c;本文介绍各选项的基本操作。 实验数据&#xff1a; 链接&#xff1a;https://pan.baidu.com/s/1ZK4_ShrQ5BsbyWfJ6fVW4A?pwdpiap 提取码&#xff1a;piap 一、导出到电子表格 工具…

OpenAl 视频生成模型 —— Sora技术报告解读

这里是陌小北&#xff0c;一个正在研究硅基生命的碳基生命。正在努力成为写代码的里面背诗最多的&#xff0c;背诗的里面最会写段子的&#xff0c;写段子的里面代码写得最好的…厨子。 写在前面 早上醒来&#xff0c;就看到OpenAl推出的视频模型Sora炸锅了&#xff0c;感觉所…

代码随想录 Leetcode452. 用最少数量的箭引爆气球

题目9&#xff1a; 代码&#xff08;首刷看解析 2024年2月17日&#xff09;&#xff1a; class Solution { private:const static bool cmp(vector<int>& a, vector<int>& b) {return a[0] < b[0];} public:int findMinArrowShots(vector<vector<…

Linux下解压tar.xz文件的命令

tar -c: 建立压缩档案-x&#xff1a;解压-t&#xff1a;查看内容-r&#xff1a;向压缩归档文件末尾追加文件-u&#xff1a;更新原压缩包中的文件 ------------------------------------------ 这五个是独立的命令&#xff0c;压缩解压都要用到其中一个&#xff0c;可以和别的…

如何用Qt实现一个无标题栏、半透明、置顶(悬浮)的窗口

在Qt框架中&#xff0c;要实现一个无标题栏、半透明、置顶&#xff08;悬浮&#xff09;的窗口&#xff0c;需要一些特定的设置和技巧。废话不多说&#xff0c;下面我将以DrawClient软件为例&#xff0c;介绍一下实现这种效果的四个要点。 要点一&#xff1a;移除标题栏&#…

定时器外部时钟

一、相较于内部时钟中断改动&#xff1a; 1.Timer.c RCC_APB2PeriphClockCmd(RCC_APB2Periph_GPIOA, ENABLE); //开启GPIOA的时钟/*GPIO初始化*/GPIO_InitTypeDef GPIO_InitStructure;GPIO_InitStructure.GPIO_Mode GPIO_Mode_IPU;GPIO_InitStructure.GPIO_Pin GPIO_Pin_…

BUGKU-WEB 头等舱

题目描述 题目截图如下&#xff1a; 进入场景看看&#xff1a; 解题思路 先看看源码再看看F12请求和响应 相关工具 略 解题步骤 查看源码&#xff0c;好家伙真的什么也没有 2. 看看F12请求和响应&#xff0c;找到了 得到Flag flag{a49c7aba1014c3673ec9982946d0545a…

fastposter v2.18.0 一分钟完成开发海报-云服务来袭

fastposter v2.18.0 一分钟完成开发海报-云服务来袭 fastposter 是一款快速开发海报的工具&#xff0c;已经服务众多电商、行业海报、分销系统、电商海报、电商主图等海报生成和制作场景。 什么是 fastposter &#x1f525;&#x1f525;&#x1f525;fastposter 是一款海报…

【数据结构】16 二叉树的定义,性质,存储结构(以及先序、后序、中序遍历)

二叉树 一个二叉树是一个有穷的结点集合。 它是由根节点和称为其左子树和右子树的两个不相交的二叉树组成的。 二叉树可具有以下5种形态。 性质 一个二叉树第i层的最大结点数为 2 i − 1 2^{i-1} 2i−1, i ≥ 1 i \geq 1 i≥1 每层最大结点可以对应完美二叉树&#xff08;…

CSDN如何获得更多勋章?

文章目录 前言一、如何找到自己的勋章&#xff1f;二、如何获得更多勋章&#xff1f;三、重点勋章、易得勋章介绍&推荐1.创作能手2.五一创作勋章3.创作纪念日IT一周年勋章4.新秀勋章5.话题达人6.128天创作纪念日&#xff08;IT博客专属&#xff09;7.GitHub绑定勋章8.其他 …

TIM(Timer)定时中断 P1

难点&#xff1a;定时器级联、主从模式 一、简介&#xff1a; 1.TIM&#xff08;Timer&#xff09;定时器 定时器可以对输入的时钟进行计数&#xff0c;并在计数值达到设定值时触发中断 补充&#xff1a; { 定时器本质上是一个计数器&#xff0c;可以工作在定时或计数模式&…

如何简单上手清华AutoGPT并搭建到本地环境

一、准备工作 安装Docker&#xff1a;确保你的本地机器上已经安装了Docker。如果还没有安装&#xff0c;请访问Docker官方网站并按照指引进行安装。--点击进入Docker官网 获取清华AutoGPT的Docker镜像&#xff1a;清华AutoGPT团队可能已经提供了一个Docker镜像&#xff0c;方便…

java8新特性——StreamAPI

说明&#xff1a; java8中有两大最为重要的改变。第一个是Lambda表达式&#xff1b;另外一个则是Stream API。 Stream API&#xff08;java.util.stream&#xff09;把真正的函数式编程风格引入java。这是目前为止对java类库最好的补充&#xff0c;因为Stream API可以极大提供j…

AI:130-基于深度学习的室内导航与定位

🚀点击这里跳转到本专栏,可查阅专栏顶置最新的指南宝典~ 🎉🎊🎉 你的技术旅程将在这里启航! 从基础到实践,深入学习。无论你是初学者还是经验丰富的老手,对于本专栏案例和项目实践都有参考学习意义。 ✨✨✨ 每一个案例都附带有在本地跑过的关键代码,详细讲解供…

红蓝对抗:网络安全领域的模拟实战演练

引言&#xff1a; 随着信息技术的快速发展&#xff0c;网络安全问题日益突出。为了应对这一挑战&#xff0c;企业和组织需要不断提升自身的安全防护能力。红蓝对抗作为一种模拟实战演练方法&#xff0c;在网络安全领域得到了广泛应用。本文将介绍红蓝对抗的概念、目的、过程和…

硬盘合并分区失败显示未格式化怎么办?看这里

在当今数字化时代&#xff0c;硬盘分区已成为计算机存储管理的重要手段。然而&#xff0c;在使用过程中&#xff0c;我们有时会遇到“合并分区失败显示未格式化”的问题&#xff0c;这不仅影响了我们的工作效率&#xff0c;还可能造成数据丢失的风险。下面将根据不同情况给予不…

【LeetCode: 107. 二叉树的层序遍历 II + BFS】

&#x1f680; 算法题 &#x1f680; &#x1f332; 算法刷题专栏 | 面试必备算法 | 面试高频算法 &#x1f340; &#x1f332; 越难的东西,越要努力坚持&#xff0c;因为它具有很高的价值&#xff0c;算法就是这样✨ &#x1f332; 作者简介&#xff1a;硕风和炜&#xff0c;…

node+vue3+mysql前后分离开发范式——实现视频文件上传并渲染

文章目录 ⭐前言⭐ 功能设计与实现💖 node上传文件写入file_map映射表💖 vue3前端上传文件回显⭐ 效果⭐结束⭐前言 大家好,我是yma16,本文分享关于 node+vue3+mysql前后分离开发范式——实现视频文件上传并渲染。 技术选型 前端:vite+vue3+antd 后端:node koa 数据库…

JVM-垃圾回收(标记算法,收集器)

申明&#xff1a;文章内容是本人学习极客时间课程所写&#xff0c;文字和图片基本来源于课程资料&#xff0c;在某些地方会插入一点自己的理解&#xff0c;未用于商业用途&#xff0c;侵删。 原资料地址&#xff1a;课程资料 垃圾回收的基本原理 1 什么是垃圾&#xff1f; 在…

阿里云ECS香港服务器性能强大_安全可靠香港免备案服务器

阿里云香港服务器中国香港数据中心网络线路类型BGP多线精品&#xff0c;中国电信CN2高速网络高质量、大规格BGP带宽&#xff0c;运营商精品公网直连中国内地&#xff0c;时延更低&#xff0c;优化海外回中国内地流量的公网线路&#xff0c;可以提高国际业务访问质量。阿里云服务…