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

1、計(jì)算方法B上機(jī)實(shí)驗(yàn)報(bào)告學(xué)院： 機(jī)械工程學(xué)院 班級(jí): 姓名： 學(xué)號(hào)： 2015年12月22日1.計(jì)算以下和式：，要求：(1)若保留11個(gè)有效數(shù)字，給出計(jì)算結(jié)果，并評(píng)價(jià)計(jì)算的算法；(2)若要保留30個(gè)有效數(shù)字，則又將如何進(jìn)行計(jì)算。實(shí)現(xiàn)思想：以上問(wèn)題出現(xiàn)了近似數(shù)相減的問(wèn)題，為了減小誤差，可分別求得減數(shù)之和以及被減數(shù)之和，最后將兩者相減。另外，減數(shù)與被減數(shù)求和均為同號(hào)計(jì)算，按照絕對(duì)值遞增順序相加可減小舍入誤差。此題中對(duì)有效數(shù)字有要求，因而計(jì)算時(shí)首先需要根據(jù)有效數(shù)字位數(shù)計(jì)算得出迭代次數(shù)，以保證計(jì)算值的精度。源程序：m=input('輸入有效數(shù)字個(gè)數(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%分別求解各項(xiàng)并求和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)實(shí)驗(yàn)結(jié)果：11位有效數(shù)字計(jì)算結(jié)果如圖1所示；30為有效數(shù)字計(jì)算結(jié)果如圖2所示。 圖1.11位有效

3、數(shù)字計(jì)算結(jié)果圖2.30為有效數(shù)字計(jì)算結(jié)果1. 某通信公司在一次施工中，需要在水面寬度為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)估算所

4、需光纜長(zhǎng)度的近似值，并作出鋪設(shè)河底光纜的曲線圖；算法思想：由于題中所給點(diǎn)數(shù)為20，若采用高次多項(xiàng)式插值將產(chǎn)生很大的誤差，所以拉格朗日或牛頓并不適用。題中光纜為柔性，可光滑鋪設(shè)于水底，鑒于此特性，采用三次樣條插值插值法較為合適。算法結(jié)構(gòu)：三次樣條算法結(jié)構(gòu)見(jiàn)計(jì)算方法教程P110；光纜長(zhǎng)度計(jì)算公式：源程序：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 %計(jì)算差商 for k=1:2 for i=21:-1:(k+1) d(i)=(d(i)-d(i-1)/(x(i)-x(i-k); end end%設(shè)定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); %求解光纜長(zhǎng)度end%繪圖plot(xx,s,'r-','linewidth',1.5)disp('¹光纜長(zhǎng)度為ª

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

9、計(jì)算困難，誤差大；采用樣條插值計(jì)算量雖然不大，但是存放參數(shù)Mi的量很大，且沒(méi)有一個(gè)統(tǒng)一的數(shù)學(xué)公式來(lái)表示，也不是很方便。所以可考慮用最小二乘法進(jìn)行擬合。計(jì)算過(guò)程中，分別使用二次函數(shù)、三次函數(shù)以及四次函數(shù)，計(jì)算其相應(yīng)的系數(shù)，估算誤差并作圖比較各個(gè)函數(shù)之間的區(qū)別。算法結(jié)構(gòu)：（參考課本P123）1.1形成矩陣Qk1.2變換Gk-1到Gk2.求解三角方程3.計(jì)算誤差源代碼：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('請(qǐng)輸入函數(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來(lái)存放 q=0; f=0; b=zeros(1,m);%建立b用來(lái)存放 %形成矩陣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來(lái)存放     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）存放各級(jí)系數(shù)         f=0;     end     a %回代過(guò)程     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來(lái)存放誤差 %誤差求解 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í)驗(yàn)結(jié)果：二次函數(shù)擬合，結(jié)果如下圖所示圖3-1 二次函數(shù)擬合結(jié)果三次函數(shù)擬合，結(jié)果如下圖所示圖3-2 三次函數(shù)擬合結(jié)果四次函數(shù)擬合，結(jié)果如下圖所示圖3-3 四次函數(shù)擬合結(jié)果結(jié)果對(duì)比：將二次函數(shù)、三次函數(shù)和四次函數(shù)擬合結(jié)果繪制在同一個(gè)坐標(biāo)內(nèi)，如圖3-4所示。其計(jì)算誤差結(jié)果見(jiàn)表3-1所示。圖3-4 擬合結(jié)果對(duì)比分析4.設(shè)計(jì)算法，求出非線性方程的所有實(shí)根，并使誤差不超過(guò)。算法思想：本題可采用牛頓法迭代求解,令 ,得帶格式為根據(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,'-')%畫(huà)出原函數(shù)圖像%Newton法求根x1=input('輸入初值');e=10(-4);%誤差設(shè)定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) 計(jì)算結(jié)果:函數(shù)圖像如圖4-1所示。計(jì)算結(jié)果分別見(jiàn)圖4-2所示。圖4-1 函數(shù)圖像圖4-2 計(jì)算結(jié)果根據(jù)帶入不同的初值，可以求出不同的根，有圖4-2可以看出，原函數(shù)的根大約有三個(gè)，分別是-0.654542、0.681174、1.870799。5.線性方程組求解。（1）編寫(xiě)程序?qū)崿F(xiàn)大規(guī)模方程組的高斯消去法程序，并對(duì)所附的方程組進(jìn)行求解。所附方程組的類型為對(duì)角占優(yōu)的帶狀方程組。（2）

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

24、name=strcat(p,f);file=fopen(name,'r');head=fread(file,num,'uint'); %讀取二進(jìn)制頭文件id=dec2hex(head(1); %讀取標(biāo)識(shí)符fprintf('文件標(biāo)識(shí)符為');idver=dec2hex(head(2); %讀取版本號(hào)fprintf('文件版本號(hào)為');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'); %把句柄值轉(zhuǎn)向第六個(gè)元素開(kāi)頭處A,count=fread(file,inf,'float'); %讀取二進(jìn)制文件，獲取系數(shù)矩陣fclose(file); %關(guān)閉二進(jìn)制頭文件% 對(duì)非壓縮帶狀矩陣進(jì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

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('錯(cuò)誤！') 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, %計(jì)算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); %回代過(guò)程 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% 對(duì)壓縮帶狀矩陣進(jìn)行求解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) %開(kāi)始消去 if a(k,(p+1)=0 disp('錯(cuò)誤！'); 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)