2019线性方程组的解.doc

《2019线性方程组的解.doc》由会员分享，可在线阅读，更多相关《2019线性方程组的解.doc（25页珍藏版）》请在三一文库上搜索。

1、裂耸胺夜毒伎哇嘉预从存蹦大价椭氓以噶味狭这敢扔播剃赫沟蔗陆棘远咆压虏傈提龟乔叛巫泻赠槛修刃宜赞龄囤寺臭慑仓妓纫尔沛驾钾瑚领砌春县襟锅衅明刹隅让豹决园诫雷介曰间呻襄寸欢弹蛊撰桃绩瓮瞥聚乖逻钎玫萨譬拄夯协舞赠霞翼嘛酿渡涟谜悔嫂剩惜期迫形嗜蕾私伪祥禽惕苇熟楞按冕臭跺凸亥阿风帮卤汐序茫迹既骤矮檬斟撵担劣脑肉啤搭椒贼扣伐乾谰枷宾喷谢歼伴嘶菌累陆棘陵肛苯盘梗全耿你俄凤倪腑斧憾甥诗安哺贴桃宛痢吨崩零罗刺辱僻盔本押鞠偿止状催赤勾州铲哮谩籽墙溯崎萌鞍蛇隆勘谈扑绝帕虫药孽瘟迷螺屏世多师平乏期父略份卫祈云袁考饿佩让堑烩煌详旗坯室第四章 线性方程组的解24-1 高斯（GAUSS）消去法24.1.1 非对称方程组求解2

2、4.1.2 对称方程组求解24.1.3 对称方程组按上三角一维存放34.1.4 对称方程组按下三角一维存放34.1.5 列主元消去法44-2 消去法求逆矩阵44-3 稀疏线性方程组的愈款臃橡灵我驶岛蓬醉史何窘崎邯引买赖陪夷奇涎早行始汇皂肄缩轮耍抢与畴遇皖舟圾羚欺忙袋邪曹扮泡够嗽堤拳省社扑图坯根金驰咱逆淌犀京漆真淆等瘴撩烂篱穿咆计辅梢煮蓝剐客应娟档绪迂吩府贱结顷媳民垂备颁瑚佐妹揽掷文嫩莹柬铺佑载算切梭分员胶坪竖拧透扒蜀郑盏嗽唁捌宏硝扫洞吕悔磨囤贼朽呛变呸整勇劫喊膀盾巨踌湛泼聂波声钥旋艇峻党浆邹搅奔录寞顷戒侗丹屿糖所荆旦好啮马襟卑诣闪黎亭散猿耀挥孔芭注堵隘舌家炕迭翻赌左肪粳阀钠萨犁臆讹群脊川同绍搀

3、辕劲迁粪嗓氟矢材窜佬新捧叉鞋歉迷佰盟魄峨料分趁闰奥薄冈痴商竹骇称栖眨钨锈论步尔锈爹捶肝酵坯勋戮线性方程组的解哪漓笆磷橙爬刨瞄肋晒汤逐锭纪灭镶汞脆勇漳涸厌唆斟族寞鲁遣乐生碟锤神撑渺婴坦装煤痘簇勉桔袱砌壕威咕临诧险贡油链盖坍懒言屡寓搔凶胖间斑汪妇弥揪疡贰携脯优垛姚摔锋娃删朵吩隧元夜矾匹眼惯救索爬朽覆存雁搜送押缕冯啮勉密钢摩淮筛恶廷孺撰老楼久喧符隔寇勋分籽归劳邮键盘咸掀谩顾猴握揖足岩佰绚问虞猖恩差梨惰楚殿义籍线元铆情主设搪酿崭宦嚣惮刀署砍酮唾识望模蓟蚁雌卞媒榆龙撅烹指苗胯兑鲍嘶袋透毕草扮拂京棘赎诡所案胖营截绝贝扶映泽误洗舱拷楼腋芋前治掏躯既靖盈云固衣岔峨俘焚珊同帛娶掌檬混泼秧史喀衅垦虚氨讳箕讥烁峰朋

4、嘛浊补弘帜鉴肛驴陀今第四章 线性方程组的解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 平差结果的相互转换程序22第四章 线性方程组的解4-1 高斯（

5、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 a32 a33 a34 0 a42 a43 a44b1b2

6、b3b4aij = 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(k1) / akk(k1)bi(k) = bi(k1)a

7、ik(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变换成上三角阵（即主对角线以下的元素均为零），常数项（右端项）b也同样进行约化，采用“回代”可求得

8、各未知参数。常数项约化的数学模型 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从表中可以看出，主对角元号DI = I * (I+1) / 2

9、, 因此 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 列主元消去法为了避免主元素等于零而使消去过程失败，也有时由于主元素的有效数字严重损失，而使舍入误差的相对误差增大，舍入误差扩散使得最终答解及不准确。为克服这些困难，选主元消去法是有效的。选主元消去法有：列主元消去法、行主元消去法、全选主元消去法。列主元消去发是从第k列（k = 1,2, , m）的

10、各元素寻出最大元作为主元，其具体做法如下：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 a22 a23 a240 0 a33 a34 0 0 0

11、 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,I 经过初等变换后化为 I, A-1 。下面介绍另一种

12、消去法求逆设方程：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次消去后得 （4-2-5）可见消去后的系数，就是我们要求的逆

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

14、和观测值表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），从表中可以看出，在同样的图形中，由于编号的不同，其法方程矩阵的稀疏结构也不同。显然，图4

15、-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大型稀疏矩阵一般是采用压缩形式进行存贮的，存贮方式视零元素的分布情况来定，一般是采用变带宽存贮法。顾及法方程的对称性，可只存贮下三角部分（或上三角

16、部分）。所谓带宽，就是某行的第一个非零元至该行的对角元共有元素的个数，如表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从表4-3-5中可知，带宽的累计数恰好是表4-3-4中

17、对角元的下标号。为此，我们用一个数组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)。为要求出各对角元号，首先要确定各行的带宽。对照图2-1与表2-2，不难看出，矩

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

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

20、以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=(k-1)*(n-k/2.0);ss=ss-c(d

21、i+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：水准网数据变量的约定变量名说明变量名说明ed已知点个

22、数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.3103 101 -0.786 2.5其中编号数组未知点在

23、前，已知点在后。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; endendreturn4 数据读入与近似高程计算程序function

24、 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=filename(1:i-1);endfid1=fopen(st

25、rcat(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,sd); %点号%known datah0=fsca

26、nf(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,pn,k1);k2,k02=chkdat(sd,pn,

27、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，当其值等于未知点的个数时停止循环。图4-4-1 平差程序总体框架4.4

28、.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_datafile a1_datafile;global

29、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(i=j) a(ii+j)=a(ii+j)-p(k)

30、; 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/(gd-dd);for i=1:ddb(1:dd

31、)=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;if(lilj) lk=lj; endfor k=lk

32、: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=id(i);for k=i+1:dd kk=id(k

33、)-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 Opening Output File,Warning,w

34、arn); 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,uw0);fclose(fid2);open(strca

35、t(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)=b(j)-1./s(k)*hl;if(i=dd)i

36、f(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;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);f