09级矩阵与数值分析上机作业

《矩阵与数值分析》课程数值实验报告矩阵05班20906284马强土木水利学院道路与铁道工程专业任课教师：董波一、解线性方程组1．问题(1)分别Jacobi迭代法和Gauss-Seidel迭代法求解教材85页例题。迭代法计算停止的条件为：6)()1(3110maxkjkjjxx;(2).用Gauss列主元消去法、QR方法求解如下方程组：123424311282006,504023900516xxxx2．程序及计算结果(1)a.源代码：clearallA=[3.4336-0.523800.67105-0.15270;-0.523803.28326-0.73051-0.26890;0.67105-0.730514.02612-0.09835;-0.15272-0.268900.018352.75702];%A-方程组系数矩阵b=[-1.0;1.5;2.5;-2.0];%b-方程组的矩阵表达形式Ax=b的右端列向量表c=0.000001;%终止条件N=length(b);x0j=zeros(N,1);%Jacobi迭代法初始值x0g=zeros(N,1);%Gauss-Seidel迭代法初始值kj=0;kg=0;%迭代次数计数D=zeros(N,N);L=zeros(N,N);U=zeros(N,N);%将A分解为D-L-U的形式fori=1:ND(i,i)=A(i,i);endfori=2:Nforj=1:i-1;L(i,j)=-A(i,j);endendfori=1:N-1forj=i+1:N;U(i,j)=-A(i,j);endend%Jacobi迭代法矩阵形式Bj=inv(D)*(L+U);fj=inv(D)*b;%Gauss-Seidel迭代法矩阵形式Bg=inv(D-L)*U;fg=inv(D-L)*b;%Jacobi迭代fori=1:100xj=Bj*x0j+fj;kj=kj+1;ifmax(xj-x0j)=cxjkjbreak;endx0j=xj;end%Gauss-Seidel迭代fori=1:100xg=Bg*x0g+fg;kg=kg+1;ifmax(xg-x0g)=cxgkgbreak;endx0g=xg;endb.计算结果:xj=，迭代步数kj=11xg=,迭代步数kg=6分析：G—S法应用了上一步计算出来的较精确的值，所以可以用较少的迭代步数满足误差要求。(2)a.源代码：A=[2431;8200;5040;9005];b=[1262316];[N,N]=size(A);x=zeros(N,1);y=zeros(N,1);z=zeros(N,1);R=1:N;forp=1:N-1[max1,j]=max(abs(A(p:N,p)));c=A(p,:);A(p,:)=A(j+p-1,:);A(j+p-1,:)=c;d=b(p);b(p)=b(j+p-1);b(j+p-1)=d;ifA(p,p)==0break;end;fork=p+1:Nmult=A(k,p)/A(p,p);b(k)=b(k)-mult*b(p);A(k,p:N)=A(k,p:N)-mult*A(p,p:N);end;end;fork=N:-1:1x(k)=(b(k)-A(k,k+1:N)*x(k+1:N))/A(k,k)endb.计算结果:x=二、非线性方程的迭代解法1、问题：（1）．用Newton迭代法、割线法求方程01xxexf在5.0x附近的正．计算停止的条件为：6110kkxx；（2）．利用迭代法求多项式354214.66213.092154.2072631.9296013.7=xxxxx的所有零点。2、程序及计算结果：(1)a.源代码：①Newton迭代法clearallx0=0.5;k=0;c=0.000001;fori=1:10000f=x0*exp(x0)-1;df=(x0+1)*exp(x0);x=x0-f/df;k=k+1;ifabs(x-x0)=cbreak;endx0=x;endxk②割线法clearallx0=0.5;k=0;c=0.000001;fori=1:10000f=x0*exp(x0)-1;df=(x0+1)*exp(x0);x=x0-f/df;k=k+1;ifabs(x-x0)=cbreak;endx0=x;endxk②割线法clearallx0=0.5;x1=0.7;k=0;c=0.000001;fori=1:10000x=x1-(x1*exp(x1)-1)*(x1-x0)/((x1*exp(x1)-1)-(x0*exp(x0)-1));;k=k+1;ifabs(x-x1)=cbreak;endx0=x1;x1=x;endxkb.计算结果:①Newton迭代法x=k=4②割线法x=k=5分析：割线法的收敛阶小于Newton迭代法，所以用了较多的迭代步数以满足误差要求。(2)a.源代码：①clearallx=-5:0.5:11fori=1:33f(i)=x(i)^5-13.7*x(i)^4+14.66*x(i)^3+213.092*x(i)^2-154.2072*x(i)-631.9296endplot(x,f)gridon②functionx1=Newton(b)x=b;c=0.000001;fori=1:10000f=x^5-13.7*x^4+14.66*x^3+213.092*x^2-154.2072*x-631.9296;df=5*x^4-54.8*x^3+43.98*x^2+426.184*x-15402072;x1=x-f/df;ifabs(x1-x)=cbreak;endx=x1;endx1b.计算结果:①-6-4-2024681012-10000-8000-6000-4000-2000020004000②由上图选取初始值：-3.5、-2、2、5、10.6由Newton迭代法得到方程的五个解：x=三、数值积分1、问题：用Romberg方法求dxex1022的近似值，要求误差不超过510.2、程序及计算结果：a.源代码：clearalla=0;b=1;h=(b-a);T(1,1)=(h/2)*1/sqrt(pi)*(exp(-(a^2)/2)+exp(-(b^2)/2));k=2;while(k=2)sum=0;forp=1:2^(k-2)sum=sum+1/sqrt(pi)*exp(-(a+(2*p-1)*h/2^(k-1))^2/2);endT(1,k)=0.5*T(1,k-1)+sum*h/2^(k-1);form=2:kT(m,k-m+1)=(4^(m-1)*T(m-1,k-m+2)-T(m-1,k-m+1))/(4^(m-1)-1);if(abs(T(m,1)-T(m-1,1))=1.0e-5)k=0;endendifk~=0k=k+1;endendT(m,1)m-1abs(T(m,1)-T(m-1,1))b.计算结果:I=，加速次数M=3，绝对误差c=分析：应用加速法可以使之前不收敛的式子收敛，有效的加快了收敛速度。四、插值与逼近1、问题：(1)．给定1,1上的函数22511xxf，请做如下的插值逼近：⑴构造等距节点分别取5n，8n，10n的Lagrange插值多项式；⑵构造分段线性取10n的Lagrange插值多项式；⑶取Chebyshev多项式xnxTnarccoscos的零点：nkxk212cos，nk,,1作插值节点构造10n的插值多项式xf和上述的插值多项式均要求画出曲线图形（用不同的线型或颜色表示不同的曲线）。(2)．已知函数值kx012345678910ky00.791.532.192.713.033.272.893.063.193.29和边界条件：(0)0.8,(10)0.2ss．求三次样条插值函数()ysx并画出其图形．(3)．分别用axbey和2cxbxay拟合下列数据，并比较两种方法得到的误差。kx12345678ky15.320.527.436.649.165.687.8117.62、程序及计算结果：(1)a.源代码：clearalla=-1;b=1;n1=5;n2=8;n3=10;x00=a:(b-a)/(n1-1):b;x01=a:(b-a)/(n2-1):b;x02=a:(b-a)/(n3-1):b;x=a:(b-a)/99:b;fori=1:100f(i)=1/(1+25*x(i)*x(i));endfori=1:n1y00(i)=1/(1+25*(x00(i))^2);endfori=1:n2y01(i)=1/(1+25*(x01(i))^2);endfori=1:n3y02(i)=1/(1+25*(x02(i))^2);endfori=1:n3y03(i)=1/(1+25*(x02(i))^2);endfori=1:n3xk(i)=cos((2*i-1)/(2*n3));endp10=lagrange(x00,y00,x);p11=lagrange(x01,y01,x);p12=lagrange(x02,y02,x);p20=interp1(x02,y02,x);p30=lagrange(x02,y02,xk);plot(x,f,'r',x,p10,'y',x,p11,'g',x,p12,'m',x,p20,'k',xk,p30)axis([-1,1,-0.5,1.2])xlabel('x,xk'),ylabel('f,p10,p11,p12,p20,p30');gtext('f'),gtext('p10'),gtext('p11'),gtext('p12'),gtext('p20'),gtext('p30')b.计算结果:对f,p10,p11,p12,p20取插值点数x=100，对p30取n＝10点的Chebyshev多项式零点为插值点，绘出各插值多项式的图形如下。分析：1.拉格朗日插值法随着插值节点的选择的增多而更加接近真实曲线，但是还是会在区间［-1，-0.4］及［0.4，1］上，产生容格现象。2.分段插值法则优于前者，因为分段插值在每个小区间上都可以最大限度的满足边界条件。3.利用切片比雪夫迭代法可以更接近真实曲线。(2)a.源代码：formatlongg;xk=[012345678910];%给定节点y=[00.791.532.192.713.033.272.893.063.193.29];%节点处的值h=zeros(1,10);%分别设置初始状态为0向量c=zeros(1,9);d=zeros(1,9);g=zeros(1,9);b=zeros(1,9);%设置初始状态均为0向量m2=zeros(9,1);m1=zeros(9,1);m0=0.8;m10=0.2;%已知的两个边界条件fori=1:1:10h(i)=xk(i+1)-xk(i);%算出hkendfori=1:1:9c(i)=h(i+1)/(h(i+1)+h(i));d(i)=h(i)/(h(i+1)+h(i));endfori=1:1:9g(i)=3*(d(i)*(y(i+2)-y(i+1))/h(i)+c(i)*(y(i+1)-y(i))/h(i));endv=[222222222];%对角元均为2A=diag(v);%生成矩阵fori=1:1:8A(i,i+1)=d(i);A(i+1,i)=c(i+1);b(i)=g(i);endb(1)=g(1)-c(1)*m0;b(9)=g(9)-d(9)*m10;[L,U]=lu(A);%将A进行LU分解，解方程，得出m向量m1=L\b';%LU法求解m2=U\m1;m(1)=m0;m(11)=m10;fori=2:1:10m(i)=m2(i-1);endfori=1:1:10;%计算S(x)的格式x=i-