計(jì)算方法上機(jī)作業(yè)集合

1、第一次&第二次上機(jī)作業(yè)上機(jī)作業(yè)：1.在Matlab上執(zhí)行：>> 5.1-5-0.1和>> 1.5-1-0.5給出執(zhí)行結(jié)果，并簡要分析一下產(chǎn)生現(xiàn)象的原因。解：執(zhí)行結(jié)果如下： 在Matlab中，小數(shù)值很難用二進(jìn)制進(jìn)行描述。由于計(jì)算精度的影響，相近兩數(shù)相減會出現(xiàn)誤差。2.（課本181頁第一題）解：（1）n=0時(shí)，積分得I0=ln6-ln5，編寫如下圖代碼從以上代碼顯示的結(jié)果可以看出，I20的近似值為0.012712966517465（2）In=01xn5+xdx,可得01xn6dx01xn5+xdx01xn5dx,得 16(n+1)In15(n+1),則有1126I2

2、01105,取I20=1105,以此逆序估算I0。代碼段及結(jié)果如下圖所示 結(jié)果是從I19逆序輸出至I0，所以得到I0的近似值為0.088392216030227。（3）從I20估計(jì)的過程更為可靠。首先根據(jù)積分得表達(dá)式是可知，被積函數(shù)隨著n的增大，其所圍面積應(yīng)當(dāng)是逐步減小的，即積分值應(yīng)是隨著n的遞增二單調(diào)減小的，（1）中輸出的值不滿足這一條件，（2）滿足。設(shè)Sn表示In的近似值，Sn-In=(-5)n(S0-I0)(根據(jù)遞推公式可以導(dǎo)出此式),可以看出，隨著n的增大，誤差也在增大，所以順序估計(jì)時(shí)，算法不穩(wěn)定性逐漸增大，逆序估計(jì)情況則剛好相反，誤差不斷減小，算法逐漸趨于穩(wěn)定。2.（課本181頁第二

3、題）（1）上機(jī)代碼如圖所示求得近似根為0.09058（2）上機(jī)代碼如圖所示得近似根為0.09064；（3）牛頓法上機(jī)代碼如下計(jì)算所得近似解為0.09091第三次上機(jī)作業(yè)上機(jī)作業(yè)181頁第四題線性方程組為1.13483.83260.53011.78751.16513.40172.53301.54353.41294.93171.23714.99988.76431.314210.67210.0147x1x2x3x4=9.53426.394118.423116.9237（1） 順序消元法A=1.1348,3.8326,1.1651,3.4017;0.5301,1.7875,2.5330,1.5435;

4、3.4129,4.9317,8.7643,1.3142;1.2371,4.9998,10.6721,0.0147;b=9.5342;6.3941;18.4231;16.9237;上機(jī)代碼（函數(shù)部分）如下function b = gaus( A,b )%用b返回方程組的解B=A,b;n=length(b);RA=rank(A);RB=rank(B);dif=RB-RA;if dif>0 disp('此方程組無解'); returnendif RA=RB if RA=n format long; disp('此方程組有唯一解'); for p=1:n-1 fo

