西安交大計(jì)算方法B上機(jī)作業(yè)

1、實(shí)用標(biāo)準(zhǔn)計(jì)算方法（B ）上機(jī)作業(yè)一、三次樣條擬合某通信公司在一次施工中，需要在水面寬度為 20米的河溝底部沿直線走向 鋪設(shè)一條溝底光纜。在鋪設(shè)光纜之前需要對(duì)溝底的地形進(jìn)行初步探測(cè)，從而估計(jì) 所需光纜的長度，為工程預(yù)算提供依據(jù)。已探測(cè)到一組等分點(diǎn)位置的深度數(shù)據(jù)（單 位:米）如下表所示：分點(diǎn)0123456深度9.018.967.967.978.029.0510.13分點(diǎn)78910111213深度11.1812.2613.2813.3212.6111.2910.22分點(diǎn)14151617181920深度9.157.907.958.869.8110.8010.93（1）請(qǐng)用合適的曲線擬合所測(cè)數(shù)據(jù)點(diǎn)；估算

2、所需光纜長度的近似值，并作出鋪設(shè)河底光纜的曲線圖;解:1、算法實(shí)現(xiàn)的思想及依據(jù)題目（1）為曲線擬合問題多項(xiàng)式插值、分段插值和最小二乘法。多項(xiàng)式插 值，隨著插值數(shù)據(jù)點(diǎn)的數(shù)目增多，誤差也會(huì)隨之增大，因此不選用。最小二乘法 適于數(shù)據(jù)點(diǎn)較多的場合，在此也不適用。故選用分段插值。分段插值又分為分段線性插值、分段二次插值、三次樣條插值及更高階的多 項(xiàng)式插值。由本題的物理背景知，光纜正常工作時(shí)各點(diǎn)應(yīng)該是平滑過渡， 因此至 少選用三次樣條插值法。對(duì)于更高階的多項(xiàng)式插值，由于“龍格現(xiàn)象”而不選用。題目（2）求光纜長度，即求擬合曲線在0到20的長度，對(duì)弧長進(jìn)行積分19即可。光纜長度的第一型線積分表達(dá)式為If (x

3、)2dx。2、算法實(shí)現(xiàn)的結(jié)構(gòu)參考教材給出的SPLINEM算法和TTS算法，在選定邊界條件和選定插值點(diǎn)等 距分布后，可以先將數(shù)據(jù)點(diǎn)的二階差商求出來并賦值給右端向量_d，再根據(jù)TSS 解法求解M。光纜長度的第一型線積分表達(dá)式為I ： 1 ； 1一f(x)2dx。k 03、程序運(yùn)行結(jié)果及分析三戲轉(zhuǎn)辛也謠mm三撿律壬由嘆m冒專一料即齊辛広屯血*制揚(yáng)氣K WEY -T沁* -1. r -T r *ffJ* ;:匕 1*J-u/.7r*/- . /-世筈曲#. gji:V .1J4?A1L1101抑.總生盤KS_!1 1E la F圖1.1三種邊界條件下三次樣條插值趙鋼亍團(tuán)口你想便用第幾抽邊界糸件丫請(qǐng)喻人

4、2、2之一： 1LODCI夭第1種迪界條件T悵度處.密陰米圖1.2光纜長度4、MATLAB 代碼:1)自己編程實(shí)現(xiàn)代碼clear;clc;I=input(你想使用第幾種邊界條件？請(qǐng)輸入 1、2、3之一：);x=0:20;y=9.01 8.96 7.96 7.97 8.02 9.05 10.13 11.18 12.26 13.28 13.32 12.61 11.29 10.229.15 7.90 7.95 8.86 9.81 10.8 10.93;plot(x,-y,k.,markersize,15)%y為深度，取負(fù)號(hào)hold on%計(jì)算一階差商y1=o nes(1,21);for i=2:1:

