《测量程序设计三ppt课件》由会员分享，可在线阅读，更多相关《测量程序设计三ppt课件（67页珍藏版）》请在金锄头文库上搜索。
1、测量程序设计曾喆中国石油大学地球科学与技术学院第三讲 法方程式的解算及系数阵的求逆主要内容：高斯约化法高斯-约旦法变量循环重新编号法求逆平方根法求法方程式系数阵的逆改正的平方根法求法方程式系数阵的逆补充内容平差法方程式条件平差法方程间接平差法方程1.高斯约化法主要步骤1、消元2、回代求解实例解方程组第一步,将第一行乘-2加与第二行相加；第一行乘-1加到第三行,得到第二步,将上面第二行乘-2/3加到第三行,得到回代:解第三行得x3,将x3代入第二行得x2,将x2,x3代入第一行得x1,得到解x*=(2,1,-1)T设方程组为：其初始的矩阵形式为：如果，我们首先保留矩阵的第一行，并利用它来消去其余
2、三行中的第一列。即 消元过程第i行减去乘以第0行，得到其中，同理，如果保留第0行和第1行基础上消去第1列对角线下方元素。即：第i行减去乘以第1行，得到其中，同理，若，还可以进一步消元同样地，第3行减去乘以第2行，得到其中，经过这个消元过程，最开始的方程组就变为以上将原始的方程组化为上式的过程为消元过程。下面将从方程组的最下一开始上推出所有x。 回代过程如果，根据下式回代得到所有的x的解小结从上面可以看出整个高斯消去法解方程的方法主要分为两步1、消去过程从k=0开始，在第k步中，定义因子然后通过下式计算k+1步的方程系数通过n-1步，即k=n-2时就可以得到原方程组系数的上三角阵2、回代过程代码
3、/消元过程doublem;for(intk=0;kn-1;k+)tryfor(inti=k+1;in;i+)m=ai,k/ak,k;for(intj=k+1;j=0;i-)doublesum=0.0;for(j=i+1;jn;j+)sum=sum+aij*xj;xi=(ci-sum)/aii;2.高斯-约旦法高斯消去法在消元时始终消去对角线下方的元素，而高斯约当消去法则同时消去对角线上方和下方的元素。第一次消元第二次消元故方程组的解为高斯约当消去法的特点：(1)消元和回代同时进行；(2)乘除法的次数要比高斯消去法大，所以通常用于同时求解系数矩阵相同的多个方程组或求逆矩阵。求逆：则不为0则上面方
4、程组第一式除之,得到,交换与并将代入余下各式,得到:反复运行,可以将x全部放到左边,y放到右边,这样y的系数阵就是C的逆阵计算公式:即3.变量循环重新编号法求逆方程式逆关系这里就有通过第一式计算出x0，再代入余下的方程式可得到下式新系数通过如下计算可得对上式按下面规则重新编号经过n次这种变换，变量均恢复原来的状态按上法进行了k步，A为正定对称，就有1、对于任意k，A11和A22均为正定2、A12=-A21TA为对称正定矩阵，用“变量循环重新编号法”经过n次变换得到A-1，计算如下：主要代码for(intk=0;kn;k+)doublea00=a0;if(a00+1.0=1.0)deletea0
5、;returnfalse;for(inti=0;in;i+)doubleai0=ai*(i+1)/2;if(i=n-k-1)a0i=-ai0/a00;elsea0i=ai0/a00;for(intj=1;j=i;j+)a(i-1)*i/2+j-1=ai*(i+1)/2+j+ai0*a0j;for(i=1;i=k=j)这样，A的逆就为平方根法步骤1、对A做分解LLT2、求R=L-13、求A-1=RRT平方根法就是正定对称矩阵的Cholesky分解代码double sum;for (i=0;in;i+) for (j=i;j=0;k-) sum -= aik*ajk; if (i = j) if
6、(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.改进的平方根法平方根法求平方是很耗时的，对于正定对称阵来说有n阶对称正定矩阵 A的 分解公式：6.补充内容数值计算的误差C#编程的问题数值计算的误差问题1、机器的精度(截断误差)2、数值计算的舍入误差截断误差浮点数的表达用一个尾数（Mantissa ），一个基数（Base），一个指数（Exponent）以及一个表示正负的符号来表达实数。机器精度浮点数在32位机子上有两种精
7、度， float占32位，double占64位。IEEE标准754（1985） 符号位 阶码 尾数 长度float 1 8 23 32double1 11 52 64 例一书中代码（C+）P17，P18+0.000001例二（C#）如下C#代码：floata=0.65f;floatb=0.6f;floatc=a-b;此时c为多少？0.05？错误！此时c为0.0499999523！小数部分为1/2,1/4,1/8,这样的没有问题截断误差原因floatvar=5.2f;二进制表示：整数部分：101小数部分：0.2*2=0.4*2=0.8*2=1.6(0.6)*2=1.2(0.2)*2=0.4*2=
8、0.8*2=1.6(0.6)*2=1.2.00110011.float是32位，后面尾数长度最大23位。5.2的二进制表示就为：101.00110011001100110011一共23位。使用科学计数法表示，结果为：1.0100110011001100110011*22最高位为1，可以省略，在最后一位补一位：.01001100110011001100110*22010000001010011001100110011001101位|8位|-23位-|符号|阶码|尾数|二进制到浮点数101.0011001100110011001105.xxxxxxxxxxxxxxxxxxxxxxxxxxxxx0.
9、0011001100110011001100+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高斯消去法，第一步消去后，得到下一步消去乘数为2.5*103，继续采用高斯消去得到最后一行：假设5位的精度，做截断圆整后得到：这样就得到：回代得到：高斯消去法中，应该选择最大的元素，做行交换后到主元素。避免不必要的圆整误差带来的解得精度降低问题。线性方程组的形态及解的误差设x*是计算得到的解（计算机的精度导致的误差）A非
10、奇异，理论解一定精度的计算解与理论解的差异例假设在十进制精度为3位的机器上计算：换元后高斯消去，乘数为消去得计算解与理论解的差异残差值均是小于10-3，残差很小。将主元素换为较大的值可以保证的仅仅是残差很小。虽然采用了换元的高斯消去，但是由于截断圆整对于这样的方程仍然得到的不是正确的解将精度扩展到6位的十进制，有解得，什么原因？范数对于方程Ax=bM/m比值称为A的条件数由定义可得由上式可以看出，条件数可以大致度量出b的相对误差对于解x的放大的倍数改变b值，得到C#编程中的一些关键问题1、委托机制2、界面设计delegate double ProcessDelegate（double para
11、m1，double param2）；static double Multiply(double param1,double param2) return param1*param2;static double Divide(double param1,double param2) return param1/param2;staticvoidMain(stringargs)ProcessDelegateprocess;Console.WriteLine(“Enter2numbersseparatedwithacomma:”);stringinput=Console.ReadLine();int
12、commaPos=input.Indexof(,);doubleparam1=Convert.ToDouble(input.Substring(0,commaPos);doubleparam2=Convert.ToDouble(input.Substring(commaPos+1,input.Length-commaPos-1);Console.WriteLine(“EnterMtomultiplyorDtodivide:”);input=Console.ReadLine();If(input=“M”)process=newprocessDelegate(Multiply);Elseprocess=newprocessDelegate(Divide);Console.WriteLine(“Result:(0)”,process(param1,param2);Console.ReadKey();
