双向解析光束法上课讲义

1、双向解析光束法双向解析光束法光束法程序有问题，在 Geteleme nt这个函数里便出现索引超限，这个问题一直 解决不了光束法的流程：1根据同名像点对对相交理论求系数阵A，系数阵B和常数阵La11=(a1f+a3x)/Z; a12=(b1f+b3x)/Z; a13=(c1f+c3x)/Z; a14=ysin(omega)-x/f(xcos(kappa)-ysi n(kappa)+fcos(kappa)cos(omega);a15=-fsi n(kappa)-x/f(xsi n(kappa)+ycos(kappa);a16=y;a21=(a2f+a3y)/Z; a22=(b2f+b3y)/Z;

2、a23=(c2f+c3y); a24=-xsin(omega)-x/f(xcos(kappa)-ysi n(kappa)-fsi n(kappa)cos(omega)a25=-fcos(kappa)-y/f(xsi n(kappa)+ycos(kappa);a26=-x;2. 求方程的改化法方程求出外方位元素和物方坐标改正数3. 判断改正数的值，如果小于限差则输出结果光束法是最严密的一种方法的原因：在一张相片中，待定点与控制点的像点与摄影中心及相应地面点均构成一条光 束，该方法是以每张相片所组成的一束光线作为平差的基本单元，已共线条件 方程作为平差的基础方程，通过各个光束在空间中的旋转和平移，

3、使模型之间 公共点的光线实现最佳交汇，并使整个区域纳入到已知的地面控制点坐标系 中，所以要建立全区域统一的误差方程，整体解求全区域内每张相片的六个外 方位元素及所有待定点坐标，光束法区域网平差是基于摄影时像点，物点和摄 站点三点共线提出来的。由单张相片构成区域，其平差的数学模型是共线条件 方程，平差单元是单个光束，像点坐标是观测值，未知数是每张相片的外方位 元素及所有待定点坐标。误差方程直接由像点坐标的观测值列出，能对像点坐 标进行系统误差改正。光束法的程序代码为：计算像片外方位元素，架设 phi=0,omega=0,kappa=0求Xs,Ys,Zs求左片的Xs,Ys,ZsXsl = (str

4、X0 + strX2) / 2;Ysl = (strY0 + strY2) / 2;L = Math.Pow(Math.Pow(strX0 - strX2, 2) + Math.Pow(strY0 - strY2, 2) + Math.Pow(strZ0 - strZ2, 2), 0.5);l = Math.Pow(Math.Pow(strXI0 - strXI2, 2) + Math.Pow(strYI0 - strYI2, 2), 0.5);H = f * L / l;Zsl = (strZ0 + strZ2) / 2 + H;计算片左的加密点的物方坐标Class1 xyz = new C

5、lass1(8, 3, strXY);Classi xyzt=xyz.Tra nspose ();phi,omega,kappa为零，故旋转矩阵为单位阵/像空间辅助坐标系中的坐标Classi UVW = xyzt.Multiply(H);求地面摄影测量坐标系中的坐标for (int i = 0; i < 8; i+)UVW.SetElement(0, i, UVW.GetElement(0, i) + Xsl);UVW.SetElement(1, i, UVW.GetElement(1, i) + Ysl);UVW.SetElement(2, i, UVW.GetElement(2, i

6、) + Zsl);求右片的Xs,Ys,ZsXsr = (strX1 + strX3) / 2;Ysr = (strY1 + strY3) / 2;L = Math.Pow(Math.Pow(strX1 - strX3, 2) + Math.Pow(strY1 - strY3, 2) + Math.Pow(strZ1 - strZ3, 2), 0.5);l = Math.Pow(Math.Pow(strXr1 - strXr3, 2) + Math.Pow(strYr1- strYr3, 2), 0.5);H = f * L /I;Zsr = (strZ1 + strZ3) / 2 + H;右片

7、加密点double strxy;strxy = new double24;Console.WriteLine("请输入右片加密点的相片坐标");for (int i = 0; i < 8; i+)/strxy0 + 3 * i = Con vert.ToDouble(Co nsole.ReadLi ne();/strxy1 + 3 * i = Con vert.ToDouble(Co nsole.ReadLi ne(); strxy2 + 3 * i = -f;for (int i = 0; i < 24; i+)strxyi = strxyi * 0.0000

8、25;计算片左的加密点的物方坐标Classi XYZ = new Class1(8, 3, strxy);Classi XYZt=XYZ .Tra nspose ();phi,omega,kappa为零，故旋转矩阵为单位阵/像空间辅助坐标系中的坐标Classi uvw = XYZt.Multiply(H);求地面摄影测量坐标系中的坐标for (int i = 0; i < 8; i+)uvw.SetEleme nt(0, i, (uvw.GetEleme nt(0, i) + Xs叶 UVW.GetEleme nt (0,i)/2);uvw.SetEleme nt(1, i, (uvw.

9、GetEleme nt(1, i) + Ys叶 UVW.GetEleme nt(1.1) )/2);uvw.SetElement(2, i, (uvw.GetElement(2, i) + Zsr+UVW .GetElement(2.1) )/2);/从此开始循环光束法进行平差/循环体do求左片旋转矩阵Rdouble mtrRI;mtrRI = new double9;mtrRI0 = Math.Cos(phil) * Math.Cos(kappal) - Math.Sin(phil) * Math.S in (omegal) * Math.Si n( kappal);mtrRl1 = -Ma

10、th.Cos(phil) * Math.Si n( kappal) - Math.Si n( phil) * Math.Si n(omegal) * Math.Cos(kappal);mtrRl2 = -Math.Sin(phil) * Math.Cos(omegal);mtrRl3 = Math.Cos(omegal) * Math.Si n(kappal);mtrRl4 = Math.Cos(omegal) * Math.Cos(kappal);mtrRl 5 = -Math.Si n(omegal);mtrRl6 = Math.Sin(phil) * Math.Cos(kappal) +

11、 Math.Cos(phil) ” Math.S in (omegal) * Math.Si n( kappal);mtrRl7 = -Math.Sin(phil) * Math.Sin(kappal) + Math.Cos(phil) ” Math.Si n(omegal) * Math.Cos(kappal);mtrRl8 = Math.Cos(phil) * Math.Cos(omegal);求右片旋转矩阵Rrdouble mtrRr;mtrRr = new double9;mtrRr0 = Math.Cos(phir) * Math.Cos(k appar) - Math.Si n(p

12、hir) Math.Sin(omegar) * Math.Sin(kappar);mtrRr1 = -Math.Cos(phir) * Math.Si n(kappar) - Math.Si n(phir)Math.Sin(omegar) * Math.Cos(kappar);mtrRr2 = -Math.Si n(phir) * Math.Cos(omegar);mtrRr3 = Math.Cos(omegar) * Math.Si n(k appar);mtrRr4 = Math.Cos(omegar) * Math.Cos(k appar);mtrRr5 = -Math.Si n(ome

13、gar);mtrRr =Math.Si n(phir) * Math.Cos(k appar) + Math.Cos(phir) Math.Sin(omegar) * Math.Sin(kappar);mtrRr7 = -Math.Sin(phir) * Math.Sin(kappar) + Math.Cos(phir) Math.Sin(omegar) * Math.Cos(kappar);mtrRr8 = Math.Cos(phir) * Math.Cos(omegar);/求系数阵Adouble mtrA;mtrA = new double576;for (int i = 0; i &l

14、t; 8; i+)mtrA0 + 48 * i = (mtrRI0 * f + mtrRI2 * xyz.GetElement(i, 0) / (mtrRI2 * (uvw.GetElement(0, i) - Xsl) + mtrRI5 * (uvw.GetElement(1, i) - Ysl) + mtrRI8 * (uvw.GetElement(2, i) - Zsl);mtrA1 + 48 * i = (mtrRI3 * f + mtrRI5 * xyz.GetElement(i, 0) / (mtrRI2 * (uvw.GetElement(0, i) - Xsl) + mtrRI

15、5 * (uvw.GetElement(1, i) - Ysl) + mtrRI8 * (uvw.GetElement(2, i) - Zsl);mtrA2 + 48 * i = (mtrRI6 * f + mtrRI8 * xyz.GetElement(i, 0) / (mtrRI2 * (uvw.GetElement(0, i) - Xsl) + mtrRI5 * (uvw.GetElement(1, i) - Ysl) + mtrRI8 * (uvw.GetElement(2, i) - Zsl);mtrA3 + 48 * i = xyz.GetElement(i, 1) * Math.

16、Sin(omegal)-(xyz.GetElement(i, 0) * (xyz.GetEIement(i, 0) * Math.Cos(kappal)- xyz.GetElement(i, 1) * Math.Sin(kappal) / f) + f * Math.Cos(kappal) * Math.Cos(omegal);mtrA4 + 48 * i = -f * Math.Sin(kappal) - xyz.GetElement(i, 0) *(xyz.GetElement(i, 0) * Math.Sin(kappal) + xyz.GetElement(i, 1) * Math.C

17、os(kappal) / f;mtrA5 + 48 * i = xy z. GetEleme nt(i, 1);for (i nt j = 0; j < 6; j+)mtrAj + 6 + 48 * i = 0;mtrA12 + 48 * i = (mtrRI1 * f + mtrRI2 * xyz.GetEIeme nt(i, 1) / (mtrRI2 * (uvw.GetElement(0, i) - Xsl) + mtrRI5 * (uvw.GetElement(1, i) - Ysl) + mtrRI8 * (uvw.GetElement(2, i) - Zsl);mtrA13

18、+ 48 * i = (mtrRI4 * f + mtrRI 5 * xyz.GetEIeme nt(i, 1) / (mtrRI2 * (uvw.GetElement(0, i) - Xsl) + mtrRI5 * (uvw.GetElement(1, i) - Ysl) + mtrRI8 * (uvw.GetElement(2, i) - Zsl);mtrA14 + 48 * i = (mtrRI7 * f + mtrRI8 * xyz.GetEIeme nt(i, 1) / (mtrRI2 * (uvw.GetElement(0, i) - Xsl) + mtrRI5 * (uvw.Ge

19、tElement(1, i) - Ysl) + mtrRI8 * (uvw.GetElement(2, i) - Zsl);mtrA15 + 48 * i = -xyz.GetEIement(i, 0) * Math.Sin(omegal)-(xyz.GetElement(i, 0) * (xyz.GetEIement(i, 0) * Math.Cos(kappal)- xyz.GetElement(i, 1) * Math.Sin(kappal) / f - f * Math.Sin(kappal) * Math.Cos(omegal);mtrA16 + 48 * i = -f * Math

20、.Cos(k appal) - xyz. GetEleme nt(i, 1) *(xyz.GetElement(i, 0) * Math.Sin(kappal) + xyz.GetElement(i, 1) * Math.Cos(kappal) / f;mtrA17 + 48 * i = -xyz.GetEIement(i, 0);for (i nt j = 0; j < 12; j+)mtrAj + 18 + 48 * i = 0;mtrA30 + 48 * i = (mtrRrO * f + mtrRr2 * XYZ.GetElement(i, 0) / (mtrRr2 * (uvw

21、.GetElement(0, i) - Xsr) + mtrRr5 * (uvw.GetElement(1, i) - Ysr) + mtrRr8 * (uvw.GetElement(2, i) - Zsr);mtrA31 + 48 * i = (mtrRr3 * f + mtrRr5 * XYZ.GetElement(i, 0) / (mtrRr2 * (uvw.GetElement(0, i) - Xsr) + mtrRr5 * (uvw.GetElement(1, i) - Ysr) + mtrRr8 * (uvw.GetElement(2, i) - Zsr);mtrA32 + 48

22、* i = (mtrRr6 * f + mtrRr8 * XYZ.GetElement(i, 0) / (mtrRr2 * (uvw.GetElement(0, i) - Xsr) + mtrRr5 * (uvw.GetElement(1, i) - Ysr) + mtrRr8 * (uvw.GetElement(2, i) - Zsr);mtrA33 + 48 * i = XYZ.GetElement(i, 1) * Math.Sin(omegar)- (XYZ.GetEleme nt(i, 0) * (XYZ.GetEleme nt(i, 0) * Math.Cos(kappar)- XY

23、Z.GetElement(i, 1) * Math.Sin(kappar) / f) + f * Math.Cos(kappar) * Math.Cos(omegar);mtrA34 + 48 * i = -f * Math.Sin(kappar) - XYZ.GetElement(i, 0) * (XYZ.GetElement(i, 0) * Math.Sin(kappar) + XYZ.GetElement(i, 1) * Math.Cos(k appar) / f;mtrA35 + 48 * i = XYZ.GetEleme nt(i, 1);for (i nt j = 0; j <

24、; 6; j+)mtrAj + 36 + 48 * i = 0;mtrA42 + 48 * i = (mtrRr1 * f + mtrRr2 * XYZ.GetEleme nt(i, 1) / (mtrRr2 * (uvw.GetElement(0, i) - Xsr) + mtrRr5 * (uvw.GetElement(1, i) - Ysr) + mtrRr8 * (uvw.GetElement(2, i) - Zsr);mtrA43 + 48 * i = (mtrRr4 * f + mtrRr5 * XYZ.GetElement(i, 1) / (mtrRr2 * (uvw.GetEl

25、ement(0, i) - Xsr) + mtrRr5 * (uvw.GetElement(1, i) - Ysr) + mtrRr8 * (uvw.GetElement(2, i) - Zsr);mtrA44 + 48 * i = (mtrRr7 * f + mtrRr8 * XYZ.GetEleme nt(i, 1) / (mtrRr2 * (uvw.GetElement(0, i) - Xsr) + mtrRr5 * (uvw.GetElement(1, i) - Ysr) + mtrRr8 * (uvw.GetElement(2, i) - Zsr);mtrA45 + 48 * i =

26、 -XYZ.GetElement(i, 0) * Math.Sin(omegar)- (XYZ.GetEleme nt(i, 0) * (XYZ.GetEleme nt(i, 0) * Math.Cos(kappar)- XYZ.GetElement(i, 1) * Math.Sin(kappar) / f - f * Math.Sin(kappar) * Math.Cos(omegar);mtrA46 + 48 * i = -f * Math.Cos(kappar) - XYZ.GetElement(i, 1) * (XYZ.GetElement(i, 0) * Math.Sin(kappa

27、r) + XYZ.GetElement(i, 1) * Math.Cos(k appar) / f;mtrA47 + 48 * i = -XYZ.GetElement(i, 0);for (int i = 8; i < 12; i+)mtrAO + 48 * i = (mtrRIO * f + mtrRI2 * strXIi - 8) / (mtrRI2 * (strXi - 8 - Xsl) + mtrRI5 * (strYi - 8 - Ysl) + mtrRI8 * (strZi - 8 - Zsl);mtrA1 + 48 * i = (mtrRI3 * f + mtrRI 5 *

28、 strXIi - 8) / (mtrRI2 * (strXi - 8 - Xsl) + mtrRI 5 * (strYi - 8 - Ysl) + mtrRI8 * (strZi - 8 - Zsl);mtrA2 + 48 * i = (mtrRI6 * f + mtrRI8 * strXIi - 8) / (mtrRI2 * (strXi - 8 - Xsl) + mtrRI 5 * (strYi - 8 - Ysl) + mtrRI8 * (strZi - 8 - Zsl);mtrA3 + 48 * i = strYIi - 8 * Math.Si n(omegal) - (strXIi

29、 - 8 * (strXIi - 8 * Math.Cos(k appal) - strYIi - 8 * Math.S in (kappal) / f) + f * Math.Cos(k appal) * Math.Cos(omegal);mtrA4 + 48 * i = -f * Math.Sin(kappal) - strXIi - 8 * (strXIi - 8 * Math.Si n( kappal) + strYIi - 8 * Math.Cos(kappal) / f;mtrA5 + 48 * i = strYIi - 8;for (i nt j = 0; j < 6; j

30、+)mtrAj + 6 + 48 * i = 0;mtrA12 + 48 * i = (mtrRI1 * f + mtrRI2 * strYIi - 8) / (mtrRI2' (strXi - 8 - Xsl) + mtrRI 5 * (strYi - 8 - Ysl) + mtrRI8 * (strZi - 8 - Zsl);mtrA13 + 48 * i = (mtrRI4 * f + mtrRI 5 * strYIi - 8) / (mtrRI2' (strXi - 8 - Xsl) + mtrRI 5 * (strYi - 8 - Ysl) + mtrRI8 * (s

31、trZi - 8 - Zsl);mtrA14 + 48 * i = (mtrRI7 * f + mtrRI8 * strYIi - 8) / (mtrRI2' (strXi - 8 - Xsl) + mtrRI 5 * (strYi - 8 - Ysl) + mtrRI8 * (strZi - 8 - Zsl);mtrA15 + 48 * i = -strXIi - 8 * Math.Sin(omegal) - (strXIi - 8 * (strXIi - 8 * Math.Cos(k appal) - strYIi - 8 * Math.Si n( kappal) / f - f

32、* Math.Si n( kappal) * Math.Cos(omegal);mtrA16 + 48 * i = -f * Math.Cos(k appal) - strYIi - 8 * (strXIi - 8 Math.Si n( kappal) + strYIi - 8 * Math.Cos(kappal) / f;mtrA17 + 48 * i = -strXIi - 8;for (i nt j = 0; j < 12; j+)mtrAj + 18 + 48 * i = 0;mtrA30 + 48 * i = (mtrRr0 * f + mtrRr2 * strXri - 8)

33、 / (mtrRr2 (strXi - 8 - Xsr) + mtrRr5 * (strYi - 8 - Ysr) + mtrRr8 * (strZi - 8 - Zsr);mtrA31 + 48 * i = (mtrRr3 * f + mtrRr5 * strXri - 8) / (mtrRr2 (strXi - 8 - Xsr) + mtrRr5 * (strYi - 8 - Ysr) + mtrRr8 * (strZi - 8 - Zsr);mtrA32 + 48 * i = (mtrRr6 * f + mtrRr8 * strXri - 8) / (mtrRr2 (strXi - 8

34、- Xsr) + mtrRr5 * (strYi - 8 - Ysr) + mtrRr8 * (strZi - 8 - Zsr);mtrA33 + 48 * i = strYri - 8 * Math.Sin(omegar) - (strXri - 8 * (strXri - 8 * Math.Cos(k appar) - strYri - 8 * Math.Si n(kappar) / f) + f * Math.Cos(k appar) * Math.Cos(omegar);mtrA34 + 48 * i = -f * Math.Si n(kap par) - strXri - 8 * (

35、strXri - 8 Math.S in (kappar) + strYri - 8 * Math.Cos(k appar) / f;mtrA35 + 48 * i = strYri - 8;for (i nt j = 0; j < 6; j+)mtrAj + 36 + 48 * i = 0;mtrA42 + 48 * i = (mtrRr1 * f + mtrRr2 * strYri - 8) / (mtrRr2 (strXi - 8 - Xsr) + mtrRr5 * (strYi - 8 - Ysr) + mtrRr8 * (strZi - 8 - Zsr);mtrA43 + 48

36、 * i = (mtrRr4 * f + mtrRr5 * strYri - 8) / (mtrRr2 (strXi - 8 - Xsr) + mtrRr5 * (strYi - 8 - Ysr) + mtrRr8 * (strZi - 8 - Zsr);mtrA44 + 48 * i = (mtrRr7 * f + mtrRr8 * strYri - 8) / (mtrRr2 (strXi - 8 - Xsr) + mtrRr5 * (strYi - 8 - Ysr) + mtrRr8 * (strZi - 8 - Zsr);mtrA45 + 48 * i = -strXri - 8 * M

37、ath.Si n(omegar) - (strXri - 8 * (strXri - 8 * Math.Cos(kappar) - strYri - 8 * Math.Sin(kappar) / f - f * Math.Si n( kappar) * Math.Cos(omegar);mtrA46 + 48 * i = -f * Math.Cos(k appar) - strYri - 8 * (strXri - 8 Math.S in (kappar) + strYri - 8 * Math.Cos(k appar) / f;mtrA47 + 48 * i = -strXri - 8; /

38、将数组A转化为矩阵Classi A = new Class1(48, 6, mtrA);/求系数阵B double mtrB;mtrB = new double1728; for (int i = 0; i < 1728; i+) mtrBi = 0;/给B矩阵赋值for (int i = 0; i < 8; i+)mtrB0 + 147 * i = -mtrA0 + 48 * i; mtrB1 + 147 * i = -mtrA1 + 48 * i; mtrB2 + 147 * i = -mtrA2 + 48 * i; mtrB36 + 147 * i = -mtrA12 + 4

39、8 * i; mtrB37 + 147 * i = -mtrA13 + 48 * i; mtrB38 + 147 * i = -mtrA14 + 48 * i; mtrB72 + 147 * i = -mtrA30 + 48 * i; mtrB73 + 147 * i = -mtrA31 + 48 * i; mtrB74 + 147 * i = -mtrA32 + 48 * i; mtrB108 + 147 * i = -mtrA42 + 48 * i; mtrB109 + 147 * i = -mtrA43 + 48 * i; mtrB110 + 147 * i = -mtrA44 + 48

40、 * i;/转化矩阵BClassi B = new Class1(48, 36, mtrB);/求常数项double mtrL;mtrL = new double48;for (int i = 0; i < 8; i+)mtrL0 + 4 * i = xyz.GetElement(i, 0) + f * (mtrRI0 *(uvw.GetElement(0, i) - Xsl) + mtrRl3 * (uvw.GetElement(1, i) - Ysl) + mtrRl6 * (uvw.GetElement(2, i) - Zsl) / (mtrRl2 * (uvw.GetElemen

41、t(0, i) - Xsl) + mtrRl5 ” (uvw.GetElement(1, i) - Ysl) + mtrRl8 * (uvw.GetElement(2, i) - Zsl);mtrL1 + 4 * i = xyz. GetEleme nt(i, 1) + f * (mtrRl1 *(uvw.GetElement(0, i) - Xsl) + mtrRl4 * (uvw.GetElement(1, i) - Ysl) + mtrRl7 * (uvw.GetElement(2, i) - Zsl) / (mtrRl2 * (uvw.GetElement(0, i) - Xsl) +

42、 mtrRl5 ” (uvw.GetElement(1, i) - Ysl) + mtrRl8 * (uvw.GetElement(2, i) - Zsl);mtrL2 + 4 * i = XYZ.GetElement(i, 0) + f * (mtrRr0 *(uvw.GetElement(0, i) - Xsr) + mtrRr3 * (uvw.GetElement(1, i) - Ysr) + mtrRr6 (uvw.GetElement(2, i) - Zsr) / (mtrRr2 * (uvw.GetElement(0, i) - Xsr) + mtrRr5 (uvw.GetElem

43、ent(1, i) - Ysr) + mtrRr8 * (uvw.GetElement(2, i) - Zsr);mtrL3 + 4 * i = XYZ.GetEleme nt(i, 1) + f * (mtrRr1 *(uvw.GetElement(0, i) - Xsr) + mtrRr4 * (uvw.GetElement(1, i) - Ysr) + mtrRr7 (uvw.GetElement(2, i) - Zsr) / (mtrRr2 * (uvw.GetElement(0, i) - Xsr) + mtrRr5 (uvw.GetElement(1, i) - Ysr) + mt

44、rRr8 * (uvw.GetElement(2, i) - Zsr);for (int i = 8; i < 12; i+)mtrL0 + 4 * i = strXli - 8 + f * (mtrRl0 * (strXi - 8 - Xsl) + mtrRl3 * (strYi - 8 - Ysl) + mtrRl6 * (strZi - 8 - Zsl) / (mtrRl2 * (strXi - 8- Xsl) + mtrRl 5 * (strYi - 8 - Ysl) + mtrRl8 * (strZi - 8 - Zsl);mtrL1 + 4 * i = strYli - 8

45、+ f * (mtrRl1 * (strXi - 8 - Xsl) + mtrRl4 * (strYi - 8 - Ysl) + mtrRl7 * (strZi - 8 - Zsl) / (mtrRl2 * (strXi - 8- Xsl) + mtrRl 5 * (strYi - 8 - Ysl) + mtrRl8 * (strZi - 8 - Zsl);mtrL2 + 4 * i = strXri - 8 + f * (mtrRr0 * (strXi - 8 - Xsr) + mtrRr3 * (strYi - 8 - Ysr) + mtrRr6 * (strZi - 8 - Zsr) /

46、 (mtrRr2 * (strXi - 8 -Xsr) + mtrRr5 * (strYi - 8 - Ysr) + mtrRr8 * (strZi - 8 - Zsr);mtrL3 + 4 * i = strYri - 8 + f * (mtrRr1 * (strXi - 8 - Xsr) + mtrRr4 * (strYi - 8 - Ysr) + mtrRr7 * (strZi - 8 - Zsr) / (mtrRr2 * (strXi - 8 -Xsr) + mtrRr5 * (strYi - 8 - Ysr) + mtrRr8 * (strZi - 8 - Zsr); /转化矩阵LClassi LL = new Class1(48, 1, mtrL);/求改化法方程的系数Classi Nil = A.Tra nspose().Multiply(A); Classi N12 = A.Tra nspose().Multiply(B);Class1 N21 = B.Tra nspose().Multiply(A); Classi N22 = B.Tra nspose().Multiply(B);Classi ui = A.