线性方程组的解.doc

《线性方程组的解.doc》由会员分享，可在线阅读，更多相关《线性方程组的解.doc（21页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、【精品文档】如有侵权，请联系网站删除，仅供学习与交流线性方程组的解.精品文档.第四章 线性方程组的解24-1 高斯（GAUSS）消去法24.1.1 非对称方程组求解24.1.2 对称方程组求解24.1.3 对称方程组按上三角一维存放34.1.4 对称方程组按下三角一维存放34.1.5 列主元消去法44-2 消去法求逆矩阵44-3 稀疏线性方程组的形式与解算64.4 改进的平方根法求逆阵84.4 水准网平差144.4.1 观测数据的组织与输入与转换144.4.2 一维压缩存储法方程平差程序174.4.3 上三角存储法方程平差程序204.4.4 利用Matlab矩阵运算平差程序214.4.5 平差

2、结果的相互转换程序22第四章 线性方程组的解4-1 高斯（GAUSS）消去法高斯消去法在“算法语言”课程中有过详细介绍，本章将抄录一些消去法程序和求逆矩阵程序。4.1.1 非对称方程组求解 设有非对称方程Ax = b，其中A是 mm阶非奇异实矩阵，b是给定的m1向量。为了说明简便，下面以4阶方程组为例介绍消去法的具体做法。（设矩阵A的主元素不等于零）：表4-1-1 消去法解非对称线性方程a11 a12 a13 a14a21 a22 a23 a24a31 a32 a33 a34a41 a42 a43 a44b1b2b3b4消去第一列a11 a12 a13 a14 0 a22 a23 a24 0

3、a32 a33 a34 0 a42 a43 a44b1b2b3b4aij = aijai1 a1j / a11bi = bi ai1 b1 / a11i , j = 2,3,4消去第二列a11 a12 a13 a140 a22 a23 a240 0 a33 a340 0 a43 a44b1b2b3b4aij=aijai2a2j / a22bi=biai2 b2 / a22i , j = 3,4根据上表的消去过程，归纳如下：1. 消去第k个未知数是从k+1行开始，使第k列所有元素为零。（k = 1,2, m1）2. 从第k+1列开始，各元素为aij(k) = aij(k1)aik(k1) akj

4、(k1) / akk(k1)bi(k) = bi(k1)aik(k1) bk(k1) / akk(k1)用赋值语句写为 aij = aijaik akj / akk(4-1-1) bi = biaik bk / akk(4-1-2)k = 1,2, m1，i = k+1, m，j = k+1, m4.1.2 对称方程组求解当矩阵A是对称正定时，可以按行逐个消去，也就是按测量平差中高斯约化法进行。其数学模型为（以赋值语句写出）： aij = aijaik akj / akk（4-1-3）其中i =1,2, m，j = i,i+1, m通过上述消去法将矩阵A变换成上三角阵（即主对角线以下的元素均为

5、零），常数项（右端项）b也同样进行约化，采用“回代”可求得各未知参数。常数项约化的数学模型 biA(DII + J)例如第2行第3列，即DI = m + 1, I=2, J=3, 即，A(m+12+3) = A(m+2), 与上表相符合。为方便起见，令 II = DII = (I1)*(MI/2)(4-1-6)则二维数组与一维存放的对照公式为A (I, J) = A (II+J)，其中II代表第i行。按（4-1-6）式计算，J为第j列。4.1.4 对称方程组按下三角一维存放下三角一维存贮的编号如下表所示1 2 3 m 列 123m行a1 a2 a3 a4 a5 a6 am(m+1) / 2从表

6、中可以看出，主对角元号DI = I * (I+1) / 2, 因此 II DII = (IL) * I / 2（4-1-7）则二维数组与下三角一维存放的对照公式为A (I, J) = A (II+J)。必须指出，下三角存贮与上三角存贮其行列恰好相反，因此在应用（4-1-3）式约化时，应将aki akj) /akk改写为aik ajk) /akk。4.1.5 列主元消去法为了避免主元素等于零而使消去过程失败，也有时由于主元素的有效数字严重损失，而使舍入误差的相对误差增大，舍入误差扩散使得最终答解及不准确。为克服这些困难，选主元消去法是有效的。选主元消去法有：列主元消去法、行主元消去法、全选主元消

