北理工数值分析大作业

1、1.1计算积分In=01xnex-1dx，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项时，所得的结果即为In=01xnex-1dx，n=9时的准确积分值。取6位有效数字可得I9=0.0916122.1.2分别将区间-10

2、.10分为100,200,400等份，利用mesh或surf命令画出二元函数z=e-|x|+cosx+y+1x2+y2+1的三维图形。程序：>> 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

3、)+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')结果截图及分析：由图可知，步长越小时，绘得的图形越精确。试用MATLAB编程实现追赶法求三对角方程组的算法，并考虑梯形电路电阻问题：电路中的电流满足下列线性方程组：设

4、，求各段电路的电流量。处理思路：观察该方程的系数矩阵可知，它是一个三对角矩阵，故可运用追赶法对其进行求解。程序:for i=1:8 a(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:8x(i)=d(i);endx结果截图及分析：在MA

5、TLAB中运行以上代码，得到结果如下图所示：图中8个值依次为的数值。试分别用（1）Jacobi迭代法；（2）Gauss-Seidel解线性方程组迭代初始向量取.3.1 Jacobi迭代法程序：>> A=10 1 2 3 4; 1 9 -1 2 -3; 2 -1 7 3 -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-DA;g=Db;y=B*x0+g;n=1;while norm(y-x0)>=1.0e-6 x0=y

6、; y=B*x0+g; n=n+1;endfprintf('%8.6fn',y);n结果截图及分析：得到此结果时迭代次数为67次，达到精度要求。3.2 Gauss-Seidel迭代法：程序：>> A=10 1 2 3 4; 1 9 -1 2 -3; 2 -1 7 3 -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)

7、>=1.0e-6 x0=y; y=M*x0+g; n=n+1;endfprintf('%8.6fn',y);结果截图及分析：Gauss-Seidel迭代法只需要迭代38次即可满足精度要求。设A=，取x(0)=（1,1,1）T，先用幂法迭代3次，得到A的按模最大特征值的近似值，取*为其整数部分，再用反幂法计算A的按模最大特征值的更精确的近似值，要求误差小于10-10.程序：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

8、=k+1;fprintf('eig1 equals %6.4fn',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.10fn',k, eig);k=k+1;end实验截图及分析：由图可知，由幂法3次迭代后得到的特征值为19.4，而由反

9、幂法得到的特征值为20.3999999999.误差小于10-10。试编写MATLAB函数实现Newton插值，要求能输出插值多项式。对函数f(x)=11+4x2在区间-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

10、(n,n);d(:,1)=y'for j=2:n for k=j:n d(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:1 c=conv(c,poly(x(k); m=length(c); c(m)=c(m)+d(k,k);endend结果及分析：ans = Columns 1 through 2 -0.000049595763049 Columns 3 through 4 0.002740165908483 0.000000000000000 Columns 5 through 6 -0

11、.051421507076720 0.000000000000000 Columns 7 through 8 0.392014985282312 0.000000000000000 Columns 9 through 10 -1.143284048351025 0.000000000000001 Column 11 1.00000000000000010次插值多项式由高到低系数为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

12、=polyval(n,x0);plot(x0,vn,'-r',x0,y0,'-b');xlabel('x');ylabel('y');实验结果截图：原函数与插值多项式的图形如上图所示，蓝色为原函数的图形，红色为插值多项式的图形。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('

13、;x');ylabel('y');实验结果截图：误差图如上图所示。炼钢厂出钢时所用的盛钢水的钢包，在使用过程中由于钢液及炉渣对包衬耐火材料的腐蚀，使其容积不断加大。经试验，钢包的容积与相应的使用次数的数据列表如下：使用次数 x 容积y使用次数 x容积y2106.4211110.593108.2612110.605109.5814110.726109.5016110.907109.8617110.769110.0019111.1010109.9320111.30选用双曲线对数据进行拟合，使用最小二乘法求出拟合函数，作出拟合曲线图。处理思路：用Y替代1/y，用X替代1/x，

14、原曲线化为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:14 k1=k1+1/x(i);endfor i=1:14 k2=k2+1/y(i);endfor i=1:14 k3=k3+1/(x(i)2;endfor i=1:14 k4=k4+

15、1/(x(i)*y(i);endb=(k1*k2-14*k4)/(k12-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实验结果截图与分析：即最小二乘法求出拟合函数为：1y=0.008973+0.0008421x拟合曲线图为：考纽螺线的形状象钟表的发条，也称回旋曲线，它在直角坐标系中的参数方程为曲线关于原点对称，取a=1，参数s的变化范围-5,5，容许误差限分别是10-6和10-10

16、。选取适当的节点个数，利用数值积分方法计算曲线上点的坐标，并画出曲线的图形。误差限为10-6时：程序：x=zeros(101,1);y=zeros(101,1);f1=inline('cos(1/2*(t.2)');f2=inline('sin(1/2*(t.2)');i=1;for s= -5:0.1:5 x(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)

17、9;);ylabel('y(s)');实验截图：误差限为10-6时得到曲线如下图所示：误差限为10-10时：程序：x=zeros(101,1);y=zeros(101,1);f1=inline('cos(1/2*(t.2)');f2=inline('sin(1/2*(t.2)');i=1;for s= -5:0.1:5 x(i,1)=quad(f1,0,s,1e-10); y(i,1)=quad(f2,0,s,1e-10); i=i+1;endplot(x,y,'r-');title('误差限-1e-10');xl

18、abel('x(s)');ylabel('y(s)');实验结果截图：求方程的非零根。程序：x=sym('x');f=sym('log(513+0.6651*x)/(513-0.6651*x)-x/(1400*0.0981)');df=diff(f,x);FX=x-f/df;Fx=inline(FX);format long;i=1;x0=760;for i=1:10 disp(x0); x0=feval(Fx,x0); endx0结果截图及分析：取初值为760，迭代9次，两个非零根为765.48及-765.48.设有常微分方程的

19、初值问题其精确解为。选取步长h使四阶Adams预测-校正算法和经典RK法均稳定，分别用这两种方法求解微分方程，将数值解与精确解进行比较，输出结果。其中多步法需要的初值由经典RK法提供。程序：f=inline('-y+2*cos(x)','x','y');n=20;Y=zeros(1,n+1);Y(1)=1;h=0.05*pi;X=0:h:pi;for i=1:20 K1=f(X(i),Y(i); K2=f(X(i)+h/2,Y(i)+(h*K1)/2); K3=f(X(i)+h/2,Y(i)+(h*K2)/2); K4=f(X(i)+h,Y(i)+h*K3); Y(i+1)=Y(i)+h*(K1+2*K2+2*K3+K4)/6;end;Y1=cos(X)+sin(X);disp(' 步长 经典RK法 精确解');disp(X' Y' Y1');plot(X,Y1,'k*',X,Y,'-r');grid on;title('经典RK法解非线性方程');legend('精确解 ','经典RK法');实验结果截图：四阶Adams预测-校正算法程序：f=inlin