goodtutorialonc++templateforbeginners.txt

给c++初学者一个用于科学计算的好例子

  l[i]=a[i]/u[i-1];
  u[i]=b[i]-(l[i])*(c[i-1]);
 }
 
 DOUBLE* y=new DOUBLE[nRow];
 y[0]=pb[0];
 
 for(i=1;i<=nRow-1;i++)
  y[i]=pb[i]-(l[i])*(y[i-1]);
 
 //DOUBLE* x=new DOUBLE[nRow];
 
 pb[nRow-1]=y[nRow-1]/u[nRow-1];
 
 
 for(i=nRow-2;i>=0;i--)
  pb[i]=(y[i]-(c[i])*(pb[i+1]))/u[i];
 
 delete [] a;
 delete [] b;
 delete [] c;
 delete [] l;
 delete [] u;
 delete [] y;
}

BOOL CMatrix::ReadMatrix(char c)
{
 //从文件以文本形式读入数据到数组；
 
 
 //检查隔离符是否正确
 if(c!=' ' && c!=',' && c!='$' && c!=' ')
 {
  AfxMessageBox("数据间隔符只能用',','$',空格和Tab键!");
  exit(0);
 }
 
 int i,j;
 typedef CArray<DOUBLE,DOUBLE> DblArray;
 DblArray p;
 
 static char BASED_CODE szFilter[] = "Ascii Files(*.MAsc)|*.MAsc|";
 //static char BASED_CODE szFilter[] = "Chart Files (*.xlc)|*.xlc|Worksheet Files (*.xls)|*.xls|Data Files (*.xlc;*.xls)|*.xlc; *.xls|All Files (*.*)|*.*||";
 
 CFileDialog dlg(true,"MAsc",NULL,OFN_HIDEREADONLY | OFN_OVERWRITEPROMPT,szFilter);
 if(dlg.DoModal()==IDOK)
 {
  CFile file(dlg.GetPathName(),CFile::modeRead);
  CArchive ar(&file,CArchive::load);
        
  CString str;
  CStringArray strA;
  ar.ReadString(str);
  while(str!="")
  {
   strA.Add((LPTSTR)(LPCTSTR)str);
   ar.ReadString(str);
  }
  
  for(i=0;i<strA.GetSize();i++)
  {
   CString str=strA[i];
   CString strL,strR;
   int k,l;
   
   str.TrimLeft();
   str.TrimRight();
   
   l=str.GetLength();
   k=str.Find(c);
   
   while(k!=-1)
   {
    strL=str.Left(k);
    strR=str.Right(l-k-1);
    
    str=strR;
    str.TrimLeft();
    str.TrimRight();
    
    //strL.TrimLeft();//
    strL.TrimRight();
    DOUBLE t=atof((LPTSTR)(LPCTSTR)strL);
    p.Add(t);
    
    strL=strR;
    strL.TrimLeft();
    strL.TrimRight();
    l=strL.GetLength();
    k=strL.Find(c);
   }
   DOUBLE t;
   strL.TrimLeft();
   strL.TrimRight();
   t=atof((LPTSTR)(LPCTSTR)strL);
   p.Add(t);
  }
  int m,n,L;
  m=strA.GetSize();
  L=p.GetSize();
  
  if((L%m)!=0)
  {
   AfxMessageBox("数据读取错误!");
   exit(0);
  }
  else
   n=L/m; 
  
  allocMemory(m,n);
  
  for(i=0;i<m;i++)
   for(j=0;j<n;j++)
   {
    Val[i][j]=p[n*i+j];
    DOUBLE tt=Val[i][j];
   }
   return true;//读取到数据；
 }
 else
  return false;//用户选择了取消；
}

 