7、去法。列主元消去发是从第k列（k = 1,2, , m）的各元素寻出最大元作为主元，其具体做法如下：1. 从第1列开始，寻找本列绝对值最大的元素，记下该元素所在的行数（i0），并把该行（i0）与第1行（k行）的所有元素互换（即把第i0行移至第k行，使主元素是本列中最大元）。2. 按非对称方程组求解方法进行，亦即按（2-1）式进行消去aij = aijaikakj /akk，k = 1,2, ,m1，i , j = k+1, , m4-2 消去法求逆矩阵方程：Ax = b，其中A是mm阶非奇异矩阵。若按（4-1-1）和（4-1-2）式进行消元，则可得出下面等价方程a11 a12 a13 a140

8、 a22 a23 a240 0 a33 a34 0 0 0 a44b1b2b3b4如果每次消元均从第一行开始，则（4-1-1）式可以写为aij = aijaikakj /akk，k = 1,2, ,m1，i = 1,2, m (ik)，j = k+1, m这样消元结果为a11 0 0 00 a22 0 00 0 a33 00 0 0 a44b1b2b3b41 0 0 00 1 0 00 0 1 00 0 0 1 x1x2x3x4可见，按照此方法消元可以不回代，就可以求得未知数xi 。如果把右端项b改为I,则所求得未知数就是逆矩阵Q(Q = A-1)，或写为AQ = I =IQ = A-1即 A

9、,I 经过初等变换后化为 I, A-1 。下面介绍另一种消去法求逆设方程：b1 = a11 x1 + a12 x2 + + a1m xmb2 = a21 x1 + a22 x2 + + a2m xm (4-2-1)bm = am1 x1 + am2 x2 + + amm xm消去x1式中： 并令： （4-2-2） 为了编写程序方便，将第1行移至最后一行，第1列移到最后一列，其它各元素的行列各上升一行一列。因此，（4-2-2）式可写为 （4-2-3）调换后，将仍存放在原单元，则为 （4-2-4）这样（4-2-4）式与（4-2-1）式形式完全相同，可以用（4-2-3）式相同方法消去，依此经m次消去

10、后得 （4-2-5）可见消去后的系数，就是我们要求的逆矩阵A-1，则，上述过程适用于非对称方程组。4-3 稀疏线性方程组的形式与解算在平差中，当测量控制网较大时，它的误差方程式和法方程式的系数大多是零元素。例如一个边长误差方程式就有个系数，其中只有4个非零元素，其余个都是零元素。其中是未知数的个数），所组成的系数矩阵是稀疏矩阵。在贮存和解算稀疏线性方程组时，应充分利用其稀疏性质，避免零元素的存储和零的运算，以降低存贮量和运算次数。求解此类问题，最常用的方法仍然是高斯消去法。这里重要的问题，就是怎样形成这个稀疏矩阵。图4-3-1 水准网稀疏矩阵中零元素的分布与控制网的图形和点号的排列次序有关。下

11、面先以水准网平差为例来叙述。图4-3-1的水准网，已知和观测值表4-3-1。 表4-3-1：水准网观测值表线号高差距离A 11 22 33 B244 55 B+1.036m-0.845-0.133+0.917+1.121-0.164-0.1693.5km1.02.03.31.22.02.5误差方程式为v1 = x1 l1 权 0.4v2= x1 + x2 l2 1.0v3= x 2+ x3 l3 0.5v4= x3 l4 0.3v5= x2 + x4 l5 0.8v6= x4 + x5 l6 0.5v7= x5 l7 0.4组成法方程系数矩阵（表4-3-2），从表中可以看出，在同样的图形中，由

