矩阵与数值分析课程数值实验题目.docx

1、矩阵与数值分析课程数值实验题目矩阵与数值分析课程数值实验题目一、设，分别编制从小到大和从大到小的顺序计算，并指出有效位数。解：程序代码如下：function A,B=first(n)format long gT=1/2*(1+1/2-1/n-1/(n+1);k=0;S=0;for i=2:n S=S+1/(i2-1);endA=S; x=abs(A-T);i=0;while 1 if x1/2*10(-i) i=i+1; else break endendN1=k+i-1;display(从大到小:),A,N1S=0;for i=n:-1:2 S=S+1/(i2-1);endB=S;x=abs

2、(B-T);i=0;while 1 if x first(100)从大到小:A = 0.740049504950495N1 =15从小到大:B = 0.740049504950495N2 =15 first(10000)从大到小:A =0.749900004999506N1 =13从小到大:B =0.7499000049995N2 =323 first(1000000)从大到小:A =0.749999000000522N1 =13从小到大:B =0.7499990000005N2 =15分析：在做加法运算时，按照从小到大计算的顺序得到的结果要比按从大到小计算得到的结果有效数字位数更多。二、解线

3、性方程组 1分别Jacobi迭代法和Gauss-Seidel迭代法求解线性方程组。迭代法计算停止的条件为：。解：（1）Jacobi迭代法程序代码：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); 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)

4、= Tj*x(:,k) + cj; disp(k); disp(x = );disp(x(:,k+1); 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

5、 met after k iterations 60x = 0.79999879906731 0.59999842795870 0.39999805685009 0.199*505（2）Gauss-Seidel迭代法程序代码：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 = zer

6、os(n,1); while k = N x(:,k+1) = Tg*x(:,k) + cg; disp(k); disp(x = );disp(x(:,k+1); 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 successfulCo

7、ndition |x(k+1) - x(k)| tol was met after k iterations 31x = 0.79999921397935 0.59999897108561 0.39999916759077 0.199*539分析：GS法应用了上一步计算出来的较精确的值，所以可以用较少的迭代步数满足误差要求。2. 用Gauss列主元消去法、QR方法求解如下方程组：解：(1)Gauss列主元消去法程序代码：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

8、(N,1);y=zeros(N,1);c=0;d=0;A(:,N+1)=b;for k=1:N-1 for i=k:4 if cabs(A(i,k) d=i;c=abs(A(i,k); end end y=A(k,:); A(k,:)=A(d,:); A(d,:)=y; for i=k+1:N c=A(i,k); for j=1:N+1 A(i,j)=A(i,j)-A(k,j)*c/A(k,k); end endend b=A(:,N+1); x(N)=b(N)/A(N,N); for k=N-1:-1:1 x(k)=(b(k)-A(k,k+1:N)*x(k+1:N)/A(k,k); end结

9、果输出 ans = 18.00000000000000 -9.57142857142857 6.00000000000000 -0.42857142857143（2）QR方法：程序代码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;Q,R=qr(A)y=Q*b;x=Ry结果输出Q = -0.31622776601684 -0.04116934847963 -0.75164602800283 0.57735026918962 -0.63245553203368 -0.49403218175557 -0.15

10、0*056 -0.57735026918963 0.63245553203368 -0.74104827263336 -0.22549380840085 -0.00000000000000 -0.31622776601684 -0.45286283327594 0.60131682240226 0.57735026918963R = -3.16227766016838 -6.00832755431992 -0.94868329805051 2.84604989415154 0 -2.42899156029822 -4.65213637819829 -4.158*272 0 0 -0.67648

11、142520255 -0.52615221960200.0414*x = 17.99999999999989 -9.57142857142851 5.99999999999997 -0.42857142857143分析：Gauss列主元消去法和QR方法结果是一致的。三、非线性方程的迭代解法1用Newton迭代法求方程的根，计算停止的条件为：；解：Newton法程序代码：clc;clear;syms xf=inline(exp(x)+2(-x)+2*cos(x)-6);df=inline(exp(x)-2(-x)*log(2)-2*sin(x);x0=-3;delta=1e-06;while 1

12、 x1=x0-feval(f,x0)/feval(df,x0); err=abs(x1-x0); x0=x1; y=feval(f,x0); if (errdelta),break,endendx0x0=0;while 1 x1=x0-feval(f,x0)/feval(df,x0); err=abs(x1-x0); x0=x1; y=feval(f,x0); if (errdelta),break,endendx0结果输出:x0 =-2.98650806938193x0 =1.82938360193385分析：Newton法收敛速度较快，但是需要一个好的初始点。2利用Newton迭代法求多项

13、式的所有实零点，注意重根的问题。解：由于不知道重根的个数，采用试探法求重根。clc;clear;syms x;delta=1e-06;f=inline(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

14、*x-4477/500)2*(12*x2-162/5*x+528/25);j=1;for i=0:1:3 x0=i; while 1 x1=x0-feval(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=

15、n结果输出I = 1.38654458753674step =53分析：该算法运算精度高，但是计算时间长、计算量大。（2）Romberg公式程序代码：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/

16、2; sum=0.0; for j=1:2:2(i-1)-1 sum=sum+feval(f,a+j*h); end; A(i,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

17、estimate: ttttt %.4e n,s,errest); s=sprintf(%s tt number of function 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给定上的函数，请做如

18、下的插值逼近：构造等距节点分别取，的Lagrange插值多项式；构造分段线性取的Lagrange插值多项式；取Chebyshev多项式的零点：， 作插值节点构造的插值多项式和上述的插值多项式均要求画出曲线图形（用不同的线型或颜色表示不同的曲线）。解：程序代码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;

19、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); 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 els

20、eif N=8 set(p,Color,r); set(p,LineStyle,-); lagrange_8=y else set(p,Color,b); set(p,LineStyle,-) lagrange_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

21、,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)=cos(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);

22、ff=vpa(ff,4); yy=ff;ff=collect(ff,x);yy=ff;lagrange_chebshev_10=yycheb=ezplot(yy,a,b); grid onYLIM(-0.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*x2l

23、agrange_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.336*x2+.5288e-3*x分析：（1）拉格朗日插值法随着插值节点增多而更加接近真实曲线，但是还是会在区间-1，-0.4及0.4，1上，由于高次振荡产生容格现象。（2）分段插值法则优于前者，因为分段插值在每个小区间上都可以最大限度的满足边界条件。（3）利用契比雪夫迭代法可以更接近真实曲线。2已知函数值0123456789102.513.304.044.705.225.545.785.405.

24、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);lan=zeros(n-2,1);mu=zeros(n-2,1);g=zeros(n-2,1);m=zeros(n,1);

25、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);for i=2:n-2A(i,i)=2;A(i,i-1)=mu(i);A(i-1,i)=lan(i);endg(1)=g(1

26、)-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)/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