用C语言实现雅可比迭代法
C语言实现jacobi迭代
用C++编程实现简单迭代法班级：09医软一班 制作人：王法金小组成员：王法金 09713038 源代码编写、测试及程序最终修改。
汪晓阳 09713037 程序设计讨论和测试 王腾飞 09713039 程序设计讨论和完善简单迭代法原理及公式：⎪⎪⎩⎪⎪⎨⎧=+++=+++=+++nn nn n n n n n n b x a x a x a b x a x a x a b x a x a x a 22112222212111212111设n 阶线性方程组()()()⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧+----=+----=+----=--n n n n n n nn n n n n n b x a x a x a a x b x a x a x a a x b x a x a x a a x 11,22112232312122211313212111111 ()n i a ii ,,2,10 =≠的系数矩阵A 非奇异，且 ，化为等价方程组一、Jacobi 迭代计算公式}{)(m x T m n m m m x x x x ),,,()()(2)(1)( =可得向量序列 ，其中 α=∞→)(lim m m x 如果，那么 α就是原方程组的解.这种求线性方程组的解的方法称为简单迭代法.或称为雅可比（Jacobi ）迭代法.()()()()()()()()()()()()()()()⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧+----=+----=+----=--+++nm n n n m n m n nn m n m n n m m m m n n m m m b x a x a x a a x b x a x a x a a x b x a x a x a a x 11,21112223121221211132121111233111 ()()()()()Tnx x x x 002010,,, =任给初始向量，由迭代公式二、Jacobi 迭代的矩阵形式⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=--00001,21323121n n n n a a a a a a L⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=--0000,122311312n n n n a a a a a a U ⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=nn a a a D2211若令()bD x U L D x 11--++=()),2,1,0( )(1)(11 =++=--+m b D x U L D x m m 则令,),(111b D g U L D B --=+=),2,1,0( 1)(1)1( =+=+m g x B x m m 则方程组Ax=b 化为等价方程组 于是迭代公式为：为简单迭代法的矩阵形式程序测试：测试例子1：《数值计算方法》 p201 例7.测试结果分析：经测试，和书上例题迭代6次结果相同！手动算了迭代前5次结果和测试结果相同，说明了本程序有了一定的正确性！迭代11次，是在精度p=0.001情况下所需要的迭代次数，迭代11次结果很接近X值，同时也说明了程序的正确性。
雅克比迭代法与高斯-赛德尔
实验三 线性方程组的迭代法班级：**计本 ** 班 姓名：** 座号： ** 时间：2010/6/2一、 实验目的（1） 熟悉VC++开发平台和开发语言。
（2） 掌握雅可比及高斯-塞德尔迭代法解方程组的迭代法，并能根据给定的精度要求计算（包括迭代过程）；比较两种方法的优劣。
（3） 培养编程和上机调试能力。
二、 实验设备一台PC 机，XP 操作系统，VC++软件三、 实验内容用雅可比迭代法和高斯-塞德尔迭代法解方程组1231231238322041133631236x x x x a x x x x x -+=⎧⎪+-=⎨⎪++=⎩，其精度为10-5。
四、 算法描述1. 雅克比迭代法公式：X[0][k+1] = ( 1 / a[0][0] )*[ b[0] –a[0][1]*X[1][k] – a[0][2]*X[2][k] ]; X[1][k+1] = ( 1 / a[1][1] )*[ b[1] –a[0][0]*X[0][k] – a[0][2]*X[2][k] ]; X[2][k+1] = ( 1 / a[2][2] )*[ b[2] –a[0][0]*X[0][k] – a[0][1]*X[1][k] ];2. 高斯-赛德尔迭代法公式：X[0][k+1] = ( 1 / a[0][0] )*[ b[0] –a[0][1]*X[1][k] – a[0][2]*X[2][k] ]; X[1][k+1] = ( 1 / a[1][1] )*[ b[1] –a[0][0]*X[0][k+1] – a[0][2]*X[2][k] ]; X[2][k+1] = ( 1 / a[2][2] )*[ b[2] –a[0][0]*X[0][k+1] – a[0][1]*X[1][k+1] ];3. 流程图雅克比流程图高斯-赛德尔流程图五、程序代码/*雅克比方法类Jacobi.h*/#ifndef A#define A#define MAX_SIZE 50class Jacobi{public:void Solve(int n,double a[MAX_SIZE][MAX_SIZE],double b[MAX_SIZE],intMAX_XunHuan,double e,int pre);};#endif/*雅克比方法类方法实现Jacobi.cpp*/#ifndef B#define B#include "Jacobi.h"#include "iostream"#include "math.h"#include "iomanip"using namespace std;double X[MAX_SIZE][MAX_SIZE]={0};double E[MAX_SIZE];void Jacobi::Solve(int n,double a[MAX_SIZE][MAX_SIZE],double b[MAX_SIZE],int MAX_XunHuan,double e,int pre){int i,j=1,z;double s=0;int flag=1;cout<<"请输入初始值X[i][0]:";for (i=0;i<n;i++){cin>>X[i][0];}for (j=1;j<MAX_XunHuan&&flag;j++){flag=0;for (i=0;i<n;i++){s=0;for (z=0;z<n;z++){if (i!=z){s+=a[i][z]*X[z][j-1];}}X[i][j]=1.0/a[i][i]*(b[i]-s);}E[j]=fabs(X[0][j]-X[0][j-1]);for (i=1;i<n;i++){if(fabs(X[i][j]-X[i][j-1])>E[j])E[j]=fabs(X[i][j]-X[i][j-1]);}if (E[j]>e){flag=1;}}if(flag==0){for (i=0;i<n;i++){cout<<" X"<<i<<"(k) ";}cout<<" 误差"<<endl;for (i=0;i<j;i++){cout<<setiosflags(ios::left)<<setiosflags(ios::fixed)<<setprecision(pre)<<X[0][i]<<" "<<X[1][i]<<" "<<X[2][i]<<" ";if(i>0){if(i!=j-1)cout<<E[i]<<endl;elsecout<<E[i]<<" < "<<e<<endl;}elsecout<<" ∞"<<endl;}}elsecout<<"没有结果,最大循环次数不够\n";cout<<"迭代次数为"<<j<<endl;}#endif/*/*高斯-赛德尔方法类Gauss.h*/#ifndef C#define C#define MAX_SIZEG 50class Gauss{public:void Solve(int n,double a[MAX_SIZEG][MAX_SIZEG],double b[MAX_SIZEG],int MAX_XunHuan,double e,int pre);};#endif/*高斯-赛德尔方法类方法实现Gauss.cpp*/#ifndef D#define D#include"Gauss.h"#include "math.h"#include "iostream"#include "iomanip"using namespace std;double XG[MAX_SIZEG][MAX_SIZEG]={0};double EG[MAX_SIZEG];void Gauss::Solve(int n,double a[MAX_SIZEG][MAX_SIZEG],double b[MAX_SIZEG],int MAX_XunHuan,double e,int pre){int i,j=1,z;double s=0;int flag=1;cout<<"请输入初始值XG[i][0]:";for (i=0;i<n;i++){cin>>XG[i][0];}for (j=1;j<MAX_XunHuan&&flag;j++){flag=0;for (i=0;i<n;i++){s=0;for (z=0;z<n;z++){if (i!=z){if(z>i)s+=a[i][z]*XG[z][j-1];elses+=a[i][z]*XG[z][j];}}XG[i][j]=1.0/a[i][i]*(b[i]-s);}EG[j]=fabs(XG[0][j]-XG[0][j-1]);for (i=1;i<n;i++){if(fabs(XG[i][j]-XG[i][j-1])>EG[j])EG[j]=fabs(XG[i][j]-XG[i][j-1]);}if (EG[j]>e){flag=1;}}if(flag==0){for (i=0;i<n;i++){cout<<" XG"<<i<<"(k) ";}cout<<" 误差"<<endl;for (i=0;i<j;i++){cout<<setiosflags(ios::left)<<setiosflags(ios::fixed)<<setprecision(pre)<<XG[0][i]<<" "<<XG[1][i]<<" "<<XG[2][i]<<" ";if(i>0){if(i!=j-1)cout<<EG[i]<<endl;elsecout<<EG[i]<<" < "<<e<<endl;}elsecout<<" ∞"<<endl;}}else cout<<"没有结果,最大循环次数不够\n";cout<<"迭代次数为"<<j<<endl;}#endif/*主函数Main*/#include "iostream"#include "Jacobi.h"#include "Jacobi.cpp"#include "Gauss.h"#include "Gauss.cpp"using namespace std;int main(){Jacobi J;Gauss G;int key;int i,j;int pre;double a[MAX_SIZE][MAX_SIZE]={0};double b[MAX_SIZE]={0};int n,MAX_XunHuan;double e;do{cout<<"请输入方程组的阶数：";cin>>n;cout<<"请输入系数矩阵A：";for (i=0;i<n;i++){for (j=0;j<n;j++){cin>>a[i][j];}}cout<<"请输入常数矩阵B：";for (i=0;i<n;i++){cin>>b[i];}cout<<"请输入最大循环次数：";cin>>MAX_XunHuan;cout<<"请输入最小误差：";cin>>e;cout<<"请输入精度输出控制（小数点后的位数）"; cin>>pre;cout<<"结束0：\n雅克比迭代法求解1：\n高斯-赛德尔求解2：";cin>>key;if (key==1){J.Solve(n,a,b,MAX_XunHuan,e,pre);}if (key==2){G.Solve(n,a,b,MAX_XunHuan,e,pre);}} while (key);return 0;}六、实验结果1.这次试验很简单，主要是公式问题，做得很顺利。
东南大学计算方法上机报告实验报告完整版
实习题11. 用两种不同的顺序计算644834.11000012≈∑=-n n，试分析其误差的变化解：从n=1开始累加，n 逐步增大，直到n=10000；从n=10000开始累加，n 逐步减小，直至1。
算法1的C 语言程序如下： #include<stdio.h> #include<math.h> void main() { float n=0.0; int i; for(i=1;i<=10000;i++) { n=n+1.0/(i*i); } printf("%-100f",n); printf("\n"); float m=0.0; int j; for(j=10000;j>=1;j--) { m=m+1.0/(j*j); } printf("%-7f",m); printf("\n"); }运行后结果如下：结论： 4.设∑=-=Nj N j S 2211，已知其精确值为)11123(21+--N N 。
1）编制按从大到小的顺序计算N S 的程序； 2）编制按从小到大的顺序计算N S 的程序；3）按2种顺序分别计算30000100001000,,S S S ，并指出有效位数。
解：1）从大到小的C语言算法如下：#include<stdio.h>#include<math.h>#include<iostream>using namespace std;void main(){float n=0.0;int i;int N;cout<<"Please input N"<<endl;cin>>N;for(i=N;i>1;i--){n=n+1.0/(i*i-1);N=N-1;}printf("%-100f",n);printf("\n");}执行后结果为：N=2时，运行结果为：N=3时，运行结果为：N=100时，运行结果为：N=4000时，运行结果为：2）从小到大的C语言算法如下：#include<stdio.h>#include<math.h>#include<iostream>using namespace std;void main(){float n=0.0;int i;int N;cout<<"Please input N"<<endl;cin>>N;for(i=2;i<=N;i++){n=n+1.0/(i*i-1);}printf("%-100f",n);printf("\n");}执行后结果为：N=2时，运行结果为：N=3时，运行结果为：N=100时，运行结果为：N=4000时，运行结果为：结论：通过比较可知：N 的值较小时两种算法的运算结果相差不大，但随着N 的逐渐增大，两种算法的运行结果相差越来越大。
计算方法C语言编程
计算方法C 语言编程第二章2已知方程043=-+x x 在区间[1，2]内有一根，试问用二分法求根，使其具有5位有效数字至少应二分多少次？ 【程序设计】#include<math.h> main(){int n=0; float x1=1.0,x2=2.0,x=1.0,x0;do{ x0=x; x=(x1+x2)/2; n++; if(x*x*x+x-4>0) x2=x; else x1=x;} while(fabs(x-x0)>0.00005); printf("N=%d\n",n);} 〖运行结果〗N=154用迭代法求02.05=--x x 的正根，要求准确到小数点后第5位。
【程序设计】#include<math.h> main(){ float x0, x=1.5,y;do{ y=(log(x+0.2))/5; x0=x; x=exp(y);} while(fabs(x-x0)>0.000005); printf(" X=%f\n",x);}- 〖运行结果〗x=1.0447639用牛顿法求方程0133=--x x在x0=2附近的根，要求准确到小数点后第3位。
【程序设计】 #include<math.h> main(){ float x=2.0,x0; do{x0=x; x=x-(x*x*x-3*x-1)/(3*x*x-3); } while(fabs(x-x0)>0.0005); printf("X=%f\n",x);} 〖运行结果〗x=1.87938511.分别用单点和双点弦截法求方程x3-x-1=0在[1,1.5]内的根。
要求|xn+1-xn|<=0.000005 【程序设计】#include<math.h> float f(float x){ float f; f=x*x*x-x-1; return f;} float g(float x1,float x2){ float g; g=(x1*f(x2)-x2*f(x1))/(f(x2)-f(x1)); return g;} main(){ float x1=1.0,x2=1.5,x,x0; x=x1;do{ x0=x; x=g(x1,x2); if(f(x)>0) x2=x; else x1=x;} while(fabs(x-x0)>0.000005); printf(" X=%f\n",x);} 〖运行结果〗x=1.324717第三章1.分别用列主元素消去法求解下列方程组.(计算取4位小数).⎪⎪⎩⎪⎪⎨⎧=+++=+++=+++=+++8471.12671.12568.02471.02368.07471.12271.02168.12071.01968.06471.11871.01768.01675.11582.05471.11490.01397.01254.01161.14321432143214321x x x x x x x x x x x x x x x x 【程序设计】#define n 4 main(){floata[n][n]={1.1161,0.1254,0.1397,0.1490,0.1582,1.1675,0.1768,0.1871, 0.1968,0.2071,1.2168,0.2271,0.2368,0.2471,0.2568,1.2671},b[n]={1.5471,1.6471,1.7471,1.8471},k,x[n+1],y[n+1]; int i,j,t; y[n]=0.0; for(i=0;i<n;i++)for(j=i+1;j<n;j++) { k=(-1)*a[j][i]/a[i][i]; a[j][i]=0; for(t=i+1;t<n;t++) a[j][t]=a[i][t]*k+a[j][t]; b[j]=b[i]*k+b[j]; } for(i=n-1;i>=0;i--){ for(j=n-2;j>=0;j--) { if(a[i][j]!=0) b[i]=b[i]-y[j+1]*a[i][j+1]; else break;} y[i]=b[i]/a[i][j+1]; x[i+1]=y[i]; }for(i=1;i<=n;i++) printf("X%d=%f\n",i,x[i]);} 〖运行结果〗x1=1.040584x2=0.986957 x3=0.935052 x4=0.8812979设有方程组 ⎪⎪⎩⎪⎪⎨⎧=+-=++--=++31032201225321321321x x x x x x x x x取初始向量=x )0()1,1,3(-T,分别用雅可比迭代法与赛德尔迭代法求解,要求x xk ik ii )()1(31max -+≤≤3≤时迭代终止.【程序设计】#include<math.h>main() /*赛德尔*/ { float x1[100],x2[100],x3[100],y; int k=0; x1[0]=-3;x2[0]=1;x3[0]=1; do{x1[k+1]=((-2)*x2[k]-x3[k]-12)/5； x2[k+1]=(x1[k+1]-2*x3[k]+20)/4; x3[k+1]=((-2)*x1[k+1]+3*x2[k+1]+3)/10;if(fabs(x3[k+1]-x3[k])>fabs(x2[k+1]-x2[k])) y=x3[k+1]-x3[k]; else if(fabs(x2[k+1]-x2[k])>fabs(x1[k+1]-x1[k])) y=x2[k+1]-x2[k]; else y=x1[k+1]-x1[k]; k++; } while(fabs(y)>0.001);printf("x1=%f\nx2=%f\nx3=%f\n",x1[k],x2[k],x3[k]); printf("K=%d\n",k); } 〖运行结果〗 x1=-3.999974x2=3.000043 x3=2.000008 k=7#include<math.h>main() /*雅可比*/{ float x1[100],x2[100],x3[100],y; int k=0; x1[0]=-3;x2[0]=1;x3[0]=1; do{ x1[k+1]=((-2)*x2[k]-x3[k]-12)/5; x2[k+1]=(x1[k+1]-2*x3[k]+20)/4;x3[k+1]=((-2)*x1[k+1]+3*x2[k+1]+3)/10;if(fabs(x3[k+1]-x3[k])>fabs(x2[k+1]-x2[k])) y=x3[k+1]-x3[k]; else if(fabs(x2[k+1]-x2[k])>fabs(x1[k+1]-x1[k])) y=x2[k+1]-x2[k]; else y=x1[k+1]-x1[k]; k++; } while(fabs(y)>0.001);printf("x1=%f\nx2=%f\nx3=%f\n",x1[k],x2[k],x3[k]); printf("K=%d\n",k);} 〖运行结果〗x1=-4.000152x2=2.999648 x3=2.000160 k=1310设有方程组 ⎪⎪⎩⎪⎪⎨⎧=++=++=-+5223122321321321x x x x x x x x x(1) 证明解此方程组的雅可比迭代法收敛,而相应的赛德尔迭代法发散. (2) 取初始向量)0,0,0()0(Tx =,用雅可比迭代法求解,要求迭代三次.【程序设计】main(){ float x1[10],x2[10],x3[10]; int k; x1[0]=0; x2[0]=0; x3[0]=0; for(k=0;k<=2;k++) { x1[k+1]=1-2*x2[k]+2*x3[k]; x2[k+1]=3-x1[k]-x3[k]; x3[k+1]=5-2*x1[k]-2*x2[k]; }printf("x1=%f\nx2=%f\nx3=%f\n",x1[k],x2[k],x3[k]);} 〖运行结果〗 x1=1.000000x2=1.000000 x3=1.00000011设有方程组 ⎪⎪⎪⎭⎫ ⎝⎛=⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎭⎫ ⎝⎛----010101021112321x x x其等价形式为 ⎪⎪⎩⎪⎪⎨⎧=+-=++-=1132123211x x x x x x x x(1) 证明解等价方程组的简单迭代法发散,而赛德尔法收敛 (2) 取初始向量)0,0,0()0(Tx =,用赛德尔迭代法求解,要求迭代四次【程序设计】main(){ float x1[10],x2[10],x3[10]; int k; x1[0]=0;x2[0]=0; x3[0]=0;for(k=0;k<=3;k++) { x1[k+1]=(x2[k]+x3[k])/2; x2[k+1]=(x1[k+1]+1)/2; x3[k+1]=x1[k+1]; } printf("x1=%f\nx2=%f\nx3=%f\n",x1[k],x2[k],x3[k]);} 〖运行结果〗x1=0.578125x2=0.789062 x3=0.578215第五章1.已知函数表：应用拉格朗日插值公式计算f(1.1300)的近似值。
matlabjacobi迭代法
matlabjacobi迭代法Jacobi迭代法是一种求解线性方程组的迭代法，其基本思想是将原方程组的系数矩阵分解为对角部分和非对角部分，对于对角矩阵使用前、后代替法求解，对于非对角部分使用迭代更新法求解。
Jacobi迭代法的基本形式如下：$\begin{cases}a_{11}x_1+a_{12}x_2+...+a_{1n}x_n=b_1 \\a_{21}x_1+a_{22}x_2+...+a_{2n}x_n=b_2 \\... \\a_{n1}x_1+a_{n2}x_2+...+a_{nn}x_n=b_n \\\end{cases}$其中，$a_{ij}$表示系数矩阵的第$i$行第$j$列的元素，$b_i$表示方程组的第$i$个方程的解。
设向量$x^{(k)}=(x_1^{(k)},x_2^{(k)},...,x_n^{(k)})$表示Jacobi迭代法的第$k$次迭代结果，则迭代公式为：$x_i^{(k+1)}=\frac{1}{a_{ii}}(b_i-\sum_{j=1,j\ne i}^n a_{ij}x_j^{(k)}),i=1,2,...,n$迭代公式的意义是，将第$i$个变量的系数$a_{ii}$看成系数矩阵的一个主对角元，将剩下的系数$a_{ij}(i\ne j)$看成非对角元，同时将当前未知量向量$x^{(k)}$看成已知量，利用这些参数求解第$i$个方程中未知量$x_i$。
Jacobi迭代法的收敛条件为原矩阵的对角线元素不为零，且矩阵的任意一行中非对角线元素绝对值之和小于对角线元素绝对值。
在Matlab中，可通过编写函数的方式实现Jacobi迭代法。
函数jacobi实现了迭代公式，并以向量形式返回迭代结果，如下所示：```function xnew = jacobi(A, b, xold)% Jacobi迭代法求解线性方程组Ax=b% A为系数矩阵，b为常数向量，xold为迭代初值% 输出迭代后的解向量xnew% 初始化迭代初值n = length(b);xnew = zeros(n,1);% 迭代更新for i = 1:nxnew(i) = (b(i) - A(i,:)*xold + A(i,i)*xold(i)) / A(i,i);endend```在主程序中可按以下步骤使用函数jacobi求解线性方程组：1.构造系数矩阵A和常数向量b；2.设定迭代初值xold；3.利用jacobi函数求解迭代结果，并对迭代过程进行循环。
雅克比迭代c++代码
一，雅可比迭代解方程C++代码#include<stdio.h>#include<math.h>#include <iostream>using namespace std;#define max_n 20#define maxrept 100#define epsilon 0.00001int main(){int n;int i,j,k;double err;static double a[max_n][max_n],b[max_n][max_n],c[max_n],g[max_n]; static double x[max_n],nx[max_n];cout<<"Input n value (dim of Ax=c):"<<endl;cin>>n;if(n>max_n){cout<<"The input n is larger than max_n,please redifine the max_n"; return 1;}if(n<=0){cout<<"please input a numbaer between 1 and "<<max_n;return 1;}cout<<"Now input the matrix a(i,j),i,j=0,..."<<n-1<<endl;for(i=0;i<n;i++)for(j=0;j<n;j++)cin>>a[i][j];cout<<"Now input the matrix c(i), i=0..."<<n-1<<endl;for(i=0;i<n;i++) cin>>c[i];for(i=0;i<n;i++)for(j=0;j<n;j++){b[i][j]=-a[i][j]/a[i][i];g[i]=c[i]/a[i][i];}for(i=0;i<maxrept;i++){for(j=0;j<n;j++)nx[j]=g[j];for(j=0;j<n;j++){for(k=0;k<n;k++){if(j==k) continue;nx[j]+=b[j][k]*x[k];}}err=0;for(j=0;j<n;j++)if(err<fabs(nx[j]-x[i])) err=fabs(nx[j]-x[j]);for(j=0;j<n;j++)x[j]=nx[j];if(err<epsilon){cout<<"Solve...x_i="<<endl;for(i=0;i<n;i++)cout<<x[i]<<endl;return 0;}}cout<<"After "<<maxrept<<" repeat,no result."<<endl; return 1;}。
研究生数值分析(11)---雅可比(Jacobi)迭代法
a x (k) n,n1 n1
bn )
取初始向量
X
(0)
(
x (0) 1
,
x (0) 2
,
,
x (0) n
)T
利用（4）反复迭代可以得到一个向量序列 {X (k)}
称式（4）为雅可比迭Jacobi代公式。
若记
a11
D
a22
0
a21 0
0 a12
0
L a31
a32
0
U
特征方程 I D1(L U ) 0
又可以写成 D1 D L U 0 因为 D1 0 ，所以 D L U 0 上式左端为将系数矩阵 A 的对角元同乘以 λ 后所得新矩阵的行列式。
例8 用雅可比迭代法求解方程组
10x1 2x2 x3 3 2x1 10x2 x3 15 x1 2x2 5x3 10
由迭代矩阵的特征方程
10 2 1 2 10 1 0 1 2 5
展开得到
(10 2)(50 2 10 3) 0
解得
1
1 5
, 2
1 10
7
, 3
1 10
7
于是 (J ) 1 7 0.3646 1
10
因而雅可比迭代公式是收敛的。
练习:考察用雅可比Jacobi迭代法解方程组 AX=b的收敛性，
解：相应的雅可比迭代公式为
x1(
k
1)
1 10
(2x2(k )
x (k) 3
3)
x2(k
1)
1 10
(2
x1(
k
)
x (k) 3
15)
x3(k
1)
1 5
(
x1(
雅克比迭代法python
雅克比迭代法python
-Jacobi迭代法是一种数值计算技术，用于求解非线性系统的迭代方法。
它将
非线性系统拆解为若干个独立的一维或二维子系统，利用迭代过程不断迭代地改进参数，最终收敛到满足约束条件的最优解。
Jacobi迭代法有利于在非线性系统中有效求解问题，它具有以下特点：
1、计算简单：Jacobi迭代法只需要计算每次迭代的细节，不需要求解任何解析表
达式，这种迭代方法可以有效地减少计算量和计算时间；
2、易编译：Jacobi迭代法只需要将等式离散化，然后对每组等式进行迭代，在程
序上比较容易编译；
3、稳定性好：Jacobi迭代法能够很快地收敛到最优解。
因此，Jacobi迭代法在不断优化的求解参数的过程中，以及优化非线性系统
的运算效率上，都具有较高的效率和准确度。
它的优势在于计算简单性和高效的稳定性，可以有效地提升计算效率，作为从大规模非线性系统中求解问题的一种技术，它在机器学习、信号处理、图像处理、线性系统控制、建模和函数优化等诸多领域都得到了广泛应用，受到学术界和实际应用界的高度重视。
雅可比迭代法解线方程组的mpi实现
雅可比迭代法解线方程组的MPI实现一、概述1.1 研究背景雅可比迭代法是一种常用的迭代法求解线性方程组的方法，通过不断迭代更新变量值，最终达到一定的精度要求。
MPI（Message Passing Interface）是一种消息传递接口，常用于并行计算中，能够充分利用多核和多机的计算资源，提高求解效率。
1.2 研究意义现代科学技术领域中，线性方程组的求解是一个常见且重要的问题。
通过研究雅可比迭代法在MPI并行环境下的实现，可以提高线性方程组求解的效率，促进科学计算领域的发展。
二、雅可比迭代法原理2.1 基本原理雅可比迭代法是一种逐次逼近的方法，通过不断迭代更新变量值，最终达到一定的精度要求。
假设方程组为Ax=b，其中A为系数矩阵，x为待求解的变量，b为常数向量。
迭代公式如下：x(k+1) = D^(-1)*(b-R*x(k))其中，D为系数矩阵A的对角线部分组成的对角矩阵，R为A的其余部分（不含对角线），k为迭代次数。
2.2 算法流程- 初始化变量x0- 重复执行以下步骤：- 计算新的变量值x(k+1)- 判断迭代是否收敛- 如果达到精度要求，停止迭代；否则，继续迭代三、MPI并行计算简介3.1 MPI概述MPI是一种消息传递接口，常用于并行计算中。
它定义了一组函数和语义，用于在多个进程之间传递消息，实现进程间的通信和同步。
3.2 并行计算模型在MPI并行计算中，通常采用主从模型，其中一个进程作为主进程，负责调度和管理其他子进程，而其他子进程负责实际的计算任务。
主进程可以将任务分发给子进程，并收集子进程的计算结果，实现并行计算。
四、雅可比迭代法在MPI中的实现4.1 算法说明在MPI并行环境下，雅可比迭代法的实现需要考虑如何将计算任务分配给不同的进程，以及如何进行进程间通信和同步。
一般可以采用分块迭代的方式，将系数矩阵A分块并分配给不同的进程，每个进程负责计算自己分配到的部分，然后进行通信和同步，最终得到整体的解。
数值分析 迭代法 C++程序
课题三解线性方程组的迭代法实验目标：分别采用Jacobi迭代法，Gauss-Seidel迭代法和SOR迭代法求解线性方程组。
Jocabi迭代法：#include<iostream>#include<math.h>using namespace std;int i,j,k; //计数器int M = 2000;int Array(double ***Arr, int n){double **p;int i;p=(double **)malloc(n*sizeof(double *));if(!p)return 0;for(i=0;i<n;i++){p[i]=(double *)malloc(n*sizeof(double));if(!p[i])return 0;}*Arr=p;return 1;}void main(){double eps ;cout<<"默认最多迭代次数为2000次"<<endl<<"迭代精度为：";cin>>eps;double **matrix;int n;cout<<"矩阵大小为：";cin>>n;double *X;X= new double[n];double *Y;Y= new double[n];double *G;G= new double[n];for(i=0;i<n;i++){Y[i]=0;}if(!Array(&matrix,n))cout<<"内存分配失败！";elsecout<<"请输入矩阵:"<<endl;for( i=0;i<n;i++){for( j=0;j<n;j++){cin>>matrix[i][j];}}cout<<"请输入右端项："<<endl;double *B;B = new double[n];for(i=0;i<n;i++){cin>>B[i];}for (i = 0 ;i< n;i++){if (fabs(matrix[i][i]) < eps){cout <<"打印失败"<<endl;return;}double T = matrix[i][i];for ( j = 0 ; j< n;j++){matrix[i][j] = -matrix[i][j]/T;}matrix[i][i] = 0;G[i] = B[i]/T;}int counter = 0;while (counter < M){for (i = 0;i < n; i++){double temp = 0;for (j = 0;j<n; j++){temp = temp + matrix[i][j]*Y[j];}X[i] = G[i] + temp;}double temp = 0;for (i = 0 ;i< n ; i++){temp = temp + fabs(X[i] - Y[i]);}if (temp <= eps)break;else{for( i = 0; i < n ;i++){Y[i] = X[i];}}counter++;}cout << "迭代次数为："<<counter<<"次。
雅可比迭代实验报告
雅可比迭代实验报告1. 实验目的通过实验了解雅可比迭代的原理和实现方法，掌握雅可比迭代程序的编写和调试能力，同时了解雅可比迭代的优缺点。
2. 实验原理雅可比迭代是一种常见的线性方程组求解方法。
通过逐次迭代使得线性方程组的解逼近真实解。
在雅可比迭代中，每一次迭代都会利用上一次迭代得到的解来更新新解的某个分量的值，并保持其他分量不变。
这样直至迭代收敛为止，得到解的精确值。
3. 实验步骤（1）建立线性方程组和初始解（2）计算出雅可比矩阵（3）利用当前解和雅可比矩阵计算新解（4）迭代步骤（3）直至收敛4. 实验要求（1）能够熟练掌握线性代数的基本概念、知识和方法，能够合理建立线性方程组和初始解。
（2）能够掌握矩阵的基本操作，计算雅可比矩阵和更新新解。
（3）能够控制迭代次数和精确度，进行合理的收敛判断。
（4）能够对比雅可比迭代和其他线性方程组求解方法的优缺点。
5. 实验结果设线性方程组为：\begin{cases}2x_1+x_2+x_3=5\\4x_1+3x_2+x_3=10\\3x_1+2x_2+4x_3=12\end{cases}初始解为：\begin{cases}x_{1(0)}=0\\x_{2(0)}=0\\x_{3(0)}=0\end{cases}设精度为10^{-4}，最大迭代次数为1000次。
利用以上参数，使用Python编写程序，计算出解为：\begin{cases}x_1=2\\x_2=1\\x_3=2\end{cases}6. 优缺点分析雅可比迭代的优点是简单易理解、易于实现、收敛速度快；缺点是迭代次数较多，尤其是在矩阵较大的情况下，计算量大，效率较低。
同时，当矩阵不满足对角占优条件时，收敛速度会变得很慢甚至不收敛。
因此，与其他线性方程组求解方法相比，雅可比迭代并不是最优解。
编程实现简单迭代
用C++编程实现简单迭代法班级：09医软一班汪晓阳 09713037雅克比（Jacobi）迭代法原理设有n阶方程组（4.1）若系数矩阵非奇异，且(i = 1, 2,…, n)，将方程组（4.1）改写成然后写成迭代格式（4.2）（4.2）式也可以简单地写为（4.3）对（4.2）或（4.3）给定一组初值后，经反复迭代可得到一向量序列，如果X (k)收敛于，则就是方程组（4.1）的解。
这一方法称为雅克比（Jacobi）迭代法或简单迭代法，（4.2）或（4.3）称为Jacobi迭代格式。
下面介绍迭代格式的矩阵表示：设D = diag (a11, a22, …, ann)，将AX = b改写为：AX = (D – (D - A)) x = bDX = (D - A) x + bX = (I – D-1A) x + D-1b记 B = I – D-1A F = D-1 b则迭代格式的向量表示为Ｂ称为雅克比迭代矩阵。
程序框图如下：Jacobi source code：//求解线性方程组——简单迭代法#include <iostream>#include <cmath>using namespace std;float *one_array_malloc(int n); //一维数组分配float **two_array_malloc(int m,int n); //二维数组分配float matrix_category(float* x,int n); int main(){const int MAX=100;//最大迭代次数int n,i,j,k;float** a;float* x_0; //初始向量float* x_k; //迭代向量float precision; //精度cout<<"输入精度e：";cin>>precision;cout<<endl<<"输入系数矩阵的阶数,N："; cin>>n;a=two_array_malloc(n,n+1);cout<<endl<<"输入增广矩阵的各值：\n"; for(i=0;i<n;i++){for(j=0;j<n+1;j++){cin>>a[i][j];}}x_0=one_array_malloc(n);cout<<endl<<"输入初始向量：\n"; for(i=0;i<n;i++){cin>>x_0[i];}x_k=one_array_malloc(n);float temp;//迭代过程for(k=0;k<MAX;k++){for(i=0;i<n;i++){temp=0;for(j=0;j<i;j++){temp=temp+a[i][j]*x_0[j];}x_k[i]=a[i][n]-temp;temp=0;for(j=i+1;j<n;j++){temp=temp+a[i][j]*x_0[j];}x_k[i]=(x_k[i]-temp)/a[i][i]; }//求两解向量的差的范数for(i=0;i<n;i++){x_0[i]=x_k[i]-x_0[i];}if(matrix_category(x_0,n)<precision) {break;}else{for(i=0;i<n;i++){x_0[i]=x_k[i];}}}//输出过程if(MAX==k){cout<<"迭代不收敛\n";}cout<<"迭代次数为："<<k+1<<endl;cout<<"解向量为：\n";for(i=0;i<n;i++){cout<<"x"<<i<<" "<<x_k[i]<<endl;}return 0;}float *one_array_malloc(int n) //一维数组分配{float *a;a=(float *)malloc(sizeof(float)*n);return a;}float **two_array_malloc(int m,int n) //二维数组分配{float **a;int i;a=(float **)malloc(m*sizeof(float *));for (i=0;i<m;i++){a[i]=(float *)malloc(n*sizeof(float));}return a;}float matrix_category(float* x,int n)//求向量范数{int i;float temp=fabs(x[0]);for(i=1;i<n;i++){if(temp<fabs(x[i]))temp=fabs(x[i]);}return temp;}测试例子：数值计算方法第205页例8，第201页例7,第212页11题。
雅可比(Jacobi)迭代法
0 a21 A a31 an1
0 a32
an2
0
ann1
a11 0
a22
0
ann
a12 0
a13 a23 0
a1n
a2n
an1n
0
记作 A = L + D + U
则 Ax b 等价于 (L D U )x b
数值计算方法
雅可比（Jacobi）迭代法
1.1雅可比迭代法算法构造
例4.2 用雅可比迭代法求解方程组
8x1 3x2 2x3 20 4x1 11x2 x3 33 6x1 3x2 12x3 36
解建：立从迭方代程公组式的三个方程中分离出 x1, x2 和 x3
x1x( k 11) xxxx32(( kk 3211))
B (I D 1 A) f D 1b
x (k 1) Bx (k ) f （k = 0,1,2…）
称为雅可比迭代公式, B称为雅可比迭代矩阵
其中
0
B
(I
D 1 A)
a21 a22
an1
ann
a12 a11
0
an2
ann
a1n a11
a2n
a22
0
在例4.2中,由迭代公式写出雅可比迭代矩阵为
i 1,2,, n
若上xia(式xkiii1称)0为ja(1解a1i1iiii方((1bb,程2ii , 组jj的,njnn1i )1Ja,aaijc分xijo(xjbk离i)j)迭)出代变i公i量式1,。21x,,2i , n , n ji
1.2 雅可比迭代法的矩阵表示
设方程组 Ax b 的系数矩阵A非奇异，且主对
雅克比高斯赛德尔迭代法
第八节 雅可比迭代法与高斯—塞德尔迭代法一 雅可比迭代法设线性方程组b Ax = (1) 的系数矩阵A 可逆且主对角元素nn a ,...,a ,a 2211均不为零,令()nna ,...,a ,a diag D 2211=并将A 分解成()D D A A +-= (2)从而(1)可写成 ()b x A D Dx +-=令11f x B x +=其中b D f ,A D I B 1111--=-=. (3) 以1B 为迭代矩阵的迭代法(公式)()()111f x B x k k +=+ (4)称为雅可比(Jacobi)迭代法(公式),用向量的分量来表示,(4)为⎩⎨⎧[],...,,k ,n ,...,i x a ba xnij j )k (j j i iii)k (i21021111==∑-=≠=+ (5)其中()()()()()Tn x ,...x ,x x 002010=为初始向量.由此看出,雅可比迭代法公式简单,每迭代一次只需计算一次矩阵和向量的乘法.在电算时需要两组存储单元,以存放()k x 及()1+k x . 例1 例1 用雅可比迭代法求解下列方程组⎪⎩⎪⎨⎧=+--=-+-=--2453821027210321321321.x x x .x x x .x x x解 将方程组按雅可比方法写成⎪⎪⎩⎪⎪⎨⎧++=++=++=8402020830201072020*******2321.x .x .x .x .x .x .x .x .x取初始值()()()()()()T T ,,,x ,x ,x x 0000302010==按迭代公式()()()()()()()()()⎪⎪⎩⎪⎪⎨⎧++=++=++=+++840202083020107202010211331123211.x .x .x .x .x .x .x .x .x k k k k k k k k k进行迭代，其计算结果如表1所示表1二 高斯—塞德尔迭代法由雅可比迭代公式可知,在迭代的每一步计算过程中是用()k x的全部分量来计算()1+k x的所有分量,显然在计算第i 个分量()1+k i x 时,已经计算出的最新分量()()1111+-+k i k x ,...,x 没有被利用，从直观上看,最新计算出的分量可能比旧的分量要好些.因此，对这些最新计算出来的第1+k 次近似()1+k x的分量()1+k jx 加以利用,就得到所谓解方程组的高斯—塞德（Gauss-Seidel ）迭代法.把矩阵A 分解成U L D A --= (6)其中()nn a ,...,a ,a diag D 2211=,U ,L --分别为A 的主对角元除外的下三角和上三角部分,于是,方程组(1)便可以写成 ()b Ux x L D +=-即 22f x B x +=其中()()b L D f ,U L D B 1212---=-= (7)以2B 为迭代矩阵构成的迭代法(公式)()()221f x B x k k +=+ (8)称为高斯—塞德尔迭代法(公式),用 量表示的形式为⎩⎨⎧[],...,,k ,n ,,i x a x a b a xi j n i j )k (j ij )k (j ij i ii)k (i21021111111==∑∑--=-=+=++ (9)由此看出,高斯—塞德尔迭代法的一个明显的优点是,在电算时,只需一组存储单元(计算出()1+k ix 后()k ix 不再使用,所以用()1+k i x 冲掉()k ix ,以便存放近似解.例2 例2 用高斯——塞德尔迭代法求解例1.解 取初始值()()()()()()TT,,,x ,x ,x x 0000302010==，按迭代公式()()()()()()()()()⎪⎩⎪⎨⎧++=++=++=++++++840202083020107202010121113311123211.x .x .x .x .x .x .x .x .x k k k k k k k k k进行迭代，其计算结果如下表2从此例看出,高斯—塞德尔迭代法比雅可比迭代法收敛快(达到同样的精度所需迭代次数少),但这个结论,在一定条件下才是对的,甚至有这样的方程组,雅可比方法收敛，而高斯—塞德尔迭代法却是发散的.三 迭代收敛的充分条件定理1 在下列任一条件下,雅可比迭代法(5)收敛.①111<∑=≠=∞nij j iiij ia a max B ;②1111<∑=≠=nij i iiij ja a max B ;③ 111<∑=-≠=∞-nji i jjij jTa a max AD I定理2 设21B B ,分别为雅可比迭代矩阵与高斯—塞德尔迭代矩阵,则∞∞≤12B B (10)从而，当111<∑=≠=∞nij j iiij ia a max B时，高斯—塞德尔迭代法(8)收敛. 证明 由21B B ,的定义,它们可表示成()U L D B +=-11()()U D L D I U L D B 11112-----=-=用e 表示n 维向量()T,...,,e 111=,则有不等式eB e B ∞≤11UD L D B 111--+=这里,记号｜·｜表示其中矩阵的元素都取绝对值,而不等式是对相应元素来考虑的,于是()()()Ie B L D I eL D B e U D ∞------≤-=111111容易验证()11==--nnL D L D所以,L D I 1--及L D I 1--可逆，且()()()1111111111-----------=++≤+++=-L D I LD ...L D I L D ...L D I LD I n n()I L D I ≥---11从而有()()((){}e I B L D I L D I eU D LD I e B ∞----------≤⋅-≤111111121{()()}eB eL D I I B I ∞--∞≤-⋅--=11111因此必有∞∞≤12B B因为已知11<∞B 所以12<∞B .即高斯—塞德尔迭代法收敛.若矩阵A 为对称，我们有定理3 若矩阵A 正定,则高斯—塞德尔迭代法收敛.证明 把实正定对称矩阵A 分解为T L L D A --=()TL U =,则D 为正定的,迭代矩阵()T L L D B 12--=设λ是2B 的任一特征值,x 为相应的特征向量,则()()x x L L D T λ=--1以L D -左乘上式两端,并由T L L D A --=有()Ax x L T λλ=-1用向量x 的共轭转置左乘上式两端,得()Ax x x L xTTT--=-λλ1 (11)求上式左右两端的共轭转置,得Ax x x L x T T ----=⎪⎭⎫ ⎝⎛-λλ1以λ--1和λ-1分别乘以上二式然后相加,得()()Axx x L L x T T T -----⎪⎭⎫ ⎝⎛-+=+⎪⎭⎫ ⎝⎛--λλλλλλ211 由TL L D A --=,得()()Axx x A D x T T -----⎪⎭⎫ ⎝⎛-+=-⎪⎭⎫ ⎝⎛--λλλλλλ211即()Ax x x L x TT---=-λλλ2211 (12)因为A 和D 都是正定的,且x 不是零向量,所以由(11)式得1≠λ,而由(12)式得012>-λ, 即1<λ,从而()12<B ρ,因而高斯—塞德尔迭代法收敛.定义1 设()nn ij a A ⨯=为n 阶矩阵.① ①如果n,...,i ,a a nij i j ij ii 21=∑>≠= (13)即A 的每一行对角元素的绝对值都严格大于同行其他元素绝对值之和，则称A 为严格对角优势矩阵.② ②如果n,...,i ,a a nij i j ij ii 21=∑≥≠=且至少有一个不等式严格成立,则称A 为弱对角优势矩阵.例如⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-310131012是严格对角优势矩阵，⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--310121011是弱对角优势矩阵. 定义2 设()nn ij a A ⨯=是n 阶矩阵，如果经过行的互换及相应列的互换可化为⎥⎦⎤⎢⎣⎡2212110A A A ,即存在n 阶排列矩阵P,使⎥⎦⎤⎢⎣⎡=2212110A A A AP P T其中2211A ,A 为方阵，则称A 是可约的,否则称A 为不可约的.A 是可约矩阵,意味着b Ax =可经过若干次行列重排,化为两个低阶方程组,事实上,b Ax =可化为 ()b P x P AP P TT T =,记()()()()⎥⎦⎤⎢⎣⎡==⎥⎦⎤⎢⎣⎡==2121d d d b P ,y y y x P TT于是，求解b Ax =化为求解()()()()()⎪⎩⎪⎨⎧=+=+22221212111d y A d y A y A可以证明,如果A 为严格对角优势矩阵或为不可约弱对角优势矩阵,则A 是非奇异的.定理4 如果A 为严格对角优势矩阵或为不可约弱对角优势矩阵,则对任意()0x ,雅可比迭代法(4)与高斯—塞德尔迭代法(8)均为收敛的.证明 下面我们以A 为不可约弱对角优势矩阵为例，证明雅可比迭代法收敛，其他证明留给读者.要证明雅可比迭代法收敛，只要证()11<B ρ,1B 是迭代矩阵.用反证法,设矩阵1B 有某个特征值μ,使得1≥μ,则()01=-B I det μ,由于A 不可约,且具有弱对角优势，所以1-D 存在，且 ()()D A D D A D I I B I -+=--=---μμμ111从而()0=-+D A D detμ另一方面，矩阵()D A D -+μ与矩阵A 的非零元素的位置是完全相同的，所以()D A D -+μ也是不可约的,又由于1≥μ,且A 弱对角优势，所以n,...,i ,a a a nij i j ij ii ii 21=∑≥≥≠=μ并且至少有一个i 使不等号严格成立.因此,矩阵()D A D -+μ弱对角优势，故()D A D -+μ为不可约弱对角优势矩阵.从而()0≠-+D A D det μ矛盾，故1B 的特征值不能大于等于1，定理得证.。
雅可比(Jacobi)迭代法
高斯赛德尔迭代矩阵BG一般不容易计算，所以实际使用 时采用分量形式的算法，参见程序 GaussSeidelit2.m
例子：p.55(p.52)例8 ，10-3的精度，迭代6 次。
3x1x12xx22
5 5
x(k 1) 1
x(k) 2 3
x(k) i
(bi
a x( k1) ij j
aij
x
( j
k
)
)
/
aii
j 1
ji
不同的 的值会影响SOR迭代的收敛性、收敛 速度。
20
例（7）SOR迭代法
8 3 2 A 4 11 1
6 3 12
取 =1.5，则迭代矩阵：
1 / 2 9 / 16
3 / 8
B 3 /11 71/ 88 15 / 44
|| B || 20, || B ||1 17, || B ||2 14.4, (B) 13
不收敛。
14
（2）简单构造迭代法-2
8x1 3x2 2x3 20 4x1 11x2 x3 33 6x1 3x2 12x3 36
2
3
4x1 20 4x1 3x2 2x3
9x2 33 4x1 2x2 x3
举例：
8 4
x1 x1
3x2 2x3 11x2 x3
20 33
6 x1 3 x2 12x3 36
精确解
3 2 1
13
（1）简单迭代法
8 3 2 7 3 2 B I A I 4 11 1 4 10 1
6 3 12 6 3 11 20 b' 33 36
SOR迭 代( 1.3545)， 17次 ， (B) 0.452847
雅可比迭代实验报告
雅可比迭代实验报告一、实验目的通过实验掌握雅可比迭代方法的基本原理及步骤，能够在Matlab软件平台上编写该迭代方法，进一步加深对线性方程组迭代算法的认识。
二、实验原理雅可比迭代法是一种求解线性方程组解的迭代算法，它是最古老的、最基本的迭代算法之一，常用于解决一些线性方程组的求解问题。
假设一个线性方程组为Ax=b，其中A为一个系数矩阵，b为一个已知项向量，x为一个未知向量，且A是对称正定矩阵。
由于方程组的解是未知的，我们先取一个近似的解向量，用它替换掉上式中的x，得到：Ax(1)=b 或 Ax1=b将式子改写为：其中D是A的对角线构成的对角矩阵。
那么，我们就可以进行一次迭代，将上一步求得的解向量代入到式中，继续进行计算，如此类推。
直到误差满足一定的条件时，则结束迭代。
三、实验步骤本次实验的具体步骤如下：1、输入待求解的线性方程组，其中系数矩阵A需是对称正定矩阵。
2、根据所求的系数矩阵A，求解对角矩阵D，以及矩阵A-D的逆矩阵。
3、根据输入的误差限和最大迭代次数，设定初始的近似解向量x0，并计算误差。
4、进行迭代，重复进行以下操作：（1）代入上一步求得的解向量，计算新的解向量。
（3）如果误差小于所设定的误差限或者迭代次数已经达到了设定的最大迭代次数，则结束迭代。
5、输出求解结果及迭代次数。
四、实验结果通过设计实验，我们在Matlab软件平台上实现了雅可比迭代的算法。
在进行MATLAB编程时，我们遵循上述步骤，利用所编写的Matlab程序求解了如下线性方程组：1 1 1 x1 52 3 5 x2 63 6 9 x3 9并且设定误差限为10^-5和最大迭代次数为100，并取初始迭代向量为[0 0 0]。
运行该程序后，得到的计算结果为：迭代次数：39解向量：[5 0.999 1]五、实验分析通过对Matlab代码的编写及运行，我们了解了雅可比迭代算法的具体实现过程，以及在实际应用中它具有的一些优点和不足。
雅可比迭代法的 C++编程实现 - search readpudncom
雅可比迭代法的C++编程实现¾ 算法的数学说明：设有方程组=b （i=1,2,…,n ）记为x a j n j ij ∑=1iAx=b, （1式） A 为非奇异矩阵且0（i=1,2,…,n ）。
将A 分裂为A=D-L-U ，其中a ij ≠D ＝0a 21 0a 31 a 320. . .. . .a n1a n2 0L ＝－U ＝－将（1式）第i 个方程（i=1,2,…,n ）用去除再移项，得到等价方程组a ijx i ＝ii a 1(b －) （i=1,2,…,n ）， （2式）i ∑≠=n i j j j ij x a 1简记为x=Bx+f ，其中B=I －D −A= D (L+U)，f= D b 。
11−1−对方程2式）使用迭代法，得到解（1式）的雅可比（Jacobi ）迭代公式 组（ x(0)=(x ，…，x )T （初始向量）； (0)1(0)n x 1)(k i +=a ii1(b i －∑)， （3式） ≠=n i j j ija 1(k)j x其中，=(x ，…，)为第k 次迭代向量，设已经算出，由（3式）可计算下一次迭代向量（k=0,1,2,…；i=1,2,…,n ）。
x (k)(k)1x x (k)n T x (k)1)(k +迭代（3式）的矩阵形式为公式x(0)x (k 初始向量， 1)+=Bx +f ， （4式）(k)其中B 为雅可比方法迭代矩阵，初始向量用高斯列主元素消去法求解。
¾ 程序源代码如下：/*-----------------------------------------------假设有如下方程组：Ax=b用Jacobi 迭代法求解方程组的解方法：将A 分裂为A=D-L-U ，等价的迭代方程组为x=Bx+f 。
-----------------------------------------------#include <iostream.h>#include <stdlib.h>#include <math.h>double* allocMem(int ); //分配内存空间函数void GaussLineMain(double*,double*,double*,int );//采用高斯列主元素消去法求解x的初始向量值void Jacobi(double*,double*,double*,double*,int,int);//利用雅可比迭代公式求解x的值void main(){matrixNum;//矩阵的行数(列数)short//矩阵A，初始系数矩阵*matrixA;double*matrixD;//矩阵D为A中的主对角阵double*matrixL; //矩阵L为A中的下三角阵double//矩阵U为A中的上三角阵double*matrixU;*B; //矩阵B为雅可比方法迭代矩阵doubledouble*f; //矩阵f为中间的过渡的矩阵*x; //x为一维数组，存放结果double*xk; //xk为一维数组，用来在迭代中使用doubledouble*b; //b为一维数组，存放方程组右边系数i,j,k;intcout<<"<<请输入矩阵的行数（列数与行数一致）>>：";cin>>matrixNum;//分别为A、D、L、U、B、f、x、b分配内存空间matrixA=allocMem(matrixNum*matrixNum);matrixD=allocMem(matrixNum*matrixNum);matrixL=allocMem(matrixNum*matrixNum);matrixU=allocMem(matrixNum*matrixNum);B=allocMem(matrixNum*matrixNum);f=allocMem(matrixNum);x=allocMem(matrixNum);xk=allocMem(matrixNum);b=allocMem(matrixNum);//输入系数矩阵各元素值cout<<endl<<endl<<endl<<"<<请输入矩阵中各元素值（为 "<<matrixNum<<"*"<<matrixNum<<"，共计 "<<matrixNum*matrixNum<<" 个元素）"<<">>："<<endl<<endl;for(i=0;i<matrixNum;i++){cout<<"请输入矩阵中第 "<<i+1<<" 行的 "<<matrixNum<<" 个元素：";for(j=0;j<matrixNum;j++)cin>>*(matrixA+i*matrixNum+j);}//输入方程组右边系数b的各元素值cout<<endl<<endl<<endl<<"<<请输入方程组右边系数各元素值，共计"<<matrixNum<<" 个"<<">>："<<endl<<endl;for(i=0;i<matrixNum;i++)cin>>*(b+i);/* 下面将A分裂为A=D-L-U *///首先将D、L、U做初始化工作for(i=0;i<matrixNum;i++)for(j=0;j<matrixNum;j++)*(matrixD+i*matrixNum+j)=*(matrixL+i*matrixNum+j)=*(matrixU+i*matrixNum+j)=0;//D、L、U分别得到A的主对角线、下三角和上三角；其中D取逆矩阵、L和U 各元素取相反数for(i=0;i<matrixNum;i++)for(j=0;j<matrixNum;j++)if(i==j&&*(matrixA+i*matrixNum+j))*(matrixD+i*matrixNum+j)=1/(*(matrixA+i*matrixNum+j));*(matrixL+i*matrixNum+j)=-*(matrixA+i*matrixNum+j);if(i>j)else*(matrixU+i*matrixNum+j)=-*(matrixA+i*matrixNum+j);else//求B矩阵中的元素for(i=0;i<matrixNum;i++)for(j=0;j<matrixNum;j++){doubletemp=0;for(k=0;k<matrixNum;k++)temp+=*(matrixD+i*matrixNum+k)*(*(matrixL+k*matrixNum+j)+*(matrixU+k*matrixNu m+j));*(B+i*matrixNum+j)=temp;}//求f中的元素for(i=0;i<matrixNum;i++){temp=0;doublefor(j=0;j<matrixNum;j++)temp+=*(matrixD+i*matrixNum+j)*(*(b+j));*(f+i)=temp;}/* 计算x的初始向量值 */GaussLineMain(matrixA,x,b,matrixNum);/* 利用雅可比迭代公式求解xk的值 */intJacobiTime;cout<<endl<<endl<<endl<<"<<雅可比迭代开始，请输入希望迭代的次数>>：";cin>>JacobiTime;while(JacobiTime<=0){cout<<"迭代次数必须大于0，请重新输入：";cin>>JacobiTime;}Jacobi(x,xk,B,f,matrixNum,JacobiTime);//输出线性方程组的解 */cout<<endl<<endl<<endl<<"<<方程组运算结果如下>>"<<endl;//设置输出精度，以此比较不同迭代次数的结果 cout.precision(20);for(i=0;i<matrixNum;i++)"<<*(xk+i)<<endl;cout<<"x"<<i+1<<"=cout<<endl<<endl<<"求解过程结束..."<<endl<<endl;//释放掉所有动态分配的内存delete [] matrixA;delete [] matrixD;delete [] matrixL;delete [] matrixU;delete [] B;delete [] f;delete [] x;delete [] xk;delete [] b;}/*----------------------分配内存空间函数----------------------*/double* allocMem(int num){double*head;double[num])==NULL)if((head=new{cout<<"内存空间分配失败，程序终止运行！"<<endl; exit(0);}head;return}/*-----------------------------------------------计算Ax=b中x的初始向量值，采用高斯列主元素消去法-----------------------------------------------*/void GaussLineMain(double* A,double* x,double* b,int num) {inti,j,k;//共处理num-1行for(i=0;i<num-1;i++){//首先每列选主元，即最大的一个lineMax=fabs(*(A+i*num+i));doublelineI=i;intfor(j=i;j<num;j++)lineI=j;if(fabs(*(A+j*num+i))>fabs(lineMax))//主元所在行和当前处理行做行交换，右系数b也随之交换for(j=i;j<num;j++){//A做交换lineMax=*(A+i*num+j);*(A+i*num+j)=*(A+lineI*num+j);*(A+lineI*num+j)=lineMax;//b中对应元素做交换lineMax=*(b+i);*(b+i)=*(b+lineI);*(b+lineI)=lineMax;}continue; //如果当前主元为0，本次循环结束 if(*(A+i*num+i)==0)//将A变为上三角矩阵，同样b也随之变换for(j=i+1;j<num;j++){temp=-*(A+j*num+i)/(*(A+i*num+i));doublefor(k=i;k<num;k++){*(A+j*num+k)+=temp*(*(A+i*num+k));}*(b+j)+=temp*(*(b+i));}}验证Ax=b是否有唯一解，就是验证A的行列式是否为0；/*如果|A|!=0，说明有唯一解*/determinantA=1;doublefor(i=0;i<num;i++)determinantA*=*(A+i*num+i);if(determinantA==0){cout<<endl<<endl<<"通过计算,矩阵A的行列式为|A|=0，即A没有唯一解。
用C语言实现雅可比迭代法
用C语言实现雅可比迭代法#include "stdlib.h"#include "stdio.h"#include "conio.h"#include "string.h"#include "math.h"#define N 100float Table(int n,float a[N][N],float b[N]) {int i,j;float c[N][N];printf("请输入矩阵的行!\n");for(i=0;i<n;i++){printf("请输入行 %d:",i+1);for(j=0;j<n;j++)scanf("%f",&a[i][j]);}printf("请输入向量b:");for(i=0;i<n;i++)scanf("%f",&b[i]);for(i=0;i<n;i++)for(j=0;j<n;j++){if(i==j){c[i][j]=0;continue;}c[i][j]=-a[i][j]/a[i][i];}printf("\n矩阵A和向量b:\n");for(i=0;i<n;i++){for(j=0;j<n;j++)printf("%10.5f",a[i][j]);printf("%10.5f",b[i]);printf("\n");}printf("\n高斯-赛德尔迭代的方案:\n"); for(i=0;i<n;i++){for(j=0;j<n;j++)printf("%10.5f",c[i][j]);printf("%10.5f",b[i]/a[i][i]);printf("\n");}}float init_vec(int n,float x[N]){int i;printf("请输入初始迭代向量x:");for(i=0;i<n;i++)scanf("%f",&x[i]);printf("\n初始迭代向量x:\n");for(i=0;i<n;i++)printf("%10.5f",x[i]);printf("\n");}float jacobi(int n,float a[N][N],float b[N],float x[N]) {int i,j,k;float tmp,x2[N];for(k=0;;k++){for(i=0;i<n;i++)x2[i]=x[i];for(i=0;i<n;i++){tmp=0.0;for(j=0;j<n;j++){if(j==i) continue;tmp+=a[i][j]*x2[j];}x[i]=(b[i]-tmp)/a[i][i];}for(i=0,j=0;i<n;i++)if(fabs(x2[i]-x[i])<0.00001) j++;if(j==n){printf("\n这个高斯-赛德尔迭代方案收敛!\n");printf("迭代次数:%d",k+1);break;}if(k==499){printf("\n这雅可比迭代计划可能不收敛!");}}printf("\n结果为：\n"); for(i=0;i<n;i++) printf("%12.7f",x[i]); }int main(){int n;float x[N],a[N][N],b[N]; printf("请输入n:"); scanf("%d",&n); Table(n,a,b);init_vec(n,x);jacobi(n,a,b,x); getch();}</n;i++)</n;i++)</n;j++)</n;i++)</n;i++)</n;i++)</n;i++)</n;j++)</n;i++)</n;j++)</n;j++) </n;i++) </n;i++) </n;j++) </n;i++)。
雅克比迭代法……
{ a[i]=(float *)malloc(n*sizeof(float)); return a; } float matrix_category(float* x,int n) { int i; float temp=0; for(i=0;i<n;i++) { temp=temp+fabs(x[i]); } return temp; }
结束
④三角分解法 #include<iostream> using namespace std; int main() { const int N=100; static double a[N][N],b[N]; int i,j,k,num,p; double m,t,q; cout<<"请输入矩阵阶数 请输入矩阵阶数： cout<<"请输入矩阵阶数："; cin>>num;
cout<<endl<<"输入初始向量 输入初始向量: cout<<endl<<"输入初始向量:\n"; for(i=0;i<n;i++) { cin>>x_0[i]; } x_k=one_array_malloc(n); cout<<"输入松弛因子 (1<w<2)： cout<<"输入松弛因子 w (1<w<2)：\n"; cin>>w; float temp; for(k=0;k<MAX;k++) for(k=0;k<MAX;k++) { for(i=0;i<n;i++) { temp=0; for(j=0;j<i;j++) { temp=temp+a[i][j]*x_k[j]; } x_k[i]=a[i][n]x_k[i]=a[i][n]-temp; temp=0; for(j=i+1;j<n;j++) { temp=temp+a[i][j]*x_0[j]; } x_k[i]=(x_k[i]x_k[i]=(x_k[i]-temp)/a[i][i]; x_k[i]=(1x_k[i]=(1-w)*x_0[i]+w*x_k[i]; } for(i=0;i<n;i++) x_0[i]=x_k[i]{ x_0[i]=x_k[i]-x_0[i]; } if(matrix_category(x_0,n)<precision) { break; } else { for(i=0;i<n;i++) { x_0[i]=x_k[i]; } } } if(MAX==k) cout<<"迭代不收敛 迭代不收敛\ { cout<<"迭代不收敛\n"; } cout<<"迭代次数为 迭代次数为： cout<<"迭代次数为："<<k<<endl; cout<<"解向量为 解向量为： cout<<"解向量为：\n"; for(i=0;i<n;i++) { cout<<"x"<<i<<": "<<x_k[i]<<endl; } return 0; } float *one_array_malloc(int n) { float *a; a=(float *)malloc(sizeof(float)*n); return a; } float **two_array_malloc(int m,int n) { float **a; int i; **)malloc(m*sizeof(float a=(float **)malloc(m*sizeof(float *)); for (i=0;i<m;i++)
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
for(i=0,j=0;i<n;i++)
if(fabs(x2[i]-x[i])<0.00001) j++;
if(j==n)
}
int main()
{
int n;
float x[N],a[N][N],b[N];
printf("请输入n:");
scanf("%d",&n);
Table(n,a,b);
{
printf("\n这个高斯-赛德尔迭代方案收敛!\n");
printf("迭代次数:%d",k+1);
break;
}
if(k==499)
init_vec(n,x);
jacobi(n,a,b,x);
getch();
}
}
printf("请输入向量b:");
for(i=0;i<n;i++)
scanf("%f",&b[i]);
for(i=0;i<n;i++)
for(j=0;j<n;j++)
}
printf("\n高斯-赛德尔迭代的方案:\n");
for(i=0;i<n;i++)
{
for(j=0;j<n;j++)
printf("%10.5f",c[i][j]);
printf("%10.5f",b[i]/a[i][i]);
for(j=0;j<n;j++)
{
if(j==i) continue;
tmp+=a[i][j]*x2[j];
}
x[i]=(b[i]-tmp)/a[i][i];
}
printf("\n矩阵A和向量b:\n");
for(i=0;i<n;i++)
{
for(j=0;j<n;j++)
printf("%10.5f",a[i][j]);
printf("%10.5f",b[i]);
printf("\n");
printf("\n");
}
}
float init_vec(int n,float x[N])
{
int i;
printf("请输入初始迭代向量x:");
for(i=0;i<n;i++)
scanf("%f",&x[i]);
#include "stdlib.h"
#include "stdio.h"
#include "conio.h"
#include "string.h"
#include "math.h"
#define N 100
float Table(int n,float a[N][N],float b[N])
{
int i,j,k;
float tmp,x2[N];
for(k=0;;k++)
{
for(i=0;i<n;i++)
x2[i]=x[i];
for(i=0;i<n;i++)
{
tmp=0.0;
{
if(i==j)
{
c[i][j]=0;
continue;
}
c[i][j]=-a[i][j]/a[i][i];
{
int i,j;
float c[N][N];
printf("请输入矩阵的行!\n");
for(i=0;i<n;i++)
{
printf("请输入行 %d:",i+1);
for(j=0;j<n;j++)
scanf("%f",&a[i][j]);
printf("\n初始迭代向量x:\n");
for(i=0;i<n;i++)
printf("%10.5f",x[i]);
printf("\n");
}
float jacobi(int n,float a[N][N],float b[N],float x雅可比迭代计划可能不收敛!");
break;
}
}
printf("\n结果为：\n");
for(i=0;i<n;i++)
printf("%12.7f",x[i]);