12、于编号的不同，其法方程矩阵的稀疏结构也不同。显然，图4-3-1的编号是由左向右方向推进，所得的稀疏矩阵的非零元素集中在主对角线附近。当水准网较大时，这个性质更为明显。 表4-3-2 法方程系数矩阵x1x2x3x4x51.41.01.02.30.51.80.50.81.81.30.50.50.9图4-3-2 水准网（重新编号）如果按图4-3-2编号，则法方程矩阵如表4-3-2。 表4-3-2 法方程系数矩阵x1x2x3x4x52.30.51.81.00.520.80.90.51.80.51.31.01.4大型稀疏矩阵一般是采用压缩形式进行存贮的，存贮方式视零元素的分布情况来定，一般是采用变带宽存

13、贮法。顾及法方程的对称性，可只存贮下三角部分（或上三角部分）。所谓带宽，就是某行的第一个非零元至该行的对角元共有元素的个数，如表2-2的法方程矩阵的带宽为k为：1 = 1，2 = 2，3 =2，4 =3，5 = 2变带宽存贮是根据每行不同的元素个数，按一维数组存放。本例存放顺序是1 = 1.4，2 = 1.0，3 = 2.3，4 = 0.5，5 = 0.8，6 = 1.8，7 = 0，8 = 1.3，9 = 0.5，10 = 0.9为了直观我们仍按二维排列于表4-3-4。12345678910 表4-3-4 一维数组存放 表4-3-5 带宽以及累计数行号带宽累计1112233254385210

14、从表4-3-5中可知，带宽的累计数恰好是表4-3-4中对角元的下标号。为此，我们用一个数组ID(i)贮存这些累计数，如ID(1) = 1, ID(2) = 3, ID(3) = 5, ID(4) = 8, ID(5) = 10等（ID(i) = ）。在这个基础上，要索取某行某列的元素是方便的。如果要从表4-3-4中取出第i行第j列元素，它就是在一维数组A中的第（ID(i)i+j）个元素，亦即A(ID(I) I+J)，例如，i = 4，j = 2, 即A（84+2）= A（6）与表4-3-4吻合。 可见，变带宽一维数组存贮的关键问题，在于如何形成对角元号的数组ID(i)。为要求出各对角元号，首先

15、要确定各行的带宽。对照图2-1与表2-2，不难看出，矩阵第i行带宽与i点相连接的各点中最小点号有关，例如第4点的相邻接点是2，5两点，最小点号是2，因此，带宽4 = 42+1（即i =i最小点号+1），也就是第4行的第一个非零元是在第2列上，所以带宽是第2、3、4三个元素。4.4 改进的平方根法求逆阵由于法方程是对称正定阵，因此可以采用改进的平方根法进行解算。平方根法是对称正定矩阵非常有效的三角分解法，由定理 设A为n阶方阵，如果其所有顺序主子式均不为零，则其存在唯一的分解式：其中，由于此处A有对称性，得，又根据A阵正定的性质，可证明D均为正数。现在设：，即 则 ，为方便记为，称为Choles

16、ky分解，即正定对称矩阵的平方根分解法。解等价于求解两个三角方程组 和。在用平方根分解法进行计算时，需要进行n次开方运算。为避免开方，可以直接采用对称正定的分解式对平方根法进行改进。从而解方程组可以按如下步骤进行：把分解成，则变成，即等价于 ，由此可解出Y和X。这称为改进的平方根法，在计算中避免了开方运算。平方根法和改进的平方根法的计算量和存贮量比消去法节约近一半，而且不需要选主元，能得到比较精确的数值解。法方程用改进平方根法解算的过程如下： 分解：其中，纯量计算公式为： 求逆： ，由，得到纯量计算公式： 通式为： 求积： 其纯量计算公式：计算R阵需要除以S的对角元，计算Q阵要乘S的对角元，做

