cpp1.cpp

简单的小程序 就是用jacobi迭代来算线性方程组的解 见笑见笑

#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++;
 }

}