桂林理工大学数值分析上机实验作业

1、数值分析上机实验报告 学校：桂林理工大学 学院：环境科学与工程学院 专业：市政工程 学号： 姓名 ：邢梦龙第2章 数值代数2P45-46 实验题2（1）（3），5（A为12阶改为5阶），8，10（n=300改为n=10）；第三章 迭代法25P71-72 实验题2，3；习题2，12第四章 数据建模34P106 实验题1，2，4；习题3第五章 数值微积分43P135 实验题1（1）-（4）4，5，8；第六章 数值分析及其MATLAB实验53P166 实验题1（1），5，6第二章 数值代数实验 2 求下列矩阵的行列式、逆、特征值、特征向量、各种范数、和条件数：（1）（1） A=4 1 -1;3 2

2、-6;1 -5 3;a=det(A),B=inv(A),V,D=eig(A),t=eig(A)%矩阵的行列式a、逆B、特征值t、特征向量Vnorm(A),norm(A,1),norm(A,inf)%分别为矩阵A的2,1，-范数cond(A),cond(A,1),cond(A,inf)%分别为矩阵A的2,1，-条件数a = -94B = 0.2553 -0.0213 0.0426 0.1596 -0.1383 -0.2234 0.1809 -0.2234 -0.0532V = 0.0185 -0.9009 -0.3066 -0.7693 -0.1240 -0.7248 -0.6386 -0.41

3、58 0.6170D = -3.0527 0 0 0 3.6760 0 0 0 8.3766t = -3.0527 3.6760 8.3766ans = 8.6089 10.0000 11.0000ans = 3.7494 5.9574 5.7340（2） A=5 7 6 5;7 10 8 7;6 8 10 9;5 7 9 10;a=det(A),B=inv(A),V,D=eig(A),t=eig(A)%矩阵的行列式a、逆B、特征值t、特征向量Vnorm(A),norm(A,1),norm(A,inf)%分别为矩阵A的2,1，-范数cond(A),cond(A,1),cond(A,inf)%分

4、别为矩阵A的2,1，-条件数a = 1.0000B = 68.0000 -41.0000 -17.0000 10.0000 -41.0000 25.0000 10.0000 -6.0000 -17.0000 10.0000 5.0000 -3.0000 10.0000 -6.0000 -3.0000 2.0000V = 0.8304 0.0933 0.3963 0.3803 -0.5016 -0.3017 0.6149 0.5286 -0.2086 0.7603 -0.2716 0.5520 0.1237 -0.5676 -0.6254 0.5209D = 0.0102 0 0 0 0 0.8

5、431 0 0 0 0 3.8581 0 0 0 0 30.2887t = 0.0102 0.8431 3.8581 30.2887ans = 30.2887 33.0000 33.0000ans = 1.0e+03 * 2.9841 4.4880 4.4880实验5 Hilbert矩阵是著名的病态矩阵，n阶Hilbert矩阵定义为A=（aij），其中aij=1/(i+j-1)。设A为5阶Hilbert矩阵，计算cond（A），A-1，norm（A-1 A-E）及|A|A-1|-1，并分析结果的精度。再比较MATLAB求解HIlbert矩阵及其逆函数hilb（5），invhilb（5）。 A=

6、hilb(5)A = 1.0000 0.5000 0.3333 0.2500 0.2000 0.5000 0.3333 0.2500 0.2000 0.1667 0.3333 0.2500 0.2000 0.1667 0.1429 0.2500 0.2000 0.1667 0.1429 0.1250 0.2000 0.1667 0.1429 0.1250 0.1111 cond(A)ans = 4.7661e+005 inv(A)ans = 1.0e+005 * 0.0002 -0.0030 0.0105 -0.0140 0.0063 -0.0030 0.0480 -0.1890 0.2688

