《应用计算方法教程》matlab作业二

作业六6-1试验目的计算特征值，实现算法试验内容：随机产生一个10阶整数矩阵，各数均在-5和5之间。（1）用MATLAB函数“eig”求矩阵全部特征值。（2）用幂法求A的主特征值及对应的特征向量。（3）用基本QR算法求全部特征值（可用MATLAB函数“qr”实现矩阵的QR分解）。原理幂法：设矩阵A的特征值为12n||||||并设A有完全的特征向量系12,,,n(它们线性无关)，则对任意一个非零向量0nVR所构造的向量序列1kkVAV有11()lim()kjkkjVV，其中()kjV表示向量的第j个分量。为避免逐次迭代向量kV不为零的分量变得很大（1||1时）或很小（1||1时），将每一步的kV按其模最大的元素进行归一化。具体过程如下：选择初始向量0V，令1max(),,,1kkkkkkkVmVUVAUkm，当k充分大时1111,max()max()kkUV。QR法求全部特征值：11111222111,1,2,3,kkkkkAAQRRQAQRkRQAQR由于此题的矩阵是10阶的，上述算法计算时间过长，考虑采用改进算法——移位加速。迭代格式如下：1kkkkkkkkAqIQRARQqI计算kA右下角的二阶矩阵()()1,11,()(),1,kknnnnkknnnnaaaa的特征值()()1,kknn，当()()1,kknn为实数时，选kq为()()1,kknn中最接近(),knna的。程序A=-5+round(10*rand(10));[V,D]=eig(A)[lamdau]=lab6_2_power(A,[1;1;1;1;1;1;1;1;1;1],10^(-5),1000)d=lab6_3_qr2(A,10^(-5))function[lamdau]=lab6_2_power(a,v,eps,N)lamda=0;err=1;k=1;while(k=N&&erreps)u=a*v;[mj]=max(abs(u));dc=abs(lamda-m);u=u/m;dv=norm(u-v);err=max(dc,dv);v=u;lamda=m;k=k+1;endfunctionD=lab6_3_qr2(A,eps)[n,n]=size(A);m=n;D=zeros(n,1);B=A;while(m1)while(abs(B(m,m-1))=eps*(abs(B(m-1,m-1))+abs(B(m,m))))S=eig(B(m-1:m,m-1:m));[j,k]=min([abs(B(m,m)-S(1)),abs(B(m,m)-S(2))]);[Q,U]=qr(B-S(k)*eye(m));B=U*Q+S(k)*eye(m);endA(1:m,1:m)=B;m=m-1;B=A(1:m,1:m);endD=diag(A);界面(1)(2)(3)作业七7-1试验目的：熟悉代数插值试验内容：已知在f(x)在7个点的函数值如下表所示，分别使用拉格朗日插值法和牛顿插值法求f(0.596)与f(0.906)的近似值。ix0.40.50.60.70.80.91.0iy11.751.962.192.442.713.00原理拉格朗日插值多项式：01110011100()()()()()()()()()()()()njjnnjjjjjjjjjnnnijjijiijxxxxxxxxxxLxyxxxxxxxxxxxxyxx牛顿插值多项式：001001201001()()[,]()[,,]()()[,,]()()nnnNxfxfxxxxfxxxxxxxfxxxxxx其中00011()[,,]()()()()kjkjjjjjjjkfxfxxxxxxxxxx。程序functiony1=lab7_1_Lagrange(x,y,x1)y1=0;[mn]=size(x);n=n-1;fork=1:n+1t=1;fori=1:n+1if(i~=k)t=t*(x1-x(i))/(x(k)-x(i));endendy1=y1+t*y(k);endfunctiony1=lab7_2_Newton(x,y,x1)[mn]=size(x);n=n-1;forj=1:nfori=n+1:-1:j+1y(i)=(y(i)-y(i-1))/(x(i)-x(i-j));endendy1=y(n+1);forj=n:-1:1y1=y(j)+(x1-x(j))*y1;end界面作业八8-1试验目的：熟悉最小二乘法拟合多项式试验内容：给定数据点(ix,iy)，X0.400.550.650.800.901.05f(x)0.410750.578150.696750.888111.026521.25386用3次最小二乘多项式拟合数据，并求平方误差。原理要作三次最小二乘拟合，令00112233()()()()(),()jjSxaxaxaxaxxx，计算法方程Gad，其中000101011101(,)(,)(,)(,)(,)(,)(,)(,)(,)nnnnnnG，01naaaa，12ndddd，00(,)()()(),(,0,1,,),()1(,)()()(),(0,1,,)mjkijikiiimkkiikiixxxjknxdfxfxxkn。其平方误差为20()[()()]miiiisxSxfx。程序x=[0.40.550.650.800.901.05];f=[0.410750.578150.696750.888111.026521.25386];G=zeros(4,4);forj=1:4fork=1:4fori=1:6G(j,k)=G(j,k)+x(i)^(j+k-2);endendendd=zeros(4,1);fork=1:4fori=1:6d(k,1)=d(k,1)+f(i)*x(i)^(k-1);endenda=G\ds=0;fori=1:6s=s+(f(i)-(a(1)+a(2)*x(i)+a(3)*x(i)^2+a(4)*x(i)^3))^2;ends界面作业九9-1试验目的：熟悉数值积分公式，掌握数值计算定积分的方法试验内容：采用不同方法数值计算积分10ln(1)xdxx编写复合梯形公式和复合Simpson公式通用子程序，分别采用4，8，16，32，64等分区间计算。原理复合梯形公式：将区间[a,b]作n等分，bahn，结点(0)ixaihin，11101[()()][()2()()]22nnniiiiihhTfxfxfafxfb复合Simpson公式：将区间[a,b]作2n等分，记2bahn，111222122212001[()4()()][()4()2()()]33nnnniiiiiiiihhSfxfxfxfafxfxfb程序functiony=lab9_f(x)y=(log(1+x))/x;functiony=lab9_1_fTrapezoid(a,b,eps,n)f0=0;h=(b-a)/n;fori=0:n-1f0=f0+lab9_f(a+i*(b-a)/n)+lab9_f(a+(i+1)*(b-a)/n);endy=f0*h/2;functiony=lab9_2_fSimpson(a,b,eps,n)f0=0;h=(b-a)/n;k=n/2;fori=0:k-1f0=f0+lab9_f(a+2*i*(b-a)/n)+4*lab9_f(a+(2*i+1)*(b-a)/n)+lab9_f(a+(2*i+2)*(b-a)/n);endy=f0*h/3;界面作业十10-1试验目的：学会用Euler法、改进Euler法、经典的4阶Runge-Kutta法求解常微分方程初值问题。试验内容：分别用1)Euler法（步长h=0.025）2)改进Euler法（步长h=0.05）3)4阶Runge-Kutta（步长h=0.1）求解下面的初值问题：(1)(3),02yyyx比较公共节点解的误差。精确解为2132(1)xye。原理Euler法：令(,())()fxyxyx，1(,)nnnnyyhfxy（banh）。改进Euler法（梯形公式）：121112(,)(,)()2nnnnnnkfxykfxhyhkhyykk4阶Runge-Kutta：121324311234(,)1(,)221(,)22(,)1(22)6nnnnnnnnnnkhfxyhkhfxykhkhfxykkhfxhykyykkkk程序y0=-2;h=0.025;n=2/h;y(1)=y0-h*(y0+1)*(y0+3);fori=2:ny(i)=y(i-1)-h*(y(i-1)+1)*(y(i-1)+3);endyforj=1:ny1(j)=-3+2/(1+exp(-4*j/n));endy1x=0.025:0.025:2;plot(x,y1,'r',x,y,'b')y0=-2;h=0.05;n=2/h;k1=-(y0+1)*(y0+3);k2=-(y0+h*k1+1)*(y0+h*k1+3);y21(1)=y0+h/2*(k1+k2);fori=2:nk1=-(y21(i-1)+1)*(y21(i-1)+3);k2=-(y21(i-1)+h*k1+1)*(y21(i-1)+h*k1+3);y21(i)=y21(i-1)+h/2*(k1+k2);endy21forj=1:ny22(j)=-3+2/(1+exp(-4*j/n));endy22x=0.05:0.05:2;plot(x,y21,'r',x,y22,'b')y0=-2;h=0.1;n=2/h;k1=-h*(y0+1)*(y0+3);k2=-h*(y0+k1/2+1)*(y0+k1/2+3);k3=-h*(y0+k2/2+1)*(y0+k2/2+3);k4=-h*(y0+k3+1)*(y0+k3+3);y31(1)=y0+(k1+2*k2+2*k3+k4)/6;fori=2:nk1=-h*(y31(i-1)+1)*(y31(i-1)+3);k2=-h*(y31(i-1)+k1/2+1)*(y31(i-1)+k1/2+3);k3=-h*(y31(i-1)+k2/2+1)*(y31(i-1)+k2/2+3);k4=-h*(y31(i-1)+k3+1)*(y31(i-1)+k3+3);y31(i)=y31(i-1)+(k1+2*k2+2*k3+k4)/6;endy31forj=1:ny32(j)=-3+2/(1+exp(-4*j/n));endy32x=0.1:0.1:2;plot(x,y31,'r',x,y32,'b')界面