   /*--------------------------------------------- 假设有如下方程组： Ax=b 用Jacobi迭代法求解方程组的解 方法：将A分裂为A=D-L-U，等价的迭代方程组为x=Bx+f。 有关算法的详细说明，请点击此地址查看。 ---------------------------------------------*/ #include <iostream.h> #include <stdlib.h> #include <math.h> /*楼竞网站www.LouJing.com 拥有该程序的版权,转载请保留该版权. 谢谢合作!*/ 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; } /*-------------------- 分配内存空间函数 www.LouJing.com --------------------*/ double* allocMem(int num) { double *head; if((head=new double[num])==NULL) { cout<<"内存空间分配失败，程序终止运行！"<<endl; exit(0); } return head; } /*--------------------------------------------- 计算Ax=b中x的初始向量值，采用高斯列主元素消去法 www.LouJing.com ---------------------------------------------*/ 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的精确值 www.LouJing.com 迭代公式为：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++; } } |
|
|