7、 -0.1260 0.0105 -0.1890 0.7938 -1.1760 0.5670 -0.0140 0.2688 -1.1760 1.7920 -0.88200.0063 -0.1260 0.5670 -0.8820 0.4410 norm(inv(A)*A-eye(5)ans = 4.8038e-012 det(A)*det(inv(A)-1ans = 2.4158e-013 det(hilb(5)ans = 3.7493e-012 det(inv(hilb(5)ans = 2.6672e+011实验8 分别用顺序Gauss消去法（程序2.1）和选列主元Gauss消去法（程序2.2）

8、求解下列方程组，验证计算结果，并分析误差产生的原因：顺序Gauss消去法程序function x=nagauss(a,b,flag)if nargin A=0.3*10(-15) 59.14 3 1;5.291 -6.13 -1 2;11.2 9 5 2;1 2 1 1;b=59.17;46.78;1;2; x=nagauss(A,b)a = 1.0e+018 * 0.0000 0.0000 0.0000 0.0000 0.0000 0 -1.0430 -0.0529 -0.0176 -1.0436 0 -2.2079 -0.1120 -0.0373 -2.2090 0 -0.1971 -0.

9、0100 -0.0033 -0.1972a = 1.0e+018 * 0.0000 0.0000 0.0000 0.0000 0.0000 0 -1.0430 -0.0529 -0.0176 -1.0436 0 0 -0.0000 -0.0000 0 0 0 0 0.0000 0a = 1.0e+018 * 0.0000 0.0000 0.0000 0.0000 0.0000 0 -1.0430 -0.0529 -0.0176 -1.0436 0 0 -0.0000 -0.0000 0 0 0 0 0.0000 0x = 23.6848 1.0005 0 0选列主元Gauss消去法程序func

10、tion x=nagauss2(a,b,flag)if nargink, t=a(k,:);a(k,:)=a(p,:);a(p,:)=t;enda(k+1):n,(k+1):(n+1)=a(k+1):n,(k+1) :(n+1)-a(k+1):n,k)/a(k,k)*a(k,(k+1):(n+1);a(k+1):n,k)=zeros(n-k,1);if flag=0,a,endendx=zeros(n,1);x(n)=a(n,n+1)/a(n,n);for k=n-1:-1:1 x(k,:)=(a(k,n+1)-a(k,(k+1):n)*x(k+1):n)/a(k,k);end计算 A=0.3

11、*10(-15) 59.14 3 1;5.291 -6.13 -1 2;11.2 9 5 2;1 2 1 1;b=59.17;46.78;1;2; x=nagauss2(A,b)a = 11.2000 9.0000 5.0000 2.0000 1.0000 0 -10.3817 -3.3621 1.0552 46.3076 0 59.1400 3.0000 1.0000 59.1700 0 1.1964 0.5536 0.8214 1.9107a = 11.2000 9.0000 5.0000 2.0000 1.0000 0 59.1400 3.0000 1.0000 59.1700 0 0

12、-2.8354 1.2307 56.6946 0 0 0.4929 0.8012 0.7137a = 11.2000 9.0000 5.0000 2.0000 1.0000 0 59.1400 3.0000 1.0000 59.1700 0 0 -2.8354 1.2307 56.6946 0 0 0 1.0151 10.5689x = 3.8457 1.6095 -15.4761 10.4113误差产生的原因由于第一列的主元素0.310-15的绝对值过于小，当它在消元过程作分母时把中间过程数值放大3.31015倍，使中间结果“吃掉了”原始数据，从而造成数值不稳定。 针对以上问题，考虑选用绝对

13、值大的数作为主元素。即选用选列主元Gauss消去法。实验10(追赶法的速度) 考虑n阶三对角方程组 ,(1) 用选列主元Gauss消去法(程序2.2)求解；(2) 编写追赶法程序并求解；(3) 用矩阵除法求解；(4) 比较三种方法的计算时间和精度.(1)用选列主元Gauss消去法(程序2.2)求解；function x=nagauss2(a,b,flag)if nargink, t=a(k,:);a(k,:)=a(p,:);a(p,:)=t;enda(k+1):n,(k+1):(n+1)=a(k+1):n,(k+1) :(n+1)-a(k+1):n,k)/a(k,k)*a(k,(k+1):(n

14、+1);a(k+1):n,k)=zeros(n-k,1);if flag=0,a,endendx=zeros(n,1);x(n)=a(n,n+1)/a(n,n);for k=n-1:-1:1 x(k,:)=(a(k,n+1)-a(k,(k+1):n)*x(k+1):n)/a(k,k);end A=2*ones(1,10);%对角线元素B=1*ones(1,9);%对角线上方的元素,个数比A少一个C=1*ones(1,9);%对角线下方的元素,个数比A少一个D=diag(A)+diag(B,1)+diag(C,-1);b=3;4*ones(8,1);3;x=nagauss2(D,b)a = 2.

15、0000 1.0000 0 0 0 0 0 0 0 0 3.0000 0 1.5000 1.0000 0 0 0 0 0 0 0 2.5000 0 1.0000 2.0000 1.0000 0 0 0 0 0 0 4.0000 0 0 1.0000 2.0000 1.0000 0 0 0 0 0 4.0000 0 0 0 1.0000 2.0000 1.0000 0 0 0 0 4.0000 0 0 0 0 1.0000 2.0000 1.0000 0 0 0 4.0000 0 0 0 0 0 1.0000 2.0000 1.0000 0 0 4.0000 0 0 0 0 0 0 1.0000

16、 2.0000 1.0000 0 4.0000 0 0 0 0 0 0 0 1.0000 2.0000 1.0000 4.0000 0 0 0 0 0 0 0 0 1.0000 2.0000 3.0000a = 2.0000 1.0000 0 0 0 0 0 0 0 0 3.0000 0 1.5000 1.0000 0 0 0 0 0 0 0 2.5000 0 0 1.3333 1.0000 0 0 0 0 0 0 2.3333 0 0 1.0000 2.0000 1.0000 0 0 0 0 0 4.0000 0 0 0 1.0000 2.0000 1.0000 0 0 0 0 4.0000

17、 0 0 0 0 1.0000 2.0000 1.0000 0 0 0 4.0000 0 0 0 0 0 1.0000 2.0000 1.0000 0 0 4.0000 0 0 0 0 0 0 1.0000 2.0000 1.0000 0 4.0000 0 0 0 0 0 0 0 1.0000 2.0000 1.0000 4.0000 0 0 0 0 0 0 0 0 1.0000 2.0000 3.0000a = 2.0000 1.0000 0 0 0 0 0 0 0 0 3.0000 0 1.5000 1.0000 0 0 0 0 0 0 0 2.5000 0 0 1.3333 1.0000

18、 0 0 0 0 0 0 2.3333 0 0 0 1.2500 1.0000 0 0 0 0 0 2.2500 0 0 0 1.0000 2.0000 1.0000 0 0 0 0 4.0000 0 0 0 0 1.0000 2.0000 1.0000 0 0 0 4.0000 0 0 0 0 0 1.0000 2.0000 1.0000 0 0 4.0000 0 0 0 0 0 0 1.0000 2.0000 1.0000 0 4.0000 0 0 0 0 0 0 0 1.0000 2.0000 1.0000 4.0000 0 0 0 0 0 0 0 0 1.0000 2.0000 3.0

19、000a = 2.0000 1.0000 0 0 0 0 0 0 0 0 3.0000 0 1.5000 1.0000 0 0 0 0 0 0 0 2.5000 0 0 1.3333 1.0000 0 0 0 0 0 0 2.3333 0 0 0 1.2500 1.0000 0 0 0 0 0 2.2500 0 0 0 0 1.2000 1.0000 0 0 0 0 2.2000 0 0 0 0 1.0000 2.0000 1.0000 0 0 0 4.0000 0 0 0 0 0 1.0000 2.0000 1.0000 0 0 4.0000 0 0 0 0 0 0 1.0000 2.000

20、0 1.0000 0 4.0000 0 0 0 0 0 0 0 1.0000 2.0000 1.0000 4.0000 0 0 0 0 0 0 0 0 1.0000 2.0000 3.0000a = 2.0000 1.0000 0 0 0 0 0 0 0 0 3.0000 0 1.5000 1.0000 0 0 0 0 0 0 0 2.5000 0 0 1.3333 1.0000 0 0 0 0 0 0 2.3333 0 0 0 1.2500 1.0000 0 0 0 0 0 2.2500 0 0 0 0 1.2000 1.0000 0 0 0 0 2.2000 0 0 0 0 0 1.166

21、7 1.0000 0 0 0 2.1667 0 0 0 0 0 1.0000 2.0000 1.0000 0 0 4.0000 0 0 0 0 0 0 1.0000 2.0000 1.0000 0 4.0000 0 0 0 0 0 0 0 1.0000 2.0000 1.0000 4.0000 0 0 0 0 0 0 0 0 1.0000 2.0000 3.0000a = 2.0000 1.0000 0 0 0 0 0 0 0 0 3.0000 0 1.5000 1.0000 0 0 0 0 0 0 0 2.5000 0 0 1.3333 1.0000 0 0 0 0 0 0 2.3333 0

22、 0 0 1.2500 1.0000 0 0 0 0 0 2.2500 0 0 0 0 1.2000 1.0000 0 0 0 0 2.2000 0 0 0 0 0 1.1667 1.0000 0 0 0 2.1667 0 0 0 0 0 0 1.1429 1.0000 0 0 2.1429 0 0 0 0 0 0 1.0000 2.0000 1.0000 0 4.0000 0 0 0 0 0 0 0 1.0000 2.0000 1.0000 4.0000 0 0 0 0 0 0 0 0 1.0000 2.0000 3.0000a = 2.0000 1.0000 0 0 0 0 0 0 0 0

23、 3.0000 0 1.5000 1.0000 0 0 0 0 0 0 0 2.5000 0 0 1.3333 1.0000 0 0 0 0 0 0 2.3333 0 0 0 1.2500 1.0000 0 0 0 0 0 2.2500 0 0 0 0 1.2000 1.0000 0 0 0 0 2.2000 0 0 0 0 0 1.1667 1.0000 0 0 0 2.1667 0 0 0 0 0 0 1.1429 1.0000 0 0 2.1429 0 0 0 0 0 0 0 1.1250 1.0000 0 2.1250 0 0 0 0 0 0 0 1.0000 2.0000 1.000

24、0 4.0000 0 0 0 0 0 0 0 0 1.0000 2.0000 3.0000a = 2.0000 1.0000 0 0 0 0 0 0 0 0 3.0000 0 1.5000 1.0000 0 0 0 0 0 0 0 2.5000 0 0 1.3333 1.0000 0 0 0 0 0 0 2.3333 0 0 0 1.2500 1.0000 0 0 0 0 0 2.2500 0 0 0 0 1.2000 1.0000 0 0 0 0 2.2000 0 0 0 0 0 1.1667 1.0000 0 0 0 2.1667 0 0 0 0 0 0 1.1429 1.0000 0 0

25、 2.1429 0 0 0 0 0 0 0 1.1250 1.0000 0 2.1250 0 0 0 0 0 0 0 0 1.1111 1.0000 2.1111 0 0 0 0 0 0 0 0 1.0000 2.0000 3.0000a = 2.0000 1.0000 0 0 0 0 0 0 0 0 3.0000 0 1.5000 1.0000 0 0 0 0 0 0 0 2.5000 0 0 1.3333 1.0000 0 0 0 0 0 0 2.3333 0 0 0 1.2500 1.0000 0 0 0 0 0 2.2500 0 0 0 0 1.2000 1.0000 0 0 0 0

26、2.2000 0 0 0 0 0 1.1667 1.0000 0 0 0 2.1667 0 0 0 0 0 0 1.1429 1.0000 0 0 2.1429 0 0 0 0 0 0 0 1.1250 1.0000 0 2.1250 0 0 0 0 0 0 0 0 1.1111 1.0000 2.1111 0 0 0 0 0 0 0 0 0 1.1000 1.1000x = 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000（2）编写追赶法程序并求解；function x,L,U=thomas(a,b,c

27、,f) n=length(b); % 对A进行分解 u(1)=b(1); for i=2:n if(u(i-1)=0) l(i-1)=a(i-1)/u(i-1); u(i)=b(i)-l(i-1)*c(i-1); else break; endendL=eye(n)+diag(l,-1);U=diag(u)+diag(c,1);x=zeros(n,1); y=x; % 求解Ly=b y(1)=f(1); for i=2:n y(i)=f(i)-l(i-1)*y(i-1);end % 求解Ux=y if(u(n)=0) x(n)=y(n)/u(n); endfor i=n-1:-1:1 x(i)=(y(i)-c(i)*