2016级矩阵与数值分析上机作业(共72页).docx

《2016级矩阵与数值分析上机作业(共72页).docx》由会员分享，可在线阅读，更多相关《2016级矩阵与数值分析上机作业(共72页).docx（72页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、精选优质文档-倾情为你奉上 2016级矩阵与数值分析上机作业 学生班级： 学生姓名： 任课教师： 所在学院：电子信息与电气学部 学生学号： 使用软件：MATLAB2016年12月10号1. 考虑计算给定向量的范数：输入向量，输出，请编制一个通用程序，并用你编制的程序计算如下向量的范数： ， 对n = 10，100，1000甚至更大的n计算其范数，你会发现什么结果？你能否修改你的程序使得计算结果相对精确呢？（1）计算范数的程序function Nor1,Nor2,Nor3=Nor(a) b=length(a);format long Nor11=0; format longNor21=0;i=1

2、;while it t=abs(a(i);end i=i+1; endNor1=Nor11; %计算1范数Nor2=sqrt(Nor21); %计算2范数Nor3=t; %计算无穷范数end(2)x与y向量的程序function a= afun(n)a=zeros(n);a(1)=1;for i=1:1:n; a(i)=1/i;endendfunction b= bfun(n)b=zeros(n);b(1)=1;for i=1:1:n; b(i)=i;endend运行： n=10; a=afun(n); b=bfun(n); f1,f2,f3=Nor(a); h1,h2,h3=Nor(b);f

3、1 =2.8254 f2 =1.5769 f3 = 1h1 =55 h2 =19.8583 h3 = 10 n=100;a=afun(n);b=bfun(n);f1,f2,f3=Nor(a); h1,h2,h3=Nor(b);f1 = 5.9621 f2 =1.3052 f3 = 1h1 = 5050 h2 = 5.1153e+02 h3 =100 n=1000;a=afun(n);b=bfun(n);f1,f2,f3=Nor(a); h1,h2,h3=Nor(b);f1 = 7.0343 f2 =1.1847 f3 = 1h1 = h2 = 1.2642e+04 h3 = 1000 上述结果

4、由于MATLAB的高精度运算较为准确，但是由于当n非常大时可能会出现了大数吃小数的现象使误差增大，而y向量的范数结果较为准确是由于y向量是按从小到大排列。改进范数计算程序使输入序列按从小到大排列再计算其范数：function Nor1,Nor2,Nor3=Nor(c) a=sort(c);b=length(a);format long Nor11=0; format longNor21=0;i=1;while it t=abs(a(i);end i=i+1; endNor1=Nor11; %计算1范数Nor2=sqrt(Nor21); %计算2范数Nor3=t; %计算无穷范数end2.考虑，

5、其中定义，此时是连续函数。用此公式计算当时的函数值，画出图像。另一方面，考虑下面算法： 用此算法计算x1015,1015时的函数值，画出图像。比较一下发生了什么？函数：function F=Piecewise1_x(x)F= log(x+1)/x*(x=0&x0);endfunction F=Piecewise2_x(x)F= log(x)/(x-1)*(x=1&x1);end运行：clcclearx=linspace(-10(-15),10(-15); d=x+1;%原来图像（红色曲线）F=Piecewise1_x(x);%计算相应函数值 plot(x,F,r);%绘制曲线 hold on;

6、 %改变算法后的图像（蓝色曲线）F1=Piecewise2_x(d);%计算相应函数值 plot(x,F1,b);%绘制曲线 hold on; plot(0*ones(1,2),ylim,g:);%画区间间隔线plot(xlim,1*ones(1,2),g:);%画区间间隔线 xlabel(变量 X)ylabel(变量 Y1 & Y2)由上图可知改变后的算法画出的曲线（蓝色）比直接画出的曲线（红色）更贴近实际值误差更小，且由于ln（x+1)j A(i,j)=-1; elseif (ij)&(jj A(i,j)=-1; elseif (ij)&(jj A(i,j)=-1; elseif (ij)

7、&(jn) A(i,j)=0; else A(i,j)=1; end endende=ones(n,n);for k=1:n b(k)=e(:,k);P,L,U=PLU(A);b1=b(k);b=P*b1;y=zeros(n,1);y(1)=b(1); for i=2:n y(i)=b(i)-sum(L(i,1:i-1).*y(1:i-1) ; end% yx(n)=y(n)/U(n,n);for i=n-1:-1:1 x(i)=(y(i)-sum(U(i,i+1:n).*x(i+1:n)/U(i,i);endB(:,k)=x;endP=A*B %验证是否正确% x%d=norm(x1-x);

8、 %计算算法误差endclearfor n=5:30Bn-4=zeros(n);Bn-4= PLU2(n)end求出的n从5到30的逆矩阵存储在B中： B = Columns 1 through 5 5x5 double 6x6 double 7x7 double 8x8 double 9x9 double Columns 6 through 9 10x10 double 11x11 double 12x12 double 13x13 double Columns 10 through 13 14x14 double 15x15 double 16x16 double 17x17 double

9、 Columns 14 through 17 18x18 double 19x19 double 20x20 double 21x21 double Columns 18 through 21 22x22 double 23x23 double 24x24 double 25x25 double Columns 22 through 25 26x26 double 27x27 double 28x28 double 29x29 double Column 26 30x30 double5. 编制计算对称正定阵的Cholesky分解的通用程序，并用你编制的程序计算Ax=b，其中。b可以由你自己取

10、定，对n从10到20验证程序的可靠性。% Cholesky分解的通用程序function L = Cholesk(A)%如果是正定矩阵，可以进行下面分解操作L(1,1)=sqrt(A(1,1);%确定第一列元素n=length(A);for i=2:n L(i,1)=A(i,1)/L(1,1);endsum1=0;for j=2:n-1 %确定第j列元素 for k=1:j-1 sum1=sum1+L(j,k)2; end L(j,j)=sqrt(A(j,j)-sum1); sum1=0; for i=j+1:n for k=1:j-1 sum1=sum1+L(i,k)*L(j,k); end

11、 L(i,j)=(A(i,j)-sum1)/L(j,j); end sum1=0;endfor k=1:n-1 %确定第n列元素 sum1=sum1+L(n,k)2;endL(n,n)=sqrt(A(n,n)-sum1);P=L*L ; %验证分解是否正确end%Cholesky分解的通用程序用于生成A,b矩阵的函数程序function x1,b,A = Ch(n) x1=zeros(n,1);for k=1:n x1(k)=1;endA=zeros(n,n);for i=1:n for j=1:n A(i,j)=1/(i+j-1); endendb=A*x1;end% 计算Cholesky分

12、解的通用程序的误差程序function e = Cho(n) x1,b,A = Ch( n );m,n=size(A);if m=n %判断输入的矩阵是不是方阵 e=inf; return;endd=eig(A);%根据方阵的特征值判定是不是正定矩阵for i=1:n if d(i) cleard=zeros(19,1);for n=2:20 d(n-1)= Cho(n);enddd = 1.0e+17 * 0.0000 0.0000 0.0000 0.0042 0.0000 0.6285 0.0000 0.1891 0.7943 1.0288 0.4207 Inf 0.2969 Inf In

13、f Inf Inf Inf Inf 由上面程序运算结果是n从220的的程序误差结果可知，低阶是程序可靠性较好，随着n的增大程序的误差总体越来越大，（运行结果后面输出为inf（程序设定）的n阶矩阵不是正定矩阵）程序的可靠性变差。6. (1)编制程序，其作用是对输入的向量x，输出单位向量u使得 。(2)编制Householder变换阵乘以的程序HA，注意，你的程序并不显式的计算出。(3)考虑矩阵 用你编制的程序计算使得HA的第一列为的形式，并将HA的结果显示。(1)function u= house(x) e=eye(length(x); w=x-norm(x)*e(:,1); u=w./norm

14、(w);end(2)function H ,HA= HA( x,A )u=house(x);H=eye(length(u)-2*u*u;HA=H*A;end(3)clearA=1,2,3,4;-1,3,2(1/2),3(1/2);-2,2,exp(1),pi;-10(1/2),2,-3,7;0,2,7,2.5 ;x=A(:,1);H,HA=HA(x,A) H = 0.2500 -0.2500 -0.5000 -0.7906 0 -0.2500 0.9167 -0.1667 -0.2635 0 -0.5000 -0.1667 0.6667 -0.5270 0 -0.7906 -0.2635 -0

15、.5270 0.1667 0 0 0 0 0 1.0000 HA = 4.0000 -2.8311 1.4090 -6.5378 -0.0000 1.3896 0.8839 -1.7805 -0.0000 -1.2208 1.6576 -3.8836 -0.0000 -3.0925 -4.6770 -4.1078 0 2.0000 7.0000 2.50007. 用Jacobi和Gauss-Seidel迭代求解下面的方程组，输出迭代每一步的误差： 假设迭代法计算停止的条件为：（1）Jacobi迭代法clc;clear;A=5,-1,-3;-1,2,4;-3,4,15 ;b=-2,1,10;de

16、lta=10(-4);%误差n=length(A);k=0;x=zeros(n,1);y=zeros(n,1);while 1 for i=1:n y(i)=b(i); for j=1:n if j=i y(i)=y(i)-A(i,j)*x(j); end end y(i)=y(i)/A(i,i); %迭代出的xendd=norm(y-x);if ddelta %迭代停止条件 break; end x=y;k=k+1;B(1,k)=k;B(2,k)=d;endBy运行结果：（B的第一行是迭代次数k，第二行是每一步迭代误差） B = Columns 1 through 8 1.0000 2.00

17、00 3.0000 4.0000 5.0000 6.0000 7.0000 8.0000 0.9244 1.6268 0.9517 1.3383 0.9542 1.1490 0.9261 1.0146 Columns 9 through 16 9.0000 10.0000 11.0000 12.0000 13.0000 14.0000 15.0000 16.0000 0.8794 0.9114 0.8242 0.8269 0.7665 0.7547 0.7097 0.6912 Columns 17 through 24 17.0000 18.0000 19.0000 20.0000 21.00

18、00 22.0000 23.0000 24.0000 0.6554 0.6342 0.6043 0.5826 0.5567 0.5356 0.5126 0.4926 Columns 25 through 32 25.0000 26.0000 27.0000 28.0000 29.0000 30.0000 31.0000 32.0000 0.4719 0.4531 0.4343 0.4168 0.3996 0.3835 0.3678 0.3528 Columns 33 through 40 33.0000 34.0000 35.0000 36.0000 37.0000 38.0000 39.00

19、00 40.0000 0.3384 0.3246 0.3114 0.2987 0.2865 0.2748 0.2636 0.2529 Columns 41 through 48 41.0000 42.0000 43.0000 44.0000 45.0000 46.0000 47.0000 48.0000 0.2426 0.2327 0.2232 0.2141 0.2054 0.1970 0.1890 0.1813 Columns 49 through 56 49.0000 50.0000 51.0000 52.0000 53.0000 54.0000 55.0000 56.0000 0.173

20、9 0.1668 0.1600 0.1535 0.1472 0.1412 0.1355 0.1299 Columns 57 through 64 57.0000 58.0000 59.0000 60.0000 61.0000 62.0000 63.0000 64.0000 0.1246 0.1196 0.1147 0.1100 0.1055 0.1012 0.0971 0.0931 Columns 65 through 72 65.0000 66.0000 67.0000 68.0000 69.0000 70.0000 71.0000 72.0000 0.0893 0.0857 0.0822

21、0.0789 0.0756 0.0726 0.0696 0.0668 Columns 73 through 80 73.0000 74.0000 75.0000 76.0000 77.0000 78.0000 79.0000 80.0000 0.0640 0.0614 0.0589 0.0565 0.0542 0.0520 0.0499 0.0479 Columns 81 through 88 81.0000 82.0000 83.0000 84.0000 85.0000 86.0000 87.0000 88.0000 0.0459 0.0440 0.0422 0.0405 0.0389 0.

22、0373 0.0358 0.0343 Columns 89 through 96 89.0000 90.0000 91.0000 92.0000 93.0000 94.0000 95.0000 96.0000 0.0329 0.0316 0.0303 0.0290 0.0279 0.0267 0.0256 0.0246 Columns 97 through 104 97.0000 98.0000 99.0000 100.0000 101.0000 102.0000 103.0000 104.0000 0.0236 0.0226 0.0217 0.0208 0.0200 0.0192 0.018

23、4 0.0176 Columns 105 through 112 105.0000 106.0000 107.0000 108.0000 109.0000 110.0000 111.0000 112.0000 0.0169 0.0162 0.0156 0.0149 0.0143 0.0137 0.0132 0.0126 Columns 113 through 120 113.0000 114.0000 115.0000 116.0000 117.0000 118.0000 119.0000 120.0000 0.0121 0.0116 0.0112 0.0107 0.0103 0.0098 0

24、.0094 0.0091 Columns 121 through 128 121.0000 122.0000 123.0000 124.0000 125.0000 126.0000 127.0000 128.0000 0.0087 0.0083 0.0080 0.0077 0.0074 0.0071 0.0068 0.0065 Columns 129 through 136 129.0000 130.0000 131.0000 132.0000 133.0000 134.0000 135.0000 136.0000 0.0062 0.0060 0.0057 0.0055 0.0053 0.00

25、51 0.0049 0.0047 Columns 137 through 144 137.0000 138.0000 139.0000 140.0000 141.0000 142.0000 143.0000 144.0000 0.0045 0.0043 0.0041 0.0039 0.0038 0.0036 0.0035 0.0033 Columns 145 through 152 145.0000 146.0000 147.0000 148.0000 149.0000 150.0000 151.0000 152.0000 0.0032 0.0031 0.0029 0.0028 0.0027