西安交大計算方法B2017大作業(yè)

1、實用標(biāo)準(zhǔn)計算方法B上機(jī)報告姓名：學(xué)號：班級：學(xué)院：任課教師：2017年12月29日題目一:1.1 題目內(nèi)容某通信公司在一次施工中，需要在水面寬度為 20米的河溝底部沿直線走向 鋪設(shè)一條溝底光纜。在鋪設(shè)光纜之前需要對溝底的地形進(jìn)行初步探測，從而估計 所需光纜的長度，為工程預(yù)算提供依據(jù)。已探測到一組等分點位置的深度數(shù)據(jù)（單 位：米）如下表所示：分點0123456深度9.018.967.967.978.029.0510.13分點78910111213深度11.1812.2613.2813.3212.6111.2910.22分點14151617181920深度9.157.907.958.869.811

2、0.8010.93（1）請用合適的曲線擬合所測數(shù)據(jù)點；（2）估算所需光纜長度的近似值，并作出鋪設(shè)河底光纜的曲線圖;1.2 實現(xiàn)題目的思想及算法依據(jù)首先在題目（1）中要實現(xiàn)的是數(shù)據(jù)的擬合，顯然用到的是我們在第三章中數(shù) 據(jù)近似的知識內(nèi)容。多項式插值時，這里有21個數(shù)據(jù)點，則是一個20次的多項 式，但是多項式插值隨著數(shù)據(jù)點的增多， 會導(dǎo)致誤差也會隨之增大，插值結(jié)果會 出現(xiàn)龍格現(xiàn)象，所以不適用于該題目中點數(shù)較多的情況。 為了避免結(jié)果出現(xiàn)大的 誤差，同時又希望盡可能多地使用所提供的數(shù)據(jù)點，提高數(shù)據(jù)點的有效使用率， 這里選擇分段插值方法進(jìn)行數(shù)據(jù)擬合。分段插值又可分為分段線性插值、分段二 次插值和三次樣條插

3、值。由于題目中所求光纜的現(xiàn)實意義，而前兩者在節(jié)點處的 光滑性較差，因此在這里選擇使用三次樣條插值。根據(jù)課本SPLINEMI法和TSS算法，采用第三種真正的自然邊界條件，在選 定邊界條件和選定插值點等距分布后，可以先將數(shù)據(jù)點的二階差商求出并賦值給 右端向量d,再根據(jù)TSS解法求解三對角線線性方程組從而解得 M值。求出M后， 對區(qū)間進(jìn)行加密，計算200個點以便于繪圖以及光纜長度計算。對于問題（2），使用以下的公式：20L= 0 f（x）ds19k =0二:f(x) 1 f'(x)2dxk 1;2-k f(x) d f (x)2dx k1.3算法結(jié)構(gòu)1. For i = 0,1,2,； n1

4、.1 yi = Mi2. For k = 1,22.1 For i = n,n -1, ,k2.1.1 (Mj -M j 二)/化 一乂鋰)二 M i3. % -Xq= h14. For i = 1,2, ,n -14.1 x羊X = h*4.2 hv/(hi +hg)二 Cj；1Ci=即2二 b4.3 6Mj j = dj5. do= Mo；dn= Mn； 0= Co2= boLn= an；2= bn6. D = RdL 17. For k=2,3，m ! 獲取M的矩陣元素個數(shù)，存入7.1 ak/k7.2 0 - Ik Ck" k7.3 dk - lk ' kd : k8.

5、 m / "m = M m9. For k 二 m -1, m -2,19.1 ( k 7 Mk J/% 二 Mk10. 1= k!獲取x的元素個數(shù)存入11. Fori =12，s-111.1 if x <xi then i= k； breakelse i 1- k文案大全12. xk -xkJ=- h； x -x= x；x-xkJ= ?3xM k j M k6$(g-M.2.2h _hk i )x，(yk-M k )刃 / h =661.4 matlab源程序n=20;x=O:n;y=9.01 8.96 7.96 7.97 8.02 9.05 10.13 11.18 12.2