5、r k=p+1:n m=B(k,p)/B(p,p); B(k,p:n+1)=B(k,p:n+1)-m*B(p,p:n+1); end end %順序消元形成上三角矩陣 b=B(1:n,n+1); A=B(1:n,1:n); b(n)=b(n)/A(n,n); for q=n-1:-1:1 b(q)=(b(q)-sum(A(q,q+1:n)*b(q+1:n)/A(q,q); end %回代求解 else disp('此方程組有無數(shù)組解'); endend上機(jī)運(yùn)行結(jié)果為>> A=1.1348,3.8326,1.1651,3.4017;0.5301,1.7875,2.53

6、30,1.5435;3.4129,4.9317,8.7643,1.3142;1.2371,4.9998,10.6721,0.0147;b=9.5342;6.3941;18.4231;16.9237;>> X=gaus(A,b)此方程組有唯一解X = 1.000000000000000 1.000000000000000 1.000000000000000 1.000000000000000>>（2） 列主元消元法A=1.1348,3.8326,1.1651,3.4017;0.5301,1.7875,2.5330,1.5435;3.4129,4.9317,8.7643,1

7、.3142;1.2371,4.9998,10.6721,0.0147;b=9.5342;6.3941;18.4231;16.9237;函數(shù)部分代碼如下function b = lieZhu(A,b )%用b返回方程組的解B=A,b;RA=rank(A);RB=rank(B);n=length(b);dif=RB-RA;format long;if dif>0 disp('該方程組無解'); returnendif dif=0 if RA=n disp('該方程組有唯一解'); c=zeros(1,n); for i=1:n-1 max=abs(A(i,i)

8、; m=i; for j=i+1:n if max<abs(A(j,i) max=abs(A(j,i); m=j; end end %求出每一次消元時(shí)絕對值最大的一行的行號 if m=i for k=i:n c(k)=A(i,k); A(i,k)=A(m,k); A(m,k)=c(k); end d1=b(i); b(i)=b(m); b(m)=d1;%函數(shù)值跟隨方程一起換位置 end for k=i+1:n for j=i+1:n A(k,j)=A(k,j)-A(i,j)*A(k,i)/A(i,i); end b(k)=b(k)-b(i)*A(k,i)/A(i,i); A(k,i)=0

9、; end end %完成消元操作，形成上三角矩陣 b(n)=b(n)/A(n,n); for i=n-1:-1:1 sum=0; for j=i+1:n sum=sum+A(i,j)*b(j); end b(i)=(b(i)-sum)/A(i,i) ;%回代求解其他未知數(shù) end endelse disp('此方程組有無數(shù)組解');endend上機(jī)運(yùn)行，結(jié)果為>> A=1.1348,3.8326,1.1651,3.4017;0.5301,1.7875,2.5330,1.5435;3.4129,4.9317,8.7643,1.3142;1.2371,4.9998,1

10、0.6721,0.0147;b=9.5342;6.3941;18.4231;16.9237;X=lieZhu(A,b)該方程組有唯一解X = 1.000000000000000 1.000000000000002 0.999999999999999 0.999999999999999>>根據(jù)兩種方法運(yùn)算結(jié)果，順序消元法得到結(jié)果比列主元消元法得到的結(jié)果精度更高。（注：matlab使用的是2015b版本，不知道是保留小數(shù)位數(shù)太少，還是程序原因，順序消元輸出結(jié)果總是等于準(zhǔn)確解，請老師指正）第四次上機(jī)作業(yè)7.分析用下列迭代法解線性方程組 4-1-1 4 0-1-1 0 0 0-1 0 0-

11、1-1 0 4-1-1 4 0 -1-1 0 0-10 0 0-1-1 0 4-1-1 4 x1x2x3x4x5x6 = 0 5-2 5-2 6 的收斂性，并求出使X(k+1)-X(k)20.0001的近似解及相應(yīng)的迭代次數(shù)。（1） 雅可比迭代法解：上機(jī)編寫的函數(shù)如下function = Jacobi(X,b)%雅可比迭代法D=diag(diag(X);%得到對角線元素組成的矩陣B=inv(D)*(D-X);f=inv(D)*b;b(:,:)=0;x1=B*b+f;num=1;while(norm(x1-b)>0.0001)%判斷當(dāng)前的解是否達(dá)到精度要求 b=x1; x1=B*b+f;

12、num=num+1;end;fprintf('求得的解為：n');disp(x1);fprintf('迭代次數(shù)：%d次n',num); end上機(jī)運(yùn)行，結(jié)果如下>> A=4,-1,0,-1,0,0;-1,4,-1,0,-1,0;0,-1,4,-1,0,-1;-1,0,-1,4,-1,0;0,-1,0,-1,4,-1;0,0,-1,0,-1,4;>> b=0;5;-2;5;-2;6;>> Jacobi(A,b)求得的解為： 0.999981765074381 1.99995018125674 0.999975090628368

13、1.99995018125674 0.999975090628368 1.99996353014876迭代次數(shù)：28次滿足要求的解如輸出結(jié)果所示，總共迭代了28次（2） 高斯-賽德爾迭代法上機(jī)程序如下所示function =Gauss_Seidel( A,b )%高斯賽德爾迭代法D=diag(diag(A);L=D-tril(A);U=D-triu(A);B=inv(D-L)*U;f=inv(D-L)*b;b(:,:)=0;x0=B*b+f;num=1;while(norm(x0-b)>0.0001) num=num+1; b=x0; x0=B*b+f;end;fprintf('

14、結(jié)果為n');disp(x0);fprintf('迭代次數(shù)為：%d次n',num);end >> A=4,-1,0,-1,0,0;-1,4,-1,0,-1,0;0,-1,4,-1,0,-1;-1,0,-1,4,-1,0;0,-1,0,-1,4,-1;0,0,-1,0,-1,4;>> b=0;5;-2;5;-2;6;>> Gauss_Seidel(A,b)結(jié)果為 0.999960143810658 1.99995676152139 0.999963508299833 1.99996662162874 0.999972527179715

15、1.99998400886989迭代次數(shù)為：15次滿足精度要求的解如上述程序打印輸出所示，迭代了15次（3） SOR迭代法（w依次取1.334,1.95,0.95）上機(jī)代碼如下function = SOR(A,b,w )%SOR迭代法¨D=diag(diag(A);L=D-tril(A);U=D-triu(A);B=inv(D-w*L)*(1-w)*D+w*U);f=w*inv(D-w*L)*b;b(:,:)=0;x0=B*b+f;num=1;while(norm(x0-b)>0.0001) num=num+1; b=x0; x0=B*b+f;end;fprintf('

16、結(jié)果為n');disp(x0);fprintf('迭代次數(shù)為%dn',num); end 上機(jī)運(yùn)行>> A=4,-1,0,-1,0,0;-1,4,-1,0,-1,0;0,-1,4,-1,0,-1;-1,0,-1,4,-1,0;0,-1,0,-1,4,-1;0,0,-1,0,-1,4;>> b=0;5;-2;5;-2;6;>> SOR(A,b,1.334)結(jié)果為 1.00001878481009 1.99998698322858 1.00001815013068 2.00000041318053 0.999991858543476 2.

17、0000068413569迭代次數(shù)為13>> SOR(A,b,1.95)結(jié)果為 0.999984952088107 2.00000960832604 0.999959115182729 2.0000168426006 1.00000443526697 1.99997885113446迭代次數(shù)為241>> SOR(A,b,0.95)結(jié)果為 0.999961518309351 1.99995706825231 0.999963054838453 1.99996580572033 0.999971141727589 1.9999827300678迭代次數(shù)為17由以上輸出得到w

18、取值不同的情況下，得到的滿足精度要求的結(jié)果，迭代次數(shù)分別如輸出所示第五次上機(jī)作業(yè)8.從函數(shù)表x0.00.10.20.30.4010.5f(x)0.398940.396950.391420.381380.368120.35206出發(fā)，用下列方法計(jì)算f(0.15),f(0.31)及f(0.47)的近似值（1） 分段線性插值（2） 分段二次插值（3） 全區(qū)間上拉格朗日插值解：（1）線性插值編寫函數(shù)如下function R = div_line( x0,y0,x )%線性插值p=length(x0);n=length(y0);m=length(x);if(p=n)%x的個(gè)數(shù)與y的個(gè)數(shù)不等 error(

19、'數(shù)據(jù)輸入有誤，請重新輸入'); return;else fprintf('線性插值n'); for t=1:m z=x(t); if(z<x0(1)|z>x0(p) fprintf('x%d不在所給自變量范圍內(nèi)，無法進(jìn)行插值',t); continue; else for i=1:p-1 if(z<x0(i+1) break; end; end; R(t)=y0(i)*(x(t)-x0(i+1)/(x0(i)-x0(i+1)+y0(i+1)*(x(t)-x0(i)/(x0(i+1)-x0(i); end; end; end;e

20、nd 上機(jī)運(yùn)行如下>> x0=0.0 0.1 0.195 0.3 0.401 0.5;>> y0=0.39894 0.39695 0.39142 0.38138 0.36812 0.35206;>> x=0.15 0.31 0.47;>> div_line(x0,y0,x)線性插值ans = 0.39404 0.38007 0.35693即結(jié)果為f(0.15)0.39404,f(0.31) 0.38007,f(0.47) 0.35693（2）分段二次插值編寫的函數(shù)如下function R = div2line(x0,y0,x)%分段二次插值p=l

21、ength(x0);m=length(y0);n=length(x);if(p=m) error('輸入錯(cuò)誤，請重新輸入數(shù)據(jù)');end;for t=1:n if(x(t)<x0(1)|x(t)>x0(p) fprintf('x%d不在所給區(qū)間上',t); continue; end; tag=2; m=abs(x(t)-x0(1)+abs(x(t)-x0(2)+abs(x(t)-x0(3); for i=3:p-1 sum=abs(x(t)-x0(i-1)+abs(x(t)-x0(i)+abs(x(t)-x0(i+1); if(sum<m)

22、m=sum; tag=i; end end; fprintf('tag=%dn',tag);R(t)=y0(tag-1)*(x(t)-x0(tag)*(x(t)-x0(tag+1)/(x0(tag-1)-x0(tag)*(x0(tag-1)-x0(tag+1)+y0(tag)*(x(t)-x0(tag-1)*(x(t)-x0(tag+1)/(x0(tag)-x0(tag-1)*(x0(tag)-x0(tag+1)+y0(tag+1)*(x(t)-x0(tag-1)*(x(t)-x0(tag)/(x0(tag+1)-x0(tag-1)*(x0(tag+1)-x0(tag);End

23、上機(jī)運(yùn)行，執(zhí)行結(jié)果為>> x0=0.0 0.1 0.195 0.3 0.401 0.5;>> y0=0.39894 0.39695 0.39142 0.38138 0.36812 0.35206;>> x=0.15 0.31 0.47;>>div2line(x0,y0,x)ans = 0.39448 0.38022 0.35725即分段二次插值方法下，f(0.15)0.39448,f(0.31) 0.38022,f(0.47) 0.35725（3）上機(jī)編寫的程序如下function R = lagrange(x0,y0,x)%全區(qū)間上拉格朗日插值

24、p=length(y0);n=length(x0);m=length(x); %計(jì)算函數(shù)表和x的長度if p = n error('數(shù)據(jù)輸入有誤，請重新輸入'); %若函數(shù)表的x與y長度不一致則輸入有誤 else fprintf('拉格朗日插值n'); for t=1:m %利用循環(huán)計(jì)算每個(gè)x的插值 s=0.0; z=x(t); for k=1:n p=1; for i=1:n if i=k p=p*(z-x0(i)/(x0(k)-x0(i); end end s=s+y0(k)*p; end %根據(jù)拉格朗日插值公式求解yR(t)=s; %輸出插值結(jié)果 ende

25、nd上機(jī)運(yùn)行結(jié)果為>> x0=0.0 0.1 0.195 0.3 0.401 0.5;>> y0=0.39894 0.39695 0.39142 0.38138 0.36812 0.35206;>> x=0.15 0.31 0.47;>> lagrange(x0,y0,x)拉格朗日插值ans = 0.39447 0.38022 0.35722即分段二次插值方法下，f(0.15)0.39447,f(0.31) 0.38022,f(0.47) 0.357229.解：上機(jī)程序如下，為方便起見，將所有操作分在四個(gè)函數(shù)中進(jìn)行入口函數(shù)function =sp

26、line( X,Y,xx,y1_0,y1_18 )%輸出自變量所對應(yīng)的函數(shù)值M=getM(X,Y,y1_0,y1_18);%先得到Ms=xx;k=length(xx);for a=1:k s(xx(a)=getExp(X,Y,M,xx(a);%計(jì)算自變量所在小區(qū)間對應(yīng)曲線的表達(dá)式，并根據(jù)表達(dá)式計(jì)算對應(yīng)的函數(shù)值 fprintf('s(%d)=%fn',xx(a),s(xx(a); %輸出打印函數(shù)值endend獲取Mfunction M = getM(X,Y,y1_0,y1_1)%得到Mn=length(X);a=0*X;b=a;c=a;h=a;f=a;b=b+2;h(2:n)=X

27、(2:n)-X(1:n-1); % h(1)不用 a(2:n-1)=h(2:n-1)./(h(2:n-1)+h(3:n);c(2:n-1)=1-a(2:n-1);a(n)=1;c(1)=1; Yf(2:n)=Y(2:n)-Y(1:n-1); f(2:n-1)=6*(Yf(3:n)./h(3:n)-Yf(2:n-1)./h(2:n-1)./(h(2:n-1)+h(3:n);f(1)=6*(Yf(2)/h(2)-y1_0)/h(2);f(n)=6*(y1_1-Yf(n)/h(n)/h(n);M=CalM(n,a,b,c,f);%計(jì)算Mend計(jì)算Mfunction f = CalM(n,a,b,c,

28、f)% 追趕法求解Meps=0.1e-15; %防止參數(shù)過小，是的計(jì)算誤差過大if abs(b(1)<eps disp('除數(shù)為0，停止計(jì)算'); returnelse c(1)=c(1)/b(1);end%追趕法：根據(jù)遞推算式計(jì)算for i=2:n-1 b(i)=b(i)-a(i)*c(i-1); if abs(b(i)<eps disp('除數(shù)為0，停止計(jì)算'); return else c(i)=c(i)/b(i); endendb(n)=b(n)-a(n)*c(n-1);%追趕法:根據(jù)遞推算式計(jì)算f(1)=f(1)/b(1);for i=2:

29、n f(i)=(f(i)-a(i)*f(i-1)/b(i);end%以下求解Ux=y, x的值存入ffor i=n-1:-1:1 f(i)=f(i)-c(i)*f(i+1);endreturnend得到自變量所在區(qū)間的表達(dá)式，并求自變量對應(yīng)的函數(shù)值function y = getExp(X,Y,M,x)%根據(jù)X、Y、M計(jì)算表達(dá)式，并根據(jù)表達(dá)式計(jì)算對應(yīng)的函數(shù)值n=length(X);h(2:n)=X(2:n)-X(1:n-1); %判斷x落在哪個(gè)小區(qū)間n1=1;n2=n;while n2=n1+1 n5=fix(n1+n2)/2); if x>X(n5) n1=n5; else n2=n5

30、; endend%計(jì)算yy=M(n1)*(X(n2)-x)3/(6*h(n2)+ M(n2)*(x-X(n1)3/(6*h(n2);y=y+(Y(n1)-M(n1)*h(n2)*h(n2)/6)*(X(n2)-x)/h(n2);y=y+(Y(n2)-M(n2)*h(n2)*h(n2)/6)*(x-X(n1)/h(n2);%結(jié)束end上機(jī)試運(yùn)行，過程如下>> X=0.52 3.1 8.0 17.95 28.65 39.62 50.65 78 104.6 156.6 208.6 260.7 312.5 364.4 416.3 468 494 507 520;>> Y=5.2

31、8794 9.4 13.84 20.2 24.9 28.44 31.1 35 36.5 36.6 34.6 31.0 26.34 20.9 14.8 7.8 3.7 1.5 0.2;>> xx=2 4 6 12 16 30 60 110 180 280 400 515;>> y1_0=1.86548;>> y1_18=-0.046115;spline(X,Y,xx,y1_0,y1_18)s(2)=7.825123s(4)=10.481311s(6)=12.363477s(12)=16.575574s(16)=19.091594s(30)=25.386402s

32、(60)=32.804283s(110)=36.647878s(180)=35.917148s(280)=29.368462s(400)=16.799784s(515)=0.542713根據(jù)上述程序運(yùn)行結(jié)果，可得所有自變量對應(yīng)的函數(shù)值，如上輸出結(jié)果所示第六次上機(jī)作業(yè)10.已知一組實(shí)驗(yàn)數(shù)據(jù)i123456789xi1345678910yi1054211234試用最小二乘法求它的多項(xiàng)式擬合曲線，并求出最低點(diǎn)位置。解：試用matlab繪圖命令，將以上個(gè)點(diǎn)繪制在坐標(biāo)圖上，如下圖所示該函數(shù)圖像形如二次函數(shù)的圖像，現(xiàn)使用二次擬合程序如下function xmin,ymin = SecondFitting(x

33、,y)%最小二乘法，二次擬合,形如a+bx+cx2x=x'y=y'm=length(x);A(:,1)=ones(m,1);A(:,2)=x;A(:,3)=x.2;cc=Ay;a=cc(1);b=cc(2);c=cc(3);fprintf('擬合的曲線為%f%fx+%fx2n',a,b,c);xx=1:0.01:10;yy=cc(1)+cc(2)*xx+cc(3)*xx.2;plot(x,y,'r*',xx,yy);ymin=min(yy);c1=a-ymin;p=c b c1;xmin=roots(p);fprintf('最低點(diǎn)的橫縱坐

34、標(biāo)分別為');end上機(jī)運(yùn)行>> x=1 3 4 5 6 7 8 9 10;>> y=10 5 4 2 1 1 2 3 4;>> SecondFitting(x,y)擬合的曲線為13.459664-3.605309x+0.267571x2最低點(diǎn)的橫縱坐標(biāo)分別為ans = 6.74 6.7342圖像如下圖所示根據(jù)程序運(yùn)行結(jié)果，得到擬合方程為y=13.459664-3.605309x+0.267571x2,最低點(diǎn)坐標(biāo)為（6.74，6.7342）。計(jì)算方法第七次上機(jī)作業(yè)15.用龍貝格算法計(jì)算橢圓x2400+y2100 =1的周長，使誤差不超過10-4.解：

35、橢圓周長L計(jì)算公式為L=402300(sin)2+100 d上機(jī)程序如下所示（注：把每次計(jì)算結(jié)果存儲于一個(gè)二維數(shù)組之中，為了輸出如書中所示的表格形式，每一次結(jié)果的下標(biāo)有一定的調(diào)整）function = Romberg( fun,a,b,e )%fun表示被積函數(shù)，a為區(qū)間下限，b為上限，e為精度要求%龍貝格算法求橢圓周長T=zeros();T(1,1)=(b-a)/2*(feval(fun,a)+feval(fun,b);T(2,1)=0.5*T(1,1)+(b-a)/2*(feval(fun,a+(b-a)/2);T(2,2)=(4*T(2,1)-T(1,1)/3; T(3,1)=0.5*T

36、(2,1)+(b-a)/4*(sum(feval(fun,a+(2*(1:2)*(b-a)/4);T(3,2)=(4*T(3,1)-T(2,1)/3;T(3,3)=(16*T(3,2)-T(2,2)/15;T(4,1)=0.5*T(3,1)+(b-a)/8*(sum(feval(fun,a+(2*(1:4)*(b-a)/8);T(4,2)=(4*T(4,1)-T(3,1)/3;T(4,3)=(16*T(4,2)-T(3,2)/15;T(4,4)=(64*T(4,3)-T(3,3)/63;T(5,1)=0.5*T(4,1)+(b-a)/16*(sum(feval(fun,a+(2*(1:8)*(b-a)/16);T(5,2)=(4*T(5,1)-T(4,1)/3;T(5,3)=(16*T(5,2)-T(4,2)/15;T(5,4)=(64*T(5,3)-T(4,3)/63;%以上為求出序列表中前五行的值k=5;while(4*abs(T(k,4)-T(k-1,4)>e) k=k+1; T(k,1)=0.5*T(k-1,1)+(b-a)/(2(k-1)*(sum(fev