数值计算方法课后实验-李维国

实验3.1Gauss消去法的数值稳定性实验实验目的：观察和理解高斯消元过程中出现小主元(即|a(k)kk|很小)时引起方程组解数值不稳定性．实验内容：求解方程组Ax=b,其中实验要求：（1）计算矩阵的条件数，判断系数矩阵是良态的还是病态的．（2）用高斯列主元消去法求得L和U及解向量x1,x2．（3）用不选主元的高斯消去法求得L和U及解向量x1,x2．（4）观察小主元并分析对计算结果的影响．解：（1）cond(A1)=68.4296cond(A2)=8.9939A1矩阵条件数远大于1，A1矩阵为病态;A2矩阵条件数没有远远大于1，A2矩阵为良态。（2）高斯列主元消去法：程序如下ClearA1=[0.3*power(10,-15),59.14,3,1;5.291,-6.130,-1,2;11.2,9,5,2;1,2,1,1];b1=[59.17;46.18;1;2];%A1=[1,2,-1,1;1,1,2,-1;3,-1,1,1;2,1,3,-1];%b1=[1;-2;6;-1];A2=[10,-7,0,1;-3,2.099999999999,6,2;5,-1,5,-1;0,1,0,2];b2=[8;5.900000000001;5;1];A=A1;b=b1;n=input('n=');fork=1:n-1[a,b3]=max(A(k:n,k));A([kk-1+b3],:)=A([k-1+b3k],:);b([kk-1+b3],:)=b([k-1+b3k],:);ifA(k,k)~=0A(k+1:n,k)=A(k+1:n,k)./A(k,k);A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n);elsestopendendA;b;L=tril(A,-1);U=triu(A,0);fori=1:1:nL(i,i)=1;endL;forj=1:n-1b(j)=b(j)/L(j,j);b(j+1:n)=b(j+1:n)-b(j)*L(j+1:n,j);endb(n)=b(n)/L(n,n);b;y=b;forj=n:-1:2y(j)=y(j)/U(j,j);y(1:j-1)=y(1:j-1)-y(j)*U(1:j-1,j);endy(1)=y(1)/U(1,1);结果如下:方程组一：L1：U1：解向量x1结果如下：方程组二结果如下：L2：U2:解向量x2结果如下：（3）不选主元的分解程序如下：ClearA1=[0.3*power(10,-15),59.14,3,1;5.291,-6.130,-1,2;11.2,9,5,2;1,2,1,1];b1=[59.17;46.18;1;2];%A1=[6,2,1,-1;2,4,1,0;1,1,4,-1;-1,0,-1,3];%b1=[6;-1;5;-5];A2=[10,-7,0,1;-3,2.099999999999,6,2;5,-1,5,-1;0,1,0,2];b2=[8;5.900000000001;5;1];A=A2;b=b2;n=input('n=');fork=1:n-1A(k+1:n,k)=A(k+1:n,k)/A(k,k);A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n);endL=tril(A,-1);U=triu(A,0);fori=1:1:nL(i,i)=1;endL;forj=1:n-1b(j)=b(j)/L(j,j);b(j+1:n)=b(j+1:n)-b(j)*L(j+1:n,j);endb(n)=b(n)/L(n,n);b;y=b;forj=n:-1:2y(j)=y(j)/U(j,j);y(1:j-1)=y(1:j-1)-y(j)*U(1:j-1,j);endy(1)=y(1)/U(1,1);方程组一结果：L1:U1:解向量x1结果如下：方程组二结果：L2:U2:解向量x2结果如下：(4)观察方程在两种不同方法下的结果可知：由于计算机字长是一定的，小主元会造成大数除以小数的结果超出字长，结果发生很大的变化。实验3.2方程组的性态和条件数实验实验目的：理解条件数的意义和方程组的性态对解向量的影响。实验内容：已知两个方程组A1x=b和A2x=b，其中实验要求：解：程序如下：Clearn=input('n=');A1=ones(n+1,n+1);A2=ones(n,n);fori=0:1:nforj=1:1:n+1A1(i+1,j)=power(1+0.1*i,j-1);endend%A1(2,2)=A1(2,2)+1e-14;%A1(6,6)=A1(6,6)+1e-14;fori=1:1:nforj=1:1:nA2(i,j)=(1/(i+j-1));endend%A2(2,2)=A2(2,2)+1e-7;A2(6,6)=A2(6,6)+1e-14;A1;A2;b1=sum(A1,2);b2=sum(A2,2);x1=A1\b1;x2=A2\b2;（1）求解第二范数下的条件数如下：N=4时，A1条件数为：A2条件数为：N=6时，A1条件数为：A2条件数为：N=8时，A1条件数为：A2条件数为：结论：随着n的增大，条件数也逐渐增大，矩阵病态更加严重；n=4,6,8时，矩阵三种条件数非常大，都为病态矩阵。（2）n=5时，解向量X1与X2结果如下：A1结果为A2结果为：（3）N=5时，对A1中a22项施加扰动1e-14，解向量X为：对A1中a66项施加扰动1e-14，解向量X为：对b元素施加扰动1e-4，解向量X为:（4）取n=6不变,对A2中a22和a66施加扰动1e-7，解向量分别为：(5)对A1,A2的元素项进行微小扰动时，对解向量无明显影响；对b元素进行微小扰动时，对解向量的影响较大;结论：对Ax=b进行扰动时，解向量更容易受b元素的扰动影响，不容易受A元素的扰动影响。（6）求~xxxN=5时，A1中a22和a66扰动对应的范数比分别为：实验3.3平方根法实验目的：掌握用平方根法、改进的平方根法求解系数矩阵正定的方程组的方法。实验内容：分别用平方根法、改进的平方根法求解方程组Ax=b，其中实验要求：（1）对于不同的n，编辑matlab程序求解方程组。（2）与数学软件自带函数求解计算结果对比。（3）与选主元的高斯消去法对比计算量。解（1）对于不同的n，编辑matlab程序求解方程组。程序如下：平方根法程序Clearn=input('n=');A=ones(n,n);fori=1:1:nforj=1:1:nA(i,j)=(1/(i+j-1));endendb=sum(A,2);A;b;%A=[4,-1,1;-1,4.25,2.75;1,2.75,3.5];%b=[6;-0.5;1.25];fork=1:nA(k,k)=sqrt(A(k,k));A(k+1:n,k)=A(k+1:n,k)/A(k,k);forj=k+1:nA(j:n,j)=A(j:n,j)-A(j:n,k)*A(j,k);endendL=tril(A,0);forj=1:n-1b(j)=b(j)/L(j,j);b(j+1:n)=b(j+1:n)-b(j)*L(j+1:n,j);endb(n)=b(n)/L(n,n);b;y=b;U=L';forj=n:-1:2y(j)=y(j)/U(j,j);y(1:j-1)=y(1:j-1)-y(j)*U(1:j-1,j);endy(1)=y(1)/U(1,1);改进的平方根法程序：ClearA=ones(n,n);n=input('n=');fori=1:1:nforj=1:1:nA(i,j)=(1/(i+j-1));endendb=sum(A,2);%A=[4,-2,4,2;-2,10,-2,-7;4,-2,8,4;2,-7,4,7];%b=[8;2;16;6];forj=1:1:nv=zeros(j-1,1);fori=1:j-1v(i)=A(j,i)*A(i,i);endv;A(j,j)=A(j,j)-A(j,1:j-1)*v(1:j-1);A(j+1:n,j)=(A(j+1:n,j)-A(j+1:n,1:j-1)*v(1:j-1))/A(j,j);endA;L=tril(A,0);D=diag(A);forj=1:nL(j,j)=1;endL;forj=1:n-1b(j)=b(j)/L(j,j);b(j+1:n)=b(j+1:n)-b(j)*L(j+1:n,j);endb(n)=b(n)/L(n,n);b;z=b;y=z./D;U=L';forj=n:-1:2y(j)=y(j)/U(j,j);y(1:j-1)=y(1:j-1)-y(j)*U(1:j-1,j);endy(1)=y(1)/U(1,1);（2）与数学软件自带函数求解计算结果对比。数学软件自带求解结果如下：X=linsolve（A,b）N=5时，N=6时，N=7时，而平方根法与改进的平方根法的结果都为与数学软件自带的求解结果不一致。（3）与选主元的高斯消去法对比计算量。全主元Gauss消去法计算量为1/3*n^3，列主元Gauss消去法计算量小于1/3*n^3。平方根法与改进平方根法计算量为1/3*n^3，平方根法避免了开方。由此看来列主元高斯消去法的计算量相对较小。