水准网平差c++代码

水准网平差c++代码
水准平差网
结果
#include<iostream>
#include<stdlib.h>
#include<iomanip>
#include<math.h>
#define max 50
Using namespace std;
Class CMatrix
{
public:
CMatrix(){row=0; column=0;}; // 默构造函数认认认认认
CMatrix(int i, int j){row=i;column=j;} // 构造函数一
CMatrix(constCMatrix& m); // 认制构造函数
~CMatrix(void){/*cout<<"认认认认认认认认认认认使用,矩所占空以放:"<<endl;*/} // 默析构函数认认认认认
CMatrix& operator=(constCMatrix& m); // 认认运算符
bool operator==(constCMatrix& m); // 比括运算符认认认认
bool operator!=(constCMatrix& m); // 比括运算符认认认认
CMatrix operator+(constCMatrix& m); // 加运算符
CMatrix operator-(constCMatrix& m); // 减运算符
CMatrix& operator+=(constCMatrix& m); // 自加运算符
CMatrix& operator-=(constCMatrix& m); // 自减运算符
CMatrix operator-();// 取数认认
CMatrix& operator*(constCMatrix& m); // 乘法运算符
void input(); //认认入矩
void outputMatrix(); // 认认认认出矩
void setValue(int row, int column, double value) { A[row-1][column-1] = value; }// 认置(i,j)的认
double getValue(int row, int column) const { return A[row-1][column-1]; }// 认认置行、列的
voidsetRow(constint i) { row = i; }
intgetRow() const { return row; }
voidsetColunm(constint i) { column = i; }
intgetColumn() const { return column; }
CMatrix& change(int i, int j);//交矩的行认认认认认
CMatrix& transpose(); // 矩置认认认
CMatrix& inverse(); // 矩求逆认认认
void find(int& f)const;// 判断矩是否可用于迭代求解认认认认认认认认认
认认认
friend void jocabi(constCMatrix& a) ; //迭代求解
void lzys(); //列主元素法求解
void solve(); //可逆性矩求解认认认认认认
void qxnh(); //曲合认认认
private:// 成量认认认
double A[max][max];
int row;// 行
int column;// 列
};
void CMatrix::input() //认入
{ cout<<"认认认认认认认始入矩:"<<endl;
int i, j;
double z;
for(i=0;i<row;i++)
{ cout<<"认认入第"<<i+1<<"行的:认认"<<endl;
for(j=0;j<column;j++)
{cin>>z;
A[i][j]=z;
}
}
cout<<endl;
};
CMatrix::CMatrix(constCMatrix& m) // 认制构造函数{ int i, j; for(i=0;i<m.row;i++)
for(j=0;j<m.column;j++)
this->A[i][j]=m.A[i][j];};
CMatrix&CMatrix::operator=(constCMatrix& m) // 认认运算符{ inti,j;
for(i=0;i<row;i++)
{
for(j=0;j<column;j++)
A[i][j]=m.A[i][j];
}
return *this;
};
boolCMatrix::operator ==(constCMatrix& m) // 比括运算符认认认认{ inti,j,k;
for(i=0;i<m.row;i++)
{ for(j=0;j<m.column;j++)
if(this->A[i][j]=m.A[i][j]) k=1;
else k=0;
}
if(k=1) return true;
else return false;
};
boolCMatrix::operator !=(constCMatrix& m) // 比括运算符认认认认{ inti,j,k;
for(i=0;i<m.row;i++)
{ for(j=0;j<m.column;j++)
if(this->A[i][j]=m.A[i][j]) k=1;
else k=0;
}
if(k=0) return true;
else return false;
};
CMatrixCMatrix::operator+(constCMatrix& m)// 加运算符{ inti,j; if((this->row==m.row)&&(this->column==m.column))
{ for(i=0;i<m.row;i++)
for(j=0;j<m.column;j++)
this->A[i][j]+=m.A[i][j];
}
else { cout<<"此两矩不能相加,:认认认认认认认认认认"<<endl; } return *this;
};
CMatrixCMatrix::operator-(constCMatrix& m)// 减运算符
{ inti,j;
if((this->row==m.row)&&(this->column==m.column))
{ for(i=0;i<m.row;i++)
for(j=0;j<m.column;j++)
this->A[i][j]-=m.A[i][j];
}
else { cout<<"此两矩不能相加,:认认认认认认认认认认"<<endl; }
return *this;
};
CMatrix&CMatrix::operator+=(constCMatrix& m) //自加运算符{ inti,j;
for(i=0;i<m.row;i++)
for(j=0;j<m.column;j++)
this->A[i][j]=2*m.A[i][j];
return *this;
};
CMatrix&CMatrix::operator-=(constCMatrix& m) //自减运算符{ inti,j;
for(i=0;i<m.row;i++)
for(j=0;j<m.column;j++)
this->A[i][j]=m.A[i][j]-m.A[i][j];
return *this;
};
voidCMatrix::find(int& f)const{ int i;
for(i=0;i<this->row;i++)
if(this->A[i][i]!=0) f=1;
else f=0;
};
CMatrixCMatrix::operator-() // 取数认认
{ inti,j;
for(i=0;i<this->row;i++)
for(j=0;j<this->column;j++)
this->A[i][j]=-this->A[i][j];
return *this;
};
CMatrix&CMatrix::operator*(constCMatrix& m) // 乘法运算符
{ inti,j,t;
CMatrix n;
if(this->column==m.row)
{
for(i=0;i<this->row;i++)
for(j=0;j<m.column;j++)
{double sum=0.0;
for(t=0;t<m.row;t++)
sum+=this->A[i][t]*m.A[t][j];
n.A[i][j]=sum;
}
}
else {cerr<<"此两矩不能相乘,:认认认认认认认认认认"<<endl; exit(1);} return n;
};
void CMatrix::outputMatrix()// 认认认认出矩
{ inti,j;
for(i=1;i<=row;i++)
{ for(j=1;j<=column;j++)
cout<<A[i-1][j-1]<<" ";
cout<<endl;
}
};
CMatrix&CMatrix::transpose()// 矩置认认认
{ inti,j;
CMatrix m(this->column,this->row);
for(i=0;i<this->row;i++)
for(j=0;j<this->column;j++)
{ m.A[j][i]=this->A[i][j];}
return m;
};
void jocabi(constCMatrix& a) //高斯迭代求解
{ int f=1;
a.find(f);
if(f==0) {cerr<<"认认认认认认认认认认认认矩不足迭代求解条件:"<<endl; exit(1); }
else
{
CMatrixx,w;
x.setColunm(1);
x.setRow(a.getColumn());
w.setColunm(1);
w.setRow(a.getColumn());int i;
double z;
for(i=1;i<=x.row;i++)
{ cout<<"认认认认认入等式右的第"<<i<<"个:认认"; cin>>z;
w.A[i-1][0]=z;
}
for(i=1;i<=x.row;i++)
{ cout<<"认认入X"<<i<<"的近似:认认";
cin>>z;
x.setValue(i, 1, z);
}
i=1;
while(i<=20)
{ int j, k;
for(j=1;j<=x.row;j++)
{ double sum=0.0;
for(k=1;k<j;k++)
{sum=sum-(a.A[j-1][k-1])*(x.A[k-1][0]);}
for(k=j+1;k<=x.row;k++)
{sum=sum-(a.A[j-1][k-1])*(x.A[k-1][0]);}
sum+=w.A[j-1][0];
x.A[j-1][0]=sum/a.A[j-1][j-1];
}
i++;
}
for(i=1;i<=x.row;i++)
cout<<"X"<<i<<" = "<<x.A[i-1][0]<<endl;
}
};
double jdz(double a) //认认认
{ if(a>0.0) return a;else return -a;
};
CMatrix&CMatrix::change(int i, int j)//交矩的行认认认认认{ int k;
double z;
for(k=1;k<=this->column;k++)
{z=A[i-1][k-1];
A[i-1][k-1]=A[j-1][k-1];
A[j-1][k-1]=z;
}
return *this;
};
void CMatrix::lzys() //列主元素法求解
{ CMatrixx,w;
x.setColunm(1);
x.setRow(getColumn());
w.setColunm(1);
w.setRow(getColumn());int i;
double z;
for(i=0;i<row;i++)
{for(int j=0;j<row; j++)
w.A[i][j]=0.0;
}
for(i=1;i<=x.row;i++)
x.setValue(i, 1, 0.0);
for(i=1;i<=x.row;i++)
{ cout<<"认认认认认入等式右的第"<<i<<"个:认认"; cin>>z;
w.A[i-1][0]=z;
}
i=0;
while(i<x.row-1)
{ intj,t,h=i;
for(j=i;j<x.row-1;j++)
{if(jdz(A[j][i])<jdz(A[j+1][i])) h=j+1; }
this->change(i+1,h+1);
w.change(i+1,1);
for(j=i+1;j<row;j++)
{ double k;
k=A[j][i]/A[i][i];
for(t=i;t<column;t++)
{
A[j][t]=A[j][t]-k*A[i][t];
}
w.A[j][0]=w.A[j][0]-k*w.A[i][0];
}
i++;
}
if(this->A[row-1][column-1]==0)
{cerr<<"此矩的方程有无解:认认认认认认认认认认认认"<<endl; exit(1);} i=column-1;
while(i>0)
{ int j, t;
for(j=i-1;j>=0;j--)
{ double k;
k=A[j][i]/A[i][i];
for(t=i;t>j;t--)
{ A[j][t]=A[j][t]-k*A[i][t]; }
w.A[j][0]=w.A[j][0]-k*w.A[i][0];
}
i--;
}
int j;
for(j=0;j<x.row;j++)
{
x.A[j][0]=w.A[j][0]/this->A[j][j];
}
for(i=1;i<=x.row;i++)
cout<<"X"<<i<<" = "<<x.A[i-1][0]<<endl;
};
CMatrix&CMatrix::inverse() //矩求逆认认认
{ if(this->row!=this->column) {cerr<<"认认认认认认认认认认矩不符合求逆条件:"<<endl; exit(1);}
else
{
CMatrix w;
w.setColunm(this->getRow());w.setRow(getColumn());
int i;
for(i=0;i<row;i++)
{
for(int j=0;j<row; j++)
{w.A[i][j]=0.0;}
w.A[i][i]=1;
}
i=0;
while(i<column-1)
{ intj,t,h=i;
for(j=i;j<row-1;j++)
{
if(jdz(A[j][i])<jdz(A[j+1][i]))
h=j+1;
}
this->change(i+1,h+1);
w.change(i+1,h+1);
for(j=i+1;j<row;j++)
{ double k;
k=A[j][i]/A[i][i];
for(t=i;t<column;t++)
{
A[j][t]=A[j][t]-k*A[i][t];
w.A[j][t]=w.A[j][t]-k*w.A[i][t];
}
}
i++;
}
if(this->A[row-1][column-1]==0)
{cerr<<"此矩求逆不成功,其所的方程有无解:认认认认认认认认认认认认认认认认认认认认"<<endl; exit(1);}
i=column-1;
while(i>0)
{ int j, t;
for(j=i-1;j>=0;j--)
{ double k;
k=A[j][i]/A[i][i];
for(t=i;t>j;t--)
{ A[j][t]=A[j][t]-k*A[i][t]; }
for(t=column-1;t>=0;t--)
{w.A[j][t]=w.A[j][t]-k*w.A[i][t];}
}
i--;
}
intj,k;
for(j=0;j<row;j++)
{for(k=0;k<w.column;k++)
{w.A[j][k]=w.A[j][k]/this->A[j][j];}
this->A[j][j]=this->A[j][j]/this->A[j][j];
}
return w;
}
};
void CMatrix::solve() //可逆性矩求解认认认认认认{
CMatrixx,w,c;
x.setColunm(1);
x.setRow(getColumn());
w.setColunm(1);
w.setRow(getColumn());
c.setColunm(getColumn());
c.setRow(getRow());
int i;
double z;
for(i=1;i<=x.row;i++)
{ cout<<"认认认认认入等式右的第"<<i<<"个:认认";
cin>>z;
w.A[i-1][0]=z;
}
c.operator=(this->inverse());x.operator=(c.operator *(w)); cout<<"求解果:认认认"<<endl;
for(i=1;i<=x.row;i++)
cout<<"X"<<i<<" = "<<x.A[i-1][0]<<endl;
};
void CMatrix::qxnh() //曲合认认认
{
cout<<"用矩曲合:认认认认认认"<<endl;
int i, j,k,t;
double x, y;
cout<<"认认认认认认认认入已知坐点的个数:";
cin>>i;
CMatrix a(i,1), g(i,i) , w(i,1),c(i,i),q(i,i),p(i,i);
for(j=1;j<=i;j++)
{ cout<<"认认入第"<<j<<"个已知点坐认(先横坐后坐认认认认认):"; cin>>x>>y;
w.A[j-1][0]=y;
g.A[j-1][0]=1.0;
for(k=1;k<i;k++)
{ g.setValue(j,k+1,1.0);
for(t=1;t<=k;t++)
g.A[j-1][k]*=x;
}
}
c=g.transpose();
q=c*g;
p=q.inverse();
w=c*w;
a=p*w;
cout<<"Y="<<a.A[0][0];
for(j=1;j<i;j++)
{ cout<<'+'<<'('<<a.A[j][0]<<')';
for(k=0;k<j;k++)
cout<<"*x";
}
cout<<endl;
cout<<" 1 .求取其它点坐认 0. 退出 "<<endl; while(1)
{ cout<<"认认认:";
cin>>k;
switch(k)
{ case 1:
{ double x, y,sum;
cout<<"认认认认入待求点的横坐，x,:";
cin>>x;
y=a.A[0][0];
for(j=1;j<i;j++)
{sum=a.A[j][0];
for(k=0;k<j;k++)
{sum*=x;}
y+=sum;
}
cout<<"认认认坐:y="<<y<<endl;break;
}
case 0:
{
cerr<<"认认使用:"<<endl; exit(1);
break;
}
default: cout<<"认认认认认认认认有,:"<<endl; break; }
}
};
structCElvDif
{
doublevalue;// 认认认
double weight;// 认重
longstartPoint;// 起始点号认认
longendPoint;// 认认认点号
};
//水准点的认认认认
structCLevelPoint
{// 高程平差=认高程认+高程改正数认认认认
long index;// 水准点号认认
double eleValue;// 高程认
double dv;// 高程改正数，初始化认认认认认认认认认0, boolisKnown;// 是否已知点认认认认
};
classCElevationNet
{
public:// 成函数认认认
CElevationNet(){numElvDif=0;numPoints=0;numKnPoint=0;};// 构造函数~CElevationNet(){};// 析构函数void input();
void output();
voidjsgc();
intgetgz(){return numElvDif;} //返回数目认认认认认
intgetzs(){return numPoints;} //返回点数认认认
intgetys(){return numKnPoint;} //返回已知点数
intgetws(){return numPoints-numKnPoint;} //返回未知点数
double geteleValue(int i){return this->lpVec[i-1].eleValue;} //认认得高程
long getindex(int i){return lpVec[i-1].index;} //返回高程点号认认void seteleValue(int i, double value){lpVec[i-1].eleValue+=value;} //修改高程认
void setdv(inti,double value){this->lpVec[i-1].dv=value;} //修改改正数
double getdv(int i){return lpVec[i-1].dv;} //返回改正数
friend void xishu(CMatrix& B, CMatrix&X,CElevationNet A); //求取系数矩和未知点高程矩认认认认认认认认认
friend void quanzhen(CMatrix& Q, CElevationNet A); //求取认认
friend void l_zhen( CMatrix& l ,CMatrix& L, CElevationNet A); //求取认认认L矩和认认l认
private:// 成量认认认
intnumElvDif;// 高差数认认
intnumPoints;// 控制网中点的数目认认认
intnumKnPoint;//控制网中已知点的数目
CElvDifedVec[max];// 认认认认数
CLevelPointlpVec[max]; // 高程数认认认
};
voidCElevationNet::input(){ inta,b,c,i;
doublez,l,g;
cout<<"认认认认认认认认认认认认认认认认认认认认认认认认认已知点和未知点行号,注意已知点的号小于未知点的号:"<<endl<<endl;
cout<<"认认认认认认认认认认入控制网中水准点数和已知点数:";
cin>>a>>b;
numPoints=a;
numKnPoint=b;
for(i=0;i<numPoints;i++)
{ lpVec[i].index=i+1;
lpVec[i].eleValue=0.0;
lpVec[i].dv=0.0;
lpVec[i].isKnown=0;
}
for(i=0;i<numKnPoint;i++)
{ cout<<"认认入第"<<i+1<<"个已知点高程:认认";
cin>>z;
lpVec[i].eleValue=z;
lpVec[i].isKnown=1;
}
cout<<"认认认认认认认认入个数:";
cin>>c;
numElvDif=c;
for(i=0;i<numElvDif;i++)
{ cout<<"认认入第"<<i+1<<"段段的高差认认认认认认认(认位米)、度认认(认位千米)、起始点号和点号:认认认认认认认认"<<endl;
cin>>g>>l>>a>>b;
edVec[i].value=g;
edVec[i].weight=l;
edVec[i].startPoint=a;
edVec[i].endPoint=b;
}
};
voidCElevationNet::output()
{
int i;
for(i=0;i<this->numPoints;i++)
cout<<"认认号"<<this->lpVec[i].index<<"的点的高程:认认认"<<this-
>lpVec[i].eleValue<<endl;
for(i=0;i<this->numElvDif;i++)
cout<<"认认段"<<i+1<<"的平差后的:认认认"<<this->edVec[i].value<<endl;
};
void CElevationNet::jsgc()
{ inti,j;
for(i=0;i<this->numElvDif;i++)
{ if((lpVec[edVec[i].startPoint-
1].eleValue!=0.0)&&(lpVec[edVec[i].endPoint-1].eleValue==0.0)) {lpVec[edVec[i].endPoint-1].eleValue=(lpVec[edVec[i].startPoint-
1].eleValue)+(edVec[i].value);}
else if(lpVec[edVec[i].startPoint-
1].eleValue==0.0&&lpVec[edVec[i].endPoint-1].eleValue!=0.0) {lpVec[edVec[i].startPoint-1].eleValue=(lpVec[edVec[i].endPoint-1].eleValue)-(edVec[i].value);}
}
for(j=numKnPoint;j<this->numPoints;j++)
{
if(this->lpVec[j].eleValue==0.0)
{
for(i=0;i<this->numElvDif;i++)
{ if(this->edVec[i].startPoint==lpVec[j].index&&this->lpVec[this->edVec[i].endPoint-1].eleValue!
=0.0)
this->lpVec[j].eleValue=this->lpVec[this->edVec[i].endPoint-
1].eleValue-this->edVec[i].value;
else if(this->edVec[i].endPoint==lpVec[j].index&&this->lpVec[this->edVec[i].startPoint-1].eleValue!
=0.0)
{this->lpVec[j].eleValue=this->lpVec[this->edVec[i].startPoint-1].eleValue+this->edVec[i].value;}
}
}
//cout<<"第"<<this->lpVec[j].index<<"个点的高程近似:认认"<<this->lpVec[j].eleValue<<endl;
}
};
voidxishu(CMatrix& B, CMatrix&X,CElevationNet A)
{
inti,j;
for(i=1;i<=B.getRow();i++)
for(j=1;j<=B.getColumn();j++)
{ B.setValue(i,j,0);}
for(i=0;i<A.numElvDif;i++)
{ if(A.lpVec[A.edVec[i].startPoint-
1].isKnown==1&&A.lpVec[A.edVec[i].endPoint-1].isKnown==0) {B.setValue(i+1,A.edVec[i].endPoint-A.numKnPoint,1);}
else if(A.lpVec[A.edVec[i].startPoint-
1].isKnown==0&&A.lpVec[A.edVec[i].endPoint-1].isKnown==1) { B.setValue(i+1,A.edVec[i].startPoint-A.numKnPoint,-1);} else if(A.lpVec[A.edVec[i].startPoint-
1].isKnown==0&&A.lpVec[A.edVec[i].endPoint-1].isKnown==0) { B.setValue(i+1,A.edVec[i].endPoint-A.numKnPoint,1);
B.setValue(i+1,A.edVec[i].startPoint-A.numKnPoint,-1);
}
}
for(i=0;i<X.getRow();i++)
X.setValue(i+1,1,A.lpVec[A.numKnPoint+i].eleValue);
};
void quanzhen(CMatrix& Q, CElevationNet A){
inti,j;
double sum=0.0;
for(i=0;i<A.numElvDif;i++)
sum+=A.edVec[i].weight;
sum=sum/A.numElvDif;
for(i=1;i<=Q.getRow();i++)
for(j=1;j<=Q.getColumn();j++)
{ Q.setValue(i,j,0);}
for(i=1;i<=Q.getRow();i++)
Q.setValue(i,i,sum/A.edVec[i-1].weight);};
void l_zhen( CMatrix& l ,CMatrix& L, CElevationNet A){
int i;
double m;
for(i=0;i<A.numElvDif;i++)
L.setValue(i+1,1,A.edVec[i].value);
for(i=0;i<l.getRow();i++)
{
m=(A.lpVec[A.edVec[i].startPoint-1].eleValue+A.edVec[i].value-
A.lpVec[A.edVec[i].endPoint-
1].eleValue)*1000;
l.setValue(i+1,1,m);
}
};
void main()
{
CElevationNet A;
A.input();
A.jsgc();
CMatrix
B(A.getgz(),A.getws()),Q(A.getgz(),A.getgz()),X(A.getws(),1),l(A.get gz(),1),b(A.getws(),A.getgz()),V(A.getgz(),1);
CMatrix
NBB(A.getws(),A.getws()),x(A.getws(),1),N(A.getws(),A.getws()),L(A.g etgz(),1),M(A.getws(),A.getgz()),R(A.getws(),A.getgz()), W(A.getws(),A.getgz());
CMatrixvz(1,A.getgz()),vp(1,A.getgz()),p(1,1);
xishu(B,X,A);
quanzhen(Q,A);
l_zhen(l,L,A);
b.operator=(B.transpose());
M.operator=(b.operator *(Q));
NBB.operator=(M.operator *(B));
N.operator=(NBB.inverse());
R.operator=(N.operator *(b));
W.operator=(R.operator *(Q));
x.operator=(W.operator *(l));
inti,j;
for(i=A.getys()+1;i<=A.getzs();i++)
{A.setdv(i,(x.getValue(i-A.getys(),1)/1000));
A.seteleValue(i,A.getdv(i));
cout<<"平差后号认认认"<<A.getindex(i)<<"的点高程:认认"<<A.geteleValue(i)<<"米"<<endl;
}
V.operator=(B.operator *(x));
V.operator-(l);
for(j=0;j<V.getRow();j++)
cout<<"第"<<j+1<<"段高差平差:认认认认认认认认认
"<<(L.getValue(j+1,1)+V.getValue(j+1,1)/1000)<<"米
"<<endl;vz.operator=(V.transpose());
vp=(vz.operator *(Q));
p=(vp.operator *(V));
cout<<"平差后位中差:认认认认认认认认
"<<"+_"<<sqrt(p.getValue(1,1)/A.getws())<<"毫米"<<endl;}。