計算方法上機報告姓名：學號：班級：目錄TOC\o"1-3"\h\u31920題目一 -4-87341.1題目內容 -4-289011.2算法思想 -4-119191.3Matlab源程序 -5-5141.4計算結果及總結 -5-775題目二 -7-275972.1題目內容 -7-51662.2算法思想 -7-164772.3Matlab源程序 -8-90972.4計算結果及總結 -9-31690題目三 -11-7493.1題目內容 -11-159393.2算法思想 -11-227063.3Matlab源程序 -13-102433.4計算結果及總結 -14-4708題目四 -15-280904.1題目內容 -15-250824.2算法思想 -15-126414.3Matlab源程序 -15-168704.4計算結果及總結 -16-19095題目五 -18-77985.1題目內容 -18-171405.2算法思想 -18-305535.3Matlab源程序 -18-117725.3.1非壓縮帶狀對角方程組 -18-87825.3.2壓縮帶狀對角方程組 -20-4035.4實驗結果及分析 -22-246445.4.1Matlab運行結果 -22-315205.4.2總結分析 -24-160675.5本專業(yè)算例 -24-26922學習感悟 -27-題目一1.1題目內容計算以下和式：，要求：

(1)若保留11個有效數字，給出計算結果，并評價計算的算法；

(2)若要保留30個有效數字，則又將如何進行計算。1.2算法思想在程序編寫中需要把握以下幾點：①隨著n值的增加，和式的項遞減速度很快，因此我們可以認為，在確定為某一精度的前提下，n達到一定的值，加下一項將不會對最終的加和產生影響，首先我們應找到n值。②根據精度要求估計所加的項數，可以使用后驗誤差估計，通項為：③為減小舍入誤差，在計算s時所采用的方法是逆序相加，其依據是：兩個數量級相差較大的數字相加減時，較小數的有效數字會被喪失，從而導致最后的運算結果失真。為避免“大數吃小數”現象的發(fā)生，采用逆序相加。④對于實現30位有效數字，則調用從工具箱中digits(位數)或vpa（變量，精度位數）即可實現。1.3Matlab源程序>>clear;>>clc;>>m=input('輸入需要求的有效數字位數m=');s=0forn=0:200%尋找滿足條件的最小ns=(1/16^n)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6));ifs<=10^(-m)%當項小于10^-m時，停止循環(huán)breakendend;fprintf('n值加至n%d\n',n-1);%需要將n值加到的數值fori=n-1:-1:0%逐項逆序相加求和s=(1/16^i)*(4/(8*i+1)-2/(8*i+4)-1/(8*i+5)-1/(8*i+6));t=t+s;ends=vpa(t,m)1.4計算結果及總結①輸入需要求的有效數字位數m=11t=0n值加至n7s=3.1415926536②t=0n值加至n22s=3.14159265358979311599796346854從上述的算法思想中可以看出，運算中不僅要滿足誤差要求，還要盡可能地減少計算量，此外還要考慮舍入誤差的影響，這時就要對所運算數據的性質進行分析，設置合適的算法，從而提高運算的精度。而逆序算法能夠很好的滿足上述要求。題目二2.1題目內容某通信公司在一次施工中，需要在水面寬度為20米的河溝底部沿直線走向鋪設一條溝底光纜。在鋪設光纜之前需要對溝底的地形進行初步探測，從而估計所需光纜的長度，為工程預算提供依據。已探測到一組等分點位置的深度數據(單位:米)如下表所示：分點0123456深度9.018.967.967.978.029.0510.13分點78910111213深度11.1812.2613.2813.3212.6111.2910.22分點14151617181920深度9.157.907.958.869.8110.8010.93(1)請用合適的曲線擬合所測數據點；(2)估算所需光纜長度的近似值，并作出鋪設河底光纜的曲線圖；2.2算法思想利用曲線擬合數據點，即利用數據點擬合差值多項式，我們可以利用Newton法進行擬合，也可以用復化Simpson求積公式、三次樣條插值來擬合，但三次樣條插值使用方程組計算增大了計算量，同時還要附加邊界條件，分段三次樣條插值對圖形的控制能力還不夠靈活。因此這里用Newton形式的差值多項式進行擬合。首先計算出各差商，然后計算出Newton差值多項式的每一項，最后將所有項相加，即可計算出Newton差值多項式，然后利用所得的差值多項式一次算出多個點的函數值。MATLAB的plot函數進行繪圖。計算長度近似值，只需將每隔兩點之間的距離算出，然后一次相加，所得的折線長度，即為長度的近似值。2.3Matlab源程序Untitled2clearclcx=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.90-7.95-8.86-9.81-10.80-10.93];%輸入給定的數據點xi=0:20;[Nx,Ni]=Newton(x,y,xi);%調用函數，建立Newton差值多項式plot(xi,Ni);%繪制擬合的曲線圖long=0;%為長度賦初值fori=1:20%將每一段折線相加算出長度的近似值long=long+sqrt(1+((y(i)-y(i+1))^2));enddisp('需要的光纜長度為')%顯示需要的光纜長度disp(long)Newton插值法function[Nx,N0]=Newton(X,Y,x0)n=size(X);%插值點個數y=Y;Nx=Y(1);N=1;fori=1:n-1%計算Newton插值多項式forj=i+1:nyi(j)=(y(j)-y(i))/(X(j)-X(i));endm(i)=yi(i+1);N=N*(x-X(i));Nx=Nx+N*m(i);y=yi;endN0=subs(Nx,'x',x0);2.4計算結果及總結針對上述Matlab程序，鋪設海底光纜的曲線圖如下圖所示：結果如下：Nl=26.4844，即為所求近似計算光纜長度。本題利用Newton法進行擬合，既簡單又使用，運行也可以用復化Simpson求積公式、三次樣條插值來擬合，但三次樣條插值使用方程組計算增大了計算量，同時還要附加邊界條件，分段三次樣條插值對圖形的控制能力還不夠靈活。題目三3.1題目內容假定某天的氣溫變化記錄如下表所示，試用數據擬合的方法找出這一天的氣溫變化的規(guī)律;試計算這一天的平均氣溫,并試估計誤差。時刻0123456789101112平均氣溫15141414141516182020232528時刻131415161718192021222324平均氣溫3134312927252422201817163.2算法思想在本題中，數據點的數目較多。當數據點的數目很多時，用“多項式插值”方法做數據近似要用較高次的多項式，這不僅給計算帶來困難，更主要的缺點是誤差很大。用“插值樣條函數”做數據近似，雖然有很好的數值性質，且計算量也不大，但存放參數的量很大，且沒有一個統(tǒng)一的數學公式來表示，也帶來了一些不便。另一方面，在有的實際問題中，用插值方法并不合適。當數據點的數目很大時，要求通過所有數據點，可能會失去原數據所表示的規(guī)律。如果數據點是由測量而來的，必然帶有誤差，插值法要求準確通過這些不準確的數據點是不合適的。在這種情況下，不用插值標準而用其他近似標準更加合理。通常情況下，是選取使最小，這就是最小二乘近似問題。在本題中，采用“最小二乘法”找出這一天的氣溫變化的規(guī)律，使用二次函數、三次函數、四次函數以及指數型函數，計算相應的系數，估算誤差，并作圖比較各種函數之間的區(qū)別。給定數據點和一組函數。求數(假定)，使函數滿足為計算方便，只需考慮最小二乘近似問題的解應使上式成立；在中，它是ak的二次函數，因此其取極小值的必要條件為，即經過推導知，若記則最小二乘問題的法方程為解此法方程，就可以得到各項系數，數據擬合之后便可知道一天中溫度變化趨勢。3.3Matlab源程序clearclcA=0:24;B=[15,14,14,14,14,15,16,18,20,20,23,25,28,31,32,41,29,27,25,24,22,20,18,17,16];%輸入給定的數據G=zeros(25,7);%給矩陣G賦初值fori=1:25;forj=1:7;G(i,j)=A(i)^(j-1);endendx1=(G'*G)\(G'*B');%按法方程法算出各項系數a=zeros(1,7);forn=1:7a(1,n)=x1(8-n,1);end%各項系倒序排列x=0:0.1:24;plot(x,polyval(a,x),'b')%繪制曲線holdonfx=polyval(a,x)sum=0fori=1:240%計算曲線下面積sum=sum+(fx(i)+fx(i+1))*0.1/2endaverage=sum/24%計算平均溫度e=0fori=1:10:241%計算誤差e=e+(fx(i)-B((i+9)/10)).^2enddisp('平均溫度為')%輸出平均溫度disp(average)disp('誤差為')%輸出誤差值disp(e)3.4計算結果及總結鋪設海底光纜的曲線圖如下圖所示：這一天平均溫度為average=21.7729估算誤差e=118.7770本算法利用最小二乘法，這里用法方程法進行求解，相比于正交化方法，此方法運算量較大，但對于本題來說，數據量并不大，因此在電腦上運行也不會很慢，這里生成的矩陣G由7列，因此有六階算法的精度。題目四4.1題目內容設計算法，求出非線性方程的所有實根，并使誤差不超過。4.2算法思想采用Newton法進行求解，牛頓法的迭代函數該式的導數則迭代函數迭代的初值為-1,0.5,1.54.3Matlab源程序主程序x1=-1;%初值為-1時whileabs(6*x1^5-45*x1^2+20)>0.0001%誤差大于0.0001時繼續(xù)循環(huán)x1=x1-(6*x1^5-45*x1^2+20)/(30*x1^4-90*x1);%Newton迭代法endx2=0.5;%初值為0.5時whileabs(6*x2^5-45*x2^2+20)>0.0001x2=x2-(6*x2^5-45*x2^2+20)/(30*x2^4-90*x2);endx3=1.5;%初值為1.5時whileabs(6*x3^5-45*x3^2+20)>0.0001x3=x3-(6*x3^5-45*x3^2+20)/(30*x3^4-90*x3);endfprintf('x1=%.6f,x2=%.6f,x3=%.6f',x1,x2,x3)%輸出三個根繪圖程序x=-2:0.1:2y=6*x.^5-45*x.^2+20plot(x,y)xlabel('X軸');ylabel('Y軸');%坐標軸表示對象標簽gridon;%顯示網格線axison;%顯示坐標軸4.4計算結果及總結x1=-0.654542x2=0.681174x3=1.870799所得曲線所給的方程形式較為簡單，導數較為好求所以用牛頓法非常方便，先作圖知道根的大致范圍后再取初值可以有效的減小迭代次數。題目五5.1題目內容線性方程組求解。（1）編寫程序實現大規(guī)模方程組的高斯消去法程序，并對所附的方程組進行求解。所附方程組的類型為對角占優(yōu)的帶狀方程組。（2）針對本專業(yè)中所碰到的實際問題，提煉一個使用方程組進行求解的例子，并對求解過程進行分析、求解。5.2算法思想消去法的中心是“降維”，即求解n元方程組的問題轉化為先解n-1元方程組，一旦這個n-1元方程組的解取得，則剩余的一個未知量自然可以求得。這樣逐步減少未知量個數的方法，是求解多元方程組的一個重要思想。Gauss消去過程中，適當交換方程的順序對保證消去過程能夠順利進行及計算解的精確度都是必要的。消去過程中產生的數稱為第k步消去的主元。交換方程的原則是使中，絕對值最大的一個換到(k,k)位置而成為第k步消去的主元，帶有這種交換的Gauss消去法為列主元Gauss消去法。5.3Matlab源程序5.3.1非壓縮帶狀對角方程組（1）dat51.dat、dat52.dat主程序:%非壓縮格式dat51.dat;dat52.datclear;clc;formatshort%讀出數據fid=fopen('E:\2015上機題目\實驗五\dat52.dat','r');D=fread(fid,3,'long',0);n=D(3)%矩陣維數offset=20;fseek(fid,offset,'bof');A=zeros(n);B=zeros(1,n);C=fread(fid,n*(n+1),'float',0)fori=1:nforj=1:nA(i,j)=C((i-1)*n+j);endendA;fori=1:nB(i)=C(n*n+i);endB;x=GAUSSPP(A,B)（2）列主元高斯消去法Matlab函數程序：%列主元高斯消去法functionx=GAUSSPP(A,B)n=length(B);%消去過程%找出列主元fork=1:n-1MAX=abs(A(k,k));Rmax=k;fori=k+1:nifabs(A(i,k))>MAXMAX=abs(A(i,k));Rmax=i;endendifA(Rmax,k)==0break;end%交換當前行和主元所在行forj=1:ntemp=A(k,j);A(k,j)=A(Rmax,j);A(Rmax,j)=temp;endtemp=B(k);B(k)=B(Rmax);B(Rmax)=temp;fori=k+1:nA(i,k)=A(i,k)/A(k,k);forj=k+1:nA(i,j)=A(i,j)-A(i,k)*A(k,j);endB(i)=B(i)-A(i,k)*B(k);endendA%回代過程x=zeros(1,n);x(n)=B(n)/A(n,n);fork=n-1:-1:1S=B(k);forj=k+1:nS=S-A(k,j)*x(j);endx(k)=S/A(k,k);end5.3.2壓縮帶狀對角方程組（1）dat53.dat主程序:%壓縮格式dat143.datclear;formatshort%讀出數據fid=fopen('E:\2015上機題目\實驗五\dat53.dat','r');D=fread(fid,5,'long',0);n=D(3);%矩陣維數p=D(5);%上帶寬q=D(4);%下帶寬A=fread(fid,[p+q+1,n],'float');B=fread(fid,[n,1],'float');C=A';D=B';%采用求解壓縮型的列主元Gauss消去法x=GAUSSPP1(C,D,n,p,q);（2）針對壓縮格式的帶狀矩陣，尤其是大型矩陣的列主元Gauss消去法子程序GAUSSPP1如下：%針對壓縮格式的帶狀方程組的列主元Gauss消去法functionx=GAUSSPP1(A,B,n,p,q)%消去過程fork=1:n-1pmax=k+p;ifpmax>npmax=n;endfori=k+1:pmaxA(i,k-i+p+1)=A(i,k-i+p+1)/A(k,p+1);qmax=k+q;ifqmax>nqmax=n;endforj=k+1:qmaxA(i,j-i+p+1)=A(i,j-i+p+1)-A(i,k-i+p+1)*A(k,j-k+p+1);endB(i)=B(i)-B(k)*A(i,k-i+p+1);endend%回代過程x(n)=B(n)/A(n,p+1);fork=n-1:-1:1sum=0;qmax=k+q;ifqmax>nqmax=n;endforj=k+1:qmaxsum=sum+A(k,j-k+p+1)*x(j);endx(k)=(B(k)-sum)/A(k,p+1);endend5.4實驗結果及分析5.4.1Matlab運行結果方程組一：文件標識符為id=F1E1D1A0文件版本號為ver=102矩陣A的階數n=15矩陣A的上帶寬q=6矩陣A的下帶寬p=6方程組的的解為：1.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000方程組二：文件標識符為id=F1E1D1A0文件版本號為ver=102矩陣A的階數n=2050矩陣A的上帶寬q=6矩陣A的下帶寬p=5方程組的的解為：1.96001.96001.96001.96001.96001.96001.96001.96001.96001.9600方程組三：文件標識符為id=F1E1D1A0文件版本號為ver=202矩陣A的階數n=43512矩陣A的上帶寬q=4矩陣A的下帶寬p=3方程組的的解為：3.14133.14133.14133.14133.14133.14133.14133.14133.14133.14133.14135.