数值分析上机实验3

1、数值分析上机实验3例1利用不选主元的高斯消元法实现矩阵的三角分解：MATLAB程序如下：A=2,3,4;3,5,2;4,3,30 n=max(size(A);for k=1:n-1 d=A(k,k); for i=k+1:n m=A(i,k)/d; for j=k+1:n A(i,j)=A(i,j)-m*A(k,j); end A(i,k)=m; endend A程序运行后，输出结果如下 2 3 4 3/2 1/2 -4 2 -6 -2 所以，矩阵A的三角分解为例2将平面向量旋转所用的数学工具是二阶正交矩阵当矩阵A作用于向量a 时，有b = Aa ，则b 是由向量a 旋转所得，旋转角为（逆时针

2、旋转为正）。设有中心在原点，边长为2的正方形B0，将正方形旋转300，并适当缩小得另一正方形B1，如此操作六次，并将前后七个正方形绘图，程序如下：xy=-1 -1;1 -1;1 1;-1 1;-1 -1;A=cos(pi/6) -sin(pi/6);sin(pi/6) cos(pi/6);x=xy(:,1);y=xy(:,2);axis offline(x,y),pause(1)for k=1:6 xy=.73*xy*A' x=xy(:,1);y=xy(:,2); line(x,y),pause(1)end例3在某一时刻T，根据GPS的卫星广播星历计算得4颗卫星坐标分别为x 1，y 1

3、，z 1=15524471.175，-16649826.222，13512272.387x 2，y2，z 2=-2304058.534，-23287906.465，11917038.105x 3，y 3，z 3=16680243.357，-3069625.561，20378551.047x 4，y 4，z 4=-14799931.395，-21425358.240，6069947.224已测得4颗卫星在同一时刻T，到GPS接收机的距离观测值分别为：r 1=22274133.8630r 2=19920232.728r 3=24566287.852r 4=21491226.260试列出方程组并计算

4、接收机的概略位置。解：设接收机的位置为（x，y，z），由条件得将等式两端平方，得其中，固定第一式，由第二、第三、第四式分别减去第一式，得线性方程组经计算，得接收机的概略位置为（x，y，z）=（-7.337313000e+005，-5.447498718e+006，3.234494646e+006）MATLAB程序如下：S1=15524471.175,-16649826.222,13512272.387;S2=-2304058.534,-23287906.465,11917038.105;S3=16680243.357,-3069625.561,20378551.047;S4=-14799931

5、.395,-21425358.240,6069947.224;Ro1=22274133.8630;Ro2=19920232.728;Ro3=24566287.852;Ro4=21491226.260;A=S2-S1;S3-S1;S4-S1;L=.5*(Ro12-Ro22;Ro32;Ro42)-(norm(S1)2-norm(S2)2;norm(S3)2;norm(S4)2);format long eP=AL程序运行结果为P = -7.337313000858261e+005 -5.447498718473805e+006 3.234494646022576e+006例4 用牛顿迭代法求解非线

6、性方程组解：显然，这两个方程中一个是抛物线，另一个是圆。令利用牛顿迭代法计算时，假设已有一个解的初始近似值为（x0，y0）利用台劳展开式可得线性方程组，其中，近似值修正为：，。将这一过程反复进行，就是牛顿迭代算法。取初值为（x0，y0）=（1，0），用牛顿迭代法计算四次数据如下kxkyk01011.06250.125021.06730.139131.06730.139241.06730.1392程序如下：x=input('input x0=');y=input('input y0=');X=x;y;for k=1:6 G=2*x,-1;2*x-4,2*y-1;

7、F=x2-y-1;(x-2)2+(y-0.5)2-1; dX=GF; X=X-dX; x=X(1);y=X(2);X'end例5 最短路问题的动态规划算法。在一个m×n棋盘内（m，n均不超过100），每个格子内有一个正整数值（不超过100）表示占据该格子应支付的费用。一个国际象棋的武士从棋盘左下角格子（1，1）开始沿着向右或向上的方向向右上角格子（m，n）行进，要求找一条行进路径使武士支付的费用最小。以下实例是一个4×4棋盘，最小费用为15。2821223185492157解：由于武士只能向上或向右走。所以，武士在格子（1，1）处的行进路径方向取决于格子（1，2）和

8、格子（2，1）到格子（4，4）支付的最小费用。而武士在这两个格子处行进的方向又取决于该格子右边和上边格子到格子（4，4）所支付的费用，。设武士占据格子（i，j）的费用为w(i，j)，记d(i，j)为格子（i，j）到格子（4，4）的最小费用。则d(i，j) = w(i，j) + min(d(i+1，j)，d(i，j+1)这种分阶段计算的方法就是动态规划算法。程序如下w=2 1 5 7;8 5 4 9;2 2 3 1;2 8 2 1;d(4,4)=w(4,4);for j=3:-1:1 d(4,j)=w(4,j)+d(4,j+1);R(4,j)=1;endfor i=3:-1:1 d(i,4)=w

9、(i,4)+d(i+1,4);R(i,4)=0;endfor j=3:-1:1 for i=3:-1:1 r=d(i,j+1);u=d(i+1,j); if r<=u d(i,j)=w(i,j)+r;R(i,j)=1; else d(i,j)=w(i,j)+u; end endendi=1;j=1;p(1,:)=i,j;for k=2:7 if R(i,j)=1 j=j+1; else i=i+1; end p(k,:)=i,j;endd(1,1),p运行后，输出结果如下ans= 15p = 1 1 1 2 2 2 3 2 3 3 3 44 4这说明最小费用为15，路线为：（1，1）&#

10、224;（1，2）à（2，2）à（3，2）à（3，3）à（3，4）à（4，4）上机实习题一、数值实验1将平面旋转方法推广到空间，考虑向量在三个坐标面内的旋转。对应于三个正交矩阵，第一个正交矩阵用于将Y-Z平面上向量以X轴为旋转轴左旋a 角；第二个正交矩阵用于将Z-X平面上向量以Y轴为旋转轴左旋b 角，第三个正交矩阵用于将X-Y平面上向量以Z轴为旋转轴左旋g 角。现有位于第一象限内的单位立方体B1。首先绕Z轴将立方体左旋900，立方体转动到第二象限为B2；然后，再绕Z轴左旋900，立方体转动到第三象限为B3；最后，将立方体再绕Y轴右旋450，且绕

11、X轴左旋450为B4，绘出旋转前后的四个立方体形状。立方体B1可用八个顶点坐标描述 0 0 0 1 0 0 1 1 0 0 1 0 0 0 1 1 0 1 1 1 1 0 1 1例如，经第一次正交变换后，得八个顶点新的坐标。B2 = 0 0 0 0.0000 1.0000 0 -1.0000 1.0000 0 -1.0000 0.0000 0 0 0 1.0000 0.0000 1.0000 1.0000 -1.0000 1.0000 1.0000绘制原始正立方体参考程序如下：B=0 0 0;1 0 0;1 1 0;0 1 0;0 0 1;1 0 1;1 1 1;0 1 1;B0=B;p1=B

12、(1,:);p2=B(2,:);p3=B(3,:);p4=B(4,:);p5=B(5,:);p6=B(6,:);p7=B(7,:);p8=B(8,:);A=p1;p2;p3;p4;p1;p5;p6;p7;p8;p5;p6;p2;p3;p7;p8;p4;plot3(A(:,1),A(:,2),A(:,3);hold on2希尔伯特矩阵的病态性问题n阶Hilbert矩阵是一个对称正定矩阵，而且条件数 Cond(Hn) 随着n的增大而迅速增加。其逆矩阵为，这里（1）取某一种范数（例如用-范数）对于 n = 4，5，6，12，计算 ln(Cond(Hn) 的值，并绘出该对数函数的曲线。猜测ln(Con

13、d(Hn)与n之间有何种关系，并给予验证。（2）设D是由Hn的对角线元素开方构成的对角矩阵，令，不难看出该矩阵仍然是对称正定矩阵，而且其对角元素全为1。将Hn 变换为的技术称为预处理（Preconditioning）。试分析函数随n变化的规律，绘制该函数的曲线。3对下面的9阶对称矩阵分别做LU分解和LDLT分解并利用MATLAB的spy(A)命令，绘制矩阵的A非零元素结构图以及分解后的矩阵的非零元素结构图4三对角方程组求解。对如下n阶三对角矩阵A，以及n维列向量b， 用追赶法求解方程组Ax=b5在迈阿密附近的弗罗里达海湾，GPS用户接收机收到八颗GPS卫星信号，根据卫星星历计算出GPS卫星的位

14、置，根据伪码测距获得伪距。数据如下IDXYZr GPS-116414028.668660383.61820932036.90724658975.31743GPS-216896800.648-18784061.365-7418318.85622964286.41228GPS-39339639.616-14514964.65820305107.16121338550.64536GPS-4-18335582.591-11640868.30515028599.07123606547.29359GPS-5-2077142.705-20987755.987-15879741.19624263298.504

15、01GPS-6-4957166.885-23306741.03912039027.09620758264.10823GPS-717977519.820-13089823.31214331151.06521847468.81689GPS-89682727.508-24060519.4853985404.53020352077.19349考虑伪距定位模型用牛顿迭代法计算接收机的位置，在站心坐标系（接收机为空间直角坐标原点）中计算上述八颗卫星的的方位角和高度角，绘制如下图形。二、方法应用与编程练习1利用部分选主元的列主元高斯消元法实现矩阵的三角分解：2 用全部选主元的消元法求解下列方程组3 折叠矩阵

16、算法起始数位于方阵左上角，然后从起始数开始递增，层层折叠展开，最后排列为方阵。如下所示为起始数为10，行数是4的折叠方阵10 11 141913 12 15 2018 17 16 2125 24 23 22设计程序，对输入的起始数a和方阵阶数n，输出折叠方阵；并查询：第x行第y列的位置上的数，以及指定的数m在方阵中的位置，如果m不在方阵中，则需输出 “Error”。4 旋转方阵算法把整数1，2，3，n2 从矩阵的第一行第一列开始，从外层至中心按顺时针方向螺旋排列所成的n×n 方阵，称为n阶旋转方阵。如4阶旋转方阵为1 2 3 412 13 14 51116 15 610 9 8 7对输入的正整数n，输出旋转方阵。5 奇数阶幻方（Magic square）算法我国古代有关幻方的论述是3阶纵横图：九子斜排，上下对易，四维挺出，载九履一，左三右七，二四为肩，六八为足。492357816关于奇数阶幻方最早的一个算法据说是1687年从泰国传到法国，其大意是：把1填入第一行中间格子，然后，垂直向下到最后一行向右走一格将2填入，以后按对角线向右上方向行走依次填写3，4，n2。如果走到第n列，则将第一列做为下一列；如果走到第一行，则将第n行做下一行。如果遇到已经填有数的格子则下移一行继续行走，直到n2结束