矩阵与数值分析

1、2013级工科硕士研究生矩阵与数值分析课程数值实验题目一、设，分别编制从小到大和从大到小的顺序程序分别计算并指出两种方法计算结果的有效位数。Matlab程序如下：function si,sd=S(N)format long; si=0;sd=0;for j=N:-1:2si=1.0e6/(j2-1)+si;endfor j=2:Nsd=1.0e6/(j2-1)+sd;endend在matlab命令窗口中输入：si,sd=S(10000)运行结果：si =7.499000049995000e+005sd =7.499000049994994e+005在matlab命令窗口中输入：si,sd=S(

2、1000000)运行结果：si =7.499990000005000e+005sd =7.499990000005200e+0051结果分析：si为从大到小的顺序求和的值，sd为从小到大的顺序求和的值。当N分别为10000和1000000时，si分别为7.499000049995000e+005和7.499990000005000e+005，可以看出这两个数的有效值均为13位；而sd分别为7.499000049994994e+005和7.499990000005200e+005，这两个数的有效值均为16位。这就出现了我们在矩阵理论课上所学的“大数吃小数”的问题。为了使结果更为精确我们必须避免在

3、四则运算中出现“大数吃小数”的情况，应该按从小到大的顺序进行求和。二、解线性方程组 1分别利用Jacobi迭代法和Gauss-Seidel迭代法求解线性方程组，其中常向量为维随机生成的列向量，系数矩阵具有如下形式，其中为阶矩阵，为阶单位矩阵，迭代法计算停止的条件为：，给出时的不同迭代步数求解系数矩阵A和向量b的Matlab程序如下：function A,b,x0=jz(n)for i=1:n-1 %此处n赋值即可，如n=100 for j=1:n-1 if(i=j) T(i,j)=2; end if(i=(j+1) T(i,j)=-1; end if(i=(j-1) T(i,j)=-1; en

4、d endendI=eye(n-1);k=1;for mm=1:(n-1) A(k:(k+n-2),k:(k+n-2)=T+2*I; k=k+n-1;endk=1;for xx=1:(n-2) A(k:(k+n-2),(k+n-1):(k+2*n-3)=-I; A(k+n-1):(k+2*n-3),k:(k+n-2)=-I; k=k+n-1;endb=randn(n-1)2,1);x0=zeros(n-1)2,1);Jacobi迭代法Matlab程序如下：function n=jacobi(A,b,x0)eps= 1.0e-10;M = 10000;D=diag(diag(A); %求A的对角

5、矩阵L=-tril(A,-1); %求A的下三角阵U=-triu(A,1); %求A的上三角阵B=D(L+U);f=Db;x=B*x0+f;n=1; %迭代次数while norm(x-x0)>=eps x0=x; x=B*x0+f; n=n+1; if(n>=M) disp('Warning: 迭代次数太多，可能不收敛！'); return; endendGauss-Seidel迭代法Matlab程序如下：function n=gauseidel(A,b,x0)eps= 1.0e-10;M = 10000; D=diag(diag(A); %求A的对角矩阵L=-t

6、ril(A,-1); %求A的下三角阵U=-triu(A,1); %求A的上三角阵G=(D-L)U;f=(D-L)b;x=G*x0+f;n=1; %迭代次数while norm(x-x0)>=eps x0=x; x=G*x0+f; n=n+1; if(n>=M) disp('Warning: 迭代次数太多，可能不收敛！'); return; endend在matlab命令窗口中输入：A,b,x0=jz(10);J10=jacobi(A,b,x0)G10=gauseidel(A,b,x0)A,b,x0=jz(20);J20=jacobi(A,b,x0)G20=gaus

7、eidel(A,b,x0)A,b,x0=jz(30);J30=jacobi(A,b,x0)G30=gauseidel(A,b,x0) 运行结果：J10 =436；G10 =214J20 =1810；G20 =913J30 =3990；G30 =2001结果分析：N=10且M=1000时，Jacobi迭代法和Gaussseidel迭代法的迭代次数分别为436和214；N=20且M=1000时，Jacobi迭代法和Gaussseidel迭代法的迭代次数分别为1810和913；N=30且M=1000时 ，Jacobi和Gauss-seidel算法的迭代次数分别为3990和2001次。从以上结果可知在

8、该问题下Jacobi算法的收敛性没有Gauss-seidel算法的收敛性好，原因在于Jacobi算法迭代过程中，对已算出的信息未加充分利用，一般来说，后面计算的值要比前面的计算值要精确些，而Gauss-seidel算法则充分利用了已求出来的信息，所以此算法的收敛性更为好一些。 2. 用Gauss列主元消去法、QR方法求解如下方程组：Gauss列主元消去法Matlab程序如下：function x,XA=GaussXQLineMain(A,b) N = size(A);n = N(1);index = 0;for i=1:(n-1) me = max(abs(A(1:n,i); %选取列主元 f

9、or k=i:n if(abs(A(k,i)=me) index = k; %保存列主元所在的行 break; end end temp = A(i,1:n); A(i,1:n) = A(index,1:n); A(index,1:n) = temp; bb = b(index); b(index)=b(i); b(i) = bb; %交换主行 for j=(i+1):n if(A(i,i)=0) disp('对角元素为0！'); return; end l = A(j,i); m = A(i,i); A(j,1:n)=A(j,1:n)-l*A(i,1:n)/m; b(j)=b

10、(j)-l*b(i)/m; %消元 endendfor i=n:-1:1 if(i<n) s=A(i,(i+1):n)*x(i+1):n,1); else s=0; end x(i,1)=(b(i)-s)/A(i,i);end XA = A;QR法Matlab程序如下：function x=qrxq(A,b) N = size(A);n = N(1);B = A; %保存系数矩阵A(1:n,1)=A(1:n,1)/norm(A(1:n,1); %将A的第一列正规化for i=2:n for j=1:(i-1) A(1:n,i)= A(1:n,i)-dot(A(1:n,i),A(1:n,j

11、)*A(1:n,j); %使A的第i列与前面所有的列正交 end A(1:n,i)=A(1:n,i)/norm(A(1:n,i); %将A的第i列正规化endQ = A; %分解出来的正交矩阵R = transpose(Q)*B;bb=transpose(Q)*b;for i=n:-1:1 if(i<n) s=R(i,(i+1):n)*x(i+1):n,1); else s=0; end x(i,1)=(bb(i)-s)/R(i,i);end 在matlab命令窗口中输入：A=1 2 -1 1;2 5 0 3;1 7 9 2;8 -1 -2 1;b=0;4;12;-8;x_G=Gauss

12、XQLineMain(A,b)x_QR=qrxq(A,b) 运行结果：x_G =-1.00000000000000 -0.00000000000000 1.00000000000000 2.00000000000000x_QR =-1.00000000000000 -0.00000000000001 1.00000000000000 2.00000000000002 三、非线性方程的迭代解法1求方程的根，迭代停止的条件为：；使用弦截法求解此方程的根，弦截法的Matlab程序如下：function root=Secant(f,a,b)%f:方程;%a:区间左端点;%b:区间右端点;%root:求

13、得的根;eps=1.0e-10;f1=subs(sym(f),findsym(sym(f),a);f2=subs(sym(f),findsym(sym(f),b);if(f1=0) root=a;endif(f2=0) root=b;endif(f1*f2>0) disp('两端点函数值乘积大于0!'); return;else tol=1; fa=subs(sym(f),findsym(sym(f),a); fb=subs(sym(f),findsym(sym(f),b); root=a-(b-a)*fa/(fb-fa);while(tol>eps) r1=roo

14、t; fx=subs(sym(f),findsym(sym(f),r1); s=fx*fa; if(s=0) root=r1; else if(s>0) root=b-(r1-b)*fb/(fx-fb); else root=a-(r1-a)*fa/(fx-fa); end end tol=abs(root-r1); endend利用matlab的绘图功能可以确定求根区间为（1,3）在matlab命令窗口中输入：root=Secant('exp(x)+2*x2+2*sin(x)-log(x)-16',1,3)运行结果：root =1.962922683164462利用Ne

15、wton迭代法求多项式的所有实零点，注意重根的问题。Newton迭代法的Matlab程序如下：x=-3;d=-1;while abs(d-x)>0.5*10(-8) %在区间（-3，-1） d=x; x=x-(x4-3*x3-3*x2+11*x-6)/(4*x3-9*x2-6*x+11);%牛顿法endx0=dx=0;d=2;while abs(d-x)>0.5*10(-8)%在区间（0，2）的重根 d=x; x=x-2*(x4-3*x3-3*x2+11*x-6)/(4*x3-9*x2-6*x+11);endx1=dx=2.5;d=4;while abs(d-x)>0.5*1

16、0(-8)%在区间（2.5，4） d=x; x=x-(x4-3*x3-3*x2+11*x-6)/(4*x3-9*x2-6*x+11);endx2=d运行结果：x0 =-2.00000000151233x1 =1.00000000065769x2 =3.00000000000002结果分析：由以上结果可知，newton迭代法所得结果与初始值的选取密切相关。不同的初始值会产生不同的结果，甚至会影响牛顿迭代法收敛性。四、数值积分用三点Gauss型求积公式计算积分n点Gauss型求积公式matlab程序如下：function I = question4(f,a,b,n)%f为求积方程;%a,b分别为积

17、分上下限;%n为Gauss点个数;XK,AK=grule(n);XK1=(b-a)/2)*XK+(a+b)/2);I=(b-a)/2)*sum(AK.*subs(f,findsym(f),XK1);function bp,wf=grule(n)% bp,wf=grule(n)bp=zeros(n,1); wf=bp; iter=2; m=fix(n+1)/2); e1=n*(n+1);mm=4*m-1; t=(pi/(4*n+2)*(3:4:mm); nn=(1-(1-1/n)/(8*n*n);xo=nn*cos(t);for j=1:iter pkm1=1; pk=xo; for k=2:n

18、 t1=xo.*pk; pkp1=t1-pkm1-(t1-pkm1)/k+t1; pkm1=pk; pk=pkp1; end den=1.-xo.*xo; d1=n*(pkm1-xo.*pk); dpn=d1./den; d2pn=(2.*xo.*dpn-e1.*pk)./den; d3pn=(4*xo.*d2pn+(2-e1).*dpn)./den; d4pn=(6*xo.*d3pn+(6-e1).*d2pn)./den; u=pk./dpn; v=d2pn./dpn; h=-u.*(1+(.5*u).*(v+u.*(v.*v-u.*d3pn./(3*dpn); p=pk+h.*(dpn+(

19、.5*h).*(d2pn+(h/3).*(d3pn+.25*h.*d4pn); dp=dpn+h.*(d2pn+(.5*h).*(d3pn+h.*d4pn/3); h=h-p./dp; xo=xo+h;endbp=-xo-h;fx=d1-h.*e1.*(pk+(h/2).*(dpn+(h/3).*(. d2pn+(h/4).*(d3pn+(.2*h).*d4pn);wf=2*(1-bp.2)./(fx.*fx);if (m+m) > n, bp(m)=0; endif (m+m) = n), m=m-1; endjj=1:m; n1j=(n+1-jj); bp(n1j)=-bp(jj);

20、 wf(n1j)=wf(jj);% end在matlab命令窗口中输入：syms x;f=exp(-x);a=0;b=1;n=3;I = question4(f,a,b,n)运行结果：I =0.632120255664068五、插值与逼近 1给定上的函数，请做如下的插值逼近：（1）构造等距节点分别取，的Lagrange插值多项式；（2）取Chebyshev多项式的零点：，作插值节点构造的插值多项式和上述的插值多项式均要求画出曲线图形（用不同的线型或颜色表示不同的曲线）。（1）Lagrange插值matlab程序如下：function f = L1(n)syms t s; q=1/(1+s2);

21、for i=1:n x(i)=-1+2*(i-1)/(n-1); y(i)=subs(q,s,x(i);endf = 0.0;for(i = 1:n) l = y(i); for(j = 1:i-1) l = l*(t-x(j)/(x(i)-x(j); end; for(j = i+1:n) l = l*(t-x(j)/(x(i)-x(j); %计算拉格朗日基函数 end; f = f + l; %计算拉格朗日插值函数 f=collect(f); %化简endX=-1:0.1:1;Y=subs(q,s,X);F=subs(f,t,X);plot(X,Y,'g',X,F,'

22、;r');f=collect(f); %展开digits(5);f=vpa(f);1、在matlab命令窗口中输入：f = L1(5)运行结果2、在matlab命令窗口中输入：f = L1(8)运行结果3、在matlab命令窗口中输入：f =L1(10)运行结果：结果分析：由图像可以看出n=5时拟合曲线与原始曲线的误差比较大，而n=8,10时拟合曲线与原始曲线几乎重合在一起。可见节点个越多，拉格朗日插值法所得的多项式曲线与原始曲线越接近，误差就越小。（2）取Chebyshev多项式的零点为插值点,此时Lagrange插值matlab程序主体基本不变，变的只是插值点，其程序如下：func

23、tion f = L2(n)syms t s; q=1/(1+s2);for i=1:n x(i)=cos(2*i-1)*pi/2/n); y(i)=subs(q,s,x(i);endf = 0.0;for(i = 1:n) l = y(i); for(j = 1:i-1) l = l*(t-x(j)/(x(i)-x(j); end; for(j = i+1:n) l = l*(t-x(j)/(x(i)-x(j); %计算拉格朗日基函数 end; f = f + l; %计算拉格朗日插值函数endX=-1:0.1:1; %绘图Y=subs(q,s,X);F=subs(f,t,X);plot(X

24、,Y,'k',X,F,'rp');f=collect(f); %展开digits(5);f=vpa(f);在matlab命令窗口中输入：f = L2(10)运行结果：结果分析：有图像比较可以看出，N=10时，chebyshev零为拉格朗日插值时，其拟合曲线与原始曲线几乎重合，精度较高。2已知函数值346602和边界条件：求三次样条插值函数并画出其图形三次样条插值matlab程序如下：function ThrSample1 (x,y,y_1, y_N)syms t;f = 0.0;if(length(x) = length(y) n = length(x);els

25、e disp('x和y的维数不相等！'); return;end %维数检查A = diag(2*ones(1,n); %求解m的系数矩阵u = zeros(n-2,1);lamda = zeros(n-1,1);c = zeros(n,1);for i=2:n-1 u(i-1) = (x(i)-x(i-1)/(x(i+1)-x(i-1); lamda(i) = (x(i+1)-x(i)/(x(i+1)-x(i-1); c(i) = 3*lamda(i)*(y(i)-y(i-1)/(x(i)-x(i-1)+ . 3*u(i-1)*(y(i+1)-y(i)/(x(i+1)-x(i

26、); A(i, i+1) = u(i-1); A(i, i-1) = lamda(i); %形成系数矩阵及向量cendc(1) = 2*y_1;c(n) = 2*y_N;m = followup(A,c); %用追赶法求解方程组j=1;while j<=n-1 h = x(j+1) - x(j); f = y(j)*(2*(t-x(j)+h)*(t-x(j+1)2/h3 + . y(j+1)*(2*(x(j+1)-t)+h)*(t-x(j)2/h3 + . m(j)*(t-x(j)*(x(j+1)-t)2/h2 - . m(j+1)*(x(j+1)-t)*(t-x(j)2/h2; f=c

27、ollect(f); fprintf('当x属于%d,%d,f=n', x(j), x(j+1); disp(f); X=x(j):0.1:x(j+1); %绘图 F=subs(f,t,X); plot(X,F,'g'); hold on; j=j+1;endfunction x=followup(A,b)n = rank(A);for(i=1:n) if(A(i,i)=0) disp('Error: 对角有元素为0！'); return; endend;d = ones(n,1);a = ones(n-1,1);c = ones(n-1);for(i=1:n-1) a(i,1)=A(i+1,i); c(i,1)=A(i,i+1); d(i,1)=A(