矩阵分析与数值分析实验报告.docx

1、矩阵分析与数值分析实验报告矩阵分析与数值分析实验报告院系： 姓名： 学号：所在班号： 任课老师： 一设，分别编制从小到大和从大到小的顺序程序计算 并指出有效位数。程序如下：function sum3j=input(请输入求和个数 j:);A=0;B=0;double B;double A;for n=2:j m=n2-1; t=1./m; A=A+t;enddisp(从小到大：)s=Afor n=j:-1:2 m=n2-1; t=1./m; B=B+t;enddisp(从大到小:)s=B运行结果： sum3请输入求和个数 j:100从小到大：s =0.740049504950495从大到小:s

2、 =0.740049504950495 sum3请输入求和个数 j:10000从小到大：s =0.749900004999506从大到小:s =0.749900004999500 sum3请输入求和个数 j:1000000从小到大：s =0.749999000000522从大到小:s =0.749999000000500二、解线性方程组 1分别Jacobi迭代法和Gauss-Seidel迭代法求解线性方程组。迭代法计算停止的条件为：。解：（1）Jacobi迭代法程序代码：function jacobi(A, b, N)clc;clear;A=-2 1 0 0;1 -2 1 0;0 1 -2 1

3、;0 0 1 -2;b=-1 0 0 0;N=100;n = size(A,1);D = diag(diag(A);L = tril(-A,-1);U = triu(-A,1); Tj = inv(D)*(L+U);cj = inv(D)*b; tol = 1e-06;k = 1;format longx = zeros(n,1); while k = N x(:,k+1) = Tj*x(:,k) + cj; disp(k); disp(x = );disp(x(:,k+1); if norm(x(:,k+1)-x(:,k) tol disp(The procedure was success

4、ful) disp(Condition |x(k+1) - x(k)| tol was met after k iterations) disp(k); disp(x = );disp(x(:,k+1); break end k = k+1;end结果输出The procedure was successfulCondition |x(k+1) - x(k)| tol was met after k iterations 60x = 0.79999879906731 0.59999842795870 0.39999805685009 0.199*505（2）Gauss-Seidel迭代法程序代

5、码：function gauss_seidel(A, b, N)clc;clear;A=-2 1 0 0;1 -2 1 0;0 1 -2 1;0 0 1 -2;b=-1 0 0 0;N=100;n = size(A,1);D = diag(diag(A);L = tril(-A,-1);U = triu(-A,1); Tg = inv(D-L)*U; cg = inv(D-L)*b; tol = 1e-06;k = 1;x = zeros(n,1); while k = N x(:,k+1) = Tg*x(:,k) + cg; disp(k); disp(x = );disp(x(:,k+1)

6、; if norm(x(:,k+1)-x(:,k) tol disp(The procedure was successful) disp(Condition |x(k+1) - x(k)| tol was met after k iterations) disp(k); disp(x = );disp(x(:,k+1); break end k = k+1;end 结果输出The procedure was successfulCondition |x(k+1) - x(k)| tol was met after k iterations 31x = 0.79999921397935 0.5

7、9999897108561 0.39999916759077 0.199*5392. 用Gauss列主元消去法、QR方法求解如下方程组：(1)Gauss列主元消去法程序代码：function x=Gaussmain(A,b)clc;clear;format longA=1 2 1 2;2 5 3 -2;-2 -2 3 5;1 3 2 3;b=4 7 -1 0;N=length(A);x=zeros(N,1);y=zeros(N,1);c=0;d=0;A(:,N+1)=b;for k=1:N-1 for i=k:4 if ce a=a0-subs(f,a0)/subs(df,a0); a1=a0

8、; a0=a; N=N+1;endfprintf(a=%0.6f,a)N运行结果： newtonplease enter your equation:exp(x)+2(-x)+2*cos(x)-6f =exp(x)+2(-x)+2*cos(x)-6please enter you x(0):2df =exp(x)-2(-x)*log(2)-2*sin(x)a=1.829384N =42利用Newton迭代法求多项式的所有实零点，注意重根的问题。由于不知道重根的个数，采用试探法求重根。function X=newton2()clc;clear;syms x;delta=1e-06;f=inlin

9、e(x4-5.4*x3+10.56*x2-8.954*x+2.7951);df=inline(4*x3-16.2*x2+21.12*x-8.954);u=inline(x4-5.4*x3+10.56*x2-8.954*x+2.7951)/(4*x3-16.2*x2+21.12*x-8.954);du=inline(1-(x4-27/5*x3+264/25*x2-4477/500*x+27951/10000)/(4*x3-81/5*x2+528/25*x-4477/500)2*(12*x2-162/5*x+528/25);j=1;for i=0:1:3 x0=i; while 1 x1=x0-f

10、eval(u,x0)/feval(du,x0); err=abs(x1-x0); x0=x1; y=feval(f,x0); if (erreps n=n+1; h=(b-a)/n; I1=I2; I2=0; for i=0:n-1 x=a+h*i; x1=x+h; I2=I2+(h/2)*(subs(sym(f),findsym(sym(f),x)+. subs(sym(f),findsym(sym(f),x1); endendI=I2step=n结果输出I = 1.38654458753674step =532）Romberg公式程序代码：function approx = romberg

11、 (f,a,b,TOL)clc;clear;f=inline(1/x);a=2;b=8;TOL=1e-05;A=zeros(20,20);A(1,1)=(b-a)*(feval(f,a)+feval(f,b)/2;h=(b-a)/2;A(2,1)= A(1,1)/2+h*feval(f,a+h);A(2,2)=(4*A(2,1)-A(1,1)/3;errest =abs(A(2,2)-A(1,1)/2);i=2;while(errestTOL) i=i+1; h =h/2; sum=0.0; for j=1:2:2(i-1)-1 sum=sum+feval(f,a+j*h); end; A(i

12、,1)=A(i-1,1)/2+h*sum; for j=2:i power=4(j-1); A(i,j)=(power*A(i,j-1)-A(i-1,j-1)/(power-1); end; errest=abs(A(i,i)-A(i-1,i-1)/2(i-1);end; if (nargout =0) s=sprintf(tt approximate value of integral: t %.12f n,A(i,i); s=sprintf(%s tt error estimate: ttttt %.4e n,s,errest); s=sprintf(%s tt number of fun

13、ction evaluations: t %d n,s,2(i-1)+1); disp(s)else approx=A(i,i); end运行结果： approximate value of integral: 1.386297441871 error estimate: 8.7527e-006 number of function evaluations: 17五、插值与逼近 1给定上的函数，请做如下的插值逼近：构造等距节点分别取，的Lagrange插值多项式；构造分段线性取的Lagrange插值多项式；取Chebyshev多项式的零点：， 作插值节点构造的插值多项式和上述的插值多项式均要求

14、画出曲线图形（用不同的线型或颜色表示不同的曲线）。程序代码function Lagrangeclc;clear;close all;for i=1:3; if i=1 N=5; elseif i=2 N=8; else N=10; endf=inline(1/(1+25*x2);x1=zeros(1,N+1);a=-1;b=1;for i=1:N+1 x1(i)=a+(i-1)*(b-a)/N; y1(i)=feval(f,x1(i);endsyms xff=0; for i=1:N+1f=1; for j=1:i-1 f=f*(x-x1(j)/(x1(i)-x1(j); end for j=

15、i+1:N+1 f=f*(x-x1(j)/(x1(i)-x1(j); end ff=f*y1(i)+ff; f=1; endff=collect(ff,x);ff=vpa(ff,4); y=ff; p=ezplot(y,a,b); grid YLIM(-0.1 0.6); if N=5 set(p,Color,black); set(p,LineStyle,-); lagrange_5=y elseif N=8 set(p,Color,r); set(p,LineStyle,-); lagrange_8=y else set(p,Color,b); set(p,LineStyle,-) lag

16、range_10=y end hold on; xlabel(x);ylabel(y); title(y=p(x);hold on;endLag_Cheb();x=-1:0.01:1;y=1./(1+25*x.2);acu=plot(x,y);grid on;hold onset(acu,Color,m);set(acu,LineStyle,-); legend(N=5,N=8,N=10,Cheb,N=10,);function Lag_Cheb() f=inline(1/(1+25*x2);N=10;x1=zeros(1,N+1);a=-1;b=1; for i=1:N+1 x1(i)=co

17、s(2*i-1)*pi/(2*N); y1(i)=feval(f,x1(i);endsyms xff=0; for i=1:N+1f=1; for j=1:i-1 f=f*(x-x1(j)/(x1(i)-x1(j); end for j=i+1:N+1 f=f*(x-x1(j)/(x1(i)-x1(j); end ff=f*y1(i)+ff; f=1; endff=collect(ff,x);ff=vpa(ff,4); yy=ff;ff=collect(ff,x);yy=ff;lagrange_chebshev_10=yycheb=ezplot(yy,a,b); grid onYLIM(-0.

18、1 0.6);set(cheb,Color,g);set(cheb,LineStyle,:);结果输出lagrange_5 =.5673+1.202*x4-1.731*x2lagrange_8 =1.+53.69*x8-102.8*x6+61.37*x4-13.20*x2lagrange_10 =1.-220.9*x10+494.9*x8-381.4*x6+123.4*x4-16.86*x2lagrange_chebshev_10 =.7413-5.359*x10-.4e-2*x9+18.96*x8-.1321*x7-25.78*x6+.10*x5+16.81*x4+.5e-2*x3-5.33

19、6*x2+.5288e-3*x2已知函数值0123456789102.513.304.044.705.225.545.785.405.575.705.80和边界条件：求三次样条插值函数并画出其图形程序代码：function Spline3_1(x,y,df0,dfn)format shortx=0 1 2 3 4 5 6 7 8 9 10;y=2.51 3.30 4.04 4.70 5.22 5.54 5.78 5.40 5.57 5.70 5.80;plot(x,y,g*,MarkerSize,15);hold on;df0=1;dfn=0;n=length(x);h=zeros(n-1,1

20、);lan=zeros(n-2,1);mu=zeros(n-2,1);g=zeros(n-2,1);m=zeros(n,1);m(1)=df0;m(n)=dfn;for i=1:n-1 h(i)=x(i+1)-x(i);endfor i=1:n-2 mu(i)=h(i)/(h(i+1)+h(i); lan(i)=h(i+1)/(h(i+1)+h(i); g(i)=3*(mu(i)*(y(i+2)-y(i+1)/h(i+1)+lan(i)*(y(i+1)-y(i)/h(i);endA=zeros(n-2,n-2);A(1,1)=2;A(1,2)=mu(1);A(n-2,n-2)=lan(n-2)

21、;for i=2:n-2A(i,i)=2;A(i,i-1)=mu(i);A(i-1,i)=lan(i);endg(1)=g(1)-lan(1)*df0;g(n-2)=g(n-2)-mu(n-2)*dfn;b=Ag;for i=2:n-1 m(i)=b(i-1);endsyms zfor i=1:n-1 xx=x(i):0.01:x(i+1); sx1=y(i)*(xx-x(i+1).2.*(h(i)+2*(xx-x(i)/h(i)3; sx2=y(i+1)*(xx-x(i).2.*(h(i)-2*(xx-x(i+1)/h(i)3; sx3=m(i)*(xx-x(i+1).2.*(xx-x(i)

22、/h(i)2; sx4=m(i+1)*(xx-x(i).2.*(xx-x(i+1)/h(i)2; sx=sx1+sx2+sx3+sx4; z1=y(i)*(z-x(i+1).2.*(h(i)+2*(z-x(i)/h(i)3; z2=y(i+1)*(z-x(i).2.*(h(i)-2*(z-x(i+1)/h(i)3; z3=m(i)*(z-x(i+1).2.*(z-x(i)/h(i)2; z4=m(i+1)*(z-x(i).2.*(z-x(i+1)/h(i)2; zz=z1+z2+z3+z4; zz= collect(zz); vpa(zz,4) plot(xx,sx,r); hold on;

23、grid on; end程序输出ans =2.510+.1379*z3-.3479*z2+zans =2.692-.4369e-1*z3+.1969*z2+.4552*zans =2.287+.6872e-2*z3-.1065*z2+1.062*z ans =3.655-.4380e-1*z3+.3495*z2-.3060*z ans =-6.080+.1083*z3-1.476*z2+6.995*zans =41.14-.2694*z3+4.190*z2-21.34*zans =-109.8+.4295*z3-8.390*z2+54.14*z ans =133.0-.2784*z3+6.475*z2-49.91*zans =-57.75+.9411e-1*