DOUBLE CMatrix::Bspline3(DOUBLE xi, DOUBLE Vx[], DOUBLE Vy[], int n)
{
 //n为给定的插值节点总数目；
 //先判断计算点xi是否在插值节点范围内，并找出所属的区间位置；
 int i;
 if(xi>Vx[n-1] || xi<Vx[0])
 {
  AfxMessageBox("要计算的点不在已给插值点之间!");
  exit(0);
 }
 
 if(fabs(xi-Vx[0])<1e-20)
  return Vy[0];
 if((fabs(xi-Vx[n-1])<1e-20))
  return Vy[n-1];
 
 int j;
 for(i=0;i<=n-2;i++)
 {
  if(xi<Vx[i+1])
  {
   j=i;
   i=n-1;
  }
 }
 

 //n为给定的插值节点的数目；
 DOUBLE* h=new DOUBLE[n-1];
 DOUBLE* f=new DOUBLE[n-1];
 DOUBLE* u=new DOUBLE[n-1]; //u[0] is useless;
 DOUBLE* ff=new DOUBLE[n-1];//ff[0] is useless;
 DOUBLE* d=new DOUBLE[n-1];// d[0] is useless;
 DOUBLE* r=new DOUBLE[n-1];
/* -------------------------用GetParam()代替--------------------------------------
 for(i=0;i<=n-2;i++)
 {
  h[i]=Vx[i+1]-Vx[i];
  f[i]=(Vy[i+1]-Vy[i])/h[i];
 }
 for(i=1;i<=n-2;i++)
 {
  u[i]=h[i-1]/(h[i-1]+h[i]);
  r[i]=1-u[i];
  ff[i]=(f[i-1]-f[i])/(Vx[i-1]-Vx[i+1]);
  d[i]=6*ff[i];
  }
 ---------------------------------------------------------------------------*/
 GetParam(Vx,Vy,n,h,f,u,r,ff,d);
 
 
 CMatrix A(n-2,n-2);
 DOUBLE* b=new DOUBLE[n-2];
 DOUBLE* M=new DOUBLE[n];

 /*--------------------用GetM()代替；---------------------------------
 
   
  A(0,0)=2;A(0,1)=r[1];
  for(i=1;i<=n-3;i++)
  {
  A(i,i)=2;
  A(i,i+1)=r[i+1];
  A(i,i-1)=u[i+1];
  }
  A(n-2,n-2)=2;A(n-2,n-3)=u[n-1];
  
 -------------------------------------------------------------------*/
 
 
 GetA(u,r,d,A);
 
 //自然边界条件；M0=Mn=0;
 
 
 b[0]=d[1];
 for(i=1;i<=n-4;i++)
  b[i]=d[i+1];
 b[n-3]=d[n-2];
 
 A.chase_EQ(b);//用追赶法求解代数方程组；
 
 M[0]=0;
 for(i=1;i<=n-2;i++)
  M[i]=b[i-1];
 M[n-1]=0;
 

 ////////////////////////////////

 DOUBLE t1,t2,t3,t33,t44,t4;
 t1=(Vx[j+1]-xi)*(Vx[j+1]-xi)*(Vx[j+1]-xi);
 t1=t1*(M[j])/(6*h[j]);
 
 t2=(xi-Vx[j])*(xi-Vx[j])*(xi-Vx[j]);
 t2=t2*(M[j+1])/(6*h[j]);
 
 t3=(Vx[j+1]-xi)/h[j];
 t33=Vy[j]-(M[j])*(h[j])*(h[j])/6;
 t3=t3*t33;
 
 t4=(xi-Vx[j])/(h[j]);
 t44=Vy[j+1]-(M[j+1])*(h[j])*(h[j])/6;
 t4=t4*t44;
 
 delete [] h;
 delete [] f;
 delete [] u;
 delete [] ff;
 delete [] d;
 delete [] r;
 delete [] b;//
 delete [] M;
 double t5=t1+t2+t3+t4;
 return t5; 
}

DOUBLE* CMatrix::Bspline3(DOUBLE xi[], int m, DOUBLE Vx[], DOUBLE Vy[], int n)
{
 DOUBLE* s=new DOUBLE[m];
 double t;
 int i;
 for(i=0;i<=m-1;i++)
 {
   t=Bspline3(xi[i],Vx,Vy,n);
   s[i]=t;
 }
 return s;
}


void CMatrix::allocMemory(int m, int n)
{
 
 //先释放Val已经分配的内存，然后按照制定的大小重新分配内存； 
 /*---------------------------------------------------------------------
 注意：实际上，不能够根据Val是否为NULL来判断Val是否已经被分配过内存；
 比如：假如CMatrix类包含一个系统提供的默认的拷贝构造函数
 (即不进行任何初始化的操作的构造函数)，那么，Val的值将是一个随机值，
 如果对象是用此默认构造函数生成，则该对象的Val!=NULL，但是此时，Val
 被没有被分配过内存；
 
   至于在本程序中，在CMatrix类的所有的构造函数中都对Val进行了初始化，
   故可以根据Val是否为空来判断Val是否已经分配过内存；
 ----------------------------------------------------------------------*/
 //先调用freeMemory()函数；
 //freeMemory()函数会根据Val是否为空来检查Val是否已经被分配了内存；
 
 freeMemory();
 nRow=m;
 nCol=n;
 Val=new DOUBLE* [nRow];
 for(int i=0;i<nRow;i++)
  Val[i]=new DOUBLE [nCol];
}

