矩阵上机题(副)

矩阵与数值分析上机实习题王安然（21509201）工具：matlab1.（1）（2）Sn1=single(0);%存储(1)问的结果Sn2=single(0);%存储(2)问的结果N=10^2;%输入值分别为10^2,10^4,10^6fori=2:NSn1=Sn1+1/(i^2-1);end;fori=N:-1:2Sn2=Sn2+1/(i^2-1);end;（3）N分别取102,103,104结果如下（4）该题说明，改变元素和的相加顺序影响极限值。计算机中存在大小数相加的舍入误差误差。方法二的方法要比方法一要好，因为避免了大数吃小数的舍入误差。2.x=23;A=[7;3;-5;11];m=size(A);count=(m-1:-1:0);%多项式次数向量f=A'*x.^count';3.%initializationA=[31,-13,0,0,0,-10,0,0,0;-13,35,-9,0,-11,0,0,0,0;0,-9,31,-10,0,0,0,0,0;0,0,-10,79,-30,0,0,0,-9;0,0,0,-30,57,-7,0,-5,0;0,0,0,0,-7,47,-30,0,0;0,0,0,0,0,-30,41,0,0;0,0,0,0,-5,0,0,27,-2;0,0,0,-9,0,0,0,-2,29];b=[-15;27;-23;0;-20;12;-7;7;10];if(det(A)==0)disp('NoUniqueSolution');return;end;m=size(A,1);%GuassianEliminationaug=[A,b];fori=1:mforj=i+1:maug(j,:)=aug(j,:)-aug(j,i)/aug(i,i)*aug(i,:);end;end;disp(aug)x_GaussElim=solveFun(aug);fprintf('xforGuassianElimination\n');disp(x_GaussElim);%GuassianColumnPCAaug1=[A,b];fori=1:mmax=aug1(i,i);flag=i;fork=i+1:mif(aug1(k,i)^2max^2)max=aug1(k,i);flag=k;end;end;temp=aug1(flag,:);aug1(flag,:)=aug1(i,:);aug1(i,:)=temp;forj=i+1:maug1(j,:)=aug1(j,:)-aug1(j,i)/aug1(i,i)*aug1(i,:);end;end;x_GaussColu=solveFun(aug1);fprintf('xforGaussianColumnPCA\n');disp(x_GaussColu)其中解上三角方程组的函数solveFun（aug）:functionx=solveFun(ang)n=size(ang,1);x=zeros(n,1);A=ang(:,1:n);Y=ang(:,n+1);x(n,1)=Y(n,1)/A(n,n);fori=n-1:-1:1x(i,1)=(Y(i,1)-A(i,i+1:n)*x(i+1:n,1))/A(i,i);end;disp(x);4.%LUfunction[L,U]=LU(A)m=size(A,1);%LUA1=A;L=eye(m);fori=1:ml=eye(m);forj=i+1:ml(j,i)=-A1(j,i)/A1(i,i);end;A1=l*A1;L=L*l^(-1);end;U=A1;fprintf('L:\n');disp(L);fprintf('U:\n');disp(U);end%A_inversefunctionA_inverse=LU_inverse(A)[L,U]=LU(A);n=size(A,1);L_inv=zeros(n,n);U_inv=zeros(n,n);%inv(L)fori=1:1:nL_inv(i,i)=1/L(i,i);endfori=2:1:nforj=1:1:i-1L_inv(i,j)=-L_inv(i,i)*(L(i,j:(i-1))*L_inv(j:(i-1),j));endend%inv(U)fori=1:1:nU_inv(i,i)=1/U(i,i);endfori=1:1:n-1k=1;forj=i+1:1:nU_inv(k,j)=-U_inv(k,k)*(U(k,(k+1):j)*U_inv((k+1):j,j));k=k+1;endendA_inverse=L_inv*U_inv;fprintf('A_inverse\n');disp(A_inverse);end%A_detfunctionA_det=LU_det(A)m=size(A);[L,U]=LU(A);U_det=1;fori=1:mU_det=U_det*U(i,i);endA_det=U_det;fprintf('A_det\n');disp(A_det);end%LUColumnPCAfunction[L,U,P]=LU_ColumnPCA(A2)m=size(A2,1);P=eye(m);L=eye(m);p_cell=cell(1,m-1);l_cell=cell(1,m-1);fori=1:ml=eye(m);p=eye(m);max=A2(i,i);flag=i;forj=i+1:m%getmaxinthecolumnif(A2(j,i)^2max^2)%compareabsolutevaluemax=A2(j,i);flag=j;end;end;temp1=p(flag,:);p(flag,:)=p(i,:);p(i,:)=temp1;A2=p*A2;forj=i+1:ml(j,i)=-A2(j,i)/A2(i,i);end;A2=l*A2;p_cell{i}=p;l_cell{i}=l;P=p*P;end;U=A2;fprintf('P=\n');disp(P);fork=1:m-1L_temp=l_cell{k};forn=k+1:m-1L_temp=p_cell{n}*L_temp*p_cell{n}^(-1);end;L=L*L_temp^(-1);end;fprintf('L=\n');disp(L);fprintf('U=\n');disp(U);5.A=[2,1,-1,1;1,5,2,7;-1,2,10,1;1,7,1,11];b=[13,-9,6,0];disp(A);ifA~=(A')error('要求为对称正定阵');end;[m,n]=size(A);u=zeros(n,n);fori=1:nu(i,i)=A(i,i);fork=1:i-1u(i,i)=u(i,i)-u(k,i)^2;endu(i,i)=sqrt(u(i,i));forj=i+1:nu(i,j)=A(i,j);fork=1:i-1u(i,j)=u(i,j)-u(k,i)*u(k,j);end;u(i,j)=u(i,j)/u(i,i);end;end;L=u';disp(L)y=zeros(n,1);fori=1:n%计算ytemp=0;fork=1:i-1temp=temp+L(i,k)*y(k);end;y(i)=(b(i)-temp)/L(i,i);end;x=zeros(n,1);fori=n:-1:1%计算Xtemp=0;fork=i+1:ntemp=temp+L(k,i)*x(k);end;x(i)=(y(i)-temp)/L(i,i);end;disp(x);6.A=[4,4,0,0;3,-2,1,0;0,2,5,3;0,0,-1,6];b=[6;2;10;5];m=size(A,1);y=zeros(m,1);x=zeros(m,1);[L,U]=LU(A);y(1)=b(1);fori=2:my(i)=b(i)-L(i,i-1)*y(i-1);end;x(m)=y(m)/U(m,m);forj=m-1:-1:1x(j)=(y(j)-A(j,j+1)*x(j+1))/U(j,j);end;fprintf('x=\n');disp(x);7.function[Q,R]=QR(A)m=size(A);I=eye(m);Q=eye(m);R=A;fori=1:m-1q=zeros(m);w=R(i:end,i)-norm(R(i:end,i))*I(i:end,i);h=I(i:end,i:end)-2/(w'*w)*(w*w');q(i:end,i:end)=h;if(i1)forj=1:i-1q(j,j)=1;end;end;R=q*R;Q=Q*q';end;fprintf('R:\n');disp(R);fprintf('Q:\n');disp(Q)8.Jabobi方法function[x,i]=Jacobi(A,b)m=size(A,1);x=zeros(m,1);I=eye(m);D=I.*A;fprintf('D:\n');disp(D)L=D-tril(A);fprintf('L:\n');disp(L);U=D-triu(A);fprintf('U:\n');disp(U);i=1;while(true)x1=D^(-1)*(L+U)*x+D^(-1)*b;if(norm(x1-x,inf)=1e-4);break;end;x=x1;i=i+1;end;fprintf('x:\n');disp(x);fprintf('iteratetimes:\n');disp(i);GS方法function[x,i]=GS(A,b)m=size(A,1);x=zeros(m,1);I=eye(m);D=I.*A;disp(D)L=D-tril(A);disp(L);U=D-triu(A);disp(U);i=1;while(true)x1=(D-L)^(-1)*U*x+(D-L)^(-1)*b;if(norm(x1-x,inf)=1e-4);break;end;x=x1;i=i+1;end;fprintf('x:\n');disp(x);fprintf('iteratetimes:\n');disp(i);9.（1）2*x^3-5*x-1＝0；Newton法function[x1,i]=newton()symsx;f(x)=2*x^3-5*x-1;%fx0=1.5;i=0;while(true)x=x0;x1=x0-f(x)/eval(diff(f));if(abs(x0-x1)=1e-6);break;end;x0=x1;i=i+1;end;fprintf('x:\n');disp(vpa(x1,6));fprintf('iteratetimes:\n');disp(i);弦割法function[x2,i]=secant(f,x0,x1)i=0;while(true)x2=x1-f(x1)/(f(x1)-f(x0))*(x1-x0);if(abs(x2-x1)=1e-6)break;end;x0=x1;x1=x2;i=i+1;end;fprintf('x:\n');disp(x2);fprintf('iteratetimes:\n');disp(i);(2)exp(1)^x*sin(x)=0;Newton法f(x)=exp(1)^x*sin(x);x0=-3.5;弦割法10.二分法求f=x*cos(x)-2＝0；定义函数如下：functionf=myfun1(x)f=x*cos(x)-2;二分法函数：functionx=bisection(f,a,b)i=0;e=b-a;while(e=1e-6)x=(a+b)/2;if(f(x)*f(a)0)%disp(f(x));a=x;elseif(f(x)*f(a)0)b=x;elsea=x;b=x;end;e=e/2;i=i+1;end;fprintf('x