計算方法上機作業(yè)

1、計算方法上機報告計算方法上機報告姓 名：學(xué) 號：班 級：上課班級：計算方法上機報告說明： 本次上機實驗使用的編程語言是Matlab語言，編譯環(huán)境為MATLAB 7.11.0，運行平臺為Windows 7。1. 對以下和式計算：，要求： 若只需保留11個有效數(shù)字，該如何進行計算； 若要保留30個有效數(shù)字，則又將如何進行計算；（1） 算法思想1、根據(jù)精度要求估計所加的項數(shù)，可以使用后驗誤差估計，通項為： ；2、為了保證計算結(jié)果的準(zhǔn)確性，寫程序時，從后向前計算；3、使用Matlab時，可以使用以下函數(shù)控制位數(shù)： digits(位數(shù))或vpa(變量，精度為數(shù))（2）算法結(jié)構(gòu)1.；2.for if en

2、d;3.for （3）Matlab源程序clear; %清除工作空間變量clc; %清除命令窗口命令m=input('請輸入有效數(shù)字的位數(shù)m='); %輸入有效數(shù)字的位數(shù)s=0; for n=0:50 t=(1/16n)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6); if t<=10(-m) %判斷通項與精度的關(guān)系 break; endend;fprintf('需要將n值加到n=%dn',n-1); %需要將n值加到的數(shù)值for i=n-1:-1:0 t=(1/16i)*(4/(8*i+1)-2/(8*i+4)-1/(

3、8*i+5)-1/(8*i+6); s=s+t; %求和運算ends=vpa(s,m) %控制s的精度 （4）結(jié)果與分析 當(dāng)保留11位有效數(shù)字時，需要將n值加到n=7， s =3.1415926536； 當(dāng)保留30位有效數(shù)字時，需要將n值加到n=22, s =3.14159265358979323846264338328。 通過上面的實驗結(jié)果可以看出，通過從后往前計算，這種算法很好的保證了計算結(jié)果要求保留的準(zhǔn)確數(shù)字位數(shù)的要求。2. 某通信公司在一次施工中，需要在水面寬度為20米的河溝底部沿直線走向鋪設(shè)一條溝底光纜。在鋪設(shè)光纜之前需要對溝底的地形進行初步探測，從而估計所需光纜的長度，為工程預(yù)算提

4、供依據(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.8110.8010.93 請用合適的曲線擬合所測數(shù)據(jù)點； 預(yù)測所需光纜長度的近似值，作出鋪設(shè)河底光纜的曲線圖；（1）算法思想 如果使用多項式差值，則由于龍格現(xiàn)象，誤差較大，因此，用相對較少的插值數(shù)據(jù)點作插值，可以避免大的誤差，但是如果又希望將所得數(shù)據(jù)點都用上，且所用數(shù)據(jù)點越多越好，可

5、以采用分段插值方式，即用分段多項式代替單個多項式作插值。分段多項式是由一些在相互連接的區(qū)間上的不同多項式連接而成的一條連續(xù)曲線，其中三次樣條插值方法是一種具有較好“光滑性”的分段插值方法。在本題中，假設(shè)所鋪設(shè)的光纜足夠柔軟，在鋪設(shè)過程中光纜觸地走勢光滑，緊貼地面，并且忽略水流對光纜的沖擊。海底光纜線的長度預(yù)測模型如下所示，光纜從A點鋪至B點，在某點處的深度為。海底光纜線的長度預(yù)測模型計算光纜長度時，用如下公式：（2）算法結(jié)構(gòu)1.For 1.1 2.For 2.1 For 2.1.1 3.4.For 4.1 4.2 4.3 5.6.7.獲取M的矩陣元素個數(shù)，存入m8.For 8.1 8.2 8.

6、3 9.10.For 10.1 11.獲取x的元素個數(shù)存入s12.13.For 13.1 if then ；break else 14.（3）Matlab源程序clear;clc; x=0:1:20; %產(chǎn)生從0到20含21個等分點的數(shù)組X=0:0.2: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; %等分點位置的深度數(shù)據(jù)n=length(x); %等分點的數(shù)目N=length(X);% 求三次樣條插值

7、函數(shù)s(x) M=y; for k=2:3; %計算二階差商并存放在M中 for i=n:-1:k; M(i)=(M(i)-M(i-1)/(x(i)-x(i-k+1); endendh(1)=x(2)-x(1); %計算三對角陣系數(shù)a,b,c及右端向量dfor i=2:n-1; h(i)=x(i+1)-x(i); c(i)=h(i)/(h(i)+h(i-1); a(i)=1-c(i); b(i)=2; d(i)=6*M(i+1);end M(1)=0; %選擇自然邊界條件M(n)=0;b(1)=2;b(n)=2;c(1)=0;a(n)=0; d(1)=0; d(n)=0; u(1)=b(1);

8、 %對三對角陣進行LU分解y1(1)=d(1);for k=2:n; l(k)=a(k)/u(k-1); u(k)=b(k)-l(k)*c(k-1); y1(k)=d(k)-l(k)*y1(k-1);endM(n)=y1(n)/u(n); %追趕法求解樣條參數(shù)M(i)for k=n-1:-1:1; M(k)=(y1(k)-c(k)*M(k+1)/u(k);ends=zeros(1,N);for m=1:N; k=1; for i=2:n-1 if X(m)<=x(i); k=i-1; break; else k=i; endend H=x(k+1)-x(k); %在各區(qū)間用三次樣條插值函

9、數(shù)計算X點處的值 x1=x(k+1)-X(m); x2=X(m)-x(k); s(m)=(M(k)*(x13)/6+M(k+1)*(x23)/6+(y(k)-(M(k)*(H2)/6)*x1+(y(k+1)-(M(k+1)*(H2)/6)*x2)/H;end% 計算所需光纜長度L=0; %計算所需光纜長度for i=2:N L=L+sqrt(X(i)-X(i-1)2+(s(i)-s(i-1)2);enddisp('所需光纜長度為 L=');disp(L);figureplot(x,y,'*',X,s,'-') %繪制鋪設(shè)河底光纜的曲線圖xlabe

10、l('位置','fontsize',16); %標(biāo)注坐標(biāo)軸含義ylabel('深度/m','fontsize',16);title('鋪設(shè)河底光纜的曲線圖','fontsize',16);grid;（4）結(jié)果與分析 鋪設(shè)海底光纜的曲線圖如下圖所示： 仿真結(jié)果表明，運用分段三次樣條插值所得的擬合曲線能較準(zhǔn)確地反映鋪設(shè)光纜的走勢圖，計算出所需光纜的長度為 L=26.4844m。3. 假定某天的氣溫變化記錄如下表所示，試用數(shù)據(jù)擬合的方法找出這一天的氣溫變化的規(guī)律；試計算這一天的平均氣溫，并試估計誤差。時刻

11、0123456789101112平均氣溫15141414141516182020232528時刻131415161718192021222324平均氣溫313431292725242220181716（1）算法思想在本題中，數(shù)據(jù)點的數(shù)目較多。當(dāng)數(shù)據(jù)點的數(shù)目很多時，用“多項式插值”方法做數(shù)據(jù)近似要用較高次的多項式，這不僅給計算帶來困難，更主要的缺點是誤差很大。用“插值樣條函數(shù)”做數(shù)據(jù)近似，雖然有很好的數(shù)值性質(zhì)，且計算量也不大，但存放參數(shù)的量很大，且沒有一個統(tǒng)一的數(shù)學(xué)公式來表示，也帶來了一些不便。另一方面，在有的實際問題中，用插值方法并不合適。當(dāng)數(shù)據(jù)點的數(shù)目很大時，要求通過所有數(shù)據(jù)點，可能會失去原

12、數(shù)據(jù)所表示的規(guī)律。如果數(shù)據(jù)點是由測量而來的，必然帶有誤差，插值法要求準(zhǔn)確通過這些不準(zhǔn)確的數(shù)據(jù)點是不合適的。在這種情況下，不用插值標(biāo)準(zhǔn)而用其他近似標(biāo)準(zhǔn)更加合理。通常情況下，是選取使最小，這就是最小二乘近似問題。在本題中，采用“最小二乘法”找出這一天的氣溫變化的規(guī)律，使用二次函數(shù)、三次函數(shù)、四次函數(shù)以及指數(shù)型函數(shù)，計算相應(yīng)的系數(shù)，估算誤差，并作圖比較各種函數(shù)之間的區(qū)別。（2）算法結(jié)構(gòu) 本算法用正交化方法求數(shù)據(jù)的最小二乘近似。假定數(shù)據(jù)以用來生成了，并將作為其最后一列（第列）存放。結(jié)果在中，是誤差。I、使用二次函數(shù)、三次函數(shù)、四次函數(shù)擬合時1.將“時刻值”存入，數(shù)據(jù)點的個數(shù)存入2.輸入擬合多項式函數(shù)的

13、最高項次數(shù)，則擬合多項式函數(shù)為 , 根據(jù)給定數(shù)據(jù)點確定For For 2.1 2.2 3.For 3.1 形成矩陣3.1.1 3.1.2 3.1.3 For 3.1.3.1 3.1.4 3.2 變換到3.2.1 For 3.2.2 3.2.3 For 3.2.3.1 4. 解三角方程 4.1 4.2For 4.2.1 5.計算誤差 II、使用指數(shù)函數(shù)擬合時現(xiàn)將指數(shù)函數(shù)進行變形：將,代入得: 對上式左右取對數(shù)得: 令則可得多項式: （3）Matlab源程序clear; %清除工作空間變量clc; %清除命令窗口命令x=0:24; %將時刻值存入數(shù)組y=15,14,14,14,14,15,16,1

14、8,20,20,23,25,28,31,34,31,29,27,25,24,22,20,18,17,16;,m=size(x); %將數(shù)據(jù)點的個數(shù)存入mT=sum(y(1:m)/m;fprintf('一天的平均氣溫為T=%fn',T); %求一天的平均氣溫 % 二次、三次、四次函數(shù)的最小二乘近似h=input('請輸入擬合多項式的最高項次數(shù)='); %根據(jù)給定數(shù)據(jù)點生成矩陣Gn=h+1;G=;for j=0:(n-1) g=x.j; %g(x)按列排列 G=vertcat(G,g); %g垂直連接GendG=G' %轉(zhuǎn)置得到矩陣Gfor i=1:m %將

15、數(shù)據(jù)y作為G的最后一列（n+1列） G(i,n+1)=y(i);endG; for k=1:n %形成矩陣Q(k) if G(k,k)>0; sgn=1; elseif G(k,k)=0; sgn=0; else sgn=-1; end sgm=-sgn*sqrt(sum(G(k:m,k).2); w=zeros(1,n); w(k)=G(k,k)-sgm; for j=k+1:m w(j)=G(j,k); end bt=sgm*w(k); G(k,k)=sgm; %變換Gk-1到Gk for j=k+1:n+1 t=sum(w(k:m)*G(k:m,j)/bt; for i=k:m;

16、G(i,j)=G(i,j)+t*w(i); end endend A (n)=G(n,n+1)/G(n,n); %解三角方程求系數(shù)Afor i=n-1:-1:1 A (i)=(G(i,n+1)-sum(G(i,i+1:n).*A (i+1:n)/G(i,i);end e=sum(G(n+1:m,n+1).2); %計算誤差efprintf('%d次函數(shù)的系數(shù)是：',h); %輸出系數(shù)a及誤差edisp(A);fprintf('使用%d次函數(shù)擬合的誤差是%f：',h,e);t=0:0.05:24;A=fliplr(A); %將系數(shù)數(shù)組左右翻轉(zhuǎn) Y=poly2sym

17、(A); %將系數(shù)數(shù)組轉(zhuǎn)化為多項式subs(Y,'x',t);Y=double(ans);figure(1)plot(x,y,'k*',t,Y,'r-'); %繪制擬合多項式函數(shù)圖形xlabel('時刻'); %標(biāo)注坐標(biāo)軸含義ylabel('平均氣溫');title(num2str(n-1),'次函數(shù)的最小二乘曲線');grid;% 指數(shù)函數(shù)的最小二乘近似yy=log(y);n=3;G=;GG=;for j=0:(n-1) g=x.j; %g(x)按列排列 G=vertcat(G,g); %g垂直連

18、接G gg=t.j; %g(x)按列排列 GG=vertcat(GG,gg); %g垂直連接GendG=G' %轉(zhuǎn)置得到矩陣Gfor i=1:m %將數(shù)據(jù)y作為G的最后一列（n+1列） G(i,n+1)=yy(i);endG; for k=1:n %形成矩陣Q(k) if G(k,k)>0; sgn=1; elseif G(k,k)=0; sgn=0; else sgn=-1; end sgm=-sgn*sqrt(sum(G(k:m,k).2); w=zeros(1,n); w(k)=G(k,k)-sgm; for j=k+1:m w(j)=G(j,k); end bt=sgm*

19、w(k); G(k,k)=sgm; %變換Gk-1到Gk for j=k+1:n+1 t=sum(w(k:m)*G(k:m,j)/bt; for i=k:m; G(i,j)=G(i,j)+t*w(i); end endend A(n)=G(n,n+1)/G(n,n); %解三角方程求系數(shù)Afor i=n-1:-1:1 A (i)=(G(i,n+1)-sum(G(i,i+1:n).*A (i+1:n)/G(i,i);end b=-A(3);c=A(2)/(2*b);a=exp(A(1)+b*(c2);G(n+1:m,n+1)=exp(sum(G(n+1:m,n+1).2);e=sum(G(n+1

20、:m,n+1).2); %計算誤差efprintf('n指數(shù)函數(shù)的系數(shù)是：a=%f,b=%f,c=%f',a,b,c); %輸出系數(shù)及誤差efprintf('n使用指數(shù)函數(shù)擬合的誤差是：%f',e);t=0:0.05:24;YY=a.*exp(-b.*(t-c).2);figure(2)plot(x,y,'k*',t,YY,'r-'); %繪制擬合指數(shù)函數(shù)圖形xlabel('時刻'); %標(biāo)注坐標(biāo)軸含義ylabel('平均氣溫');title('指數(shù)函數(shù)的最小二乘曲線');grid;

21、（4）結(jié)果與分析a、二次函數(shù)：一天的平均氣溫為： 21.20002次函數(shù)的系數(shù): 8.3063 2.6064 -0.0938使用2次函數(shù)擬合的誤差是：280.339547二次函數(shù)的最小二乘曲線如下圖所示：b、三次函數(shù)：一天的平均氣溫為： 21.20003次函數(shù)的系數(shù): 13.3880 -0.2273 0.2075 -0.0084使用3次函數(shù)擬合的誤差是： 131.061822 三次函數(shù)的最小二乘曲線如下圖所示：c、四次函數(shù)：一天的平均氣溫為： 21.20004次函數(shù)的系數(shù): 16.7939 -3.7050 0.8909 -0.0532 0.0009使用4次函數(shù)擬合的誤差是：59.04118四次

22、函數(shù)的最小二乘曲線如下圖所示：d、指數(shù)函數(shù)：一天的平均氣溫為： 21.2000指數(shù)函數(shù)的系數(shù)是： a=26.160286，b=0.004442，c=14.081900使用指數(shù)函數(shù)擬合的誤差是： 57.034644指數(shù)函數(shù)的最小二乘曲線如下圖所示：通過上述幾種擬合可以發(fā)現(xiàn)，多項式的次數(shù)越高，計算擬合的效果越好，誤差越小，說明結(jié)果越準(zhǔn)確；同時，指數(shù)多項式擬合的次數(shù)雖然不高，但誤差最小，說明結(jié)果最準(zhǔn)確。4.設(shè)計算法，求出非線性方程的所有根，并使誤差不超過。（1）算法思想 首先，研究函數(shù)的形態(tài)，確定根的范圍；通過剖分區(qū)間的方法確定根的位置，然后利用二分法的基本原理進行求解，找到滿足精度要求的解。二分法

23、是產(chǎn)生一串區(qū)間，使新區(qū)間是舊區(qū)間的一個子區(qū)間，其長度是的一半，且有一個端點是的一個端點。由區(qū)間確定區(qū)間的方法是計算區(qū)間的中點若，則取，否則取，重復(fù)這一過程即可。顯然，每次迭代使區(qū)間長度減小一半，故二分法總是收斂的。（2）算法結(jié)構(gòu)1.；2.If then stop3.If then輸出作為根； stop4.If then輸出作為根； stop5.6. If then輸出作為根； stop7. 8. If then輸出作為根； 9. Ifthen 9.1 ；else9.2 ；10.go to 5 （3）Matlab源程序x=-100:100;y=6*(x.5)-45*(x.2)+20; %非線性方

24、程組的表達式g=; for i=-100:1:100 %確定根所在的區(qū)間 k=i+1; if (y(x=i).*y(x=k)<eps) %區(qū)間長度為1 g=g i; endendsyms x;f=6*x5-45*x2+20;n=length(g); %確定根的個數(shù) for j=1:n x0=g(j); %求根區(qū)間左端點 x1=g(j)+1; %求根區(qū)間右端點while (x1-x0)>=10(-4) if subs(f,x,x0)*subs(f,x,(x0+x1)/2)>eps x0=(x0+x1)/2; else x1=(x0+x1)/2; endendroot=x0 %輸

25、出方程的根 end（4）結(jié)果與分析該非線性方程組有三個實根，分別為1.8708，0.6812，-0.6545，且滿足誤差要求。5. 編寫程序?qū)崿F(xiàn)大規(guī)模方程組的列主元高斯消去法程序，并對所附的方程組進行求解。針對本專業(yè)中所碰到的實際問題，提煉一個使用方程組進行求解的例子，并對求解過程進行分析、求解。（1）算法思想 高斯消去法是利用現(xiàn)行方程組初等變換中的一種變換，即用一個不為零的數(shù)乘一個方程后加只另一個方程，使方程組變成同解的上三角方程組，然后再自下而上對上三角方程組求解。列主元消去法是當(dāng)高斯消元到第步時，從列的以下（包括）的各元素中選出絕對值最大的，然后通過行交換將其交換到的位置上。交換系數(shù)矩陣

26、中的兩行（包括常數(shù)項），只相當(dāng)于兩個方程的位置交換了，因此，列選主元不影響求解的結(jié)果。 程序的核心就是高斯列主元消去法。根據(jù)教材提供的算法，編寫列主元消去法的子函數(shù)與適應(yīng)于超大規(guī)模超出系統(tǒng)內(nèi)存的方程組的改編程序。同時，在Gauss消去過程中，適當(dāng)交換方程的順序?qū)ΡＷC消去過程能順利進行及計算解的精確度都是有必要的，交換方程的原則是使中，絕對值最大的一個換到（k,k)位置而成為第k步消去的主元，這就是列主元Gauss消去法。（2） 算法結(jié)構(gòu)1、數(shù)據(jù)文件的文件名為：文件名+.dat2、數(shù)據(jù)文件中的數(shù)據(jù)為二進制記錄結(jié)構(gòu)，分為以下四個部分：(1) 文件頭部分，其結(jié)構(gòu)： typedef struct lo

27、ng int id; long int ver; long int n; 其中：id:為該數(shù)據(jù)文件的標(biāo)識，值為0xF1E1D1A0，即為：十六進制的F1E1D1A0 ver:為數(shù)據(jù)文件的版本號，值為16進制數(shù)據(jù),版本號說明0x101系數(shù)矩陣為非壓縮格式稀疏矩陣0x102系數(shù)矩陣為非壓縮格式帶狀對角陣0x201系數(shù)矩陣為壓縮格式稀疏矩陣0x202系數(shù)矩陣為壓縮格式帶狀對角陣 n:表示方程的階數(shù) (2) 文件頭2:此部分說明為條狀矩陣的上下帶寬,結(jié)構(gòu)： typedef struct long int q; / 為上帶寬 long int p; / 為下帶寬 (3) 系數(shù)矩陣 a.如存貯格式非為壓縮

28、方式，則按行方式存貯系數(shù)矩陣中的每一個元素，個數(shù)為n*n，類型為float型； b.如果存貯格式是壓縮方式,則按行方式存貯,每行中只存放上下帶寬內(nèi)的非零元素,即，每行中存貯的最多元素為p+q+1個。(4) 右端系數(shù) 按順序存貯右端系數(shù)的每個元素，個數(shù)為n個，類型為float型3、二進制文件的讀?。篺=fopen('fun003.dat','r'); %打開文件，.dat文件放在 m 文件同一目錄下，a=fread(f,3,'uint') %讀取頭文件，3-讀取前 3 個，若讀取壓縮格式的，頭文件為 5 個b=fread(f,inf,'fl

29、oat'); %讀取剩下的文件，float 型id=dec2hex(a(1);ver=dec2hex(a(2); %這兩句是進行進制轉(zhuǎn)換，讀取 id 與ver1. A的階數(shù)2. For 2.1找滿足 2.2For 2.2.1 2.3 2.4 For 2.4.1 2.4.2 For 2.4.2.1 2.4.3 For （3）Matlab源程序clear; %清除工作空間變量clc; %清除命令窗口命令 % 讀取系數(shù)矩陣f,p=uigetfile('*.dat','選擇數(shù)據(jù)文件'); %讀取數(shù)據(jù)文件num=5; %輸入系數(shù)矩陣文件頭的個數(shù)name=strca

30、t(p,f);file=fopen(name,'r');head=fread(file,num,'uint'); %讀取二進制頭文件id=dec2hex(head(1); %讀取標(biāo)識符fprintf('文件標(biāo)識符為');idver=dec2hex(head(2); %讀取版本號fprintf('文件版本號為');vern=head(3); %讀取階數(shù)fprintf('矩陣A的階數(shù)');nq=head(4); %上帶寬fprintf('矩陣A的上帶寬');q p=head(5); %下帶寬fprint

31、f('矩陣A的下帶寬');p dist=4*num;fseek(file,dist,'bof'); %把句柄值轉(zhuǎn)向第六個元素開頭處A,count=fread(file,inf,'float'); %讀取二進制文件，獲取系數(shù)矩陣fclose(file); %關(guān)閉二進制頭文件% 對非壓縮帶狀矩陣進行求解if ver='102', a=zeros(n,n); for i=1:n, for j=1:n, a(i,j)=A(i-1)*n+j); %求系數(shù)矩陣a(i,j) end end b=zeros(n,1); for i=1:n, b(

32、i)=A(n*n+i); end for k=1:n-1, %列主元高斯消去法 m=k; for i=k+1:n, %尋找主元 if abs(a(m,k)<abs(a(i,k) m=i; end end if a(m,k)=0 %遇到條件終止 disp('錯誤！') return end for j=1:n, %交換元素位置得主元 t=a(k,j); a(k,j)=a(m,j); a(m,j)=t; t=b(k); b(k)=b(m); b(m)=t; end for i=k+1:n, %計算l(i,k)并將其放到a(i,k)中 a(i,k)=a(i,k)/a(k,k);

33、 for j=k+1:n a(i,j)=a(i,j)-a(i,k)*a(k,j); end b(i)=b(i)-a(i,k)*b(k); end end x=zeros(n,1); %回代過程 x(n)=b(n)/a(n,n); for k=n-1:-1:1, x(k)=(b(k)-sum(a(k,k+1:n)*x(k+1:n)/a(k,k); endend% 對壓縮帶狀矩陣進行求解if ver='202', %高斯消去法 m=p+q+1; a=zeros(n,m); for i=1:1:n for j=1:1:m a(i,j)=A(i-1)*m+j); %求a(i,j) en

34、d end b=zeros(n,1); for i=1:1:n b(i)=A(n*m+i); %求b(i) end for k=1:1:(n-1) %開始消去過程 if a(k,(p+1)=0 disp('錯誤！'); break; end st1=n; if (k+p)<n st1=k+p; end for i=(k+1):1:st1 a(i,(k+p-i+1)=a(i,(k+p-i+1)/a(k,(p+1); for j=(k+1):1:(k+q) a(i,j+p-i+1)=a(i,j+p-i+1)-a(i,k+p-i+1)*a(k,j+p-k+1); end b(i)=b(i)-a(i,k+p-i+1)*b(k); end end x=zeros(n,1); %回代過程 x(n)=b(n)/a(n,p+1); sum=0; for k=(n-1):-1:1 sum=b(k); st2=n; if (k+q)<n st2=k+q; end for j=(k+1):1:st2 sum=sum-a(k,j+p-k+1)*x(j); end x(k)=sum/a(k,p+1)