重点是方程组类Linequ的方法 Solve(),解法一是清华教材《C++语言程序设计》郑莉等编的算法,好像有错(反正测试计算出来试错的main()为测试码),且不怎么明白!于是自己写了个,自己命名为逐行消去法!请大家看看!呵呵写得不好请帮忙改进!
// Linqu.h
//Begin of file Linequ.h
#ifndef LINQU_H
#define LINQU_H
#include <iostream.h>
#include <math.h>
class Matrix //The definition of basic class
{
public: //Interface for outside
Matrix(int dim=2); //Constructor
~Matrix(); //Destructor
void setMatrix(double* rmatr); //Initialize Matrix
void printM(); //Display Matrix
protected:
int index; //Dimation of Matrix
double* MatrixA;
};
class Linequ:public Matrix //Definition of public externs class
{
public:
Linequ(int dim=2);
~Linequ();
void setLinequ(double* a,double* b); //Initialize the equation
void printL(); //Display the equation
int Solve(); //Solution of gaosi equation
void showX(); //Display the result
private:
double* sum;
double* solu; //result of equation
};
#endif
//End of file Linequ.h
// Linequ.cpp
//Begin of file Linequ.cpp
#include "Linqu.h"
//Realization of basic class
void Matrix::setMatrix(double* rmatr) //set Matrix
{
for(int i=0;i<index*index;i++)
*(MatrixA+i)=rmatr[i];
}
Matrix::Matrix(int dim) //Constructor
{
index=dim; //Give protected member value
MatrixA=new double[index*index]; //Distribution of dynamic storage
}
Matrix::~Matrix()
{
delete[] MatrixA; //Release storage
}
void Matrix::printM() //Display members
{
cout<<"The Matrix is:"<<endl;
for(int i=0;i<index;i++)
{
for(int j;j<index;j++)
cout<<*(MatrixA+i*index+j)<<" ";
cout<<endl;
}
}
//Realization of derive class
Linequ::Linequ(int dim):Matrix(dim) //Constructor of class Linequ
{
sum=new double[dim]; //Dynamic destribution
solu=new double[dim];
}
Linequ::~Linequ()
{
delete[] sum;
delete[] solu;
}
void Linequ::setLinequ(double* a,double* b) //Set equation
{
setMatrix(a); //Call basic function
for(int i=0;i<index;i++)
sum[i]=b[i];
}
void Linequ::printL() //Display equation
{
cout<<"The Line equation is:"<<endl;
for(int i=0;i<index;i++)
{
for(int j=0;j<index;j++)
cout<<*(MatrixA+i*index+j)<<" ";
cout<<" "<<sum[i]<<endl;
}
}
void Linequ::showX() //display result
{
cout<<"The Result is:"<<endl;
for(int i=0;i<index;i++)
cout<<"X["<<i<<"]="<<solu[i]<<endl;
}
//解法 一 :全选主元高斯消去法
/*int Linequ::Solve() //Solution of Line equation
{
int *js,l,k,i,j,is,p,q;
double d,t;
js=new int[index];
l=1;
for(k=0;k<=index-2;k++) //The process of expunction
{
d=0.0;
js[k] = k;
is = k;
for(i=k;i<=index-1;i++)
for(j=k;j<=index-1;j++)
{
t=fabs(MatrixA[i*index+j]);
if(t>d){d=t;js[k]=j;is=i;}
}
if(d+1.0==1.0) l=0;
else
{
if(js[k]!=k)
for(i=0;i<=index-1;i++)
{
p=i*index+k; q=i*index+js[k];
t=MatrixA[p];
MatrixA[p]=MatrixA[q];
MatrixA[q]=t;
}
if(is!=k)
{
for(j=k;j<=index-1;j++)
{
p=k*index+j; q=is*index+j;
t=MatrixA[p];
MatrixA[p]=MatrixA[q];
MatrixA[q]=t;
}
t=sum[k]; sum[k]=sum[is]; sum[is]=t;
}
}
if(l==0)
{
delete[] js; cout<<"fail"<<endl;
return 0;
}
d=MatrixA[k*index+k];
for(j=k+1;j<=index-1;j++) //****************
{
p=k*index+j;
if(d!=0)MatrixA[p]=MatrixA[p]/d;
}
sum[k]=sum[k]/d;
for(i=k+1;i<=index-1;i++) //****************
{
for(j=k+1;j<=index-1;j++)
{
p=i*index+j;
MatrixA[p]=MatrixA[p]-MatrixA[i*index+k]*MatrixA[k*index+j];
}
sum[i]=sum[i]-MatrixA[i*index+k]*sum[k];
}
}
d=MatrixA[(index-1)*index+index-1];
if(fabs(d)+1.0==1.0)
{
delete[] js;
cout<<"fail"<<endl;
return 0;
}
solu[index-1]=sum[index-1]/d;
for(i=index-2;i>=0;i--)
{
t=0.0;
for(j=i+1;j<=index-1;j++)
t=t+MatrixA[i*index+j]*solu[j];
solu[i]=solu[i]-t;
}
js[index-1]=index-1;
for(k=index-1;k>=0;k--)
if(js[k]!=k)
{
t=solu[k]; solu[k]=solu[js[k]]; solu[js[k]]=t;
}
delete[] js;
return 1;
}*/
//解法 二 逐行消去法 自创
int Linequ::Solve()
{
int *js,k,i,j,l,flag,p,q;
double d,t;
js=new int[index];
for(i=0;i<index;i++)
js[i]=0;
for(i=0;i<index;i++) //纪录系数矩阵各行行首连零个数
for(j=0;j<index;j++)
{
if(MatrixA[i*index+j]==0) js[i]++;
else break;
}
do{
flag=0; //设置循环和判断标志
for(i=0;i<index-1;i++)
{
if(js[i] > js[i+1]) //如果前一行的零个数比后一行多,便进行交换
{
for(j=0;j<index;j++)
{ //整行交换
p=i*index+j; q=(i+1)*index+j;
t=MatrixA[p]; MatrixA[p]=MatrixA[q]; MatrixA[q]=t;
}
t=sum[i]; sum[i]=sum[i+1]; sum[i+1]=t;//对应和向量交换
l = js[i]; js[i] = js[i+1]; js[i+1] = l;//零计数 交换
}
}
for(i=0;i<index;i++) //各行元素分别除以第一非零元素
for(j=0;j<index;j++)
if(MatrixA[i*index+j]!=0)
{
d=MatrixA[i*index+j]; //找到第一个非零元素
for(k=j;k<index;k++)
MatrixA[i*index+k]=MatrixA[i*index+k]/d;
sum[i]=sum[i]/d;
break; //跳出该行
}
for(i=0;i<index-1;i++)
if(js[i]==js[i+1]) //若存在非零元素个数相等的行(肯定是相零行)flag=1
{
flag=1;
break;
}
if(flag == 1) //有非零元素个数相等的行便进行相减
{
for(i=index-1;i>0;i--)
if(js[i] == js[i-1])
{
for(j=0;j<index;j++)
MatrixA[i*index+j] -= MatrixA[(i-1)*index+j];
sum[i] -= sum[i-1];
js[i] = 0;
for(j=0;j<index;j++)
if(MatrixA[i*index+j] == 0) js[i]++; //重新纪录行零个数
else break;
}
}
else if(js[index-1] >= index-1) //若不存在零数相等的行,
flag = 2; //那么看最后一行是否只有一个非零或全零,
} //是这样就结束while循环
while(flag!=2);
if(js[index-1]==index) //若最后一行全零,则无解或有无穷组解
{
cout<<"There infinit solutions or no solution!"<<endl;
delete[] js;
return 0;
}
else{ //若最后一行只有一个非零,则有唯一解
solu[index-1]=sum[index-1]; //因为最后一行一列元素为1,解等于对应的和
for(i=index-2;i>=0;i--)
{
solu[i]=sum[i]; //重倒数第二行开始逐行回代
for(j=i+1;j<index;j++)
solu[i] -= MatrixA[i*index+j]*solu[j];
}
delete[] js;
return 1;
}
}
//End of file Linequ.cpp
/*************************************************************************
/ test.cpp
**************************************************************************/
#include "Linqu.h"
void main()
{
double a[]=
{0.2368,0.2471,0.2568,1.2671,
0.1968,0.2071,1.2168,0.2271,
0.1581,1.1675,0.1768,0.1871,
1.1161,0.1254,0.1397,0.1490};
double b[]={1.8471,1.7471,1.6471,1.5471};
Linequ equl(4);
equl.setLinequ(a,b);
equl.printL();
if(equl.Solve())
equl.showX();
else
cout<<"Fail"<<endl;
int t;
cin>>t;
}
/////////////////////////////////////////////////////////////////////////