第二章 数值代数—————————————————————2 P45-46 实验题2(1)(3),5(A 为12阶改为5阶),8,10(n=300改为n=10);第三章迭代法——————————————————————25 P71-72 实验题2,3;习题2,12第四章数据建模—————————————————————34 P106 实验题1,2,4,8;习题3第五章数值微积分————————————————————43 P135 实验题1(1)-(4)4,5,8;第六章数值分析及其MATLAB 实验—————————————53 P166 实验题1(1),5,6第二章 数值代数实验 2 求下列矩阵的行列式、逆、特征值、特征向量、各种范数、和条件数:(1)⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡35-16-231-14 (1)>> A=[4 1 -1;3 2 -6;1 -5 3];a=det(A),B=inv(A),[V,D]=eig(A),t=eig(A)%矩阵的行列式a 、逆B 、特征值t 、特征向量V[norm(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.04260.1596 -0.1383 -0.22340.1809 -0.2234 -0.0532V =0.0185 -0.9009 -0.3066-0.7693 -0.1240 -0.7248-0.6386 -0.4158 0.6170D =-3.0527 0 00 3.6760 00 0 8.3766t =-3.05273.67608.3766ans =8.6089 10.0000 11.0000ans =3.7494 5.9574 5.7340(2)⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡1097591086781075675>> 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 、特征向量V[norm(A),norm(A,1),norm(A,inf)]%分别为矩阵A 的2,1,∞-范数[cond(A),cond(A,1),cond(A,inf)]%分别为矩阵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.000010.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.55200.1237 -0.5676 -0.6254 0.5209D =0.0102 0 0 00 0.8431 0 00 0 3.8581 00 0 0 30.2887t =0.01020.84313.858130.2887ans =30.2887 33.0000 33.0000ans =1.0e+03 *2.9841 4.4880 4.4880实验5Hilbert矩阵是著名的病态矩阵,n阶Hilbert矩阵定义为A=(a ij),其中a ij=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=hilb(5)A =1.0000 0.5000 0.3333 0.2500 0.20000.5000 0.3333 0.2500 0.2000 0.16670.3333 0.2500 0.2000 0.1667 0.14290.2500 0.2000 0.1667 0.1429 0.12500.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 -0.12600.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)求解下列方程组,验证计算结果,并分析误差产生的原因:⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡⨯2178.4617.59x x 11212592.1121-13.6-291.51314.59103.0432115-x x 顺序Gauss 消去法计算 >> 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.00000 -1.0430 -0.0529 -0.0176 -1.04360 -2.2079 -0.1120 -0.0373 -2.20900 -0.1971 -0.0100 -0.0033 -0.1972a =1.0e+018 *0.0000 0.0000 0.0000 0.0000 0.00000 -1.0430 -0.0529 -0.0176 -1.04360 0 -0.0000 -0.0000 00 0 0 0.0000 0a =1.0e+018 *0.0000 0.0000 0.0000 0.0000 0.00000 -1.0430 -0.0529 -0.0176 -1.04360 0 -0.0000 -0.0000 00 0 0 0.0000 0x =23.68481.0005选列主元Gauss 消去法计算>> 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=nagauss2(A,b)a =11.2000 9.0000 5.0000 2.0000 1.00000 -10.3817 -3.3621 1.0552 46.30760 59.1400 3.0000 1.0000 59.17000 1.1964 0.5536 0.8214 1.9107a =11.2000 9.0000 5.0000 2.0000 1.00000 59.1400 3.0000 1.0000 59.17000 0 -2.8354 1.2307 56.69460 0 0.4929 0.8012 0.7137a =11.2000 9.0000 5.0000 2.0000 1.00000 59.1400 3.0000 1.0000 59.17000 0 -2.8354 1.2307 56.69460 0 0 1.0151 10.5689x =3.84571.6095-15.476110.4113误差产生的原因由于第一列的主元素0.3×10-15的绝对值过于小,当它在消元过程作分母时把中间过程数值放大3.3×1015倍,使中间结果“吃掉了”原始数据,从而造成数值不稳定。
针对以上问题,考虑选用绝对值大的数作为主元素。
即选用选列主元Gauss 消去法。
实验10(追赶法的速度) 考虑n 阶三对角方程组10n 3443x 2112112112=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡, ,(1) 用选列主元Gauss 消去法(程序2.2)求解;(2) 编写追赶法程序并求解;(3) 用矩阵除法求解;(4) 比较三种方法的计算时间和精度.(1)用选列主元Gauss 消去法(程序2.2)求解;>> 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)第三章 迭代法实验2 用二分法(程序3.1)和Newton 迭代法(程序3.2)求下列方程的正跟: 05.01)1ln(22=---+-x x x x x二分法function x=nabisect(fname,a,b,e)if nargin<4,e=1e-4;end ;fa=feval(fname,a);fb=feval(fname,b);if fa*fb>0,error;endx=(a+b)/2while (b-a)>(2*e)fx=feval(fname,x);if fa*fx<0,b=x;fb=fx;else a=x;fa=fx;endx=(a+b)/2end>> fun=inline('x*log((x^2-1)^(1/2)+x)-(x^2-1)^(1/2)-0.5*x');x=nabisect(fun,2,3,0.05)x =2.5000x =2.2500x =2.1250x =2.0625x =2.0938x =2.0938Newton 迭代法function x=nanewton(fname,dfname,x0,e,N)if nargin<5,N=500;endif nargin<4,e=1e-4;endx=x0;x0=x+2*e;k=0;while abs(x0-x)>e&k<N,k=k+1;x0=x;x=x0-feval(fname,x0)/feval(dfname,x0);disp(x)endif k==N,warning;end>> syms x y dyy=x*log((x^2-1)^(1/2)+x)-(x^2-1)^(1/2)-0.5*x;dy=diff(y,x)dy =log(x + (x^2 - 1)^(1/2)) - x/(x^2 - 1)^(1/2) + (x*(x/(x^2 - 1)^(1/2) + 1))/(x + (x^2 - 1)^(1/2)) - 1/2>> fun=inline('x*log((x^2-1)^(1/2)+x)-(x^2-1)^(1/2)-0.5*x');dfun=inline('log(x + (x^2 - 1)^(1/2)) - x/(x^2 - 1)^(1/2) + (x*(x/(x^2 - 1)^(1/2) + 1))/(x + (x^2 - 1)^(1/2)) - 1/2');nanewton(fun,dfun,2,0.5e-3)2.12012.11552.1155ans =2.1155验3已知函数xx x f 2)(4-=在(-2,2)内有两个根。