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

1、第一次 滬a；b; f=Q(x)(x) +10; c=(a+b)/2: while abs(b-a)0.3i0 (-3) if 1(c) *f(b) fprintf C n x = of. f(x )=. of ti , x, f (s); x - 0. 0$0b8t(z)-0. 00067 Al 求得近似根為 0.09058 (2) 上機(jī)代碼如圖所示 x0=0 xO = 0 f=(x) (2-exp Ck) )/10 . riiiLa a.bU0-f (kO)0. 5*10* -3; sG=f (stO): x=sO： end fprintf ( Xq x5f,f (s)5f n 4x,f

2、tx): x= Ck 09064, f(x)=0. 09061 A 得近似根為0.09064; (3) 牛頓法上機(jī)代碼如下 (z) fixp(x)+10*x-2 : E=esp i x +10;冬 ill數(shù)f 的導(dǎo)函數(shù) y-e(i)x-f(x/g(x): hFQ; i=0 ; vhile abs(yCxO-xO)D.5*10(-3 sO=y (xO); 5=50 : end; fprintf C? n s = . 5W : x - 0, 09091 計(jì)算所得近似解為0.09091 第三次上機(jī)作業(yè) 上機(jī)作業(yè)181頁(yè)第四題 線性方程組為 rl.1348 3.83261.16513.40171 xl

3、- r 9.5342 053011.7S752.533015435 工2 6.3941 3.4129 4.93178.76431.3142 x3 18.4231 1.23714.999810.67210.0147 x4 16.9237 I1 (1)順序消元法 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,10.6721,0.0147; b=9.5342;6.3941;18.4231;16.9237; 上機(jī)代碼(函數(shù)部分)如下 fun cti

4、 on b = gaus( A,b )%用 b返回方程組的解 B=A,b; n=len gth(b); RA=ra nk(A); RB=ra nk(B); dif=RB-RA; if dif0 disp(此方程組無解); return end if RA=RB if RA=n format long; disp(此方程組有唯一解); for p=1: n-1 for 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)

5、; 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ù)組解); end end 上機(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,10.6721,0.0147; b二9.5342;6.3941;18.4231;16.9237; X二gaus(A,b) 此方程組有唯一解

6、 1.0000 1.0000 1.0000 1.0000 (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.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=ra nk(A); RB=ra nk(B); n=len gth(b); dif=RB-RA; format

7、 Io ng; if dif0 disp(該方程組無解); return end if dif=0 if RA=n disp(該方程組有唯一解); c=zeros(1, n); for i=1: n-1 max=abs(A(i,i); m=i; for j=i+1: n if max A二1.1348,3.8326,1.1651,3.4017;0.5301,1.7875,2.53 30,1.5435;3.4129,4.9317,8.7643,1.3142;1.2371,4.99 98,10.6721,0.0147; b=9.5342;6.3941;18.4231;16.9237; X=lieZ

8、hu(A,b) 該方程組有唯一解 X = 1.0000 1.0002 0.999999999999999 0.999999999999999 根據(jù)兩種方法運(yùn)算結(jié)果，順序消元法得到結(jié)果比列主元消元法得到的 結(jié)果精度更高。 (注：matlab使用的是2015b版本，不知道是保留小數(shù)位數(shù)太少，還 是程序原因，順序消元輸出結(jié)果總是等于準(zhǔn)確解，請(qǐng)老師指正) 第四次上機(jī)作業(yè) 7.分析用下列迭代法解線性方程組 4 -1 0 -1 0 -1 4 -1 0 -1 0 -1 4 -1 0 -1 0 -1 4 -1 0 -1 0 -1 4 0 0 -1 0 -1 -1 4 0 -1 0 X3 的收斂性，并求出使o.

9、oool的近似解及相應(yīng) 的迭代次數(shù)。 (1) 雅可比迭代法 解：上機(jī)編寫的函數(shù)如下 function = Jacobi(X,b) %雅可比迭代法 D=diag(diag(X);%得到對(duì)角線元素組成的矩陣 B=in v(D)*(D-X); f=in v(D)*b; b(：,：)=0; x仁 B*b+f; num=1; while (n orm(x1-b)0.0001)%判斷當(dāng)前的解是否達(dá)到精度要求 b=x1; x1=B*b+f; num=nu m+1; end ; fprintf(求得的解為：n); disp(xl); fprintf(迭代次數(shù)：(次n ,num); end 上機(jī)運(yùn)行，結(jié)果如下

10、A=4,-1,O,-1,O,0；-1,4,-1,O,-1,O;O,-1,4,-1,O,-1；-1,O,-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 1.99995018125674 0.999975090628368 1.99996353014876 迭代次數(shù)：28次 滿足要求的解如輸出結(jié)果所示，總共迭代了28次 (2) 高斯-賽德爾迭代法 上機(jī)程序如下所示 fun ctio n =

11、Gauss_Seidel( A,b ) %高斯賽德爾迭代法 D=diag(diag(A); L=D-tril(A); U=D-triu(A); B=i nv(D-L)*U; f=i nv(D-L)*b; b(：,：)=0； xO=B*b+f; num=1; while (norm(x0-b)0.0001) num=nu m+1; b=x0; x0=B*b+f; end ; fprintf(結(jié)果為 n); disp(x0); fprintf(迭代次數(shù)為：(次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,-

12、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 1.99998400886989 15次 迭代次數(shù)為：15次 滿足精度要求的解如上述程序打印輸出所示，迭代了 (3)SOF迭代法(w依次取 1.334,1.95,0.95) 上機(jī)代碼如下 fun ctio n = SOR(A,b,w ) %S0!迭代法 D=diag(diag

13、(A); L=D-tril(A); U=D-triu(A); B=i nv(D-w*L)*(1-w)*D+w*U); f=w*i nv(D-w*L)*b; b(:,:)=0; x0=B*b+f; num=1; while (norm(x0-b)0.0001) num=nu m+1; b=x0; x0=B*b+f; end ; fprintf(結(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

14、;0,0,-1,0,-1,4; b=0;5;-2;5;-2;6; SOR(A,b,1.334) 結(jié)果為 1.009 1.99998698322858 1.068 2.053 0.999991858543476 2.69 迭代次數(shù)為13 SOR(A,b,1.95) 結(jié)果為 0.999984952088107 2.604 0.999959115182729 2.06 1.697 1.99997885113446 迭代次數(shù)為241 SOR(A,b,0.95) 結(jié)果為 0.999961518309351 1.99995706825231 0.999963054838453 1.999965805720

15、33 0.999971141727589 1.9999827300678 迭代次數(shù)為17 由以上輸出得到w取值不同的情況下，得到的滿足精度要求的結(jié)果, 迭代次數(shù)分別如輸出所示 第五次上機(jī)作業(yè) 8.從函數(shù)表 x 0.0 0.1 0.2 0.3 0.401 0.5 f(x) 0.39894 0.39695 0.39142 0.38138 0.36812 0.35206 出發(fā)，用下列方法計(jì)算f(0.15),f(0.31) 及f(0.47)的近似值 (1) 分段線性插值 (2) 分段二次插值 (3) 全區(qū)間上拉格朗日插值 解：(1)線性插值 編寫函數(shù)如下 fun ctio n R = div_li n

16、e( x0,y0,x) %線性插值 p=le ngth(xO); n=len gth(y0); m=le ngth(x); if(p=n)%x的個(gè)數(shù)與y的個(gè)數(shù)不等 error(數(shù)據(jù)輸入有誤，請(qǐng)重新輸入); return; else fprintf(線性插值 n); for t=1:m z=x(t); if(zx0(p) fpri ntf(x%d不在所給自變量范圍內(nèi)，無法進(jìn)行插值,t); con ti nue; else for i=1:p-1 if(z x0=0.0 0.1 0.195 0.3 0.401 0.5; y0=0.39894 0.39695 0.39142 0.38138 0.36

17、812 0.35206; x=0.15 0.31 0.47; div_line(x0,y0,x) 線性插值 ans = 0.394040.380070.35693 即結(jié)果為 f(0.15)0.39404,f(0.31)0.38007,f(0.47) 乂 0.35693 (2)分段二次插值 編寫的函數(shù)如下 fun ctio n R = div2li ne(xO,yO,x) %分段二次插值 p=le ngth(xO); m=le ngth(yO); n=len gth(x); if(p=m) error(輸入錯(cuò)誤，請(qǐng)重新輸入數(shù)據(jù) ); end; for t=1: n if(x(t)x0(p) fp

18、rintf(x%d不在所給區(qū)間上,t); con ti nue; en d; tag=2; m=abs(x(t)-x0(1)+abs(x(t)-x0(2)+abs(x(t)-x 0(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 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; div2li ne(x0,y0,x) ans = 0.394

19、48 0.38022 0.35725 即分段二次插值方法下， 0.35725 f(0.15) | = 0.39448,f(0.31)= 0.38022,f(0.47) (3) 上機(jī)編寫的程序如下 fun cti on R = lagra nge(x0,y0,x) %全區(qū)間上拉格朗日插值 p=le ngth(yO); n=le ngth(xO);m=le ngth(x); %計(jì)算函數(shù)表和x的長(zhǎng)度 if p = n error(數(shù)據(jù)輸入有誤，請(qǐng)重新輸入)； %若函數(shù)表的x與y長(zhǎng)度不一致則輸入有誤 else fprintf(拉格朗日插值n); for t=1:m %利用循環(huán)計(jì)算每個(gè)x的插值 s=0.

20、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ù)拉格朗日插值公式求解y R(t)=s; %輸岀插值結(jié)果 end end 上機(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; lagra nge(x0,y0,x) 拉格朗日插值 ans = 0.394470.380220.35722

21、即分段二次插值方法下， f(0.15)0.39447,f(0.31)0.38022,f(0.47) 怎 0.35722 9.解：上機(jī)程序如下，為方便起見，將所有操作分在四個(gè)函數(shù)中進(jìn)行 入口函數(shù) function =spline( X,Y,xx,y1_0,y1_18 ) %輸出自變量所對(duì)應(yīng)的函數(shù)值 M=getM(X,Y,y1_0,y1_18);%得到 M s=xx; k=le ngth(xx); for a=1:k s(xx(a)=getExp(X,Y,M,xx(a)i%算自變量所在小區(qū)間對(duì)應(yīng)曲線的表達(dá)式，并 根據(jù)表達(dá)式計(jì)算對(duì)應(yīng)的函數(shù)值 fprin tf(s(%d)=%fn,xx(a),s(xx

22、(a); %輸出打印函數(shù)值 end end 獲取M fun ctio n M = getM(X,Y,y1_0,y1_1) %得到M n=le ngth(X); a=0*X;b=a;c=a;h=a;f=a; b=b+2; h(2: n) =X(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:

23、 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ì)算 M end 計(jì)算M fun ctio n f = CalM( n,a,b,c,f) %追趕法求解M eps=0.1e-15; %防止參數(shù)過小，是的計(jì)算誤差過大 if abs(b(1)eps disp(除數(shù)為0,停止計(jì)算); return else c(1)=c(1)/b(1); end %追趕法：根據(jù)遞推算式計(jì)算B for i=2: n-1 b(i)=b(i)-a(i)

24、*c(i-1); if abs(b(i)X(n 5) n1=n5; else n2=n5; end end % %計(jì)算y y=M( n1)*(X( n2)-xF3/(6*h( n2)+ M( n2)*(x-X( n1)F3/(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

25、 156.6 208.6 260.7 312.5 364.4 416.3 468 494 507 520; Y=5.28794 9.413.84 20.2 24.9 28.44 31.135 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;spl in e(X, Y,xx,y1_0,y1_18) s(2)=7.825123 s(4)=10.481311 s(6)=12.363477 s(12)=

26、16.575574 s(16)=19.091594 s(30)=25.386402 s(60)=32.804283 s(110)=36.647878 s(180)=35.917148 s(280)=29.368462 s(400)=16.799784 s(515)=0.542713 根據(jù)上述程序運(yùn)行結(jié)果，可得所有自變量對(duì)應(yīng)的函數(shù)值，如上輸出結(jié) 果所示 第六次上機(jī)作業(yè) 10.已知一組實(shí)驗(yàn)數(shù)據(jù) i i 2 3 4 C 6 7 8 9 J X i 3 4 5 6 7 8 C 1 0 y -i 5 4 2 1 1 2 3 4 0 試用最小二乘法求它的多項(xiàng)式擬合曲線，并求出最低點(diǎn)位置 解：試用matla

27、b繪圖命令，將以上個(gè)點(diǎn)繪制在坐標(biāo)圖上，如下圖所 10f 10 該函數(shù)圖像形如二次函數(shù)的圖像，現(xiàn)使用二次擬合 程序如下 function xmin,ymin = SecondFitting(x,y) %最小二乘法，二次擬合，形如a+bx+cxA2 x=x; y=y; m=le ngth(x); A(:,1)=on es(m,1); A(:,2)=x; A(:,3)=x.a2; cc=Ay; a=cc(1); b=cc(2); c=cc(3); fprintf(擬合的曲線為 %f%fx+%fxA2n,a,b,c); xx=1:0.01:10; yy二cc(1)+cc(2)*xx+cc(3)*xx.

28、A2; plot(x,y,r*,xx,yy); ymi n=min(yy); c1=a-y min; p=c b c1; xmi n=roots(p); fprintf(最低點(diǎn)的橫縱坐標(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; Secon dFitt in g(x,y) 擬合的曲線為 13.459664-3.605309x+0.267571xA2 最低點(diǎn)的橫縱坐標(biāo)分別為ans = 6.74 6.7342 圖像如下圖所示 11 1234567B910 根據(jù)程序運(yùn)行結(jié)果，得到擬合方程為 y=13.459664-3.6053

29、09x+0.267571x9,最低點(diǎn)坐標(biāo)為(6.74 , 6.7342)。 計(jì)算方法第七次上機(jī)作業(yè) 2 2 龍 y 15.用龍貝格算法計(jì)算橢圓400+100=1的周長(zhǎng)，使誤差不超過 解：橢圓周長(zhǎng)L計(jì)算公式為 L=4 上機(jī)程序如下所示 (注：把每次計(jì)算結(jié)果存儲(chǔ)于一個(gè)二維數(shù)組之中，為了輸出如書中所 示的表格形式，每一次結(jié)果的下標(biāo)有一定的調(diào)整) %fun表示被積函數(shù)，a為區(qū)間下限，b為上限, fun cti on = Romberg( fun, a,b,e ) e為精度要求 %龍貝格算法求橢圓周長(zhǎng) T=zeros(); T(1,1)=(b-a)/2*(feval(fu n,a)+feval(fu n

30、,b); T(2,1)=0.5*T(1,1)+(b-a)/2*(feval(fu n, a+(b-a)/2); T(2,2)=(4*T(2,1)-T(1,1)/3; T(3,1)=0.5*T(2,1)+(b-a)/4*(sum(feval(fu n,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(fu n,a+(2*(1:4)*(b-a)/8); T(4,2)=(4*T(4,1)-T(3,1)/3; T(4,3)=(

31、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(fu n,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)/(2A(k-1)*(sum(f

32、eval(fu n,a+(2*(1:(2A(k-2) )*(b-a)/0(k-1); T(k,2)=(4*T(k,1)-T(k-1,1)/3; T(k,3)=(16*T(k,2)-T(k-1,2)/15; T(k,4)=(64*T(k,3)-T(k-1,3)/63; end disp(4*T);%區(qū)間為四分之一圓周，所以乘以4輸岀為橢圓周長(zhǎng) 上機(jī)運(yùn)行，結(jié)果如下(復(fù)制數(shù)字結(jié)果，排版容易亂，截圖如下) 運(yùn)行命令 fun 二(x)power(300*(si n(x).A2)+100),0.5); Romberg(fu n, 0,pi/2,0.0001) f L3i=E,Cx pciFer i , 3M，x- ;s in (x)2J+10C : , Cl 61 ; RcnbcrE( fun，0, p i/2, (L (0 01) 9匚2右廠9別：旳詐 0 0 0 S-.制飾149“陽(yáng)眸 0 0 W W0IL27E03 1D7.3IB906M4793 107.