《测量程序设计三ppt课件.ppt》由会员分享，可在线阅读，更多相关《测量程序设计三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步，
3、即k=n-2时就可以得到原方程组系数的上三角阵,2、回代过程,代码,/消元过程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;jn;j+)ai,j=ai,j-ak,j*m;ci=ci-ck*m;catch(DivideByZeroException e)Console.WriteLine(e.Message);,/*回代过程*/xn-1=cn-1/an-1n-1;for(i=n-2;i=0;i-)double sum=0.0;for(j=i+1;jn;j+)sum=sum+aij*xj
4、;xi=(ci-sum)/aii;,2.高斯-约旦法,高斯消去法在消元时始终消去对角线下方的元素，而高斯约当消去法则同时消去对角线上方和下方的元素。,第一次消元,第二次消元,故方程组的解为,高斯约当消去法的特点：,(1)消元和回代同时进行；,(2)乘除法的次数要比高斯消去法大，所以通常用于同时求解系数矩阵相同的多个方程组或求逆矩阵。,求逆：,则,不为0则上面方程组第一式除之,得到,交换 与 并将 代入余下各式,得到:,反复运行,可以将x全部放到左边,y放到右边,这样y的系数阵就是C的逆阵,计算公式:,即,3.变量循环重新编号法求逆,方程式,逆关系,这里就有,通过第一式计算出x0，再代入余下的方
5、程式可得到下式,新系数通过如下计算可得,对上式按下面规则重新编号,经过n次这种变换，变量均恢复原来的状态,按上法进行了k步，A为正定对称，就有 1、对于任意k，A11和A22均为正定 2、A12=-A21T,A为对称正定矩阵，用“变量循环重新编号法”经过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;f
6、or(int j=1;j=i;j+)a(i-1)*i/2+j-1=ai*(i+1)/2+j+ai0*a0j;for(i=1;i n;i+)a(n-1)*n/2+i-1=a0i;an*(n+1)/2-1=1.0/a00;,4.平方根法,A是正定对称阵，则有 例：根据矩阵乘法有：,由此可得：,计算公式为,即,下三角阵LR也为下三角阵,(i=k=j),这样，A的逆就为,平方根法步骤 1、对A做分解LLT 2、求R=L-1 3、求A-1=RRT,平方根法就是正定对称矩阵的Cholesky分解,代码,double sum;for(i=0;i=0;k-)sum-=aik*ajk;if(i=j)if(sum
7、=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位机子上有两种精度，float
8、占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;二进制表示：整数部分：101小数部分：0.2*2=0.4*2=0.8*2=1.6(0.6)*2=1.2(0.2)*2=0.4*2=
9、0.8*2=1.6(0.6)*2=1.2.0 0 1 1 0 0 1 1.,float是32位，后面尾数长度最大23位。5.2的二进制表示就为：101.00110011001100110011 一共23位。使用科学计数法表示，结果为：1.0100110011001100110011*22最高位为1，可以省略，在最后一位补一位：.01001100110011001100110*22 0 1000 0001 01001100110011001100110 1位|8位|-23位-|符号|阶码|尾数|,二进制到浮点数1 01.0 0110 0110 0110 0110 011 0 5.xxxxxxxx
10、xxxxxxxxxxxxxxxxxxxxx0.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,高斯消去法，第一步消去后，得到 下一步消去乘数为2.5*103，继续采用高斯消去得到最后一行：,假设5位的精度，做截断圆整后得到：这样就得到：回代得到：,高斯消去法中，应该选择最大的元素，做行交换后到主元素。避免不必要的圆整误差带来的解得精度降低问题。
11、,线性方程组的形态及解的误差,设x*是计算得到的解（计算机的精度导致的误差）A非奇异，理论解 一定精度的计算解与理论解的差异,例,假设在十进制精度为3位的机器上计算：换元后高斯消去，乘数为消去得,计算解与理论解的差异,残差值均是小于10-3，残差很小。将主元素换为较大的值可以保证的仅仅是残差很小。虽然采用了换元的高斯消去，但是由于截断圆整对于这样的方程仍然得到的不是正确的解,将精度扩展到6位的十进制，有 解得，什么原因？,范数对于方程Ax=bM/m比值称为A的条件数,由定义可得,由上式可以看出，条件数可以大致度量出b的相对误差对于解x的放大的倍数,改变b值，得到,C#编程中的一些关键问题,1、
12、委托机制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 numbe
13、rs 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();,
