矩阵与数值分析作业.docx

1、矩阵与数值分析作业2011级工科硕士研究生矩阵与数值分析课程数值实验题目一、 对于数列，有如下两种生成方式1、首项为，递推公式为；2、前两项为，递推公式为；给出利用上述两种递推公式生成的序列的第50项。解：matlab编程1、c1(1)=1;for t1=2:50 c1(t1)=c1(t1-1)/3;endc1(50)运算结果：ans = 4.1788e-0242、c2(1)=1;c2(2)=1/3;for t2=3:50 c2(t2)=c2(t2-1)*10/3-c2(t2-2);endc2(50)运算结果：ans =-4.9661e+006二、 利用迭代格式及Aitken加速后的新迭代格式

2、求方程在内的根解：matlab编程k=1;x(1)=1.3;f=1;while (abs(f)1.0e-6) x(k+1)=sqrt(10/(x(k)+4); k=k+1; f=x(k)3+4*x(k)2-10;endx,fcleark=1;x(1)=1.3;f=1;while (abs(f)1.0e-6) y(k)=sqrt(10/(x(k)+4); z(k)=sqrt(10/(y(k)+4); x(k+1)=z(k)-(z(k)-y(k)2/(z(k)-2*y(k)+x(k); k=k+1; f=x(k)3+4*x(k)2-10;endx,f运算结果：x =1.3000 1.3735 1.

3、3641 1.3653 1.3652 1.3651 1.3651 1.3651f =5.8603e-007x =1.3000 1.3651 1.3651f =2.5046e-012分析：原迭代格式运算8次，达到精度要求；而是用Aitken加速后，只用3次迭代，即Aitken加速可以明显加快迭代速度。三、解线性方程组 1分别Jacobi迭代法和Gauss-Seidel迭代法求解线性方程组迭代法计算停止的条件为：解：matlab编程1Xj=0;0;0;0;Xg=0;0;0;0;A=6 2 1 -2; 2 5 0 -2; -2 0 8 5; 1 3 2 7;b=4;7;-1;0;D=diag(dia

4、g(A);L=-tril(A)+D;U=-triu(A)+D;Bj=inv(D)*(L+U);f=inv(D)*b;k=1;flag=1;Xj(:,k+1)=Bj*Xj(:,k)+f;while flag=1 if abs(Xj(:,k+1)-Xj(:,k)1e-6 flag=0; else k=k+1; Xj(:,k+1)=Bj*Xj(:,k)+f; endendXjBg=inv(D-L)*U;fg=inv(D-L)*b;k=1;flag=1;Xg(:,k+1)=Bg*Xg(:,k)+fg;while flag=1 if abs(Xg(:,k+1)-Xg(:,k)1e-6 flag=0; e

5、lse k=k+1; Xg(:,k+1)=Bg*Xg(:,k)+fg; endendXg运算结果：Xj =Columns 1 through 11 0 0.6668 0.2208 0.0622 0.0838 0.0567 0.0537 0.0535 0.0522 0.0522 0.05210 1.4000 1.1333 1.0478 1.1635 1.1442 1.1478 1.1515 1.1505 1.1509 1.15090 -0.1250 0.0416 0.3424 0.2212 0.2433 0.2478 0.2432 0.2448 0.2447 0.24460 0 -0.6593

6、-0.5292 -0.5557 -0.5738 -0.5681 -0.5704 -0.5706 -0.5705 -0.5706Columns 12 through 16 0.0521 0.0521 0.0521 0.0520 0.05201.1509 1.1509 1.1509 1.1509 1.15090.2445 0.2445 0.2445 0.2445 0.2445-0.5707 -0.5707 -0.5707 -0.5707 -0.5707Xg =0 0.6667 0.0843 0.0551 0.0523 0.0521 0.0521 0.0520 0.05200 1.1332 1.12

7、91 1.1490 1.1508 1.1508 1.1509 1.1509 1.15090 0.0417 0.2666 0.2464 0.2448 0.2446 0.2446 0.2446 0.24460 -0.5928 -0.5721 -0.5707 -0.5706 -0.5706 -0.5706 -0.5706 -0.57062. 用Gauss列主元消去法、QR方法求解如下方程组：2、A=2 2 1 2; 4 1 3 -1; -4 -2 0 1; 2 3 2 3;b=1;2;1;0;Ab=A b;for j=1:3 k=j; max(j)=Ab(j,j); for i=(j+1):4 if

8、 max(j)Ab(i,j) k=i; max(j)=Ab(i,j); end end if k=j t=Ab(j,:); Ab(j,:)=Ab(k,:); Ab(k,:)=t; end for i=(j+1):4 Ab(i,:)=Ab(i,:)-Ab(j,:)*Ab(i,j)/Ab(j,j); endendfor i=4:-1:1 Abt=Ab(i,5); for j=4:-1:(i+1) Abt=Abt-x(j)*Ab(i,j); end x(i)=Abt/Ab(i,i);endXg=xQ,R=qr(A);Xqr=RQ*b运算结果：Xg =1.5417-2.75000.08321.6667

9、Xqr =1.5417-2.75000.08321.6667分析：G-S迭代法与Jacob迭代法相比结果相同，但减少很多迭代步数。四、已知一组数据点，编写一程序求解三次样条插值函数满足 并针对下面一组具体实验数据0.250.30.390.450.530.50000.54770.62450.67080.7280求解，其中边界条件为.解：n=5;X=0.25 0.3 0.39 0.45 0.53;Y=0.5000 0.5477 0.6245 0.6708 0.7280;f00=0;fnn=0;for k=1:(n-1) h(k)=X(k+1)-X(k);endg(1)=3*(Y(2)-Y(1)/h

10、(1)-f00*h(1)/2;g(n)=3*(Y(n)-Y(n-1)/h(n-1)+fnn*h(n-1)/2;for k=2:(n-1) v(k)=h(k)/(h(k)+h(k-1); u(k)=h(k-1)/(h(k)+h(k-1); g(k)=3*u(k)*(Y(k+1)-Y(k)/h(k)+3*v(k)*(Y(k)-Y(k-1)/h(k-1);endu(1)=1;v(1)=;v(n-1)=1;A=2*eye(n)+diag(u,1)+diag(v,-1);m=Ag;syms x Sx;for k=1:(n-1) fprintf(%gx%g,X(k),X(k+1) Sx=(h(k)+2*(

11、x-X(k)*(x-X(k+1)2*Y(k)/h(k)3+. (h(k)-2*(x-X(k+1)*(x-X(k)2*Y(k+1)/h(k)3+. (x-X(k)*(x-X(k+1)2*m(k)/h(k)2+. (x-X(k+1)*(x-X(k)2*m(k+1)/h(k)2;Sx=vpa(expand(Sx),4)End运算结果：0.25x0.3 Sx = 4.699*x2-0.2051*x+0.3555-6.265*x3 0.3x0.39 Sx = -2.633*x2+1.995*x+0.1355+1.881*x3 0.39x0.45 Sx = 0.1064*x2+0.9261*x+0.274

12、4-0.4600*x3 0.45x0.53 Sx = -3.409*x2+2.508*x+0.3710e-1+2.144*x3分析：三次样条函数具有良好的光滑性和几何光顺性，同时收敛性和稳定性也较好。五、编写程序构造区间上的以等分结点为插值结点的Newton插值公式，假设结点数为（包括两个端点），给定相应的函数值，插值区间和等分的份数，该程序能快速计算出相应的插值公式。以，为例计算其对应的插值公式，分别取不同的值并画出原函数的图像以及插值函数的图像，观察当增大时的逼近效果.解：matlab编程function Lagrangeclc;clear;close all;for i=1:3; if

13、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=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);f

14、f=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,-) lagrange_10=y end hold on; xlabel(x);ylabel(y); title(y=p(x);hold on;endLag_Cheb();

15、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)=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

16、-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.1 0.6);set(cheb,Color,g);set(cheb,LineStyle,:);结果输出lagrange_5 =0.5673+1.202*x4-

17、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 =0.7413-5.359*x10-0.4e-2*x9+18.96*x8-0.1321*x7-25.78*x6+0.10*x5+16.81*x4+0.5e-2*x3-5.336*x2+0.5288e-3*x分析：（1）拉格朗日插值法随着插值节点增多而更加接近真实曲线，但是还是会在区间-1，-0.4及0.4，1上，由于高次振荡产生Runge现象。（2）分段插值法则优于前者，因为分段插值在每个小区间上都可以最大限度的满足边界条件，但光滑度较低。