水准网平差结果#include<iostream.h>#include<stdlib.h>#include<iomanip.h>#include<math.h>#define max 50class CMatrix{public:CMatrix(){row=0; column=0;}; // 默认构造函数CMatrix(int i, int j){row=i;column=j;} // 构造函数一CMatrix(const CMatrix& m); // 复制构造函数~CMatrix(void){/*cout<<"谢谢使用，矩阵所占空间以释放！"<<endl;*/} // 默认析构函数CMatrix& operator=(const CMatrix& m); // 赋值运算符bool operator==(const CMatrix& m); // 比括较运算符bool operator!=(const CMatrix& m); // 比括较运算符CMatrix operator+(const CMatrix& m); // 加运算符CMatrix operator-(const CMatrix& m); // 减运算符CMatrix& operator+=(const CMatrix& m); // 自加运算符CMatrix& operator-=(const CMatrix& m); // 自减运算符CMatrix operator-();// 取负数CMatrix& operator*(const CMatrix& 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]; }// 设置行、列的值void setRow(const int i) { row = i; }int getRow() const { return row; }void setColunm(const int i) { column = i; }int getColumn() const { return column; }CMatrix& change(int i, int j);//交换矩阵的行CMatrix& transpose(); // 矩阵转置CMatrix& inverse(); // 矩阵求逆void find(int& f)const;// 判断该矩阵是否可用于迭代求解friend void jocabi(const CMatrix& 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(const CMatrix& 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=(const CMatrix& m) // 赋值运算符{int i,j;for(i=0;i<row;i++){for(j=0;j<column;j++)A[i][j]=m.A[i][j];}return *this;};bool CMatrix::operator ==(const CMatrix& m) // 比括较运算符{ int i,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;};bool CMatrix::operator !=(const CMatrix& m) // 比括较运算符{ int i,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;};CMatrix CMatrix::operator+(const CMatrix& m)// 加运算符{ int i,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-(const CMatrix& m)// 减运算符{ int i,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+=(const CMatrix& m) //自加运算符{ int i,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-=(const CMatrix& 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]-m.A[i][j];return *this;};void CMatrix::find(int& f)const{ int i;for(i=0;i<this->row;i++)if(this->A[i][i]!=0) f=1;else f=0;};CMatrix CMatrix::operator-() // 取负数{ int i,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*(const CMatrix& m)// 乘法运算符{ int i,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()// 输出该矩阵{ int i,j;for(i=1;i<=row;i++){ for(j=1;j<=column;j++)cout<<A[i-1][j-1]<<" ";cout<<endl;}};CMatrix& CMatrix::transpose()// 矩阵转置{ int i,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(const CMatrix& a) //高斯迭代求解{ int f=1;a.find(f);if(f==0) {cerr<<"该矩阵不满足迭代求解条件！"<<endl; exit(1); }else{CMatrix x,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() //列主元素法求解{ CMatrix x,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){ int j,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){ int j,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--;}int j,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() //可逆线性矩阵求解{CMatrix x,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;}}};struct CElvDif{double value; // 观测值double weight; // 权重long startPoint; // 起始点编号long endPoint; // 终点编号};//水准点类的设计struct CLevelPoint{ // 高程平差值=高程值+高程值改正数long index; // 水准点编号double eleValue; // 高程值double dv; // 高程值改正数（初始化为0）bool isKnown; // 是否为已知点};class CElevationNet{public: // 成员函数CElevationNet(){numElvDif=0;numPoints=0 ;numKnPoint=0;}; // 构造函数~CElevationNet(){}; // 析构函数 void input();void output();void jsgc();int getgz(){return numElvDif;} //返回观测值数目int getzs(){return numPoints;} //返回总点数int getys(){return numKnPoint;} //返回已知点数int getws(){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(int i,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: // 成员变量int numElvDif; // 高差总数int numPoints; // 控制网中点的总数目int numKnPoint; //控制网中已知点的数目CElvDif edVec[max]; // 观测值数组CLevelPoint lpVec[max]; // 高程值数组};void CElevationNet::input(){ int a,b,c,i;double z,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;}};void CElevationNet::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(){ int i,j;for(i=0;i<this->numElvDif;i++){ if((lpVec[edVec[i].startPoint-1].eleValu e!=0.0)&&(lpVec[edVec[i].endPoint-1].eleVal ue==0.0)){lpVec[edVec[i].endPoint-1].eleValue=(l pVec[edVec[i].startPoint-1].eleValue)+(edVe c[i].value);}elseif(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)-(edVe c[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].endPoi nt-1].eleValue!=0.0)this->lpVec[j].eleValue=this->lpVec[thi s->edVec[i].endPoint-1].eleValue-this->edVe c[i].value;elseif(this->edVec[i].endPoint==lpVec[j].index& &this->lpVec[this->edVec[i].startPoint-1].e leValue!=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;}};void xishu(CMatrix& B, CMatrix& X,CElevationNet A){int i,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].isKn own==1&&A.lpVec[A.edVec[i].endPoint-1].isKn own==0){B.setValue(i+1,A.edVec[i].endPoint-A.numKn Point,1);}elseif(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);}elseif(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) {int i,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].ele Value+A.edVec[i].value-A.lpVec[A.edVec[i].e ndPoint-1].eleValue)*1000;l.setValue(i+1,1,m);}};void main(){CElevationNet A;A.input();A.jsgc();CMatrixB(A.getgz(),A.getws()),Q(A.getgz(),A.getgz( )),X(A.getws(),1),l(A.getgz(),1),b(A.getws( ),A.getgz()),V(A.getgz(),1);CMatrixNBB(A.getws(),A.getws()),x(A.getws(),1) ,N(A.getws(),A.getws()),L(A.getgz(),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));int i,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)/100 0)<<"米"<<endl;vz.operator =(V.transpose());vp=(vz.operator *(Q));p=(vp.operator *(V));cout<<"平差后单位权中误差为："<<"+_"<<sqrt(p.getValue(1,1)/A.getws())<<"毫米"<<endl; }。