北理工數(shù)值分析大作業(yè)

1、1.1計(jì)算積分In=01xnex-1dx，n=9。（要求計(jì)算結(jié)果具有6位有效數(shù)字）程序：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)結(jié)果截圖及分析：在MATLAB中運(yùn)行以上代碼，得到結(jié)果如下圖所示：當(dāng)計(jì)算到數(shù)列的第10項(xiàng)時(shí)，所得的結(jié)果即為In=01xnex-1dx，n=9時(shí)的準(zhǔn)確積分值。取6位有效數(shù)字可得I9=0.0916122.1.2分別將區(qū)間-10

2、.10分為100,200,400等份，利用mesh或surf命令畫出二元函數(shù)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')結(jié)果截圖及分析：由圖可知，步長越小時(shí)，繪得的圖形越精確。試用MATLAB編程實(shí)現(xiàn)追趕法求三對角方程組的算法，并考慮梯形電路電阻問題：電路中的電流滿足下列線性方程組：設(shè)

4、，求各段電路的電流量。處理思路：觀察該方程的系數(shù)矩陣可知，它是一個(gè)三對角矩陣，故可運(yùn)用追趕法對其進(jìn)行求解。程序: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結(jié)果截圖及分析：在MA

5、TLAB中運(yùn)行以上代碼，得到結(jié)果如下圖所示：圖中8個(gè)值依次為的數(shù)值。試分別用（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結(jié)果截圖及分析：得到此結(jié)果時(shí)迭代次數(shù)為67次，達(dá)到精度要求。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);結(jié)果截圖及分析：Gauss-Seidel迭代法只需要迭代38次即可滿足精度要求。設(shè)A=，取x(0)=（1,1,1）T，先用冪法迭代3次，得到A的按模最大特征值的近似值，取*為其整數(shù)部分，再用反冪法計(jì)算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實(shí)驗(yàn)截圖及分析：由圖可知，由冪法3次迭代后得到的特征值為19.4，而由反

9、冪法得到的特征值為20.3999999999.誤差小于10-10。試編寫MATLAB函數(shù)實(shí)現(xiàn)Newton插值，要求能輸出插值多項(xiàng)式。對函數(shù)f(x)=11+4x2在區(qū)間-5,5上實(shí)現(xiàn)10次多項(xiàng)式插值。要求：（1）輸出插值多項(xiàng)式。（2）在區(qū)間-5,5內(nèi)均勻插入99個(gè)節(jié)點(diǎn)，計(jì)算這些節(jié)點(diǎn)上函數(shù)f（x）的近似值，并在同一張圖上畫出原函數(shù)和插值多項(xiàng)式的圖形。（3）觀察龍格現(xiàn)象，計(jì)算插值函數(shù)在各節(jié)點(diǎn)處的誤差，并畫出誤差圖。5.1輸出插值多項(xiàng)式程序：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結(jié)果及分析：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次插值多項(xiàng)式由高到低系數(shù)為Columns 1至Column 115.2原函數(shù)與插值多項(xiàng)式的圖形程序：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');實(shí)驗(yàn)結(jié)果截圖：原函數(shù)與插值多項(xiàng)式的圖形如上圖所示，藍(lán)色為原函數(shù)的圖形，紅色為插值多項(xiàng)式的圖形。5.3各節(jié)點(diǎn)的誤差及誤差圖程序：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');實(shí)驗(yàn)結(jié)果截圖：誤差圖如上圖所示。煉鋼廠出鋼時(shí)所用的盛鋼水的鋼包，在使用過程中由于鋼液及爐渣對包襯耐火材料的腐蝕，使其容積不斷加大。經(jīng)試驗(yàn)，鋼包的容積與相應(yīng)的使用次數(shù)的數(shù)據(jù)列表如下：使用次數(shù) x 容積y使用次數(shù) x容積y2106.4211110.593108.2612110.605109.5814110.726109.5016110.907109.8617110.769110.0019111.1010109.9320111.30選用雙曲線對數(shù)據(jù)進(jìn)行擬合，使用最小二乘法求出擬合函數(shù)，作出擬合曲線圖。處理思路：用Y替代1/y，用X替代1/x，

14、原曲線化為Y=a+bx，雙曲線轉(zhuǎn)化為一次線性方程，使用最小二乘法求出該一次方程的系數(shù)。程序：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實(shí)驗(yàn)結(jié)果截圖與分析：即最小二乘法求出擬合函數(shù)為：1y=0.008973+0.0008421x擬合曲線圖為：考紐螺線的形狀象鐘表的發(fā)條，也稱回旋曲線，它在直角坐標(biāo)系中的參數(shù)方程為曲線關(guān)于原點(diǎn)對稱，取a=1，參數(shù)s的變化范圍-5,5，容許誤差限分別是10-6和10-10

16、。選取適當(dāng)?shù)墓?jié)點(diǎn)個(gè)數(shù)，利用數(shù)值積分方法計(jì)算曲線上點(diǎn)的坐標(biāo)，并畫出曲線的圖形。誤差限為10-6時(shí)：程序：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)');實(shí)驗(yàn)截圖：誤差限為10-6時(shí)得到曲線如下圖所示：誤差限為10-10時(shí)：程序：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)');實(shí)驗(yàn)結(jié)果截圖：求方程的非零根。程序：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結(jié)果截圖及分析：取初值為760，迭代9次，兩個(gè)非零根為765.48及-765.48.設(shè)有常微分方程的

19、初值問題其精確解為。選取步長h使四階Adams預(yù)測-校正算法和經(jīng)典RK法均穩(wěn)定，分別用這兩種方法求解微分方程，將數(shù)值解與精確解進(jìn)行比較，輸出結(jié)果。其中多步法需要的初值由經(jīng)典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(' 步長 經(jīng)典RK法 精確解');disp(X' Y' Y1');plot(X,Y1,'k*',X,Y,'-r');grid on;title('經(jīng)典RK法解非線性方程');legend('精確解 ','經(jīng)典RK法');實(shí)驗(yàn)結(jié)果截圖：四階Adams預(yù)測-校正算法程序：f=inlin