西安交通大學(xué)計(jì)算方法上機(jī)作業(yè).doc

計(jì)算方法上機(jī)作業(yè) 1.對(duì)以下和式計(jì)算：，要求：(1)若只需保留11個(gè)有效數(shù)字，該如何進(jìn)行計(jì)算；(2)若要保留30個(gè)有效數(shù)字，則又將如何進(jìn)行計(jì)算；（1）解題思想和算法實(shí)現(xiàn)：根據(jù)保留有效位數(shù)的要求，可以由公式得出計(jì)算精度要求。只需要很少內(nèi)存，時(shí)間復(fù)雜度和d呈線性，不需要高浮點(diǎn)支持。先根據(jù)while語(yǔ)句求出符合精度要求的n值的大小，然后利用for語(yǔ)句對(duì)這n項(xiàng)進(jìn)行求和，輸出計(jì)算結(jié)果及n值大小即可。（2）matlab源程序：保留11位有效數(shù)字時(shí)；clearclcformat longn=0;sum=1/(16n)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6);while sum=5*10(-11);n=n+1;sum=1/(16n)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6);endfor i=0:n-1; sum=sum+1/(16i)*(4/(8*i+1)-2/(8*i+4)-1/(8*i+5)-1/(8*i+6);endvpa(sum,11)n保留30位有效數(shù)字時(shí)；clearclcformat longn=0;sum=1/(16n)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6);while sum=5*10(-30);n=n+1;sum=1/(16n)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6);endfor i=0:n-1; sum=sum+1/(16i)*(4/(8*i+1)-2/(8*i+4)-1/(8*i+5)-1/(8*i+6);endvpa(sum,30)n（3）實(shí)驗(yàn)結(jié)果分析 圖1.1 保留11位有效數(shù)字的n值及計(jì)算結(jié)果圖 圖1.2 保留30位有效數(shù)字的n值及計(jì)算結(jié)果圖由計(jì)算結(jié)果可知，通過(guò)合理的誤差控制，分別通過(guò)7次和22次循環(huán)，可以實(shí)現(xiàn)題目所要求的精確度。2.某通信公司在一次施工中，需要在水面寬度為20米的河溝底部沿直線走向鋪設(shè)一條溝底光纜。在鋪設(shè)光纜之前需要對(duì)溝底的地形進(jìn)行初步探測(cè)，從而估計(jì)所需光纜的長(zhǎng)度，為工程預(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)預(yù)測(cè)所需光纜長(zhǎng)度的近似值，并作出鋪設(shè)河底光纜的曲線圖；（1）解題思想和算法原理給定區(qū)間a, b一個(gè)分劃 ：a=x0x1xN=b 若函數(shù)S(x)滿足下列條件：1) S(x)在每個(gè)區(qū)間xi, xj上是不高于3次的多項(xiàng)式。2) S(x)及其2階導(dǎo)數(shù)在a, b上連續(xù)。則稱(chēng)S(x)使關(guān)于分劃的三次樣條函數(shù)。（2）matlab源程序：clc,clearx=0:1: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.22 9.15 7.9 7.95 8.86 9.81 10.80 10.93;n=length(x);l(1)=0;m(n)=0;h=diff(x);df=diff(y)./diff(x);d(1)=0;d(n)=0; for j=2:n-1 l(j)=h(j)/(h(j-1)+h(j); m(j)=h(j-1)/(h(j-1)+h(j); d(j)=6*(df(j)-df(j-1)/(h(j-1)+h(j);endm=m(2:end);u=diag(m,-1);r=diag(l,1);a=diag(2*ones(1,n);A=u+r+a; M=inv(A)*d; syms gfor j=1:n-1s(j)=M(j)*(x(j+1)-g)3/(6*h(j)+M(j+1)*(g-x(j)3/(6*h(j)+(y(j)-M(j)*h(j)2/6)*(x(j+1)-g)/h(j)+(y(j+1)-M(j+1)*h(j)2/6)*(g-x(j)/h(j);endsr=0;for j=1:n-1 df=diff(s(j),g); warning off all; q=int(sqrt(1+df.2),g,j-1,j); r=r+q;endL=vpa(r,8);disp(the length of the label is L=);disp(L);for j=1:n-1 S(j,:)=sym2poly(s(j);end for j=1:n-1 x1=x(j):0.1:x(j+1); y1=polyval(S(j,:),x1); if j=1 y2=y1; else for i=1:11 k=(j-1)*10+i; y2(k)=y1(i); end endendx2=x(1):0.1:x(n); plot(x,y,o)gridhold onplot(x2,y2,r)（3）實(shí)驗(yàn)結(jié)果分析圖2.1 鋪設(shè)河底電纜長(zhǎng)度圖2.2 鋪設(shè)河底光纜的曲線圖由三次樣條插值得出的函數(shù)曲線的長(zhǎng)度和即鋪設(shè)河底電纜的長(zhǎng)度為26.498514。為了提高插值精度，用三次樣條插值可以增加插值節(jié)點(diǎn)的辦法來(lái)滿足要求，而且在給定節(jié)點(diǎn)數(shù)的條件下，三次樣條插值的精度要優(yōu)于多項(xiàng)式插值以及線性分段插值，雖然舍棄了降低誤差這個(gè)優(yōu)點(diǎn)，但是其曲線的光滑性要好一些。3.假定某天的氣溫變化記錄如下表所示，試用數(shù)據(jù)擬合的方法找出這一天的氣溫變化的規(guī)律;試計(jì)算這一天的平均氣溫,并試估計(jì)誤差。時(shí)刻0123456789101112平均氣溫15141414141516182020232528時(shí)刻131415161718192021222324平均氣溫313431292725242220181716（1）解題思想和數(shù)學(xué)原理：對(duì)于具體實(shí)驗(yàn)時(shí)，通常不是先給出函數(shù)的解析式，再進(jìn)行實(shí)驗(yàn)，而是通過(guò)實(shí)驗(yàn)的觀察和測(cè)量給出離散的一些點(diǎn)，再來(lái)求出具體的函數(shù)解析式。又因?yàn)闇y(cè)量誤差的存在，實(shí)際真實(shí)的解析式曲線并不一定通過(guò)測(cè)量給出的所有點(diǎn)。最小二乘法，形成法方程是求解這一問(wèn)題的很好的方法，本實(shí)驗(yàn)運(yùn)用這一方法實(shí)現(xiàn)對(duì)給定數(shù)據(jù)的擬合。對(duì)于給定的測(cè)量數(shù)據(jù)(xi,fi)(i=1,2,，n)，設(shè)函數(shù)分布為特別的，取為多項(xiàng)式 (j=0, 1,，m)則根據(jù)最小二乘法原理，可以構(gòu)造泛函令 (k=0, 1,，m)則可以得到法方程求該解方程組，則可以得到解，因此可得到數(shù)據(jù)的最小二乘解（2）matlab源程序：x=0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24; %給出題目數(shù)據(jù)（時(shí)間）y=15 14 14 14 14 15 16 18 20 20 23 25 28 31 32 31 29 27 25 24 22 20 18 17 16; %給出題目數(shù)據(jù)（溫度）plot(x,y, m*) %畫(huà)出各個(gè)離散數(shù)據(jù)點(diǎn)hold onfor n=2:4; %2、3、4代表擬合函數(shù)最高項(xiàng)次數(shù)alltemp=25; % alltemp代表數(shù)據(jù)點(diǎn)總共有25個(gè)A=zeros(n+1,n+1); %定義初始正規(guī)方程組的系數(shù)矩陣AC=ones(n+1,1); %定義初始正規(guī)方程組的系數(shù)向量CD=zeros(n+1,1); %定義初始正規(guī)方程組的向量Dfor i=1:n+1 for j=1:n+1 for k1=1: alltemp A(i,j)=A(i,j)+(x(k1).(i-1+j-1); end end for k2=1: alltemp D(i,1)=D(i,1)+(x(k2).(i-1).*(y(k2); endend %以上為計(jì)算出正規(guī)方程組矩陣A、D的所有元素的程序tol=1.0e-12;maxit=1000;C=bicg(A,D,tol,maxit); %使用bicg迭代算出正規(guī)方程組的系數(shù)向量Cp=0; %誤差分量E=0; %誤差總量if n=2 b=0:24; f=C(1)+C(2).*b+C(3).*(b.2); p=y(b+1)-f; for v=1:25 E=E+(p(v).2; end plot(b,f, r-)end %以上是對(duì)2階擬合函數(shù)的圖形處理與誤差計(jì)算if n=3 b=0:24; f=C(1)+C(2).*b+C(3).*(b.2)+C(4).*(b.3); p=y(b+1)-f; for v=1:25 E=E+(p(v).2; end plot(b,f, g-)end %以上是對(duì)3階擬合函數(shù)的圖形處理與誤差計(jì)算if n=4 b=0:24; f=C(1)+C(2).*b+C(3).*(b.2)+C(4).*(b.3)+C(5).*(b.4); p=y(b+1)-f; for v=1:25 E=E+(p(v).2; end plot(b,f, b-)end %以上是對(duì)4階擬合函數(shù)的圖形處理與誤差計(jì)算C,Eendn=2; %重新對(duì)n賦值，進(jìn)行指數(shù)函數(shù)擬合A=zeros(n+1,n+1); %重新對(duì)A矩陣賦初值C=zeros(n+1,1); %重新對(duì)C向量賦初值D=zeros(n+1,1); %重新對(duì)D向量賦初值for i=1:n+1 for j=1:n+1 for k=1: alltemp A(i,j)=A(i,j)+(x(k).(i-1+j-1); end end for l=1: alltemp D(i,1)=D(i,1)+(x(l).(i-1).*(log(y(l); endend %計(jì)算出A矩陣、D向量各元素?cái)?shù)值C=bicg(A,D,tol,maxit); %利用bicg迭代求解系數(shù)b=0:24;p=0;E=0;f=exp(C(1)+C(2).*b+C(3).*(b.2); p=y(b+1)-f; for v=1:25 E=E+(p(v).2; endplot(b,f, c-) %對(duì)指數(shù)擬合函數(shù)進(jìn)行圖形處理和誤差計(jì)算b=-C(3);c=C(2)/(2*b);a=exp(b*(c2)+C(1); %算出題設(shè)要求的指數(shù)擬合函數(shù)的各個(gè)系數(shù)a,b,c,Egrid onlegend(測(cè)量數(shù)據(jù),二次函數(shù),三次函數(shù),四次函數(shù),指數(shù)擬合,Location,NorthWest)hold off %hold on與hold off聯(lián)合使用，表示將各個(gè)曲線畫(huà)在同一個(gè)圖中圖3.1 二次曲線擬合系數(shù)與2范數(shù)誤差圖3.2 三次曲線擬合系數(shù)與2范數(shù)誤差圖3.3 四次曲線擬合系數(shù)與2范數(shù)誤差圖3.4 指數(shù)曲線擬合系數(shù)與2范數(shù)誤差圖3.5 數(shù)據(jù)原始點(diǎn)與擬合曲線對(duì)比圖（3） 實(shí)驗(yàn)結(jié)果分析：根據(jù)國(guó)家有關(guān)規(guī)定，用的是02、08、14、20時(shí)四個(gè)觀測(cè)時(shí)次的數(shù)據(jù)做平均，最有代表性。從圖中可以看出并不是多項(xiàng)式次數(shù)越高越好，隨著次數(shù)的增高，曲線所呈現(xiàn)出的給定點(diǎn)處和實(shí)際的吻合度越好，但對(duì)于其他地方的吻合度降低了。4.設(shè)計(jì)算法，求出非線性方程的所有根，并使誤差不超過(guò)。解：（1）解題思想和算法實(shí)現(xiàn)：對(duì)于一個(gè)非線性方程的數(shù)值解法很多。在此介紹兩種最常見(jiàn)的方法：二分法和Newton法。首先要研究函數(shù)的形態(tài)，確定根的數(shù)量和大致區(qū)間的位置。對(duì)于二分法，其數(shù)學(xué)實(shí)質(zhì)就是說(shuō)對(duì)于給定的待求解的方程f(x)，其在a,b上連續(xù)，f(a)f(b)0，且f(x)在a,b內(nèi)僅有一個(gè)實(shí)根x*，取區(qū)間中點(diǎn)c，若，則c恰為其根，否則根據(jù)f(a)f(c)e i=i+1; x0=x; x=x0-(6*x0.5-45*x0.2+20)/(30*x0.4-90*x0);%迭代enddisplay(方程的根)xdisplay(迭代的次數(shù))i（3） 實(shí)驗(yàn)結(jié)果分析：圖4.2 運(yùn)行結(jié)果對(duì)于Newton迭代法，三個(gè)初值x0都使得迭代收斂，這是非常重要的?？紤]Newton法迭代的收斂性條件：（1）存在一個(gè)區(qū)間，滿足。由曲線和所選的三個(gè)區(qū)間可知這一條件滿足。（2）f(x)是a,b上的單調(diào)函數(shù)，即對(duì)一切不變號(hào)。經(jīng)驗(yàn)證所選的三個(gè)區(qū)間滿足這一條件。（3）f(x)的凹向在a,b上不變，即在a,b上不改變符號(hào)。經(jīng)驗(yàn)證所選的三個(gè)區(qū)間滿足這一條件。（4）-111.5000 另外，可以看出Newton迭代法收斂速度也很快，且很快達(dá)到很高的精度，源于它一般是超線性收斂的。5. 編寫(xiě)程序?qū)崿F(xiàn)大規(guī)模方程組的列主元高斯消去法程序，并對(duì)所附的方程組進(jìn)行求解。針對(duì)本專(zhuān)業(yè)中所碰到的實(shí)際問(wèn)題，提煉一個(gè)使用方程組進(jìn)行求解的例子，并對(duì)求解過(guò)程進(jìn)行分析、求解。解:（1）算法原理由于一般線性方程在使用Gauss消去法求解時(shí)，從求解的過(guò)程中可以看到，若=0，則必須進(jìn)行行交換，才能使消去過(guò)程進(jìn)行下去。有的時(shí)候即使0，但是其絕對(duì)值非常小，由于機(jī)器舍入誤差的影響，消去過(guò)程也會(huì)出現(xiàn)不穩(wěn)定得現(xiàn)象，導(dǎo)致結(jié)果不正確。因此有必要進(jìn)行列主元技術(shù)，以最大可能的消除這種現(xiàn)象。這一技術(shù)要尋找行r，使得并將第r行和第k行的元素進(jìn)行交換，以使得當(dāng)前的的數(shù)值比0要大的多。這種列主元的消去法的主要步驟如下：1. 消元過(guò)程對(duì)k=1,2,n-1,進(jìn)行如下步驟。1) 選主元，記很小，這說(shuō)明方程的系數(shù)矩陣嚴(yán)重病態(tài)，給出警告，提示結(jié)果可能不對(duì)。2) 交換增廣陣A的r，k兩行的元素。 (j=k,n+1) 3）計(jì)算消元(i=k+1,n; j=k+1,n+1)2. 回代過(guò)程對(duì)k= n, n-1,1,進(jìn)行如下計(jì)算至此，完成了整個(gè)方程組的求解。（2） matlab源程序：%非壓縮,dat51.dat、dat53.datclear;clcfp=fopen(dat53.dat,rb);id=fread(fp,1,int32);ver=fread(fp,1,int32);N=fread(fp,1,int32);q=fread(fp,1,int32);p=fread(fp,1,int32);for i=1:N A(i,:)=fread(fp,N,float);endfor i=1:N d(i)=fread(fp,1,float);end%正向消去過(guò)程for i=1:N-q for k=1:p ll=A(i+k,i)/A(i,i); for j=i:i+q A(i+k,j)=A(i+k,j)-ll*A(i,j); end d(i+k)=d(i+k)-ll*d(i); endendfor i=N-q+1:N for k=1:N-i ll=A(i+k,i)/A(i,i); for j=i:N A(i+k,j)=A(i+k,j)-ll*A(i,j); end d(i+k)=d(i+k)-ll*d(i); endend%回代過(guò)程x(N)=d(N)/A(N,N);for i=N-1:-1:1 S=0; if i+qN cv=N;%cv-critical value else cv=i+q; end for j=i+1:cv S=S+A(i,j)*x(j); end x(i)=(d(i)-S)/A(i,i);endx%壓縮,dat52.dat、dat54.datclear;clcfp=fopen(dat54.dat,rb);id=fread(fp,1,int32);ver=fread(fp,1,int32);N=fread(fp,1,int32);q=fread(fp,1,int32);p=fread(fp,1,int32);for i=1:N A(i,:)=fread(fp,p+q+1,float);endfor i=1:N d(i)=fread(fp,1,float);end%正向消去過(guò)程for i=1:N if i+pN cv=p;%cv-critical value else cv=N-i; end for k=1:cv ll=A(i+k,p+1-k)/A(i,p+1); for j=p+1:p+q+1 A(i+k,j-k)=A(i+k,j-k)-ll*A(i,j); end d(i+k)=d(i+k)-ll*d(i); endend%回代過(guò)程x(N)=d(N)/A(N,p+1);for i=N-1:-1:1; S=0; ii=i+1; jj=p+2; while (ii=N)&(jj=p+q+1) S=S+A(i,jj)*x(ii); ii=ii+1; jj=jj+1; end x(i)=(d(i)-S)/A(i,p+1);end x（3） 實(shí)驗(yàn)結(jié)果：非壓縮矩陣求解結(jié)果（部分）壓縮矩陣求解結(jié)果（部分）（4） 分析心得：采用Gauss消去法時(shí)，如果在消元時(shí)對(duì)角線上的元素始終較大，那么本方法不需要進(jìn)行列主元計(jì)算，計(jì)算結(jié)果一般就可以達(dá)到要求，否則必須進(jìn)行列主元這一步，以減少機(jī)器誤差帶來(lái)的影響，使方法得出的結(jié)果正確。（5） 實(shí)例化學(xué)反應(yīng)方程式嚴(yán)格遵守質(zhì)量守恒定律，書(shū)寫(xiě)化學(xué)反應(yīng)