西安交大計算方法b大作業(yè)課件

1、計算方法B上機實驗報告學院： 機械工程學院 班級: 姓名： 學號： 2015年12月22日1.計算以下和式：，要求：(1)若保留11個有效數(shù)字，給出計算結果，并評價計算的算法；(2)若要保留30個有效數(shù)字，則又將如何進行計算。實現(xiàn)思想：以上問題出現(xiàn)了近似數(shù)相減的問題，為了減小誤差，可分別求得減數(shù)之和以及被減數(shù)之和，最后將兩者相減。另外，減數(shù)與被減數(shù)求和均為同號計算，按照絕對值遞增順序相加可減小舍入誤差。此題中對有效數(shù)字有要求，因而計算時首先需要根據(jù)有效數(shù)字位數(shù)計算得出迭代次數(shù)，以保證計算值的精度。源程序：m=input('輸入有效數(shù)字個數(shù)m=');s0=1;s1=0;s2=0;

2、n=0;%判斷迭代次數(shù)while s0>=0.5*10-(m-1)s0=4/(16n*(8*n+1)-2/(16n*(8*n+4)-1/(16n*(8*n+5)-1/(16n*(8*n+6); n=n+1;end%分別求解各項并求和for k=n-1:-1:0 a1=4/(16k*(8*k+1); a2=2/(16k*(8*k+4); a3=1/(16k*(8*k+5); a4=1/(16k*(8*k+6); s1=a1+s1; s2=a4+a3+a2+s2;endS=vpa(s1-s2,m)實驗結果：11位有效數(shù)字計算結果如圖1所示；30為有效數(shù)字計算結果如圖2所示。 圖1.11位有效

3、數(shù)字計算結果圖2.30為有效數(shù)字計算結果1. 某通信公司在一次施工中，需要在水面寬度為20米的河溝底部沿直線走向鋪設一條溝底光纜。在鋪設光纜之前需要對溝底的地形進行初步探測，從而估計所需光纜的長度，為工程預算提供依據(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 (1)請用合適的曲線擬合所測數(shù)據(jù)點；(2)估算所

4、需光纜長度的近似值，并作出鋪設河底光纜的曲線圖；算法思想：由于題中所給點數(shù)為20，若采用高次多項式插值將產生很大的誤差，所以拉格朗日或牛頓并不適用。題中光纜為柔性，可光滑鋪設于水底，鑒于此特性，采用三次樣條插值插值法較為合適。算法結構：三次樣條算法結構見計算方法教程P110；光纜長度計算公式：源程序：clear;clc; 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.22 9.15 7.90 7.95 8.86 9.81 10.80 10.93; d=y; plot(x,

5、y,'k.','markersize',15) hold on %計算差商 for k=1:2 for i=21:-1:(k+1) d(i)=(d(i)-d(i-1)/(x(i)-x(i-k); end end%設定d的邊界條件 for i=2:20 d(i)=6*d(i+1); end d(1)=0; d(21)=0; %帶狀矩陣求解（追趕法）a=0.5*ones(1,21);b=2*ones(1,21);c=0.5*ones(1,21);a(1)=0;c(21)=0;u=ones(1,21);u(1)=b(1);r=c;yy(1)=d(1); %追 for

6、k=2:21 l(k)=a(k)/u(k-1); u(k)=b(k)-l(k)*r(k-1); yy(k)=d(k)-l(k)*yy(k-1); end %趕 m(21)=yy(21)/u(21); for k=20:-1:1 m(k)=(yy(k)-r(k)*m(k+1)/u(k); end %繪制曲線 k=1; nn=100; xx=linspace(0,20,nn); l=0; for j=1:nn for i=2:20 if xx(j)<=x(i) k=i; break; else k=i+1; endendh=1;xbar=x(k)-xx(j);xmao=xx(j)-x(k-1

7、);s(j)=(m(k-1)*xbar3/6+m(k)*xmao3/6+(y(k-1)-m(k-1)*h2/6)*xbar+(y(k)-m(k)*h2/6)*xmao)/h;sp(j)=-m(k-1)*(x(k)-xx(j)2/(2*h)+m(k)*(xx(j)-x(k-1)2/(2*h)+(y(k)-y(k-1)/h-(m(k)-m(k-1)*h/6;l(j+1)=(1+sp(j)2)0.5*(20/nn)+l(j); %求解光纜長度end%繪圖plot(xx,s,'r-','linewidth',1.5)disp('¹光纜長度為ª

8、',num2str(l(nn+1),'Ã×')曲線圖如圖2-1所示，計算光纜長度如圖2-2所示。圖2-1光纜插值曲線圖圖2-1光纜計算長度顯示3.假定某天的氣溫變化記錄如下表所示，試用數(shù)據(jù)擬合的方法找出這一天的氣溫變化的規(guī)律;試計算這一天的平均氣溫,并試估計誤差。時刻0123456789101112平均氣溫15141414141516182020232528時刻131415161718192021222324平均氣溫313431292725242220181716實現(xiàn)思想：此題中所給數(shù)據(jù)點數(shù)目較多，采用拉格朗日插值法或者牛頓插值法需要很高次的多項式，

9、計算困難，誤差大；采用樣條插值計算量雖然不大，但是存放參數(shù)Mi的量很大，且沒有一個統(tǒng)一的數(shù)學公式來表示，也不是很方便。所以可考慮用最小二乘法進行擬合。計算過程中，分別使用二次函數(shù)、三次函數(shù)以及四次函數(shù)，計算其相應的系數(shù)，估算誤差并作圖比較各個函數(shù)之間的區(qū)別。算法結構：（參考課本P123）1.1形成矩陣Qk1.2變換Gk-1到Gk2.求解三角方程3.計算誤差源代碼：clear;clc; x=0:24; y=15 14 14 14 14 15 16 18 20 20 23

10、0;25 28 31 34 31 29 27 25 24 22 20 18 17 16; m=length(x); n=input('請輸入函數(shù)的次數(shù) '); plot(x,y,'k.',x,y,'-') grid; hold on; n=n+1; G=zeros(m,n+1); G(:,n+1)=y' 

11、;c=zeros(1,n);%建立c來存放 q=0; f=0; b=zeros(1,m);%建立b用來存放 %形成矩陣G for j=1:n     for i=1:m         G(i,j)=x(1,i)(j-1);     end end %建立矩陣Qk for k=1:n &#

12、160;   for i=k:m     c(k)=G(i,k)2+c(k);     end     c(k)=-sign(G(k,k)*(c(k)0.5);     w(k)=G(k,k)-c(k);%建立w來存放     for j=k+1:m    &

13、#160;   w(j)=G(j,k);     end     b(k)=c(k)*w(k);  %變換矩陣Gk-1到Gk     G(k,k)=c(k);      for j=k+1:n+1          q=0;

14、60;        for i=k:m             q=w(i)*G(i,j)+q;         end         s=q/b(k);   

15、60;              for i=k:m         G(i,j)=s*w(i)+G(i,j);         end      end end %求解三角方程 

16、Rx=h1 a(n)=G(n,n+1)/G(n,n);     for i=n-1:(-1):1         for j=i+1:n             f=G(i,j)*a(j)+f;         

17、;end         a(i)=(G(i,n+1)-f)/G(i,i); %a（i）存放各級系數(shù)         f=0;     end     a %回代過程     p=zeros(1,m);    

18、 for j=1:m         for i=1:n         p(j)=p(j)+a(i)*x(j)(i-1);         end     end     plot(x,p,'

19、;r*',x,p,'-');     E2=0;%用E2來存放誤差 %誤差求解 for i=n+1:m     E2=G(i,n+1)2+E2; end E2=E20.5; disp('誤差為'); disp(E2); t=0;for i=1:m     t=t+p(i); end t=t/m;&#

20、160;%平均溫度 disp('平均溫度為',num2str(t),'')實驗結果：二次函數(shù)擬合，結果如下圖所示圖3-1 二次函數(shù)擬合結果三次函數(shù)擬合，結果如下圖所示圖3-2 三次函數(shù)擬合結果四次函數(shù)擬合，結果如下圖所示圖3-3 四次函數(shù)擬合結果結果對比：將二次函數(shù)、三次函數(shù)和四次函數(shù)擬合結果繪制在同一個坐標內，如圖3-4所示。其計算誤差結果見表3-1所示。圖3-4 擬合結果對比分析4.設計算法，求出非線性方程的所有實根，并使誤差不超過。算法思想：本題可采用牛頓法迭代求解,令 ,得帶格式為根據(jù)函數(shù)圖像可以找出根的大致分布區(qū)間,帶入不同的初值即可解出不同

21、的根.源代碼:function y=f2(x)y=6*x.5-45*x.2+20;%定義原函數(shù)function y=f3(x)y=30*x4-90*x;%定義原函數(shù)倒數(shù)i=-5:0.1:5;y=f2(i);plot(i,y)hold onplot(i,0,'-')%畫出原函數(shù)圖像%Newton法求根x1=input('輸入初值');e=10(-4);%誤差設定Nmax=1000;%迭代最大次數(shù)限定for n=1:Nmax f0=f2(x1); if abs(f2(x1)<e fprintf('輸出的f(x)已經(jīng)足夠小'); x=x1; br

22、eak else F0=f3(x1); x=x1-f0/F0; if abs(x-x1)<e break else x1=x; end endendfprintf('輸出方程的根x=%2f',x) 計算結果:函數(shù)圖像如圖4-1所示。計算結果分別見圖4-2所示。圖4-1 函數(shù)圖像圖4-2 計算結果根據(jù)帶入不同的初值，可以求出不同的根，有圖4-2可以看出，原函數(shù)的根大約有三個，分別是-0.654542、0.681174、1.870799。5.線性方程組求解。（1）編寫程序實現(xiàn)大規(guī)模方程組的高斯消去法程序，并對所附的方程組進行求解。所附方程組的類型為對角占優(yōu)的帶狀方程組。（2）

23、針對本專業(yè)中所碰到的實際問題，提煉一個使用方程組進行求解的例子，并對求解過程進行分析、求解。算法思想:高斯消去法是利用現(xiàn)行方程組初等變換中的一種變換,將一個不為零的數(shù)乘到一個方程后加到另一個方程,使方程組變成同解的上三角方程組,然后再自下而上對上三角方程組求解。算法結構：1. 讀取二進制文件，存入計算矩陣2. 對矩陣進行初等變換，然后求解(詳見計算方法教程第2版高斯消去法以及列主元高斯消去法算法)源代碼：clear; clc; % 讀取系數(shù)矩陣f,p=uigetfile('*.dat','選擇數(shù)據(jù)文件'); %讀取數(shù)據(jù)文件num=5; %輸入系數(shù)矩陣文件頭的個數(shù)

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

25、%下帶寬fprintf('矩陣A的下帶寬');p dist=4*num;fseek(file,dist,'bof'); %把句柄值轉向第六個元素開頭處A,count=fread(file,inf,'float'); %讀取二進制文件，獲取系數(shù)矩陣fclose(file); %關閉二進制頭文件% 對非壓縮帶狀矩陣進行求解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

26、 i=1:n, b(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,

27、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); 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); e

28、nd 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)