测量程序设计三

测量程序设计曾喆中国石油大学地球科学与技术学院第三讲法方程式的解算及系数阵的求逆主要内容：高斯约化法高斯-约旦法变量循环重新编号法求逆平方根法求法方程式系数阵的逆改正的平方根法求法方程式系数阵的逆补充内容平差法方程式条件平差法方程间接平差法方程,1,0,,,,0,()aarrrTaarnnrnnrrNKWNAQAWALA,1,,ˆW0,bbtttTTbbttNxNBPBWBPl1.高斯约化法主要步骤1、消元2、回代求解实例解方程组第一步,将第一行乘-2加与第二行相加；第一行乘-1加到第三行,得到1231231232312756493xxxxxxxxx123232323134264xxxxxxx第二步,将上面第二行乘-2/3加到第三行,得到回代:解第三行得x3,将x3代入第二行得x2,将x2,x3代入第一行得x1,得到解x*=(2,1,-1)T12323323134202033xxxxxx000102030101112131(0)202122232303132333aaaacaaaacAaaaacaaaac0000110220330100111122133120021122223323003113223333axaxaxaxcaxaxaxaxcaxaxaxaxcaxaxaxaxc设方程组为：其初始的矩阵形式为：如果，我们首先保留矩阵的第一行，并利用它来消去其余三行中的第一列。000a(0)101000(0)202000(0)303000amaamaama即(0)0000(1,2,3)iiamia消元过程(1)(0)(0)00(1)(0)(0)(0)00(1,2,3;1,2,3)ijijijiiiaamajiccmc(0)第i行减去乘以第0行，得到其中，(0)0im000102030(1)(1)(1)111213(1)(1)(1)(1)(1)2122232(1)(1)(1)(1)3132333aaaacaaacAaaacaaac(1)1同理，如果保留第0行和第1行基础上消去第1列对角线下方元素。即：(1)110a(1)(1)11(1)11iiama(1)(1)2121(1)11(1)(1)3131(1)11amaama第i行减去乘以第1行，得到其中，(1)1im000102030(1)(1)(1)(1)1112131(2)(2)(2)(2)22232(2)(2)(2)32333aaaacaaacAaacaac(2)(1)(1)(1)11(2)(1)(1)(1)11,(2,3;2,3)ijijijiiiaamajiccmc同理，若，还可以进一步消元同样地，第3行减去乘以第2行，得到(2)220a(2)(2)22(2)22,(3)iiamia(2)2im000102030(1)(1)(1)(1)1112131(3)(2)(2)(2)22232(3)(3)333aaaacaaacAaacac其中，经过这个消元过程，最开始的方程组就变为以上将原始的方程组化为上式的过程为消元过程。下面将从方程组的最下一开始上推出所有x。(3)(2)(2)(2)22(3)(2)(2)(2)22,(,3)ijijijiiiaamaijccmc0000110220330(1)(1)(1)(1)1111221331(2)(2)(2)2222332(3)(3)3333axaxaxaxcaxaxaxcaxaxcaxc回代过程(3)3(3)33(2)(2)3223(2)22(1)(1)(1)2311213(1)110011022033003()2()1()0cacaxacaxaxacaxaxaxaxxxx如果，根据下式回代得到所有的x的解(3)330a小结从上面可以看出整个高斯消去法解方程的方法主要分为两步1、消去过程从k=0开始，在第k步中，定义因子然后通过下式计算k+1步的方程系数通过n-1步，即k=n-2时就可以得到原方程组系数的上三角阵()0kkka()()(),(1,...,1)kkikikkkkamikna(1)()()()(1)()()(),(,1,...,1),(1,...,1)kkkkijijikkjkkkkiiikkaamaijknccmcikn2、回代过程()()1()()1()()(1)(2,...,0)iiiiiniiiijjjiiiiciacaxiaxinxin代码//消元过程doublem;for(intk=0;kn-1;k++){try{for(inti=k+1;in;i++){m=a[i,k]/a[k,k];for(intj=k+1;jn;j++)a[i,j]=a[i,j]-a[k,j]*m;c[i]=c[i]-c[k]*m;}}catch(DivideByZeroExceptione){Console.WriteLine(e.Message);}}()()(),(1,...,1)kkikikkkkamikna(1)()()(),(,1,...,1)kkkkijijikkjaamaijkn(1)()()(),(1,...,1)kkkkiiikkccmcikn/*回代过程*/x[n-1]=c[n-1]/a[n-1][n-1];for(i=n-2;i=0;i--){doublesum=0.0;for(j=i+1;jn;j++)sum=sum+a[i][j]*x[j];x[i]=(c[i]-sum)/a[i][i];}()()(1)iiiiiciaxin1()()1()()(2,...,0)niiiijjjiiiicaxiaxin2.高斯-约旦法高斯消去法在消元时始终消去对角线下方的元素，而高斯——约当消去法则同时消去对角线上方和下方的元素。)1()1(2)1(121)1()1(2)1(1)1(2)1(22)1(21)1(1)1(12)1(11............nnnnnnnnbbbxxxaaaaaaaaa)2()2(2)2(121)2()2(2)2(2)2(22)2(1)2(12...0......0...1nnnnnnnbbbxxxaaaaaa0)1(11a第一次消元0)2(22a第二次消元)3()3(2)3(121)3()3(3)3(3)3(33)3(2)3(1...00............10...01nnnnnnnnbbbxxxaaaaaa)()(2)(1211...000...100...01nnnnnbbbxxx故方程组的解为T)()(2)(1T21......nnnnnbbbxxx高斯——约当消去法的特点：(1)消元和回代同时进行；(2)乘除法的次数要比高斯消去法大，所以通常用于同时求解系数矩阵相同的多个方程组或求逆矩阵。ByxCxy,求逆：则1BC11111221221122221122nnnnnnnnnnyCxCxCxyCxCxCxyCxCxCx(1)(1)(1)11111221(1)(1)(1)22112222(1)(1)(1)1122nnnnnnnnnnxCyCxCxyCyCxCxyCyCxCx11C不为0则上面方程组第一式除之,得到,交换与并将代入余下各式,得到:1x1x1y1x反复运行,可以将x全部放到左边,y放到右边,这样y的系数阵就是C的逆阵'''''''1,(,1,2,...,)kkkkkjkkkjikikkkijijikkjCCijkCCCCCCijnCCCC计算公式:()(1)()()(1)()(1)()()(1)()1,(,1,2,...,)kkkkkkkkkkjkkkjkkkikikkkkkkijijikkjCCijkCCCCCCijnCCCC即3.变量循环重新编号法求逆yAxxBy1BA方程式逆关系这里就有00000110,1111001111,1111,001,111,11nnnnnnnnnnyaxaxaxyaxaxaxyaxaxax0000010010,1001(1)()()nnxayaaxaax1100001110010011,1100,100111,00001,11,0010011,11,00,1001()()()()()()nnnnnnnnnnnnyaayaaaaxaaaaxyaayaaaaxaaaax通过第一式计算出x0，再代入余下的方程式可得到下式'0000'0000'0000'00001,(1,...,1)/,(1,...,1)(,1,...,1)jjiiijijijaaaaajnaaainaaaaaijn'''00000110,11'''11001111,11'''11,001,11111nnnnnnnnnnxayaxaxyayaxaxyayaxax新系数通过如下计算可得'''11111,11100'''11,111111,00'''00110,11000nnnnnnnnnnyaxaxayyaxaxayxaxaxay012110132012110132knnnknnknnnkkkxxxxxxxxxxyyyyyyyyyy对上式按下面规则重新编号经过n次这种变换，变量均恢复原来的状态按上法进行了k步，A为正定对称，就有1、对于任意k，A11和A22均为正定2、A12=-A21T11122122AAnkAAknkk'1,10,00,0,0'1,10,0,0,0,0,0,0'1,1,0,0,0,01//(1)/(1)/(1)/(1)nnjnjjijijijijijaaaajnkaaajnkaaaajnkaaaaajnkA为对称正定矩阵，用“变量循环重新编号法”经过n次变换得到A-1，计算如下：主要代码for(intk=0;kn;k++){doublea00=a[0];if(a00+1.0==1.0){delete[]a0;returnfalse;}for(inti=0;in;i++){doubleai0=a[i*(i+1)/2];if(i=n-k-1)a0[i]=-ai0/a00;elsea0[i]=ai0/a00;for(intj=1;j=i;j++){a[(i-1)*i/2+j-1]=a[i*(i+1)/2+j]+ai0*a0[j];}for(i=1;in;i++){a[(n-1)*n/2+i-1]=a0[i];}