6、6 13.28 13.32 12.61 11.29 10.229.15 7.90 7.95 8.86 9.81 10.80 10.93;M=y;h=zeros(1, n+1); c=zeros(1, n+1); d=zeros(1, n+1); a=zeros(1, n+1); b=2*o nes(1, n+1); h(2)=x(2)-x(1);%用于存放差商，此時為零階差商for i=2: n%書本110頁算法SPLINEMh(i+1)=x(i+1)-x(i);c(i)=h(i+1)/(h(i)+h(i+1);a(i)=1-c(i);enda(n +1)=-2;%十算邊界條件c(0),a(

7、n+1),采用的是第三類邊界條件c(1)=-2;for k=1:3 燦算k階差商for i=n+1:-1:k+1M(i)=(M(i)-M(i-1)/(x(i)-x(i-k);endif(k=2)d(2: n)=6*M(3: n+1);end%計算2階差商%給d賦值if(k=3) d(1)=(-12)*h(2)*M(4);d( n+1)=12*h( n+1)*M( n+1);%計算邊界條件d(0),d(n),采用的是第三類邊界條件endend l=zeros(1, n+1); r=zeros(1, n+1); u=zeros(1, n+1); q=zeros(1, n+1);u(1)=b(1);

8、 r(1)=c(1);q(1)=d(1);for k=2:n+1%利用書本49頁算法TSS求解三對角線性方程組r(k)=c(k);l(k)=a(k)/u(k-1);u(k)=b(k)-l(k)*r(k-1);q(k)=d(k)-l(k)*q(k-1);endp( n+1)=q( n+1)/u( n+1);for k=n :-1:1p(k)=(q(k)-r(k)*p(k+1)/u(k);endfprintf('三對角線性方程組的解為：');disp(p);%求擬合曲線x1=0:0.1:20;%首先對區(qū)間進(jìn)行加密，增加插值點n1=10* n;x2=zeros(1, n1+1);x3

9、=zeros(1, n1+1);s=zeros(1, n1+1);for i=1: n1+1for j=1: nif x1(i)>=x(j)&&x1(i)<=x(j+1)用書本 111 頁算法 EVASPLINE求解擬合曲線 s(x)h(j+1)=x(j+1)-x(j);x2(i)=x(j+1)-x1(i); x3(i)=x1(i)-x(j);s(i)=(p(j).*(x2(i).A3/6+p(j+1).*(x3(i).A3/6+(y(j)-p(j).*(h(j+1).A2/6).*x2(i)+.(y(j+1)-p(j+1).*(h(j+1).A2/6).*x3(i

10、)/h(j+1);endendendplot(x,-y,'x')%畫出插值點hold onplot(x1,-s)%畫出三次樣條插值擬合曲線hold ontitle（'三次樣條插值法擬合電纜曲線'）;xlabel('河流寬度 /m'); ylabel('河流深度 /m'); Len gth=0;for i=1: n1L=sqrt(x1(i+1)-x1(i)A2+(s(i+1)-s(i)A2);%計算電纜長度Len gth=Le ngth+L;endfprintf(' 電纜長度(m)='); disp(Le ngth)

11、;1.5結(jié)果與說明鋪設(shè)海底光纜的曲線如圖1.1所示-:次樣條插伉法血合沖底光纜曲WQ 246810121416 1B 20河流寬圉m圖1. 1三次樣條插值法擬合海底光纜曲線由上圖可以看出，所得到的曲線光滑，能夠較好得反映實際的河溝底部地勢 形貌。電纜長度計算結(jié)果為26.6656m (圖1.2 )。0.0410-0.76091.92250+ 8708-0.24=70+ 3&20-0. &224-1.8224電纜()=26. 6656圖1.2海底光纜長度結(jié)果題目二2.1題目內(nèi)容：已知非線性方程 丄；cos( X Sin旳dv -0在2,3中有根。請設(shè)計算法，求出該根,7T并使求出的

12、根的誤差不超過10*。2.2實現(xiàn)題目的思想及算法依據(jù)對于該題的非線性方程，可以將其分解成兩個部分：(1)求解數(shù)值積分；(2)求解非線性方程。首先求解數(shù)值積分，令 g(“ =cos(xsin"，則利用最簡單的梯形公式可以得到n h-f(x)(gU) gU)，其中h二一,耳二ih。于是就有了 f (x) =0形式的非i壬2兀n線性方程，這里選擇二分法進(jìn)行求解。算法參考課本BISECTION法。2.3算法結(jié)構(gòu)1. f (x(° )戶 f0f (x(1)戶 f2. iff0 f1 0 then stop3. iff0£% then輸出xC)作為根；stop4. if

13、63;v % then輸出x°)作為根；stop5. (xC +xC)= x26. If xOI-xvx。)then 輸出 x 作為根；stop7. f x = f8. if f £名2 then輸出x作為根；9. iff1 f ：0 then9.1 x= x0 ; f - f0Else9.2 x= x1 ; f f110. go to 52.4 Matlab源程序*function y = theta( i )y=i*h;end*function y = g( x,theta )%為了計算關(guān)于theta的數(shù)值積分，先令 g=cos(x*s in (theta) y=cos

14、(x*s in( theta);end*fun ctio n f = hsz(x) % n=10000;h=pi/n; % 將區(qū)間分成 f=0;for i=1: n計算數(shù)值積分n份f=f+h/(2*pi)*(g(x,i+1)+g(x,i); end end*error=1e-4; a=2; b=3; f0=hsz(a); f1=hsz(b); if f0*f1>0%誤差允許值%初始區(qū)間%判斷方程是否有解disp('該方程在a,b上無解');elseif f0=0x=a;elseif f1=0 x=b;else%判斷方程解是否在區(qū)間兩邊界上%二分法求解方程得解aO=a;b0

15、=b;while abs(b0-a0)/2)>=error half=(a0+b0)/2;fa=hsz(aO);fb=hsz(b0); fhalf=hsz(half);if fhalf=0%計算中點處的函數(shù)值，用以判斷解的位置x=half; break;elseif fa*fhalf<0b0=half; elsea0=half;endend%定義新區(qū)間，為原區(qū)間的一半%定義新區(qū)間x=(b0+a0)/2;%方程組的解endfprintf(' 方程組的解為：')disp(x)2.5結(jié)果與說明由于是利用梯形公式來求解，需要將區(qū)間劃分為n個區(qū)間，而n的取值是很關(guān)鍵 的，需要

16、取得適當(dāng)?shù)膎值才能滿足誤差精度。這里將區(qū)間分成10000份，可以得到結(jié)果圖2.1，x=2.4048。» ti2方程的網(wǎng)泊：2. 404E圖2. 1 n=10000 方程的解由于變量的嵌套很多及自身水平的不足，所以我將很多變量通過調(diào)用子程序的方 式實現(xiàn)，分別為 theta.m : g.m; hsz.m 。題目二3.1題目內(nèi)容：編寫程序?qū)崿F(xiàn)大規(guī)模方程組的高斯消去法程序，并對所附的方程組進(jìn)行求解。 所附方程組的類型為對角占優(yōu)的帶狀方程組。數(shù)據(jù)說明：dat20171.dat為非壓縮帶狀方程組，階數(shù)為10階，該方程組供調(diào)試程序使用，該方程組的根都為1;(2)dat20172.dat為壓縮帶狀方

17、程組，階數(shù)為20階，該方程組供調(diào)試程序使用,該方程組的根都為1;dat20173.dat為非壓縮帶狀方程組，階數(shù)在 3000階左右；dat20174.dat為壓縮帶狀方程組，階數(shù)在 50000階左右； 3.2實現(xiàn)題目的思想及算法依據(jù)高斯消去法主要分為消去和回代過程， 消去過程形成上三角矩陣，回代過程 就是自下而上實現(xiàn)方程組求解。題目中系數(shù)矩陣是嚴(yán)格對角占優(yōu)矩陣， 所以無需 進(jìn)行列主元，極大的提高了程序的效率。對于n階、上帶寬q、下帶寬p的非壓縮格式帶狀矩陣，通過比較 k+p、k+q及n 的大小，對于行消去，選擇k+p和n中的最小值作為循環(huán)終點。在用第k行進(jìn)行 消去時，只需要計算第k列的ak+1

18、,k到ak+p，k，每一行只需計算ak，k+1到ak，k+q ； 從第n-c+1到n行，第k列計算ak+1, k到知k，每一行只需計算ak，k+1到ak，n，于 是減小運算量。壓縮格式忽略了零元素，存儲方式與非壓縮格式有所不同，所以需要建立兩種格 式中元素之間的對應(yīng)關(guān)系，即壓縮 Ak，p+1 =非壓縮人,k，根據(jù)此關(guān)系進(jìn)行程序編 寫。3.3算法結(jié)構(gòu)1. A的階數(shù)=n2. Fori 二 k 1,k 2,na.2.1-=- ik ik:kk2.2 For j = k 1,k 2,，n2.2.1Ot ij 0( ikG kj 二 Of ij2.3:i -ik ： k =l-'i'n

19、/ _：nn = xn3. Fork =n -1,n -2,1n-kj xj ) / -：%k = xk j 士 13.4 Matlab 源程序fname='dat20171.dat： fprintf('目標(biāo)文件為：');disp(fname);fid=fope n( fname,'r'); fhead=fread(fid,5,'l on g'); 息 bbh=dec2hex(fhead(2);%首先計算非壓縮帶狀方程組%讀取文件，可改為 dat20173.dat%二進(jìn)制模式讀取數(shù)據(jù)%以long int形式讀取文件前5個數(shù)據(jù)，獲取系數(shù)矩陣

20、信%將版本號轉(zhuǎn)換為16進(jìn)制%判斷所讀取的文件是否為目標(biāo)文件if strcmp(bbh,'102')%比較版本號，相同則說明為非壓縮帶狀矩陣%讀取矩陣階數(shù)n%讀取上帶寬q%讀取下帶寬p%非壓縮格式，需要讀取nA2個浮點數(shù)，按行方式順序存儲%讀取右端系數(shù)disp（'目標(biāo)文件系數(shù)矩陣為非壓縮帶狀矩陣，讀取正確'） elsedisp(' 讀取錯誤')endn=fhead (3);q=fhead(4);p=fhead(5);fprintf('系數(shù)矩陣階數(shù)為:');disp( n)fprintf('上帶寬為:');disp(q

21、)fprintf('下帶寬為:');disp(q)d=fread(fid, nA2,'float'); b=fread(fid ,n ,'float');%生成系數(shù)矩陣A=zeros (n,n);k=1;for i=1:nfor j=1: nA(i,j)=d(k);%因為系數(shù)矩陣按行方式順序存儲在d中，所以依次形成A(i,j )k=k+1;endendfclose(fid);%讀取結(jié)束n=len gth(b);%高斯消去法，書本 33頁到34頁算法GAUSSPP%肖去過程for k=1: n-1fori=k+1: min( k+p ,n)nA(i

22、,k)=A(i,k)/A(k,k);for j=k+1: min( k+q ,n)A(i,j)=A(i,j)-A(i,k)*A(k,j);endb(i)=b(i)-A(i,k)*b(k);endend%回代過程x( n)=b( n)/A( n,n);for k=n-1:-1:1sum=0;forj=k+1:n sum=sum+A(k,j)*x(j);endx(k)=(b(k)-sum)/A(k,k);endm=10;fprintf('方程的前 d個根:n',m);for j=1:mfprintf(1,'%0.5f ',x(j);end%各方程出的解輸出為txt文

23、件if strcmp(fname,'dat20173.dat') file=strcat('f:',fname,'和文件名%僅對不為零的區(qū)域作高斯消去，邊界時k+p有可能大于%用第k行消去%只需處理右邊k+q個數(shù)據(jù)，邊界時k+q可能大于n%更新各列%更新右端系數(shù)b列%輸出前m個根%小數(shù)點后保留5位小數(shù)矩陣方程組的解','.txt');%生成要創(chuàng)建文件的文件路徑fid=fope n(file,'a+');%以讀寫方式打開指定文件，將文件指針指向文件末尾fprintf(fid,'%s的所有解 n',fn

24、ame);fprin tf(fid,'%0.5f ',x);%向指定文件寫入數(shù)據(jù)，保留5位小數(shù)fclose(fid);endfprin tf('n');%第二步處理壓縮帶狀方程組，即數(shù)據(jù)文件（2）和（4）, （2）供調(diào)試使用fname='dat20172.dat： fprintf('目標(biāo)文件為:');disp(fname);fid=fope n( fname,'r'); fhead=fread(fid,5,'l on g'); 息 bbh=dec2hex(fhead(2);%讀取文件，可改為 dat2017

25、4.dat%二進(jìn)制模式讀取數(shù)據(jù)%以long int形式讀取文件前5個數(shù)據(jù)，獲取系數(shù)矩陣信if strcmp(bbh,'2O2')%比較版本號，相同則說明文件數(shù)據(jù)是壓縮帶狀矩陣%將版本號轉(zhuǎn)換為16進(jìn)制 %判斷所讀取的文件是否為目標(biāo)文件disp（'目標(biāo)文件系數(shù)矩陣為壓縮帶狀矩陣，讀取正確'）elsedisp('讀取錯誤')endn=fhead (3);q=fhead(4);p=fhead(5);fprintf('系數(shù)矩陣階數(shù)為:');disp( n)fprintf('上帶寬為:');disp(q)fprintf(

26、9;下帶寬為:');disp(q)d=fread(fid, n*(p+q+1),'float'); b=fread(fid ,n ,'float');%生成系數(shù)矩陣A=zeros (n, p+q+1);k=1;for i=1: nfor j=1:p+q+1A(i,j)=d(k);(i,j )k=k+1;endendfclose(fid);%讀取矩陣階數(shù)n%讀取上帶寬q%讀取下帶寬p%壓縮格式，以n*（p+q+1）矩陣存儲，同樣按行存儲 %讀取右端系數(shù)%因為系數(shù)矩陣按行方式順序存儲在d中，所以依次形成A%讀取結(jié)束n=len gth(b);%高斯消去法，書本

27、 33頁到34頁算法GAUSSPP%肖去過程for k=1: n-1for i=1:min(p,n-k)%用第 k 行消去A(k+i,p+1-i)=A(k+i,p+1-i)/A(k,p+1);%壓縮格式 p+1 列即為原來第 k 列for j=1:min(q,n-k)%只需處理右邊q個數(shù)據(jù)，但邊界時q可能大于(n-k )實用標(biāo)準(zhǔn)A(k+i,p+1-i+j)=A(k+i,p+1-i+j)-A(k+i,p+1-i)*A(k,p+1+j);end%更新各列b(k+i)=b(k+i)-b(k)*A(k+i,p+1-i);endend徊代過程%更新右端系數(shù)b列x( n)=b( n)/A( n,p+1);

28、for i=n-1:-1:1sum=0;for j=1:mi n(q, n-i)sum=sum+x(i+j)*A(i,p+1+j);endx(i)=(b(i)-sum)/A(i,p+1);endm=20;fprintf('方程的前 d個根:n',m);%俞出前m個根for j=1:mfprin tf(1,'%0.5f ',x(j);end%小數(shù)點后保留5位小數(shù)%各方程出的解輸出為txt文件if strcmp(fname,'dat20174.dat')file=strcat('f:',fname,'和文件名矩陣方程組的解&#

29、39;,'.txt');%生成要創(chuàng)建文件的文件路徑fid=fope n( file,'a+');%以讀寫方式打開指定文件，將文件指針指向文件末尾fprintf(fid,'%s的所有解 n',fname);fprin tf(fid,'%0.5f ',x);%向指定文件寫入數(shù)據(jù)，保留5位小數(shù)fclose(fid);end3.5結(jié)果與說明首先使用dat20171和dat20172驗證程序的正確性，如圖 3.1目標(biāo)文件為dat20171.dat目標(biāo)文件系數(shù)拒陸為非壓縮帶狀矩陣，讀取正確秦數(shù)矩陣前敎為：10下帯竟為：3方程的前艮1. 000

30、001.00000 L 000001.000001.00000 L 00000 k 000001. 000001.000001+ 00000目標(biāo)文件為：dat2Q172.dat目標(biāo)文件系數(shù)拒肚芮壓縮帶狀柜陣，讀取正確菜數(shù)葩陣盼數(shù)為：20下帶寬為：3方程的前20個視1. 000001,00000 1.000001.00000 1.000001.00000 U 000001,000001.000001. 000001,00000實用標(biāo)準(zhǔn)文案大全圖3. 1 驗證dat20171和dat20172的正確性由上圖可知dat20171.dat中為上帶寬為3,下帶寬為3,階數(shù)為10的非壓縮矩陣，解都為1 d

31、at20172.dat中為上帶寬為3,下帶寬為3,階數(shù)為20的壓縮矩陣，解都為1 與題目中所給的條件相符。接下來求解dat20172和dat20174數(shù)據(jù)，如圖3.2目標(biāo)文件為.dat20173. dat目標(biāo)文件系藪拒陣擁非壓貓帶找拒陣，讀取正確系數(shù)拒陣齡埶為：2380上帶寬為：4下帝寬為.4方程的前10個根：1.01900 L 019001.019001.0190Q1.019001. 019001.019001.01900 L 019QQ1.0X900目標(biāo)文件為;dat 20174, dat目標(biāo)文件系數(shù)拒陣為壓縮帶狀矩陣，讀取正確系埶矩陸階數(shù)為：51296上帚甕為：5下帶竟為：5方程的前20

32、個根：2.017002, 017002,01700 2.01700 2.017002.017002.01700 2.017002. 017002,01700圖3. 2 dat20173 和dat20174的方程組解除此之外，還將程序運行的全部結(jié)果輸出在txt文件中分別給出部分截圖如圖da吃的所有解 1.01900 L 01900 1.01900 L0L 01900 L 019001.01900 1. 01900L 01900 1.019001. 01900 1. 01900L 01900 1. 01900L 01900 1.019001.01900 L 01900L 01900 L 01900

33、1.01900 1.01900L 01900 1. 019001.01900 1.01900L 01900 1.01900L01900 L 019001.01900 L 01900L 01900 1,019001. 01900 1. 019001.01900 1.01900L 01900 1. 01900L 01900 L 01900L 01900 1. 01900L01900 L 019001. 01900 L 019001.01900 1. 019001. 01900 1. 019001.01900 L 01900L 01900 1. 01900L01900 L 01900L01900 1

34、.019001. 01900 L 019001. 01900 1.019001. 01900 1. 019001. 01900 L 019001, 0100 L 01900L01900 1.019001.01900 L 01900L 01900 1.01900L 01900 1. 019001. 01900 1. 01900L 01900 1. 01900 L 01900 L 01900 1.01600 1. 01900 1.01500 1. 01900 L 01900 1. 01900 1.01900 1. 01900 L 01900 1. 01900 1. 01600 1. 01900 L

35、 01900 1. 01900 1.01900 1. 01900 1.01900 1. 01900 1.01500 1. 01900 L 01900 L 01900 1.01900 1. 01900 1. 01900 1. 01900 1.01900 1. 01900 1.01900 1. 01900 1.01900 1, 01900 1.01900 L 01900hat2O174. dat的所有解 2 01700 2. 01700 2. 01700 2. C2. 01700 2.01700 2. 01700 2. 01700 2. 01700 2. 0170C2. 01700 2.01700

36、 2.017002. 01700 2. D1700 2. 017002. 01700 2. 01700 2.017002. 01700 2. 01700 2.017002. 01700 2.01700 2.017002. 01700 2.01700 2.017002. 01700 2.01700 2.017002. 01700 2.01700 2.017002. 01700 2.01700 2.017002, 01700 2.01700 2.017002. 01700 2.01700 2.017002. 01700 2.01700 2.017002. 01700 2.01700 2.01700

37、2. 01700 2.01700 2.017002. 01700 2.01700 2.017002. 01700 2.01700 2.017003. 01700 2.01700 2,017002. 01700 2,01700 2.017002.01700 2.01700 2. 0170C2. 01700 2.01700 2. 017DC2. 01700 2. 01700 2. 0170C2. 01700 2. 01700 2. 0170CL 01700 2. 01700 2. 0170C2. 01700 2. 01700 2. 0170C£.01700 2.01700 2. 0170C2. 01700 2. 01700 2. 0170C2.01700 2.01700 2. 0170C2. 01700 2. 01700 2. 0170C2. 01700 2.01700 2. 0170C2、01700 2. 01700 2. 0170C2.01700 2. 01700 2. 0170(2. 01700 2. 01700 2. 0170C2.0170Q 2. 01700 2. 0170C2. 01700 2, 01700 2. 0170C£ 01700 2. 0