17、计算上的改进。令 ，则： 分解： 求逆： 乘积： 所以function c=invsqr(c,n)ss=0.0;for i=1:n di=(i-1)*(n-i/2.0); for j=i:n ss=c(di+j); for k=1:i-1 dk=(k-1)*(n-k/2.0); ss=ss-c(dk+i)*c(dk+j)/c(dk+k); end if(j=i) c(di+j)=1/ss; else c(di+j)=ss*c(di+i); end endendfor i=1:n-1di=(i-1)*(n-i/2.0);for j=i+1:nss=-c(di+j);for k=i+1:j-1dk

18、=(k-1)*(n-k/2.0);ss=ss-c(di+k)*c(dk+j);endc(di+j)=ss;endendfor i=1:n-1di=(i-1)*(n-i/2.0);for j=i:ndj=(j-1)*(n-j/2.0);if(j=i)ss=c(di+j);elsess=c(di+j)*c(dj+j);endfor k=j+1:ndk=(k-1)*(n-k/2.0);ss=ss+c(di+k)*c(dj+k)*c(dk+k); endc(di+j)=ss; endendreturn4.4 水准网平差程序设计4.4.1 观测数据的组织与输入与转换1 水准网数据变量的约定表4-4-1：

19、水准网数据变量的约定变量名说明变量名说明ed已知点个数h0已知点高程dd未知点个数be起点点号sd总点数en终点点号gd观测值个数h1高差观测值pn点号s距离观测值2 数据文件的组织下面给出一个水准网输入数据文件的例子：3 3 6（已知点个数、未知点个数、观测值个数）101 102 103 104 105 106（点号）34.788 35.259 37.825（已知点高程）104 101 1.652 4.5（起点点号、终点点号、高差观测值、距离观测值）101 102 -0.418 3.1105 102 0.714 3.4102 103 1.243 3.8106 103 -0.577 4.310

