数值分析上机作业

1、数值实验数值实验综述：线性代数方程组的解法是一切科学计算的基础与核心问题。求解方法大致可分为直接法和迭代法两大类。直接法指在没有舍入误差的情况下经过有限次运算可求得方程组的精确解的方法，因此也称为精确法。当系数矩阵是方的、稠密的、无任何特殊结构的中小规模线性方程组时，Gauss消去法是目前最基本和常用的方法。如若系数矩阵具有某种特殊形式，则为了尽可能地减少计算量与存储量，需采用其他专门的方法来求解。Gauss消去等同于矩阵的三角分解，但它存在潜在的不稳定性，故需要选主元素。对正定对称矩阵，采用平方根方法无需选主元。方程组的性态与方程组的条件数有关，对于病态的方程组必须采用特殊的方法进行求解。实

2、验一一、 实验名称：矩阵的LU分解.二、 实验目的：用不选主元的LU分解和列主元LU分解求解线性方程组 Ax=b, 并比较这两种方法.三、 实验内容与要求（1）用所熟悉的计算机语言将不选主元和列主元LU分解编成通用的子程序，然后用编写的程序求解下面的84阶方程组将计算结果与方程组的精确解进行比较，并就此谈谈你对Gauss消去法的看法。（2）写出追赶法求解三对角方程组的过程，并编写程序求该实验中的方程组（1）不选主元高斯消去法求解方程组function L,U=gauss1(A,b)n=length(A);for k=1:n-1 A(k+1:n,k)=A(k+1:n,k)/A(k,k); A(k

3、+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:n),-1)+eye(n);U=triu(A);主程序为：Clear；clc； a=ones(84,1)*6; b=ones(83,1); c=ones(83,1)*8;A=diag(a)+diag(b,1)+diag(c,-1); d=ones(82,1)*15;b=7;d;14;L U=gauss1(A,b);n=length(A);y=ones(n,1);y=Lb;x=ones(n,1);x=Uy解为：x=11.0001.0001.0000.99981.000

4、0.99931.0010.99721.0060.98861.0230.95451.0910.81811.3640.27251.4550.08981.8200.35921.2820.43671.1270.74691.5060.98741.0250.94971.1010.79861.4030.19461.6110.77831.4440.11301.7740.45201.0960.80811.3840.23231.5350.92911.1410.71871.5570.90881.0920.18071.1830.45351.811-0.49614.486-6.94716.79-30.1864.76-1