5、21y1(i)=(y(i)-y(i-1)/(x(i)-x(i-1);end%計(jì)算二階差商y2=o nes(1,21);for i=3:1:21y2(i)=(y1(i)-y1(i-1)/(x(i)-x(i-2);end%計(jì)算三階差商y3=o nes(1,21);for i=4:1:21y3(i)=(y2(i)-y2(i-1)/(x(i)-x(i-3);end%選擇邊界條件(I)個(gè)點(diǎn)的二階差商為 0if I=1d(1)=0;d(21)=0;a(21)=0;c(1)=0;%第一個(gè)點(diǎn)和最后endif I=2d(1)=6*y1(1);d(21)=-6*y1(21);a(1)=1;c(1)=1;endif

6、 I=3d(1)=-12*y3(1);d(21)=-12*y3(21);a(21)=-2;c(1)=-2;%endfor i=2:20d(i)=6*y2(i+1);end%構(gòu)造帶狀矩陣求解(追趕法)b=2*o nes(1,21);a=0.5*o nes(1,21);%a(21)=-2;c=0.5*o nes(1,21);%c(1)=-2;u(1)=b(1);r(1)=c(1);%追yz(1)=d(1);for i=2:21l(i)=a(i)/u(i-1);u(i)=b(i)-l(i)*r(i-1);r(i)=c(i);yz(i)=d(i)-l(i)*yz(i-1);end%趕xg(21)=yz

7、(21)/u(21);for i=20:-1:1xg(i)=(yz(i)-r(i)*xg(i+1)/u(i);endM=xg;%所有點(diǎn)的二階導(dǎo)數(shù)值%求函數(shù)表達(dá)式并積分t=1;h=1;N=1000x1=0:20/(N-1):20;len gth=0;for i=1:Nfor j=2:20if x1(i)=x(j)t=j；break;elset=j+1；endendf1=x(t)-x1(i);f2=x1(i)-x(t-1);S(i)=(M(t-1)*f1A3/6/h+M(t)*f2A3/6/h+(y(t-1)-M(t-1)*hA2/6)*f1+(y(t)-M(t)*hA2/6)*f2)/h;Sp(

8、i)=-M(t-1)*f1A2/2/h+M(t)*f2A2/2/h+(y(t)-y(t-1)/h-(M(t)-M(t-1)*h/6;len gth(i+1)=sqrt(1+Sp(i)A2)*(20/(N-1)+le ngth(i);%第一類線積分endfigure(1);plot(x1,-S,r-)%深度曲線griddisp(第,num2str(I),種邊界條件下長度,num2str(length(N+1),米)axis fill;xlabel(測(cè)點(diǎn) / 米);ylabel(深度/ 米);title(三次樣條曲線擬合);legend（數(shù)據(jù)點(diǎn),擬合曲線,3）;二、最小二乘近似假定某天的氣溫變化記

9、錄如下表所示，試用數(shù)據(jù)擬合的方法找出這一天的氣 溫變化的規(guī)律；試計(jì)算這一天的平均氣溫，并試估計(jì)誤差。時(shí)刻01234567891011121111111122222平均氣溫5444456800358111111122222時(shí)刻345678901234333222222111平均氣溫141975420876解：1、算法實(shí)現(xiàn)的思想及依據(jù)本題共有25個(gè)數(shù)據(jù)點(diǎn)，數(shù)據(jù)量較大，因此選用多項(xiàng)式插值會(huì)導(dǎo)致較大的誤 差，插值樣條函數(shù)雖然有較好的數(shù)值性質(zhì)，但是 Mi的計(jì)算量不小，而且沒有統(tǒng) 一的公式來表示。另外，本題要求找出溫度變化規(guī)律，插值方法要求p（x）必須通 已知數(shù)據(jù)點(diǎn)，只會(huì)導(dǎo)致擬合曲線失去本有的數(shù)據(jù)所表示的

10、規(guī)律；從表中的數(shù)據(jù)點(diǎn) 可以看出，溫度并不準(zhǔn)確（溫度只測(cè)量到整數(shù)位），因此沒有必要要求擬合曲線 通過數(shù)據(jù)點(diǎn)。因此，選取最小二乘近似。2、算法實(shí)現(xiàn)的結(jié)構(gòu)算法參考課本LSS算法。利用已有數(shù)據(jù)來生成了 G，將y作為其最后一列 （第n+1列）存放。先形成矩陣Qk，再變換Gk-1到Gk;求解三角方程Rx=h1 ； 最后根據(jù)定義計(jì)算誤差。平均溫度采用數(shù)據(jù)點(diǎn)相加求和求平均，避免數(shù)值積分的 繁瑣。文案大全3、程序運(yùn)行結(jié)果及分析克淫悄 劃曲吉磯平均溫寒:為：1.200000直輸入多頂式近似列數(shù)最鬲項(xiàng)的洪執(zhí)I： 2 各頃批倉扇頻為：CtW!0.3Ofi31次項(xiàng)系甑2. 5064G幻頁系甑譏093794決基犬小対W.

11、 ,433* +飆會(huì)也1lli09!特2J1151*甲均*為：21. 2000(10諸怕九鑒肝垃返井超埶昴左氏的衣埶h： 3 各項(xiàng)擦臺(tái)系數(shù)洵；電次頂鬲霖1盂388I次項(xiàng)事熬7. 22733項(xiàng)20745:災(zāi)啖聚數(shù)-0. OOS2SJ1逞差K小為11.44S2卅二第筑戦合的丸昭痙比蜩 平均喘摩T擁：21.200000語輸人步頃式近似超計(jì)晟高項(xiàng)的戾數(shù)十4 備項(xiàng)斟合乖塾商；腰頂系踰：屛丈頃樂數(shù)T.T帖2丈貞釆埶0一國怕：M次頃系熱7.坊岀技4盲両喬豹n. 00093485試越大小藥7.48384、MATLAB 代碼%最小二乘法擬合氣溫變化規(guī)律 clear;clc;x=0: 24;y=15,14,14,

12、14,14,15,16,18,20,20,23,25,28,31,34,31,29,27,25,24,22,20,18,17,16;N=le ngth(y);sum=O;%求一天的平均溫度 averagefor i=1:Nsum=sum+y(i);endaverage=sum/N;fprintf(平均溫度 T 為：%fn,average);h=input(請(qǐng)輸入多項(xiàng)式近似函數(shù)最高項(xiàng)的次數(shù)h :);n=h+1;%參課本LSS算法G=zeros(N,n+1);%建立一個(gè) N行，n+1列的矩陣 G%生成矩陣G各個(gè)元素for i=1:n%矩陣G中前N列元素for j=1:NG(j,i)=x(1,j)F

13、-1);endendfor i=1:1:N%將y作為矩陣 G中第(N+1)列元素G(i, n+1)=y(i);end%形成矩陣Qka=zeros(1,n);%建立存放的矩陣ab=zeros(1,N); % 建立存放B的矩陣bfor k=1: nfor i=k:Na(k)=G(i,k)A2+a(k);enda(k)=-sig n(G(k,k)*(sqrt(a(k);w(k)=G(k,k)-a(k);%建立存放3的矩陣for j=k+1:1:Nw(j)=G( j,k);endb(k)=a(k)*w(k);%變換Gk-1到GkG(k,k)=a(k);for j=k+1:1: n+1sum_wg=0;

14、for i=k:Nsum_wg=w(i)*G(i,j)+sum_wg;endt=sum_wg/b(k);for i=k:1:NG(i,j)=t*w(i)+G(i,j);endendend%解三角方程Rx=h1X( n)=G( n,n+1)/G( n,n);for i=( n-1):-1:1for j=i+1: n sum_gx=G(i,j)*X(j)+sum_gx;endX(i)=(G(i, n+1)-sum_gx)/G(i,i);sum_gx=0;end%參數(shù)X(i)的回代p=zeros(1,N);%建立一維最小二乘近似數(shù)組p(x)for j=1:Nfor i=1: nP(j)=P(j)+X

15、(i)*x(j)A(i-1);endend%顯示擬合函數(shù)各階系數(shù)disp(各項(xiàng)擬合系數(shù)為：);for i=1: ndisp(num2str(i-1),次項(xiàng)系數(shù),num2str(X(i);end%誤差E2求解E2=0;%多項(xiàng)式計(jì)算誤差for i=n+1:NE2=G(i, n+1)A2+E2;endE2=sqrt(E2);disp(誤差大小為);disp（E2）;%繪制曲線，顯示結(jié)果plot（x,y,r*,x,p,k-,Li neWidth,1.5）;xlabel（時(shí)刻 /h）;ylabel（平均氣溫 / C）;title（最小二乘法擬合的氣溫變化曲線圖）xlim（0 24）;legend（數(shù)據(jù)點(diǎn)

16、,擬合曲線）三、線性方程組求解。（1）編寫程序?qū)崿F(xiàn)大規(guī)模方程組的高斯消去法程序，并對(duì)所附的方程組進(jìn)行求解。所附方程組的類型為對(duì)角占優(yōu)的帶狀方程組。（2）針對(duì)本專業(yè)中所碰到的實(shí)際問題，提煉一個(gè)使用方程組進(jìn)行求解的例子， 并對(duì)求解過程進(jìn)行分析、求解。解：1、算法實(shí)現(xiàn)的思想及依據(jù)和算法實(shí)現(xiàn)的結(jié)構(gòu)高斯消去法主要分為兩大步驟，即消去過程和回帶過程。算法借鑒課本GAUSSPP算法。由于題目中給出的系數(shù)矩陣是對(duì)角占優(yōu)的矩陣，因此可以不選 主元直接進(jìn)行高斯消去；另外在非壓縮格式帶狀矩陣中，存在著大量的非零元素， 非零元素的運(yùn)算毫無意義并且占用大量機(jī)器資源和時(shí)間，因此對(duì)課本中給出的GAUSSPP算法進(jìn)行優(yōu)化。對(duì)

17、于n階、上帶寬q、下帶寬p的帶狀矩陣，選取p與q較大者（賦值給c）， 在第1行到第n-c行，第k列，只需要計(jì)算ak 1,k到ak c,k，每一行也只需從ai,k 1 到ai,k q（i從k+1到k+q）；在第n-c+1到n行，第k列，計(jì)算ak 1,k到an,k，每實(shí)用標(biāo)準(zhǔn)一行只需要從ai,k 1計(jì)算到ai,n (i從k+1到n)，這樣可以減小運(yùn)算量。對(duì)于壓縮帶狀矩陣，其消去過程和非壓縮帶狀矩陣基本一致，不同之處在于：壓縮格式忽略了零元素，因此需要建立壓縮格式帶狀矩陣與非壓縮格式帶狀矩陣 的對(duì)應(yīng)關(guān)系。主元素對(duì)應(yīng)關(guān)系：B(k,p+1)=A(k,k)( B是壓縮格式帶狀矩陣，A是非壓縮格式帶狀矩陣)

18、，編寫程序時(shí)需要根據(jù)此關(guān)系建立其元素間的對(duì)應(yīng)關(guān)系。2、程序運(yùn)行結(jié)果對(duì)比及分析命令行窗口dat 1 dat腐非壓縮拒陣方程的前3個(gè)根-1.00000 1. 00000 1. 00000 L 00000 L 00000 稗庫運(yùn)行耗時(shí)/S：0. C058圖3.1求解dat61文件結(jié)果命令鈕口dat62. dat力壓蕭矩陸方程的前亍個(gè)根：1.000001. OODOO 1. OOOOC 1.00000 1.0000QSyiste* time: 1B-Dsc-2O10 10:17:35程時(shí) 0.0058573 sI圖3.2求解dat62文件結(jié)果dat63r dat左程的前5；卜糧：1. 00000 1.

19、 00000 1.010C0 1.00000 t. toooo 程序運(yùn)行料時(shí)滄：245. 1247圖3.3求解dat63文件結(jié)果I文案大全實(shí)用標(biāo)準(zhǔn)文案大全命令再國口dat64. dat為壓縮拒陣方程的前3個(gè)根；2.01600 2.01600 2.01600 2.01600 2. 01600System tiae.18-Dtc-2016 10:19 25 程序運(yùn)行耗時(shí)鼻4674 SA圖3.4求解dat64文件結(jié)果3、Matlab 代碼%改進(jìn)前的程序clc;cleardata_fname=dat64.dat：%文件名file_id = fope n( data_fname, rb);%Id = f

20、read(file_id, 1, ui nt32);%以二進(jìn)制格式打開源文件數(shù)據(jù)文件標(biāo)示id=dec2hex(ld);ver = fread(file_id, 1, in t32);%版本號(hào)n = fread(file_id,1, in t32);%方程組的階數(shù)p= fread(file_id,1,i nt32);%帶狀矩陣上帶寬q =fread(file_id,1,i nt32);%帶狀矩陣下帶寬%確定數(shù)據(jù)格式if strcmp(dec2hex(ver),1O2)%比較字符串，確定數(shù)據(jù)格式str= data_fname 為非壓縮矩陣;個(gè)浮點(diǎn)數(shù)，以一維格式存p的系數(shù)矩陣A中disp(str);

21、d = fread(file_id, nRoat); %非壓縮格式，需要讀取n2儲(chǔ)到中間變量db = fread(file_id,n,float); %再讀取右端向量的n個(gè)元素%將讀取到的數(shù)據(jù)放到階數(shù)為n，上帶寬為q，下帶寬為k = 1;for i = 1:nforj = 1:nA(i,j) = d(k);k = k+1;endendfclose(file_id);%消去過程for k=1: n-1for i=k+1: nA(i,k)=A(i,k)/A(k,k);for j=k+1: n;A(i,j)=A(i,j)-A(i,k)*A(k,j);endb(i)=b(i)-A(i,k)*b(k);

22、endend%回帶過程求解方程組的根x( n)=b( n)/A( n,n);for k=n-1:-1:1S=b(k);for i=k+1: nS=S-A(k,i)*x(i);endx(k)=S/A(k,k);endend%若為壓縮矩陣if strcmp(dec2hex(ver),2O2)str= data_fname 為壓縮矩陣;disp(str);d = fread(file_id, n*(p+q+1),float);% 壓縮格式一共要讀取n*(p+q+1)個(gè)數(shù)b = fread(file_id,n,float);%再讀取右端向量的n個(gè)元素%將讀取到的數(shù)據(jù)放到n行、p+q+1列的系數(shù)矩陣中k

23、=1;for i = 1:nfor j = 1:(q+p+1)A(i,j) = d(k);k = k+1;endendfclose(file_id);%消去過程for k=1: n-qfor i=1:pA(k+i,p+1-i)=A(k+i,p+1-i)/A(k,p+1);for j=1:qA(k+i,p+1-i+j)=A(k+i,p+1-i+j)-A(k+i,p+1-i)*A(k,p+1+j);endb(k+i)=b(k+i)-b(k)*A(k+i,p+1-i);endendfor k = n-q+1: n-1for i = 1:n-kA(k+i,p+1-i)=A(k+i,p+1-i)/A(k,p+1);forj = 1: n-kA(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)-A(k+i,p+1-i)*b(