20、3 101 -0.786 2.5其中编号数组未知点在前，已知点在后。3 编号转换函数在形成误差方程及法方程时需要将点号进行转化，即将文件中的点号101 102 103 104 105 106变为计算的序号1 2 3 4 5 6。编号转换相应的程序如下：function n1,k=chkdat(sd,pn,n1)n=length(n1);k=0;for i=1:n i1=0; for j=1:sd if(n1(i)=pn(j) i1=1; n1(i)=j; break; end end if(i1=0)% fprintf(fit2,%5d %5dn,i,n1(i); k=1; endendret

21、urn4 数据读入与近似高程计算程序function ed,dd,sd,gd,pn,h0,k1,k2,h1,s=readlevelnetdataglobal pathname net_name s_datafile b_datafile;global ed dd sd pn gd h0 k1 k2 h1 s k11 k12;k1=;k2=;h=;s=;if(isempty(pathname)|isempty(net_name) filename,pathname=uigetfile(*.txt,Input filename);i=find(.=filename);net_name=filena

22、me(1:i-1);endfid1=fopen(strcat(pathname,net_name,s_datafile),rt);if(fid1=-1) msgbox(Input File or Path is not correct,Warning,warn); return;end %open a file to read%open a file to readed=fscanf(fid1,%f,1); %已知点个数dd=fscanf(fid1,%f,1); %未知点个数sd=ed+dd; %总点数gd=fscanf(fid1,%f,1); %观测值个数pn=fscanf(fid1,%f,

23、sd); %点号%known datah0=fscanf(fid1,%f,ed); %已知点高程h0(dd+1:ed+dd)=h0(1:ed);heightdiff=fscanf(fid1,%f,4,gd);heightdiff=heightdiff;k1=heightdiff(:,1); %起点k2=heightdiff(:,2); %终点k11=heightdiff(:,1); %起点k12=heightdiff(:,2); %终点h1=heightdiff(:,3); %高差s=heightdiff(:,4); %距离fclose(all);%点号转换k1,k01=chkdat(sd,p

24、n,k1);k2,k02=chkdat(sd,pn,k2);h0(1:dd)=20000.;ie=0;while(1)%计算近似高程 for k=1:gdi=k1(k);j=k2(k);if(h0(i)1e4)h0(j)=h0(i)+h1(k);ie=ie+1;endif(h0(i)1e4&h0(j)1e4)h0(i)=h0(j)-h1(k);ie=ie+1;end endif(ie=dd) break; endendh0=reshape(h0,length(h0),1);return其中包含了近似高程的计算，即进行循环，程序中设置了一个变量ie，每计算出一点的高程其值加1，当其值等于未知点的

25、个数时停止循环。图4-4-1 平差程序总体框架4.4.2 一维压缩存储法方程平差程序1 变带宽紧缩存储中对角元号数组ID（I）的形成function id=bm1(gd,dd,k1,k2)%计算一维压缩存放的数组idid=;for i=1:ddk=i;for j=1:gd i1=k1(j);i2=k2(j);if(i1=i&i2k) k=i2; endif(i2=i&i1k) k=i1; endend id(i)=k;endfor i=2:dd id(i)=id(i-1)+i-id(i)+1;endreturn2 一维压缩存储法方程平差global pathname net_name s_da

26、tafile a1_datafile;global ed dd sd pn gd h0 k1 k2 h1 s dh;p=1.0./s;id=bm1(gd,dd,k1,k2);mm=id(dd);a(1:mm)=0.0;b(1:dd )=0.0;for k=1:gd %形成法方程i=k1(k);j=k2(k);hl=h1(k)+h0(i)-h0(j);if(i=dd )ii=id(i)-i;a(ii+i)=a(ii+i)+p(k);b(i)=b(i)-hl*p(k);endif(j=dd )jj=id(j)-j;a(jj+j)=a(jj+j)+p(k);b(j)=b(j)+hl*p(k); if

27、(i=j) a(ii+j)=a(ii+j)-p(k); else a(jj+i)=a(jj+i)-p(k); end end endenda=gs5(dd,a,id);%变带宽下三角紧缩存储高斯消元法dh=cy6(a,b,id,dd,1);%常数项约化与回代子程序dh(dd+1:ed+dd)=0.0;hm(dd+1:ed+dd)=0.0;for i=1:sd h(i)=h0(i)+dh(i);endvv=0.0;for i=1:gd L(i)=h(k2(i)-h(k1(i);v(i)=h(k2(i)-h(k1(i)-h1(i);vv=vv+v(i)*v(i)/s(i);endu=sqrt(vv

28、/(gd-dd);for i=1:ddb(1:dd )=0.0;b(i)=1.0;b=cy6(a,b,id,dd,i);hm(i)=sqrt(b(i)*u;endwritelevelnetdata(pn,k1,k2,h1,v,L,h0,dh,h,hm,u);function a=gs5(dd,a,id)%变带宽高斯消去法for i=1:ddii=id(i)-i; if(i-1=0) li=1-ii; else li=id(i-1)-ii+1; endfor j=li:ijj=id(j)-j; if(j-1)=0 lj=1-jj; else lj=id(j-1)-jj+1; endlk=li;i

29、f(lilj) lk=lj; endfor k=lk:j-1kk=id(k); a(ii+j)=a(ii+j)-a(ii+k)/a(kk)*a(jj+k); endendendreturnfunction b=cy6(a,b,id,dd ,k1)%常数项约化与回代子程序for i=k1:dd ii=id(i)-i; if i=1 nd=id(i); else nd=id(i)-id(i-1); ende=0.0;for k=1:i-1if(i-k)nd) e=e+a(ii+k)*b(k); end end b(i)=(b(i)-e)/a(ii+i);endfor i=dd-1:-1:k1ii=

30、id(i);for k=i+1:dd kk=id(k)-k;nk=id(k)-id(k-1);if(k-ink) b(i)=b(i)-a(kk+i)/a(ii)*b(k); endendendreturn3 输出数据函数function writelevelnetdata(pn,be,en,h1,v,L,h0,dh,h,hm,uw0)global pathname net_name s_datafile a1_datafile;fid2=fopen(strcat(pathname,net_name,a1_datafile),wt);if(fid2=-1) msgbox(Error by Ope

31、ning Output File,Warning,warn); break;end fprintf(fid2,待定点高程平差值及中误差：n pn h0(m) dh(m) h(m) hmn);fprintf(fid2,%5d %10.4f %10.4f %10.4f %10.4fn,pn,h0,dh,h,hm);fprintf(fid2,高差观测值平差值：n pn pn h0(m) V(m) H(m)n);fprintf(fid2,%5d %5d %10.4f %10.4f %10.4fn,pn(be),pn(en),h1,v,L);fprintf(fid2,单位权中误差：%10.4fmn,uw

32、0);fclose(fid2);open(strcat(pathname,net_name,a1_datafile);return4.4.3 上三角存储法方程平差程序mm=(dd+1)*dd/2;a(1:mm)=0.0;b(1:dd)=0.0;for k=1:gdi=k1(k);j=k2(k);hl=h1(k)+h0(i)-h0(j);if(i=dd)ii=(i-1)*(dd-i/2.0);a(ii+i)=a(ii+i)+1/s(k);b(i)=b(i)+1./s(k)*hl; endif(j=dd)jj=(j-1)*(dd-j/2.0);a(jj+j)=a(jj+j)+1/s(k);b(j)

33、=b(j)-1./s(k)*hl;if(i=dd)if(ij) a(ii+j)=a(ii+j)-1/s(k);else a(jj+i)=a(jj+i)-1/s(k); end end endenda=invsqr(a,dd);for i=1:dd dh(i)=0.0; di=(i-1)*(dd-i/2.0); for j=1:dd dj=(j-1)*(dd-j/2.0); if(ji) dh(i)=dh(i)-a(dj+i)*b(j); else dh(i)=dh(i)-a(di+j)*b(j); end endenddh(dd+1:ed+dd)=0.0;hm(dd+1:ed+dd)=0.0;

34、for i=1:sd h(i)=h0(i)+dh(i);endvv=0.0;for i=1:gd L(i)=h(k2(i)-h(k1(i);v(i)=h(k2(i)-h(k1(i)-h1(i);vv=vv+v(i)*v(i)/s(i);enduw0=sqrt(vv/(gd-dd);for i=1:ddii=(i-1)*(dd-i/2.0);hm(i)=sqrt(a(ii+i)*uw0;endwritelevelnetdata(pn,k1,k2,h1,v,L,h0,dh,h,hm,uw0);4.4.4 利用Matlab矩阵运算平差程序function level3%ed,dd,sd,gd,pn,

35、h0,k1,k2,h1,s=readlevelnetdata;global pathname net_name s_datafile a1_datafile;global ed dd sd pn gd h0 k1 k2 h1 s dh;dh,h,V,L,uw0,uwh,uwl=calculatelevelnet(ed,dd,sd,gd,pn,h0,k1,k2,h1,s);writelevelnetdata(pn,k1,k2,h1,V,L,h0,dh,h,uwh,uw0); %输出水准网解算结果returnfunction dh,h,V,L,uw0,uwh,uwl=calculateleveln

36、et(ed,dd,sd,gd,pn,h0,be,en,hd,distance)% 水准网平差A=sparse(zeros(sd,gd); %求解系数阵b=(0:(gd-1)*sd;A(be+b)=-1;A(en+b)=1;A=A;A=A(:,1:dd);l=zeros(gd,1); %求解常数项 l=h0(be)-h0(en)+hd;p=diag(1./distance); %权阵dh=inv(A*p*A)*A*p*l; %高程改正数h00=h0(dd+1:sd);h0=h0(1:dd); %待定点高程近似值h=h0+dh; %待定点高程平差值V=A*dh-l; %高差观测值改正数L=hd+V; %高差观测值平差值%精度评定uw0=sqrt(V*p*V/(gd-dd); %单位权中误差Qxx=inv(A*p*A);uwh=uw0*sqrt(diag(Qxx); %待定点高程平差值中误差