5、26.1256.4-510.61024.-2046.4096.-8190.16383.-32764.65531.-.-.-.-.-.-.-.可见，这是一个病态方程，从56个跟开始发散。可见|A|=7.422e+50不等于零，所以是非奇异矩阵，这主要是一个小主元问题引起的，需要更高精度的计算机来计算。选主元的高斯消去法function L,U,P=gauss2(A,b)n,l=size(A);P=eye(n);for k=1:n-1 m,p=max(abs(A(:,k); A(k,p,:)=A(p,k,:); P(k,p,:)=P(p,k,:); if A(k,k)=0 A(k+1:n,k)=A

6、(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); endendL=tril(A(:,1:n),-1)+eye(n);U=triu(A);P;主程序为：clearclc a=ones(84,1)*6; b=ones(83,1); c=ones(83,1)*8;A=diag(a)+diag(b,1)+diag(c,-1); d=ones(82,1)*15;b=7;d;14;L U P=gauss2(A,b);n=length(A);y=ones(n,1);y=L(P*b);x=ones(n,1);x=Uy解

7、为：1.0001.0001.0001.0001.0001.0001.0001.0001.0001.0001.0001.0001.0001.0001.0001.0001.0001.0001.0001.0001.0001.0001.0001.0001.0001.0001.0001.0001.0001.0001.0001.0001.0001.0001.0001.0001.0001.0001.0001.0001.0001.0001.0001.0001.0001.0001.0001.0001.0001.0001.0000.99991.0000.99951.0010.99791.0040.99171.01

8、70.96671.0670.86661.2670.46641.0670.86571.2690.46291.0740.85261.2930.42271.1210.89091.6850.76491.9330.28241.8410.69361.1090.79611.3500.5339这个病态方程的解没有发散，用选主元的方法可以有效的提高解的精度。（2）追赶法clearclc a=ones(84,1)*6; e=ones(84,1); c=ones(84,1)*8; d=ones(82,1)*15;b=7;d;14;n,l=size(b);tic;U(1)=a(1);for i=2:n L(i)=c(

9、i)/U(i-1); U(i)=a(i)-L(i)*e(i-1);endfor i=2:n y(1)=b(1); y(i)=b(i)-L(i)*y(i-1);endfor i=n-1:-1:1 x(n)=y(n)/U(n); x(i)=(y(i)-e(i)*x(i+1)/U(i);endt3=tocx=x解为：11.0001.0001.0000.99981.0000.99931.0010.99721.0060.98861.0230.95451.0910.81811.3640.27251.4550.08981.8200.35921.2820.43671.1270.74691.5060.98741

10、.0250.94971.1010.79861.4030.19461.6110.77831.4440.11301.7740.45201.0960.80811.3840.23231.5350.92911.1410.71871.5570.90881.0920.18071.1830.45351.811-0.49614.486-6.94716.79-30.1864.76-126.1256.4-510.61024.-2046.4096.-8190.16383.-32764.65531.-.-.-.-.-.-.-.实验二一、 实验名称：实对称正定矩阵的的Cholesky分解.二、 实验目的：用平方根法和改进

11、的平方根方法求解线性方程组 Ax=b.三、 实验内容与要求用所熟悉的计算机语言将Cholesky分解和改进的Cholesky分解编成通用的子程序，然后用编写的程序求解对称正定方程组Ax=b，其中（1） b随机的选取，系数矩阵为100阶矩阵（2） 系数矩阵为40阶Hilbert矩阵，即系数矩阵A的第i行第j列元素为，向量b的第i个分量为（3） 用实验一的程序求解这两个方程组，并比较所有的计算结果，然后评价各个方法的优劣。平方根方法求解function L=cholesky1(A) n=length(A);tic;for k=1:n A(k,k)=sqrt(A(k,k); A(k+1:n,k)=A

12、(k+1:n,k)/A(k,k); for j=k+1:n A(j:n,j)=A(j:n,j)-A(j:n,k)*A(j,k); endendt1=tocL=tril(A);改进的平方根方法function L,D=cholesky2(A)n=length(A);v=ones(n,1);for j=1:n for i=1:j v(i)=A(j,i)*A(i,i); end A(j,j)=A(j,j)-A(j,1:j-1)*v(1:j-1,1); A(j+1:n,j)=(A(j+1:n,j)-A(j+1:n,1:j-1)*v(1:j-1,1)/A(j,j);endD=diag(diag(A);L

13、=tril(A(:,1:n),-1)+eye(n);使用两种方法求解对称正定方程（1）系数矩阵为100阶矩阵clearclc a=ones(100,1)*10; b=ones(99,1); c=ones(99,1);A=diag(a)+diag(b,1)+diag(c,-1);n=length(A);b=rand(n,1);L=cholesky1(A);y=ones(n,1);y=Lb;x=ones(n,1);x=Ly% L D=cholesky2(A);% z=ones(n,1);% z=Lb;% y=ones(n,1);% y=Dz;% x=ones(n,1);% x=Ly（2）系数矩阵为

14、40阶Hilbert矩阵clearclcn=40;for i=1:n for j=1:n A(i,j)=1/(i+j-1); endendfor i=1:n b(i,1)=sum(A(i,1:i);endL=cholesky1(A);y=ones(n,1);y=Lb;x=ones(n,1);x=Ly% L D=cholesky2(A);% z=ones(n,1);% z=Lb;% y=ones(n,1);% y=Dz;% x=ones(n,1);% x=Ly（3）使用两种方法求解对称正定方程实验一的程序求解这两个方程组系数矩阵为100阶矩阵，b随机的选取clearclc a=ones(100,

15、1)*10; b=ones(99,1); c=ones(99,1);A=diag(a)+diag(b,1)+diag(c,-1);n=length(A);b=rand(n,1);L1,U1=gauss1(A,b);y1=ones(n,1);y1=L1b;x1=ones(n,1);x1=U1y1; L2,U2,P=gauss2(A,b);y2=ones(n,1);y2=L2(P*b);x2=ones(n,1);x2=U2y2;x=x1,x2解为：0.538770.538770.513010.513010.819580.819580.957500.957500.128230.128230.3441

16、80.344180.172070.172070.652920.65292-0.74234-0.742340.464980.464980.395010.395010.984790.984790.640140.640140.458030.458030.801430.801430.722940.722940.573320.573320.182870.182870.798440.798440.851450.851450.090680.090680.781950.781950.737390.737390.441670.44167-0.-0.0.893840.893840.427210.427210.56

17、8290.568290.612010.612010.843850.843850.310140.310140.644340.644340.786240.786240.912060.91206-0.-0.0.0.0.507620.507620.449680.449680.205870.205870.0.0.707580.707580.683680.68368-0.-0.0.595830.595830.858040.85804-0.-0.0.058190.058190.021270.021270.0.0.207750.207750.753940.753940.741830.74183-0.-0.0.

18、481230.481230.086810.086810.693710.693710.510190.510190.637640.637640.372370.372370.749940.749940.0.0.732890.732890.740490.740490.943840.943840.994490.994490.380980.380980.327060.327060.168490.168490.263840.263840.0.0.143800.14380-0.-0.0.0.0.208110.208110.0.0.813640.813640.732450.73245-0.-0.0.037610

19、.037610.467230.467230.117180.11718-0.-0.0.965120.965120.052020.052020.635990.635990.685270.685270.898000.898000.837850.837850.0.0.938470.938470.718580.71858-0.-0.0.048370.048370.798310.798310.607800.607800.659380.659380.415990.415990.065990.06599-0.46017-0.460170.839360.83936两种方法的解是相同的。系数矩阵为40阶Hilbe

20、rt矩阵clearclcn=40;for i=1:n for j=1:n A(i,j)=1/(i+j-1); endendfor i=1:n b(i,1)=sum(A(i,1:i);endL1,U1=gauss1(A,b);y1=ones(n,1);y1=L1b;x1=ones(n,1);x1=U1y1; L2,U2,P=gauss2(A,b);y2=ones(n,1);y2=L2(P*b);x2=ones(n,1);x2=U2y2;x=x1,x2解为：-9.4510.98611682.0.4432-69474.0.4446.1.103-.1.429.1.066-.1.715.0.9974.0

21、.8529-.83320-0.0783.77148-1.300-.-2.476-.-3.836.-3.895.-3.102-.61393-3.619.24443-2.934-.44930-2.393.63981-1.022-.54634-0.09144.1.064-.2.992.3.377-.4.550.5.708-.6.321.499107.855-.781517.768-.7.818.7.017.209967.186-.200716.599-.5.959.684114.395-.311453.214.188562.374-.1.563-.0.78695-.-2949.2970.两种方法的解完全不同，用选主元的方法的解更加收敛。实验三实验名称：直接法的时间复杂性试验。实验目的：分别用三种不同方法求解线性方程组 Ax=b，不同工作量得出不同时间。实验内容与要求：生成方程组中矩阵和右端项，分别用，和三角分解法计算，并分别记录所花费的CPU时间，进行分析比较。实验要求：（1）取，随机生成的一条主对角线元素和二条次对