北理工数值分析大作业

数值分析上机作业1第1章1.1计算积分，n=9。（要求计算结果具有6位有效数字）程序：n=1:19;I=zeros(1,19);I(19)=1/2*((exp(-1)/20)+(1/20));I(18)=1/2*((exp(-1)/19)+(1/19));fori=2:10I(19-i)=1/(20-i)*(1-I(20-i));endformatlongdisp(I(1:19))结果截图及分析：在MATLAB中运行以上代码，得到结果如下图所示：当计算到数列的第10项时，所得的结果即为n=9时的准确积分值。取6位有效数字可得.21.2分别将区间[-10.10]分为100,200,400等份，利用mesh或surf命令画出二元函数z=的三维图形。程序：x=-10:0.1:10;y=-10:0.1:10;[X,Y]=meshgrid(x,y);Z=exp(-abs(X))+cos(X+Y)+1./(X.^2+Y.^2+1);subplot(2,2,1);mesh(X,Y,Z);title('步长0.1')x=-10:0.2:10;y=-10:0.2:10;[X,Y]=meshgrid(x,y);Z=exp(-abs(X))+cos(X+Y)+1./(X.^2+Y.^2+1);subplot(2,2,1);mesh(X,Y,Z);title('步长0.2')3x=-10:0.05:10;y=-10:0.05:10;[X,Y]=meshgrid(x,y);Z=exp(-abs(X))+cos(X+Y)+1./(X.^2+Y.^2+1);subplot(2,2,1);mesh(X,Y,Z);title('步长0.05')结果截图及分析：由图可知，步长越小时，绘得的图形越精确。4第2章试用MATLAB编程实现追赶法求三对角方程组的算法，并考虑梯形电路电阻问题：电路中的电流128{,,,}iii满足下列线性方程组：121232343454565676787822/252025202520252025202520250iiVRiiiiiiiiiiiiiiiiiiii设220,27VVR，求各段电路的电流量。处理思路：观察该方程的系数矩阵可知，它是一个三对角矩阵，故可运用追赶法对其进行求解。程序:fori=1:8a(i)=-2;b(i)=5;c(i)=-2;d(i)=0;enda(1)=0;b(1)=2;c(8)=0;d(1)=220/27;fori=2:8a(i)=a(i)/b(i-1);b(i)=b(i)-c(i-1)*a(i);d(i)=d(i)-a(i)*d(i-1);endd(8)=d(8)/b(8);fori=7:-1:1d(i)=(d(i)-c(i)*d(i+1))/b(i);endfori=1:8x(i)=d(i);end1x结果截图及分析：在MATLAB中运行以上代码，得到结果如下图所示：图中8个值依次为128{,,,}iii的数值。2第3章试分别用（1）Jacobi迭代法；（2）Gauss-Seidel解线性方程组1234510123412191232721735143231211743511512xxxxx迭代初始向量取(0)(0,0,0,0,0)Tx.3.1Jacobi迭代法程序：A=[101234;19-12-3;2-173-5;32312-1;4-3-5-115];b=[12;-27;14;-17;12];x0=[0;0;0;0;0];D=diag(diag(A));I=eye(5);L=-tril(A,-1);B=I-D\A;g=D\b;y=B*x0+g;n=1;whilenorm(y-x0)=1.0e-6x0=y;y=B*x0+g;n=n+1;endfprintf('%8.6f\n',y);n3结果截图及分析：得到此结果时迭代次数为67次，达到精度要求。3.2Gauss-Seidel迭代法：程序：A=[101234;19-12-3;2-173-5;32312-1;4-3-5-115];b=[12;-27;14;-17;12];x0=[0;0;0;0;0];D=diag(diag(A));U=-triu(A,1);L=-tril(A,-1);M=(D-L)\U;g=(D-L)\b;y=M*x0+g;n=1;whilenorm(y-x0)=1.0e-6x0=y;y=M*x0+g;n=n+1;endfprintf('%8.6f\n',y);4结果截图及分析：Gauss-Seidel迭代法只需要迭代38次即可满足精度要求。5第4章设A=162621666612，取先用幂法迭代3次，得到A的按模最大特征值的近似值，取为其整数部分，再用反幂法计算A的按模最大特征值的更精确的近似值，要求误差小于.程序：A=[126-6;6162;-6216];x0=[1;1;1];y=x0;b=max(abs(x0));k=1;while(k4)x=A*y;b=max(abs(x));y=x./b;k=k+1;fprintf('eig1equals%6.4f\n',b);endbb0=fix(b);I=eye(3,3);x0=[1;1;1];y=x0;l=0;bb=max(abs(x0));k=1;while(abs(bb-l)=1.0e-10)l=bb;x=(A-bb0*I)\y;bb=max(abs(x));y=x./bb;eig=l+b;fprintf('eig2(%d)equals%12.10f\n',k,eig);k=k+1;end6实验截图及分析：由图可知，由幂法3次迭代后得到的特征值为19.4，而由反幂法得到的特征值为20.3999999999.误差小于7第5章试编写MATLAB函数实现Newton插值，要求能输出插值多项式。对函数f(x)=在区间[-5,5]上实现10次多项式插值。要求：（1）输出插值多项式。（2）在区间[-5,5]内均匀插入99个节点，计算这些节点上函数f（x）的近似值，并在同一张图上画出原函数和插值多项式的图形。（3）观察龙格现象，计算插值函数在各节点处的误差，并画出误差图。5.1输出插值多项式程序：x=-5:1:5;y=1./(1+4*(x.^2));newpoly(x,y)function[c,d]=newpoly(x,y)n=length(x);d=zeros(n,n);d(:,1)=y';forj=2:nfork=j:nd(k,j)=(d(k,j-1)-d(k-1,j-1))/(x(k)-x(k-j+1));endendc=d(n,n);fork=(n-1):-1:1c=conv(c,poly(x(k)));m=length(c);c(m)=c(m)+d(k,k);end8end结果及分析：ans=Columns1through2-0.000049595763049Columns3through40.0027401659084830.000000000000000Columns5through6-0.0514215070767200.000000000000000Columns7through80.3920149852823120.000000000000000Columns9through10-1.1432840483510250.000000000000001Column111.00000000000000010次插值多项式由高到低系数为Columns1至Column115.2原函数与插值多项式的图形程序：x=-5:1:5;y=1./(1+4*(x.^2));n=newpoly(x,y);x0=-5:0.1:5;y0=1./(1+4*(x0.^2));9vn=polyval(n,x0);plot(x0,vn,'-r',x0,y0,'--b');xlabel('x');ylabel('y');实验结果截图：-5-4-3-2-1012345-1012345xy原函数与插值多项式的图形如上图所示，蓝色为原函数的图形，红色为插值多项式的图形。5.3各节点的误差及误差图程序：formatlong;x=-5:1:5;y=1./(1+4*(x.^2));n=newpoly(x,y);x0=-5:0.1:5;y0=1./(1+4*(x0.^2));vn=polyval(n,x0);plot(x0,y0-vn,'-r');xlabel('x');ylabel('y');10实验结果截图：-5-4-3-2-1012345-5-4-3-2-101xy误差图如上图所示。11第6章炼钢厂出钢时所用的盛钢水的钢包，在使用过程中由于钢液及炉渣对包衬耐火材料的腐蚀，使其容积不断加大。经试验，钢包的容积与相应的使用次数的数据列表如下：选用双曲线xbay11对数据进行拟合，使用最小二乘法求出拟合函数，作出拟合曲线图。处理思路：用Y替代1/y，用X替代1/x，原曲线化为Y=a+bx，双曲线转化为一次线性方程，使用最小二乘法求出该一次方程的系数。使用次数x容积y使用次数x容积y2106.4211110.593108.2612110.605109.5814110.726109.5016110.907109.8617110.769110.0019111.1010109.9320111.3012程序：x=[2356791011121416171920];y=[106.42108.26109.58109.5109.86110109.93110.59110.60110.72110.9110.76111.1111.3];k1=0;k2=0;k3=0;k4=0;fori=1:14k1=k1+1/x(i);endfori=1:14k2=k2+1/y(i);endfori=1:14k3=k3+1/(x(i))^2;endfori=1:14k4=k4+1/(x(i)*y(i));endb=(k1*k2-14*k4)/(k1^2-14*k3)a=k2/14-k1*b/14plot(x,y,'r*')holdonx=2:0.01:20;y=1./(a+b./x);plot(x,y)xlabel('x')ylabel('y')gridon实验结果截图与分析：13即最小二乘法求出拟合函数为：=0.008973+0.000842拟合曲线图为：14第7章考纽螺线的形状象钟表的发条，也称回旋曲线，它在直角坐标系中的参数方程为ssdtatsydtatsx020221sin)(21cos)(曲线关于原点对称，取a=1，参数s的变化范围[-5,5]，容许误差限分别是和。选取适当的节点个数，利用数值积分方法计算曲线上点的坐标，并画出曲线的图形。误差限为时：程序：x=zeros(101,1);y=zeros(101,1);f1=inline('cos(1/2*(t.^2))');f2=inline('sin(1/2*(t.^2))');i=1;fors=-5:0.1:5x(i,1)=quad(f1,0,s,1e-6);y(i,1)=quad(f2,0,s,1e-6);i=i+1;endplot(x,y,'r-');title('误差限-1e-6');xlabel('x(s)');ylabel('y(s)');15实验截图：误差限为时得到曲线如下图所示：误差限为时：程序：x=zeros(101,1);y=zeros(101,1);f1=inline('cos(1/2*(t.^2))');f2=inline('sin(1/2*(t.^2))');i=1;fors=-5:0.1:5x(i,1)=quad(f1,0,s,1e-10);y(i,1)=quad(f2,0,s,1e-10);i=i+1;end