数值分析上机作业第 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));for i=2:10I(19-i)=1/(20-i)*(1-I(20-i));endformat longdisp(I(1:19))结果截图及分析：在MATLAB中运行以上代码，得到结果如下图所示：当计算到数列的第10项时，所得的结果即为n=9时的准确积分值。

取6位有效数字可得.1.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')>>x = -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')结果截图及分析：由图可知，步长越小时，绘得的图形越精确。

第 2 章试用MATLAB 编程实现追赶法求三对角方程组的算法，并考虑梯形电路电阻问题：电路中的电流128{,,,}i i i 满足下列线性方程组：121232343454565676787822/252025202520252025202520250i i V R i i i i i i i i i i i i i i i i i i i i -=-+-=-+-=-+-=-+-=-+-=-+-=-+=设220,27V V R ==Ω，求各段电路的电流量。

处理思路：观察该方程的系数矩阵可知，它是一个三对角矩阵，故可运用追赶法对其进行求解。

程序:for i=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; for i=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); for i=7:-1:1d(i)=( d(i)-c(i)*d(i+1) )/b(i); endfor i=1:8 x(i)=d(i); endx结果截图及分析：在MATLAB 中运行以上代码，得到结果如下图所示：图中8个值依次为128{,,,}i i i 的数值。

第 3 章试分别用（1）Jacobi 迭代法；（2）Gauss-Seidel 解线性方程组1234510123412191232721735143231211743511512x x x x x ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥---⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥--=⎢⎥⎢⎥⎢⎥--⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥---⎣⎦⎣⎦⎣⎦迭代初始向量取(0)(0,0,0,0,0)T x =. 3.1 Jacobi 迭代法 程序：>> A=[10 1 2 3 4;1 9 -12 -3; 2 -1 73 -5; 3 2 3 12 -1;4 -3 -5 -1 15]; 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;while norm(y-x0)>=1.0e-6 x0=y; y=B*x0+g; n=n+1; endfprintf('%8.6f\n',y); n结果截图及分析：得到此结果时迭代次数为67次，达到精度要求。

3.2 Gauss-Seidel迭代法：程序：>> A=[10 1 2 3 4;1 9 -12 -3;2 -1 73 -5;3 2 3 12 -1;4 -3 -5 -1 15];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;while norm(y-x0)>=1.0e-6x0=y;y=M*x0+g;n=n+1;endfprintf('%8.6f\n',y);结果截图及分析：Gauss-Seidel迭代法只需要迭代38次即可满足精度要求。

第 4 章设A=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--162621666612，取先用幂法迭代3次，得到A 的按模最大特征值的近似值，取为其整数部分，再用反幂法计算A 的按模最大特征值的更精确的近似值，要求误差小于.程序：A=[12 6 -6;6 16 2; -6 2 16]; x0=[1;1;1];y=x0;b=max(abs(x0));k=1; while ( k<4 )x=A*y;b=max(abs(x));y=x./b; k=k+1;fprintf('eig1 equals %6.4f\n',b); end>> bb0=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;end实验截图及分析：由图可知，由幂法3次迭代后得到的特征值为19.4，而由反幂法得到的特征值为20.3999999999.误差小于第 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';for j=2:nfor k=j:nd(k,j)=(d(k,j-1)-d(k-1,j-1))/(x(k)-x(k-j+1));endendc=d(n,n);for k=(n-1):-1:1c=conv(c,poly(x(k)));m=length(c);c(m)=c(m)+d(k,k);endend结果及分析：ans =Columns 1 through 2-0.3049Columns 3 through 40.8483 0.0000Columns 5 through 6-0.6720 0.0000Columns 7 through 80.2312 0.0000Columns 9 through 10-1.1025 0.0001Column 111.000010次插值多项式由高到低系数为Columns 1至Column 115.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));vn=polyval(n,x0);plot(x0,vn,'-r',x0,y0,'--b');xlabel('x');ylabel('y');实验结果截图：yx原函数与插值多项式的图形如上图所示，蓝色为原函数的图形，红色为插值多项式的图形。

5.3各节点的误差及误差图程序：format long;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');实验结果截图：yx 误差图如上图所示。

第 6 章炼钢厂出钢时所用的盛钢水的钢包，在使用过程中由于钢液及炉渣对包衬耐火材料的腐蚀，使其容积不断加大。

经试验，钢包的容积与相应的使用次数的数据列表如下：选用双曲线xb a y11+=对数据进行拟合，使用最小二乘法求出拟合函数，作出拟合曲线图。

处理思路：用Y 替代1/y ，用X 替代1/x ，原曲线化为Y=a+bx ，双曲线转化为一次线性方程，使用最小二乘法求出该一次方程的系数。

程序：x=[2 3 5 6 7 9 10 11 12 14 16 17 19 20];y=[106.42 108.26 109.58 109.5 109.86 110 109.93 110.59 110.60 110.72 110.9 110.76 111.1 111.3];k1=0; k2=0; k3=0; k4=0;for i=1:14k1=k1+1/x(i);endfor i=1:14k2=k2+1/y(i);endfor i=1:14k3=k3+1/(x(i))^2;endfor i=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*')hold onx=2:0.01:20;y=1./(a+b./x);plot(x,y)xlabel('x')ylabel('y')grid on实验结果截图与分析：即最小二乘法求出拟合函数为：=0.008973+0.000842拟合曲线图为：第 7 章考纽螺线的形状象钟表的发条，也称回旋曲线，它在直角坐标系中的参数方程为⎪⎪⎩⎪⎪⎨⎧==⎰⎰ssdt at s y dt at s x 02221sin )(21cos)( 曲线关于原点对称，取a=1，参数s 的变化围[-5,5]，容许误差限分别是和。