BOOL CMatrix::ReadMatrix(CString strFilePath, char c)
{
 
 //从文件以文本形式读入数据到数组；
 
 
 //检查隔离符是否正确
 if(c!=' ' && c!=',' && c!='$' && c!=' ')
 {
  AfxMessageBox("数据间隔符只能用',','$',空格和Tab键!");
  exit(0);
 }
 
 int i,j;
 typedef CArray<DOUBLE,DOUBLE> DblArray;
 DblArray p;
 
 
 CFile file(strFilePath,CFile::modeRead);
 CArchive ar(&file,CArchive::load);
 
 CString str;
 CStringArray strA;
 if((ar.ReadString(str))==false)
 {
  AfxMessageBox("文件读取出错!检查文件名和文件路径!");
  exit(0);
 }
 while(str!="")
 {
  strA.Add((LPTSTR)(LPCTSTR)str);
  ar.ReadString(str);
 }
 
 for(i=0;i<strA.GetSize();i++)
 {
  CString str=strA[i];
  CString strL,strR;
  int k,l;
  
  str.TrimLeft();
  str.TrimRight();
  
  l=str.GetLength();
  k=str.Find(c);
  
  while(k!=-1)
  {
   strL=str.Left(k);
   strR=str.Right(l-k-1);
   
   str=strR;
   str.TrimLeft();
   str.TrimRight();
   
   //strL.TrimLeft();//
   strL.TrimRight();
   DOUBLE t=atof((LPTSTR)(LPCTSTR)strL);
   p.Add(t);
   
   strL=strR;
   strL.TrimLeft();
   strL.TrimRight();
   l=strL.GetLength();
   k=strL.Find(c);
  }
  DOUBLE t;
  strL.TrimLeft();
  strL.TrimRight();
  t=atof((LPTSTR)(LPCTSTR)strL);
  p.Add(t);
 }
 int m,n,L;
 m=strA.GetSize();
 L=p.GetSize();
 
 if((L%m)!=0)
 {
  AfxMessageBox("数据读取错误!");
  exit(0);
 }
 else
  n=L/m; 
 
 allocMemory(m,n);
 
 for(i=0;i<m;i++)
  for(j=0;j<n;j++)
  {
   Val[i][j]=p[n*i+j];
   DOUBLE tt=Val[i][j];
  }
  return true;//读取到数据；
}

void CMatrix::GSElimin_EQ(DOUBLE *pb,DOUBLE TOL)
{
/*-------------------------------------------------------
高斯消去法解线性代数方程组；
Ax=b; pb为数组b;;
A:方阵；
x为代求的解，是与b同维数的数组；
返回值保存在参数pb中；
 --------------------------------------------------------*/
 if(nRow!=nCol)
 {
  AfxMessageBox("高斯消元法只适用于方阵!");
  exit(0);
 }
 
 /*---------------------------------开始消元-----------------------*/
 int i;
 int iM;
 for(i=0;i<=nRow-2;i++)
 {
  iM=findMainElem(i);
  if(iM==i && fabs(Val[i][i]-0)<TOL)
  {
   AfxMessageBox("矩阵为奇异矩阵，无法求得唯一解!");
   exit(0);
  }
  else if(iM!=i)
  {
   exchangeRow(i,iM);
   DOUBLE bt;
   bt=pb[i];
   pb[i]=pb[iM];
   pb[iM]=bt;
  }
  
  int m,n;
  DOUBLE l;
  for(m=i+1;m<=nRow-1;m++)
  {
   l=-Val[m][i]/Val[i][i];
   Val[m][i]=0;
   for(n=i+1;n<=nCol-1;n++)
   {
    Val[m][n]=Val[m][n]+l*Val[i][n];
   }
   pb[m]=pb[m]+l*pb[i];
  } 
 }
 
 if(fabs(Val[nRow-1][nCol-1])<TOL)
 {
  AfxMessageBox("矩阵奇异，无法求得唯一解 !");
  exit(0);
 }
 /*-------------------------消元完毕----------------------*/
 
 /********************开始求解********************************/
 
 DOUBLE* x=new DOUBLE[nRow];
 
 x[nRow-1]=pb[nRow-1]/Val[nRow-1][nRow-1];
 
 int kk;
 for(kk=nRow-2;kk>=0;kk--)
 {
  DOUBLE sum=0;
  for(int jj=kk+1;jj<=nRow-1;jj++)
  {
   sum=sum+Val[kk][jj]*x[jj];
  }
  sum=pb[kk]-sum;
  x[kk]=sum/Val[kk][kk];
 }
 /********************求解完毕********************************/
 
 for(i=0;i<=nRow-1;i++)
  pb[i]=x[i];
 delete [] x;
}


