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

2、b源程序clear; %清除工作空間變量clc; %清除命令窗口命令m=input('請輸入有效數字的位數m='); %輸入有效數字的位數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) %判斷通項與精度的關系break;endend;fprintf('需要將n值加到n=%dn',n-1); %需要將n值加到的數值for i=n-1:-1:0 t=(1/16i)*(4/(8*i+1)-2/(8*i+4)-1/(8*i+5)-1/(8*i+6); s

3、=s+t; %求和運算ends=vpa(s,m) %控制s的精度 （4）結果與分析當保留11位有效數字時，需要將n值加到n=7， s =3.1415926536；當保留30位有效數字時，需要將n值加到n=22,s =3.8。通過上面的實驗結果可以看出，通過從后往前計算，這種算法很好的保證了計算結果要求保留的準確數字位數的要求。2. 某通信公司在一次施工中，需要在水面寬度為20米的河溝底部沿直線走向鋪設一條溝底光纜。在鋪設光纜之前需要對溝底的地形進行初步探測，從而估計所需光纜的長度，為工程預算提供依據。已探測到一組等分點位置的深度數據(單位:米)如下表所示：分點0123456深度9.018.96

4、7.967.978.029.0510.13分點78910111213深度11.1812.2613.2813.3212.6111.2910.22分點14151617181920深度9.157.907.958.869.8110.8010.93請用合適的曲線擬合所測數據點；預測所需光纜長度的近似值，作出鋪設河底光纜的曲線圖；（1）算法思想如果使用多項式差值，則由于龍格現象，誤差較大，因此，用相對較少的插值數據點作插值，可以避免大的誤差，但是如果又希望將所得數據點都用上，且所用數據點越多越好，可以采用分段插值方式，即用分段多項式代替單個多項式作插值。分段多項式是由一些在相互連接的區間上的不同多項式連接

5、而成的一條連續曲線，其中三次樣條插值方法是一種具有較好“光滑性”的分段插值方法。在本題中，假設所鋪設的光纜足夠柔軟，在鋪設過程中光纜觸地走勢光滑，緊貼地面，并且忽略水流對光纜的沖擊。海底光纜線的長度預測模型如下所示，光纜從A點鋪至B點，在某點處的深度為。海底光纜線的長度預測模型計算光纜長度時，用如下公式：（2）算法結構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的矩陣元素個數，存入m8.For 8.1 8.2 8.3 9.10.For 10.1 11.獲取x的元素個數存入s12.13.For 13.1 if then ；

6、breakelse 14.（3）Matlab源程序clear;clc; x=0:1:20; %產生從0到20含21個等分點的數組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; %等分點位置的深度數據n=length(x); %等分點的數目N=length(X);% 求三次樣條插值函數s(x) M=y; for k=2:3; %計算二階差商并存放在M中for i=n:-1:k; M(i)=

7、(M(i)-M(i-1)/(x(i)-x(i-k+1); endendh(1)=x(2)-x(1); %計算三對角陣系數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);endM(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); %對三對角陣進行LU分解y1(1)=d(1);for k=2:n; l(k)=a(k)/u(k-1);u(k)=b(k)

8、-l(k)*c(k-1);y1(k)=d(k)-l(k)*y1(k-1);endM(n)=y1(n)/u(n);%追趕法求解樣條參數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-1if X(m)<=x(i);k=i-1; break;else k=i;endend H=x(k+1)-x(k); %在各區間用三次樣條插值函數計算X點處的值 x1=x(k+1)-X(m); x2=X(m)-x(k); s(m)=(M(k)*(x13)/6+M(k+1)*(x23

9、)/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,'-') %繪制鋪設河底光纜的曲線圖xlabel('位置','fontsize',16); %標注坐標軸含義ylabel('深度/m',

10、'fontsize',16);title('鋪設河底光纜的曲線圖','fontsize',16);grid;（4）結果與分析 鋪設海底光纜的曲線圖如下圖所示：仿真結果表明，運用分段三次樣條插值所得的擬合曲線能較準確地反映鋪設光纜的走勢圖，計算出所需光纜的長度為 L=26.4844m。3.假定某天的氣溫變化記錄如下表所示，試用數據擬合的方法找出這一天的氣溫變化的規律；試計算這一天的平均氣溫，并試估計誤差。時刻0123456789101112平均氣溫15141414141516182020232528時刻13141516171819202122232

11、4平均氣溫313431292725242220181716（1）算法思想在本題中，數據點的數目較多。當數據點的數目很多時，用“多項式插值”方法做數據近似要用較高次的多項式，這不僅給計算帶來困難，更主要的缺點是誤差很大。用“插值樣條函數”做數據近似，雖然有很好的數值性質，且計算量也不大，但存放參數的量很大，且沒有一個統一的數學公式來表示，也帶來了一些不便。另一方面，在有的實際問題中，用插值方法并不合適。當數據點的數目很大時，要求通過所有數據點，可能會失去原數據所表示的規律。如果數據點是由測量而來的，必然帶有誤差，插值法要求準確通過這些不準確的數據點是不合適的。在這種情況下，不用插值標準而用其他近

12、似標準更加合理。通常情況下，是選取使最小，這就是最小二乘近似問題。在本題中，采用“最小二乘法”找出這一天的氣溫變化的規律，使用二次函數、三次函數、四次函數以及指數型函數，計算相應的系數，估算誤差，并作圖比較各種函數之間的區別。（2）算法結構本算法用正交化方法求數據的最小二乘近似。假定數據以用來生成了，并將作為其最后一列（第列）存放。結果在中，是誤差。I、使用二次函數、三次函數、四次函數擬合時1.將“時刻值”存入，數據點的個數存入2.輸入擬合多項式函數的最高項次數，則擬合多項式函數為,根據給定數據點確定For For 2.1 2.23.For 3.1 形成矩陣3.1.1 3.1.23.1.3 F

13、or 3.1.3.1 3.1.43.2 變換到3.2.1 For 3.2.23.2.3 For3.2.3.1 4. 解三角方程4.1 4.2For 4.2.1 5.計算誤差II、使用指數函數擬合時現將指數函數進行變形：將,代入得:對上式左右取對數得:令則可得多項式:（3）Matlab源程序clear; %清除工作空間變量clc; %清除命令窗口命令x=0:24; %將時刻值存入數組y=15,14,14,14,14,15,16,18,20,20,23,25,28,31,34,31,29,27,25,24,22,20,18,17,16;,m=size(x); %將數據點的個數存入mT=sum(y(

14、1:m)/m;fprintf('一天的平均氣溫為T=%fn',T); %求一天的平均氣溫 % 二次、三次、四次函數的最小二乘近似h=input('請輸入擬合多項式的最高項次數='); %根據給定數據點生成矩陣Gn=h+1;G=;for j=0:(n-1) g=x.j; %g(x)按列排列 G=vertcat(G,g); %g垂直連接GendG=G' %轉置得到矩陣Gfor i=1:m %將數據y作為G的最后一列（n+1列） G(i,n+1)=y(i);endG; for k=1:n %形成矩陣Q(k)if G(k,k)>0; sgn=1;elsei

15、f 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到Gkfor 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);endendendA (n)=G(n,n+1)/G(n,n); %解三角方程求系數Afor i=n-1:-1:1A (i)=(G(i,n+1)-

16、sum(G(i,i+1:n).*A (i+1:n)/G(i,i);ende=sum(G(n+1:m,n+1).2); %計算誤差efprintf('%d次函數的系數是：',h); %輸出系數a及誤差edisp(A);fprintf('使用%d次函數擬合的誤差是%f：',h,e);t=0:0.05:24;A=fliplr(A); %將系數數組左右翻轉 Y=poly2sym(A); %將系數數組轉化為多項式subs(Y,'x',t);Y=double(ans);figure(1)plot(x,y,'k*',t,Y,'r-

17、9;); %繪制擬合多項式函數圖形xlabel('時刻'); %標注坐標軸含義ylabel('平均氣溫');title(num2str(n-1),'次函數的最小二乘曲線');grid;% 指數函數的最小二乘近似yy=log(y);n=3;G=;GG=;for j=0:(n-1) g=x.j; %g(x)按列排列 G=vertcat(G,g); %g垂直連接G gg=t.j; %g(x)按列排列 GG=vertcat(GG,gg); %g垂直連接GendG=G' %轉置得到矩陣Gfor i=1:m %將數據y作為G的最后一列（n+1列） G

18、(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*w(k); G(k,k)=sgm; %變換Gk-1到Gkfor 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);endend

19、endA(n)=G(n,n+1)/G(n,n); %解三角方程求系數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);endb=-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:m,n+1).2); %計算誤差efprintf('n指數函數的系數是：a=%f,b=%f,c=%f',a,b,c);%輸出系數及誤差efprintf('n使用指數函數擬合的誤差是：%f'

20、,e);t=0:0.05:24;YY=a.*exp(-b.*(t-c).2);figure(2)plot(x,y,'k*',t,YY,'r-'); %繪制擬合指數函數圖形xlabel('時刻'); %標注坐標軸含義ylabel('平均氣溫');title('指數函數的最小二乘曲線');grid;（4）結果與分析a、二次函數：一天的平均氣溫為： 21.20002次函數的系數: 8.3063 2.6064 -0.0938使用2次函數擬合的誤差是：280.339547二次函數的最小二乘曲線如下圖所示：b、三次函數：一天的

21、平均氣溫為： 21.20003次函數的系數: 13.3880 -0.2273 0.2075 -0.0084使用3次函數擬合的誤差是： 131.061822 三次函數的最小二乘曲線如下圖所示：c、四次函數：一天的平均氣溫為： 21.20004次函數的系數: 16.7939 -3.7050 0.8909 -0.0532 0.0009使用4次函數擬合的誤差是：59.04118四次函數的最小二乘曲線如下圖所示：d、指數函數：一天的平均氣溫為： 21.2000指數函數的系數是：a=26.160286，b=0.004442，c=14.081900使用指數函數擬合的誤差是：57.034644指數函數的最小二

22、乘曲線如下圖所示：通過上述幾種擬合可以發現，多項式的次數越高，計算擬合的效果越好，誤差越小，說明結果越準確；同時，指數多項式擬合的次數雖然不高，但誤差最小，說明結果最準確。4.設計算法，求出非線性方程的所有根，并使誤差不超過。（1）算法思想首先，研究函數的形態，確定根的圍；通過剖分區間的方法確定根的位置，然后利用二分法的基本原理進行求解，找到滿足精度要求的解。二分法是產生一串區間，使新區間是舊區間的一個子區間，其長度是的一半，且有一個端點是的一個端點。由區間確定區間的方法是計算區間的中點若，則取，否則取，重復這一過程即可。顯然，每次迭代使區間長度減小一半，故二分法總是收斂的。（2）算法結構1.

23、；2.Ifthen 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; %非線性方程組的表達式g=; for i=-100:1:100 %確定根所在的區間 k=i+1;if(y(x=i).*y(x=k)<eps) %區間長度為1 g=g i;endendsyms x;f=6*x5-45*x2+20;n=length(g)

24、; %確定根的個數 for j=1:n x0=g(j); %求根區間左端點x1=g(j)+1; %求根區間右端點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 %輸出方程的根 end（4）結果與分析該非線性方程組有三個實根，分別為1.8708，0.6812，-0.6545，且滿足誤差要求。5.編寫程序實現大規模方程組的列主元高斯消去法程序，并對所附的方程組進行求解。針對本專業中所碰到的實際問題，提煉一個使用方程組進行求

25、解的例子，并對求解過程進行分析、求解。（1）算法思想高斯消去法是利用現行方程組初等變換中的一種變換，即用一個不為零的數乘一個方程后加只另一個方程，使方程組變成同解的上三角方程組，然后再自下而上對上三角方程組求解。列主元消去法是當高斯消元到第步時，從列的以下（包括）的各元素中選出絕對值最大的，然后通過行交換將其交換到的位置上。交換系數矩陣中的兩行（包括常數項），只相當于兩個方程的位置交換了，因此，列選主元不影響求解的結果。程序的核心就是高斯列主元消去法。根據教材提供的算法，編寫列主元消去法的子函數與適應于超大規模超出系統存的方程組的改編程序。同時，在Gauss消去過程中，適當交換方程的順序對保證

26、消去過程能順利進行及計算解的精確度都是有必要的，交換方程的原則是使中，絕對值最大的一個換到（k,k)位置而成為第k步消去的主元，這就是列主元Gauss消去法。（2） 算法結構1、數據文件的文件名為：文件名+.dat2、數據文件中的數據為二進制記錄結構，分為以下四個部分：(1)文件頭部分，其結構：typedef struct long int id; long int ver; long int n; 其中：id:為該數據文件的標識，值為0xF1E1D1A0，即為：十六進制的F1E1D1A0 ver:為數據文件的版本號，值為16進制數據,版本號說明0x101系數矩陣為非壓縮格式稀疏矩陣0x102

