测量程序设计三ppt课件

测量程序设计曾喆中国石油大学地球科学与技术学院第三讲 法方程式的解算及系数阵的求逆主要内容：高斯约化法高斯-约旦法变量循环重新编号法求逆平方根法求法方程式系数阵的逆改正的平方根法求法方程式系数阵的逆补充内容平差法方程式条件平差法方程间接平差法方程,1,0,0,()aarr rTaar nn rn nr rNK WNAQ A WALA,1,W0,bbtt tTTbbt tNxNB PB WB Pl1.高斯约化法主要步骤1、消元2、回代求解实例解方程组第一步,将第一行乘-2加与第二行相加；第一行乘-1加到第三行,得到1231231232312756493xxxxxxxxx 123232323134264xxxxxxx 第二步,将上面第二行乘-2/3加到第三行,得到回代:解第三行得x3,将x3代入第二行得x2,将x2,x3代入第一行得x1,得到解 x*=(2,1,-1)T12323323134202033xxxxxx 000102030101112131(0)202122232303132333aaaacaaaacAaaaacaaaac00001 1022033010011 1122133120021 1222233230031 13223333a xa xa xa xca xa xa xa xca xa xa xa xca xa xa xa xc设方程组为：其初始的矩阵形式为：如果 ，我们首先保留矩阵的第一行，并利用它来消去其余三行中的第一列。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)1 im000102030(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)ijijijiiiaamajiccm c同理，若 ，还可以进一步消元同样地，第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)ijijijiiiaamai jccmc00001 10220330(1)(1)(1)(1)1111221331(2)(2)(2)2222332(3)(3)3333a xa xa xa xca xa xa xcaxaxcaxc 回代过程(3)3(3)33(2)(2)3223(2)22(1)(1)(1)2311213(1)11001 102 203 3003()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)kkkkijijikkjkkkkiiikkaamai jknccmcikn2、回代过程()()1()()1()()(1)(2,.,0)iiiiiniiiijjj iiiiciacaxiaxinxin 代码/消元过程double m;for(int k=0;kn-1;k+)try for(int i=k+1;in;i+)m=ai,k/ak,k;for(int j=k+1;j=0;i-)double sum=0.0;for(j=i+1;jn;j+)sum=sum+aij*xj;xi=(ci-sum)/aii;()()(1)iiiiiciaxin1()()1()()(2,.,0)niiiijjj iiiicaxiaxin 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.1nnnnnnnbbbxxxaaaaaa0)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,求逆：则1BC111 11221221 122221 122nnnnnnnnnnyC xC xC xyC xC xC xyC xC xC x(1)(1)(1)11111221(1)(1)(1)22112222(1)(1)(1)1122nnnnnnnnnnxCyCxCxyCyCxCxyCyCxCx11C不为0则上面方程组第一式除之,得到 ,交换 与 并将 代入余下各式,得到:1x1x1y1x反复运行,可以将x全部放到左边,y放到右边,这样y的系数阵就是C的逆阵1,(,1,2,.,)kkkkkjkkkjikikkkijijikkjCCi jkCCCCCCi jnCCC C 计算公式:()(1)()()(1)()(1)()()(1)()1,(,1,2,.,)kkkkkkkkkkjkkkjkkkikikkkkkkijijikkjCCi jkCCCCCCi jnCCC C 即3.变量循环重新编号法求逆yA xxB y1BA方程式逆关系这里就有000001 10,11110011 11,1111,001,1 11,11nnnnnnnnnnya xa xaxya xa xaxyaxaxax0000010010,1001(1)()()nnxayaaxaax 1100001110 010011,110 0,100111,00001,11,0 010011,11,0 0,1001()()()()()()nnnnnnnnnnnnyaayaa aaxaa aaxyaayaaaaxaaaax通过第一式计算出x0，再代入余下的方程式可得到下式00000000000000001,(1,.,1)/,(1,.,1)(,1,.,1)jjiiijijijaaaaajnaaainaaa aai jn 000001 10,11110011 11,1111,001,1 1111nnnnnnnnnnxa ya xaxya ya xaxyayaxax新系数通过如下计算可得111 11,1110011,1 11111,00001 10,11000nnnnnnnnnnya xaxa yyaxaxayxa xaxa y012110132012110132knnnknnknnnkkkxxxxxxxxxxyyyyyyyyyy对上式按下面规则重新编号经过n次这种变换，变量均恢复原来的状态按上法进行了k步，A为正定对称，就有 1、对于任意k，A11和A22均为正定 2、A12=-A21T11122122AAnkAAknkk1,10,00,0,01,10,0,0,0,0,0,01,1,0,0,0,01/(1)/(1)/(1)/(1)nnjnjji jijiji jijaaaajnkaaajnkaa aajnkaaa aajnkA为对称正定矩阵，用“变量循环重新编号法”经过n次变换得到A-1，计算如下：主要代码for(int k=0;k n;k+)double a00=a0;if(a00+1.0=1.0)delete a0;return false;for(int i=0;in;i+)double ai0=ai*(i+1)/2;if(i=n-k-1)a0i=-ai0/a00;else a0i=ai0/a00;for(int j=1;j=i;j+)a(i-1)*i/2+j-1=ai*(i+1)/2+j+ai0*a0j;for(i=1;i=k=j)1110101()0(,01)1nik kjknik kiiiiiijiiik kjkjjkiik kjkjiii ijik kjkjl rijl rij ijlrrl rl rlrll r 这样，A的逆就为1TARR0iijikjkkar r平方根法步骤 1、对A做分解LLT 2、求R=L-1 3、求A-1=RRT平方根法就是正定对称矩阵的Cholesky分解代码double sum;for(i=0;in;i+)for(j=i;j=0;k-)sum-=aik*ajk;if(i=j)if(sum=0.0)throw(Cholesky failed);aii=sqrt(sum);else aji=sum/aii;for(i=0;in;i+)for(j=0;ji;j+)aji=0.;5.改进的平方根法TALDL平方根法求平方是很耗时的，对于正定对称阵来说有 nnnnnnaaaaaaaaaA212222111211令令 1111112121212121nnnnnllldddlll nnndddlll212121111n阶对称正定矩阵 A的 分解公式：，111)1(ad ，jjkjkkikijijdldlal)()111 。112)2ikkikiiidladTLLD)(6.补充内容数值计算的误差C#编程的问题数值计算的误差问题1、机器的精度(截断误差)2、数值计算的舍入误差截断误差浮点数的表达 用一个尾数（Mantissa），一个基数（Base），一个指数（Exponent）以及一个表示正负的符号来表达实数。机器精度 浮点数在32位机子上有两种精度，float占32位，double占64位。IEEE标准754（1985）符号位 阶码 尾数 长度float 1 8 23 32double1 11 52 64 例一书中代码（C+）P17，P18+0.000001例二（C#）如下C#代码：float a=0.65f;float b=0.6f;float c=a-b;此时c为多少？0.05？错误！此时c为0.0499999523！小数部分为1/2,1/4,1/8,这样的没有问题截断误差原因float var=5.2f;二进制表示：整数部分：101101小数部分：0.2*2=0.4*2=0.8*2=1.6(0.6)*2=1.2(0.2)*2=0.4*2=0.8*2=1.6(0.6)*2=1.2.0 0 1 0 0 1 1 0 0 1 1 1 0 0 1 1 .float是32位，后面尾数长度最大23位。5.2的二进制表示就为：101.00110011001100110011 101.00110011001100110011 一共23位。使用科学计数法表示，结果为：1.0100110011001100110011 1.0100110011001100110011*2 22 2最高位为1，可以省略，在最后一位补一位：.0100110011001100110011.01001100110011001100110 0*2 22 2 0 1000 0001 01001100110011001100110 0 1000 0001 01001100110011001100110 1位|8位|-23位-|符号|阶码|尾数|二进制到浮点数1 1 01.0 0110 0110 0110 0110 011 0 01.0 0110 0110 0110 0110 011 0 5.xxxxxxxxxxxxxxxxxxxxxxxxxxxxx0.0 0110 0110 0110 0110 011 0 0.0 0110 0110 0110 0110 011 0 0+0*2-1+0*2-2+1*2-3+1*2-4+0*2-5+0*2-6+1*2-7+.+0*2-215.1999998数值计算的舍入误差高斯消去中的主元选择线性方程组的形态及解的误差高斯消去中的主元选择方程其解为:(0,-1,1)T 1231070732.09963.9015156xxx高斯消去法，第一步消去后，得到 下一步消去乘数为2.5*103，继续采用高斯消去得到最后一行：12310707063.90102.0.001556xxx333(5(2.5 10)(6)(2.5(2.5 10)(6.001)x443(1.5000 105)(2.51.50025 10)x假设5位的精度，做截断圆整后得到：这样就得到：回代得到：30.9 9 9 9 3x4431.5005 101.5004 10 x20.001(6)(0.99993)6.0001x21.5x 10.35x 高斯消去法中，应该选择最大的元素，做行交换后到主元素。避免不必要的圆整误差带来的解得精度降低问题。线性方程组的形态及解的误差设x*是计算得到的解（计算机的精度导致的误差）A非奇异，理论解 一定精度的计算解与理论解的差异1xA b*exxrbAx例假设在十进制精度为3位的机器上计算：换元后高斯消去，乘数为消去得 120.7800.5630.2170.9130.6590.254xx0.7800.8(32540.963)139.120.9130.6590.25400.0010.001xx221*0.0011.00,0.0010.2540.6590.4430.9130.4431.00 xxxx 计算解与理论解的差异*0.0004600.000541rbAx残差值均是小于10-3，残差很小。将主元素换为较大的值可以保证的仅仅是残差很小。虽然采用了换元的高斯消去，但是由于截断圆整对于这样的方程仍然得到的不是正确的解将精度扩展到6位的十进制，有 解得，什么原因？120.9130000.6590000.25400000.0000010.000001xx121.000001.00000 xx范数对于方程Ax=bM/m比值称为A的条件数 1/1pnpipixxmaxminAxMxAxmx()MAm()A xxbbA xb,()bMxbm xxbAxb由定义可得由上式可以看出，条件数可以大致度量出b的相对误差对于解x的放大的倍数4.12.84.11,9.76.69.70Abx 1113.8,1bx4.119.700.340.97bx1.63,0.00072460.01,1.63xxxxbx改变b值，得到1.63()2249.40.0007246AC#编程中的一些关键问题1、委托机制2、界面设计delegate double ProcessDelegate（double param1，double param2）；static double Multiply(double param1,double param2)return param1*param2;static double Divide(double param1,double param2)return param1/param2;static void Main(string args)ProcessDelegate process;Console.WriteLine(“Enter 2 numbers separated with a comma:”);string input=Console.ReadLine();int commaPos=input.Indexof(,);double param1=Convert.ToDouble(input.Substring(0,commaPos);double param2 =Convert.ToDouble(input.Substring(commaPos+1,input.Length-commaPos-1);Console.WriteLine(“Enter M to multiply or D to divide:”);input=Console.ReadLine();If(input=“M”)process=new processDelegate(Multiply);Else process=new processDelegate(Divide);Console.WriteLine(“Result:(0)”,process(param1,param2);Console.ReadKey();