void CMatrix::exchangeRow(int i, int j)
{
 //仅仅用于高斯消元；
 int k;
 DOUBLE t;
 for(k=i;k<=nCol-1;k++)
 {
  t=Val[i][k];
  Val[i][k]=Val[j][k];
  Val[j][k]=t;
 }
}

int CMatrix::findMainElem(int row)
{
 //仅仅用于高斯消元；
 
 int i;
 DOUBLE max;
 int n;
 max=Val[row][row];
 n=row;
 for(i=row+1;i<=nRow-1;i++)
 {
  
  if(fabs(Val[i][row])>fabs(max))
  {
   max=Val[i][row];
   n=i;
  } 
 }
 return n;
}

BOOL CMatrix::SaveMatrix(BOOL bAscii)
{ 
 static char BASED_CODE szFilter[]="Ascii Files(*.MAsc)|*.MAsc|Binary Files(*.MBin)|*.MBin|";
 CFileDialog dlg(false,"MatrixFile.MAsc",NULL,OFN_HIDEREADONLY | OFN_OVERWRITEPROMPT,szFilter);
 if(dlg.DoModal()==IDOK)
 {
  CString str=dlg.GetPathName();
  CFile file(str,CFile::modeCreate|CFile::modeReadWrite);
  CArchive ar(&file,CArchive::store);
  
  int i,j;
  CString s;
  CString s1;
  if(bAscii)
  {
   for(i=0;i<=nRow-1;i++)
   {
    s="";
    for(j=0;j<=nCol-1;j++)
    {
     s1.Format("%.8lf  ",Val[i][j]);
     s=s+s1; 
    }
    ar.WriteString(s);
    ar.WriteString("\r\n");//回车换行符号；
   }   
   return true;
  }
  if(!bAscii)
  {
   Serialize(ar);
   return true;
  }
 }
 return false;
}

void CMatrix::GetSize(int &m, int &n)
{
 m=nRow;
 n=nCol;
}

BOOL CMatrix::SaveMatrix(CString strPath, BOOL bAscii)
{
 
 CFile file(strPath,CFile::modeCreate|CFile::modeReadWrite);
 CArchive ar(&file,CArchive::store);
 if(bAscii)
 {
  int i,j;
  CString s;
  CString s1;
  for(i=0;i<=nRow-1;i++)
  {
   s="";
   for(j=0;j<=nCol-1;j++)
   {
    s1.Format("%.8lf  ",Val[i][j]);
    s=s+s1; 
   }
   ar.WriteString(s);
   ar.WriteString("\r\n");//回车换行符号；
  }    
  return true;
 }
 else
 {
 
 Serialize(ar);
 return true;
 }
 return false;
}

int CMatrix::GetRowSize()
{
 return nRow;
}

int CMatrix::GetColSize()
{
 return nCol;
}

BOOL CMatrix::LoadMatrix(CString strPath)
{ 
  CFile file(strPath,CFile::modeRead);
  CArchive ar(&file,CArchive::load);

  Serialize(ar);
  return true;
}

BOOL CMatrix::LoadMatrix()
{
 static char BASED_CODE szFilter[] = "Ascii Files(*.MAsc)|*.MAsc|Binary Files(*.MBin)|*.MBin|";
    CFileDialog dlg(true,NULL,"2进制数据文件.MBin",OFN_HIDEREADONLY | OFN_OVERWRITEPROMPT,szFilter);
 
 if(dlg.DoModal()==IDOK)
 {
  CString str=dlg.GetPathName();
  CFile file(str,CFile::modeRead);
  CArchive ar(&file,CArchive::load);

  Serialize(ar);
  return true;
 }
 return false;
}


还有很多算法，比如矩阵求逆，矩阵求秩，矩阵分解，子矩阵的提取，高斯－赛德尔迭代法等暂时还没有添加进去，有空慢慢写。


Powered by BlogDriver 2.1