27、系數矩陣為非壓縮格式帶狀對角陣0x201系數矩陣為壓縮格式稀疏矩陣0x202系數矩陣為壓縮格式帶狀對角陣n:表示方程的階數(2)文件頭2:此部分說明為條狀矩陣的上下帶寬,結構：typedef struct long int q; / 為上帶寬 long int p; / 為下帶寬 (3)系數矩陣 a.如存貯格式非為壓縮方式，則按行方式存貯系數矩陣中的每一個元素，個數為n*n，類型為float型；b.如果存貯格式是壓縮方式,則按行方式存貯,每行中只存放上下帶寬的非零元素,即，每行中存貯的最多元素為p+q+1個。(4)右端系數 按順序存貯右端系數的每個元素，個數為n個，類型為float型3、二進制

28、文件的讀取：f=fopen('fun003.dat','r'); %打開文件，.dat文件放在 m 文件同一目錄下，a=fread(f,3,'uint') %讀取頭文件，3-讀取前 3 個，若讀取壓縮格式的，頭文件為 5 個b=fread(f,inf,'float'); %讀取剩下的文件，float 型id=dec2hex(a(1);ver=dec2hex(a(2); %這兩句是進行進制轉換，讀取 id 與ver1. A的階數2. For 2.1找滿足2.2For 2.2.1 2.32.4 For 2.4.1 2.4.2 For

29、2.4.2.1 2.4.3 For （3）Matlab源程序clear; %清除工作空間變量clc; %清除命令窗口命令 % 讀取系數矩陣f,p=uigetfile('*.dat','選擇數據文件'); %讀取數據文件num=5; %輸入系數矩陣文件頭的個數name=strcat(p,f);file=fopen(name,'r');head=fread(file,num,'uint'); %讀取二進制頭文件id=dec2hex(head(1); %讀取標識符fprintf('文件標識符為');idver=dec2h

30、ex(head(2); %讀取版本號fprintf('文件版本號為');vern=head(3); %讀取階數fprintf('矩陣A的階數');nq=head(4); %上帶寬fprintf('矩陣A的上帶寬');q p=head(5); %下帶寬fprintf('矩陣A的下帶寬');p dist=4*num;fseek(file,dist,'bof'); %把句柄值轉向第六個元素開頭處A,count=fread(file,inf,'float'); %讀取二進制文件，獲取系數矩陣fclose(f

31、ile); %關閉二進制頭文件% 對非壓縮帶狀矩陣進行求解if ver='102', a=zeros(n,n);for i=1:n,for j=1:n, a(i,j)=A(i-1)*n+j); %求系數矩陣a(i,j)endend b=zeros(n,1);for i=1:n, b(i)=A(n*n+i);endfor k=1:n-1, %列主元高斯消去法 m=k;for i=k+1:n, %尋找主元if abs(a(m,k)<abs(a(i,k) m=i;endendif a(m,k)=0 %遇到條件終止 disp('錯誤！')returnendfor

32、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;endfor i=k+1:n, %計算l(i,k)并將其放到a(i,k)中 a(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);end b(i)=b(i)-a(i,k)*b(k);endend 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

33、(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)endend b=zeros(n,1);for i=1:1:n b(i)=A(n*m+i); %求b(i)endfor k=1:1:(n-1) %開始消去過程if a(k,(p+1)=0 disp('錯誤！');break;end st1=n;if (k+p)<n st1=k+p;endfor i=(k+1):1:st1

34、 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);endend 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;endfor j=(k+1):1:st2 sum=sum-a(k,j+p-k+1)*x(j);end x(k)=sum/a(k,p+1); sum=0