雅可比(Jacobi)迭代算法的C++实现

admin · 发表于 2024-1-18 21:43:30

/*
　　-----------------------------------------------
　　假设有如下方程组：
　　Ax=b
　　用Jacobi迭代法求解方程组的解
　　方法：将A分裂为A=D-L-U，等价的迭代方程组为x=Bx+f。
　　有关算法的详细说明，参看http://www.loujing.com/mywork/c++/project/Jacobi.pdf
　　-----------------------------------------------
　　*/
　　#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()
　　{
　　short matrixNum; //矩阵的行数(列数)
　　double *matrixA; //矩阵A，初始系数矩阵
　　double *matrixD; //矩阵D为A中的主对角阵
　　double *matrixL; //矩阵L为A中的下三角阵
　　double *matrixU; //矩阵U为A中的上三角阵
　　double *B; //矩阵B为雅可比方法迭代矩阵
　　double *f; //矩阵f为中间的过渡的矩阵
　　double *x; //x为一维数组，存放结果
　　double *xk; //xk为一维数组，用来在迭代中使用
　　double *b; //b为一维数组，存放方程组右边系数
　　int i,j,k;
　　cout<<"<<请输入矩阵的行数（列数与行数一致）>>：";
　　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));
　　 else if(i>j) *(matrixL+i*matrixNum+j)=-*(matrixA+i*matrixNum+j);
　　 else *(matrixU+i*matrixNum+j)=-*(matrixA+i*matrixNum+j);
　　//求B矩阵中的元素
　　for(i=0;i<matrixNum;i++)
　　 for(j=0;j<matrixNum;j++)
　　 {
　　 double temp=0;
　　 for(k=0;k<matrixNum;k++)
　　 temp+=*(matrixD+i*matrixNum+k)*(*(matrixL+k*matrixNum+j)+*(matrixU+k*matrixNum+j));
　　 *(B+i*matrixNum+j)=temp;
　　 }
　　//求f中的元素
　　for(i=0;i<matrixNum;i++)
　　{
　　 double temp=0;
　　 for(j=0;j<matrixNum;j++)
　　 temp+=*(matrixD+i*matrixNum+j)*(*(b+j));
　　 *(f+i)=temp;
　　}
　　/* 计算x的初始向量值 */
　　GaussLineMain(matrixA,x,b,matrixNum);
　　/* 利用雅可比迭代公式求解xk的值 */
　　int JacobiTime;
　　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++)
　　 cout<<"x"<<i+1<<" = "<<*(xk+i)<<endl;
　　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;
　　if((head=new double[num])==NULL)
　　{
　　 cout<<"内存空间分配失败，程序终止运行！"<<endl;
　　 exit(0);
　　}
　　return head;
　　}
　　/*
　　-----------------------------------------------
　　计算Ax=b中x的初始向量值，采用高斯列主元素消去法
　　-----------------------------------------------
　　*/
　　void GaussLineMain(double* A,double* x,double* b,int num)
　　{
　　int i,j,k;
　　//共处理num-1行
　　for(i=0;i<num-1;i++)
　　{
　　 //首先每列选主元，即最大的一个
　　 double lineMax=fabs(*(A+i*num+i));
　　 int lineI=i;
　　 for(j=i;j<num;j++)
　　 if(fabs(*(A+j*num+i))>fabs(lineMax)) lineI=j;
　　 //主元所在行和当前处理行做行交换，右系数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;
　　 }
　　 if(*(A+i*num+i)==0) continue; //如果当前主元为0，本次循环结束
　　 //将A变为上三角矩阵，同样b也随之变换
　　 for(j=i+1;j<num;j++)
　　 {
　　 double temp=-*(A+j*num+i)/(*(A+i*num+i));
　　 for(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，说明有唯一解
　　*/
　　double determinantA=1;
　　for(i=0;i<num;i++)
　　 determinantA*=*(A+i*num+i);
　　if(determinantA==0)
　　{
　　 cout<<endl<<endl<<"通过计算,矩阵A的行列式为|A|=0，即A没有唯一解。\n程序退出..."<<endl<<endl;
　　 exit(0);
　　}
　　/* 从最后一行开始，回代求解x的初始向量值 */
　　for(i=num-1;i>=0;i--)
　　{
　　 for(j=num-1;j>i;j--)
　　 *(b+i)-=*(A+i*num+j)*(*(x+j));
　　 *(x+i)=*(b+i)/(*(A+i*num+i));
　　}
　　}
　　/*
　　--------------------------------------
　　利用雅可比迭代公式求解x的精确值
　　迭代公式为：xk=Bx+f
　　--------------------------------------
　　*/
　　void Jacobi(double* x,double* xk,double* B,double* f,int num,int time)
　　{
　　int t=1,i,j;
　　while(t<=time)
　　{
　　 for(i=0;i<num;i++)
　　 {
　　 double temp=0;
　　 for(j=0;j<num;j++)
　　 temp+=*(B+i*num+j)*(*(x+j));
　　 *(xk+i)=temp+*(f+i);
　　 }
　　 //将xk赋值给x，准备下一次迭代
　　 for(i=0;i<num;i++)
　　 *(x+i)=*(xk+i);
　　 t++;
　　}
　　}
　　//程序编写结束

复制代码

		自动登录	找回密码
密码			立即注册