西安交通大學(xué)計算方法上機(jī)課后復(fù)習(xí)

^算方法上機(jī)作業(yè).對以下和式計算：S—Xf上一上一，一，］，要求:n=016n18n+18n+48n+58〃+6^⑴若只需保留11個有效數(shù)字，該如何進(jìn)行計算；⑵若要保留30個有效數(shù)字，則又將如何進(jìn)行計算；(1)解題思想和算法實(shí)現(xiàn):根據(jù)保留有效位數(shù)的要求，可以由公式：三12E一得出計算精度要求。只需要很少內(nèi)存，時間復(fù)雜度和d呈線性，不需要高浮點(diǎn)支持。先根據(jù)while語句求出符合精度要求的n值的大小，然后利用for語句對這n項進(jìn)行求和，輸出計算結(jié)果及n值大小即可。2)matlab源程序:保留11位有效數(shù)字時;clearclcformatlongn=0;sum=1/(16An)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6));whilesum>=5*10A(-11)；n=n+1;sum=1/(16An)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6));endfori=0:n-1;sum=sum+1/(16Ai)*(4/(8*i+1)-2/(8*i+4)-1/(8*i+5)-1/(8*i+6));endvpa(sum,11)n保留30位有效數(shù)字時；clearclcformatlongn=0;sum=1/(16An)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6));whilesum>=5*10A(-30);n=n+1;sum=1/(16An)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6));endfori=0:n-1;sum=sum+1/(16Ai)*(4/(8*i+1)-2/(8*i+4)-1/(8*i+5)-1/(8*i+6));endvpa(sum,30)n(3)實(shí)驗結(jié)果分析anw=she= 3.L4159265358379323846204r33832S3A?I士?I圖1.1保留11位有效數(shù)字的n值及計算結(jié)果圖 圖1.2保留30位有效數(shù)字的n值及計算結(jié)果圖由計算結(jié)果可知，通過合理的誤差控制，分別通過7次和22次循環(huán)，可以實(shí)現(xiàn)題目所要求的精確度。2.某通信公司在一次施工中，需要在水面寬度為20米的河溝底部沿直線走向鋪設(shè)一條溝底光纜。在鋪設(shè)光纜之前需要對溝底的地形進(jìn)行初步探測，從而估計所需光纜的長度，為工程預(yù)算提供依據(jù)。已探測到一組等分點(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.810.9⑵預(yù)測所需光纜長度的近似值，并作出鋪設(shè)河底光纜的曲線圖;(1)解題思想和算法原理給定區(qū)間［a,b］一個分劃/:a=x0<x1<_<xN=b若函數(shù)S(x)滿足下列條件:S(X)在每個區(qū)間［XyXj］上是不高于3次的多項式。5僅汲其2階導(dǎo)數(shù)在［a,b］上連續(xù)。則稱S(x)使關(guān)于分劃/的三次樣條函數(shù)。(2)matlab源程序:clc,clearx=0:1:20;y=[9.018.967.967.978.029.0510.1311.1812.2613.2813.3212.6111.2910.229.157.97.958.869.8110.8010.93];n=length(x);l(1)=0;m(n)=0;h=diff(x);df=diff(y)./diff(x);d(1)=0;d(n)=0;forj=2:n-1l(j)=h(j)/(h(j-1)+h(j));m(j)=h(j-1)/(h(j-1)+h(j));d(j)=6*(df(j)-df(j-1))/(h(j-1)+h(j));endm=m(2:end);u=diag(m,-1);r=diag(l,1);a=diag(2*ones(1,n));A=u+r+a;M=inv(A)*d';symsgforj=1:n-1s(j)=M(j)*(x(j+1)-g)八3/(6*h(j))+M(j+1)*((g-x(j))八3/(6*h(j)))+(y(j)-M(j)*h(j)八2/6)*(x(j+1)-g)/h(j)+(y(j+1)-M(j+1)*h(j)八2/6)*(g-x(j))/h(j);endsr=0;forj=1:n-1df=diff(s(j),g);warningoffall;q=int(sqrt(1+df.A2),g,j-1,j);r=r+q;endL=vpa(r,8);disp('thelengthofthelabelisL=');disp(L);forj=1:n-1S(j,:)=sym2poly(s(j));endforj=1:n-1x1=x(j):0.1:x(j+1);y1=polyval(S(j,:),x1);ifj==1y2=y1;elsefori=1:11k=(j-1)*10+i;y2(k)=y1(i);endendendx2=x(1):0.1:x(n);plot(x,y,'o')gridholdonplot(x2,y2,'r')(3)實(shí)驗結(jié)果分析

[(3958520385365547*e)/14073748835532800-(186483313085687^"3)/562919953421312■+901/100.(5067619118220011*(g-thelengthafthelabelisL=26.498:514圖2.1鋪設(shè)河底電纜長度圖2.2鋪設(shè)河底光纜的曲線圖由三次樣條插值得出的函數(shù)曲線的長度和即鋪設(shè)河底電纜的長度為26.498514。為了提高插值精度，用三次樣條插值可以增加插值節(jié)點(diǎn)的辦法來滿足要求，而且在給定節(jié)點(diǎn)數(shù)的條件下，三次樣條插值的精度要優(yōu)于多項式插值以及線性分段插值，雖然舍棄了降低誤差這個優(yōu)點(diǎn)，但是其曲線的光滑性要好一些。3.假定某天的氣溫變化記錄如下表所示，試用數(shù)據(jù)擬合的方法找出這一天的氣溫變化的規(guī)律;試計算這一天的平均氣溫，并試估計誤差。時刻012345678911111111222平均氣溫5144144156180203258111122時刻314516718920122324332221平均氣溫134129725422018716(1)解題思想和數(shù)學(xué)原理：對于具體實(shí)驗時，通常不是先給出函數(shù)的解析式，再進(jìn)行實(shí)驗，而是通過實(shí)驗的觀察和測量給出離散的一些點(diǎn)，再來求出具體的函數(shù)解析式。又因為測量誤差的存在，實(shí)際真實(shí)的解析式曲線并不一定通過測量給出的所有點(diǎn)。最小二乘法，形成法方程是求解這一問題的很好的方法，本實(shí)驗運(yùn)用這一方法實(shí)現(xiàn)對給定數(shù)據(jù)的擬合。對于給定的測量數(shù)據(jù)(xfi)(i=12…，n)，設(shè)函數(shù)分布為y(x)=Ea①(x)j=o特別的，取甲j(x)為多項式①C)=xj(j=0,1,...，m)則根據(jù)最小二乘法原理，可以構(gòu)造泛函i=1ji=1(k=0,1,..,m)(f，P(f，P0)一30,50)

(P4)

0:1■(P4)0m求該解方程組，則可以得到解a0,a1，…,am，因此可得到數(shù)據(jù)的最小二乘解f(x)xEa①(x)j=0(2)matlab源程序：x=[0123456789101112131415161718192021222324];%給出題目數(shù)據(jù)(時間)y=[15141414141516182020232528313231292725242220181716]；%給出題目數(shù)據(jù)(溫度)plot(x,y,'m*')%畫出各個離散數(shù)據(jù)點(diǎn)holdonforn=2:4;%2、3、4代表擬合函數(shù)最高項次數(shù)alltemp=25;%alltemp代表數(shù)據(jù)點(diǎn)總共有25個A=zeros(n+1,n+1);%定義初始正規(guī)方程組的系數(shù)矩陣AC=ones(n+1,1);%定義初始正規(guī)方程組的系數(shù)向量CD=zeros(n+1,1);%定義初始正規(guī)方程組的向量Dfori=1:n+1forj=1:n+1fork1=1:alltempA(i,j)=A(i,j)+(x(k1).人(i-1+j-1));endendfork2=1:alltempD(i,1)=D(i,1)+(x(k2).人(i-1)).*(y(k2));endend%以上為計算出正規(guī)方程組矩陣A、D的所有元素的程序tol=1.0e-12;maxit=1000;C=bicg(A,D,tol,maxit);%使用bicg迭代算出正規(guī)方程組的系數(shù)向量Cp=0;%誤差分量E=0;%誤差總量ifn==2b=0:24;f=C(1)+C(2).*b+C(3).*(b.人2);p=y(b+1)-f;forv=1:25E=E+(p(v)).人2;endplot(b,f,'r-')end%以上是對2階擬合函數(shù)的圖形處理與誤差計算ifn==3b=0:24;f=C(1)+C(2).*b+C(3).*(b.人2)+C(4).*(b.人3);p=y(b+1)-f;forv=1:25E=E+(p(v)).人2;endplot(b,f,'g-')end%以上是對3階擬合函數(shù)的圖形處理與誤差計算ifn==4b=0:24;f=C(1)+C(2).*b+C(3).*(b.A2)+C(4).*(b.A3)+C(5).*(b.A4);p=y(b+1)-f;forv=1:25E=E+(p(v)).A2;endplot(b,f,'b-')end%以上是對4階擬合函數(shù)的圖形處理與誤差計算C,Eendn=2;%重新對n賦值，進(jìn)行指數(shù)函數(shù)擬合A=zeros(n+1,n+1);%重新對A矩陣賦初值C=zeros(n+1,1);%重新對C向量賦初值D=zeros(n+1,1);%重新對D向量賦初值fori=1:n+1forj=1:n+1fork=1:alltempA(i,j)=A(i,j)+(x(k).人(i-1+j-1));endendforl=1:alltempD(i,1)=D(i,1)+((x(l).人(i-1)).*(log(y(l))));endend%計算出A矩陣、D向量各元素數(shù)值C=bicg(A,D,tol,maxit);%利用bicg迭代求解系數(shù)b=0:24;p=0;E=0;f=exp(C(1)+C(2).*b+C(3).*(b.人2));p=y(b+1)-f;forv=1:25E=E+(p(v)).人2;endplot(b,f,'c-')%對指數(shù)擬合函數(shù)進(jìn)行圖形處理和誤差計算b=-C⑶；c-C(2)/(2*b);a=exp(b*(c^2)+C(1));%算出題設(shè)要求的指數(shù)擬合函數(shù)的各個系數(shù)a,b,c,Egridonlegend('測量數(shù)據(jù)',‘二次函數(shù)‘,‘三次函數(shù)‘,‘四次函數(shù)’,’指數(shù)擬合,'Location','NorthWest')holdoff%holdon與holdoff聯(lián)合使用，表示將各個曲線畫在同一個圖中8.42742.5606-0.0920E=253.6476圖3.1二次曲線擬合系數(shù)與2范數(shù)誤差C=13.4072-32161CL2032-0.0082E=110.2954圖3.2三次曲線擬合系數(shù)與2范數(shù)誤差-3.55470.S592-0.0E130.0009F=43.9299圖3.3四次曲線擬合系數(shù)與2范數(shù)誤差25.024Tb=0.004414.09139E=191a0771圖3.4指數(shù)曲線擬合系數(shù)與2范數(shù)誤差圖3.5數(shù)據(jù)原始點(diǎn)與擬合曲線對比圖(3)實(shí)驗結(jié)果分析：根據(jù)國家有關(guān)規(guī)定，用的是02、08、14、20時四個觀測時次的數(shù)據(jù)做平均，最有代表性。從圖中可以看出并不是多項式次數(shù)越高越好，隨著次數(shù)的增高，曲線所呈現(xiàn)出的給定點(diǎn)處和實(shí)際的吻合度越好，但對于其他地方的吻合度降低了。4.設(shè)計算法，求出非線性方程6x5_45x2+20-0的所有根，并使誤差不超過10-4。解：(1)解題思想和算法實(shí)現(xiàn)：對于一個非線性方程的數(shù)值解法很多。在此介紹兩種最常見的方法：二分法和Newton法。首先要研究函數(shù)的形態(tài)，確定根的數(shù)量和大致區(qū)間的位置。對于二分法，其數(shù)學(xué)實(shí)質(zhì)就是說對于給定的待求解的方程f(x)，其在［4切上連續(xù)，f(a)f(b)<0,且f(x)在［a,b］內(nèi)僅有一個實(shí)根x*，取區(qū)間中點(diǎn)c，若，則c恰為其根，否則根據(jù)f(a)f(c)<0是否成立判斷根在區(qū)間［a,c］和［c,b］中的哪一個，

從而得出新區(qū)間，仍稱為［a,b］。重復(fù)運(yùn)行計算，直至滿足精度為止。這就是二分法的計算思想。xk+1Newton法通常預(yù)先要給出一個猜測初值xk+1f1(x)

k產(chǎn)生逼近解X*的迭代數(shù)列｛xk｝，這就是Newton法的思想。當(dāng)x0接近x*時收斂很快，但是當(dāng)x0選擇不好時，可能會發(fā)散，因此初值的選取很重要。另外，若將該迭代公式改進(jìn)為__f(x)*k+1 Xkrf-(x)k其中r為要求的方程的根的重數(shù)，這就是改進(jìn)的Newton法，當(dāng)求解已知重數(shù)的方程的根時，在同種條件下其收斂速度要比Newton法快的多。本題中用matlab畫曲線y本題中用matlab畫曲線y=6x5_45x2+20如下:圖4.1y=6*x.A5-45*x.A2+20的曲線由曲線可以看出，該方程有三個實(shí)根，根所在區(qū)間依次為：[-1,-0.5] [0.5,1] [1.5,2]采用Newton迭代法，依據(jù)迭代式：(8-1)(8-2)X(k+1)-X(k)_〃x(k))，k=0,12…(8-1)(8-2)X(K+)—X(KJf'(xck))該方程的迭代式為：6X(k)5一45X(k)2+20X(k+1)=X(k)一 30x(k)4一90x(k)分別選用迭代初值X(0)二一1、x(0)=1、x(0)=2進(jìn)行迭代，分別求得滿足精度的三個實(shí)根。(2)matlab源程序：clc;clearx=-2:0.1:2;y=6*x.人5-45*x.人2+20;plot(x,y,'g')%畫相應(yīng)的曲線gridtitle('y=6*x.人5-45*x.人2+20曲線')formatlongroots([600-45020])%根的真實(shí)解clearx0=input('請輸入迭代初值x0=');formatlong

i=1；e=10人-4;%精度x=x0-(6*x0.人5-45*x0.人2+20)/(30*x0.人4-90*x0);whileabs(x-x0)>ei=i+1;x0=x;x=x0-(6*x0.人5-45*x0.人2+20)/(30*x0.人4-90**0);%迭代enddisplay('方程的根’)display(迭代的次數(shù)’)(3)實(shí)驗結(jié)果分析?Q4_2清輸;<迭代初值?Q4_2清輸;<迭代初值xQ=-L方程的根?Q4_2清輸；(迭代初值如=L方程的根>>Q4_2請輸；C迭代初值知=2方程的根-0.6545423853973500.-0.6545423853973500.6811741073263151.870799017267090迭代的次數(shù)迭代的次數(shù)迭代的次數(shù)迭代的次數(shù)迭代的次數(shù)迭代的次數(shù)圖4.2運(yùn)行結(jié)果對于Newton迭代法，三個初值x0都使得迭代收斂，這是非常重要的。考慮Newton法迭代的收斂性條件：(1)存在一個區(qū)間，滿足。由曲線和所選的三個區(qū)間可知這一條件滿足。（2）f（x）是［a,b］上的單調(diào)函數(shù)，即對一切不變號。經(jīng)驗證所選的三個區(qū)間滿足這一條件。（3）f（x）的凹向在［a,b］上不變，即在［a,b］上不改變符號。經(jīng)驗證所選的三個區(qū)間滿足這一條件。（4）-111.5>0>0>0另外，可以看出Newton迭代法收斂速度也很快，且很快達(dá)到很高的精度，源于它一般是超線性收斂的。5.編寫程序?qū)崿F(xiàn)大規(guī)模方程組的列主元高斯消去法程序，并對所附的方程組進(jìn)行求解。針對本專業(yè)中所碰到的實(shí)際問題，提煉一個使用方程組進(jìn)行求解的例子，并對求解過程進(jìn)行分析、求解。解：（1）算法原理由于一般線性方程在使用Gauss消去法求解時從求解的過程中可以看到，若ak廠=0廁必須進(jìn)行行交換，才能使消去過程進(jìn)行下去。有的時候即使a『）牛0，但是其絕對值非常小，由于機(jī)器舍入誤差的影響，消去過程也會出現(xiàn)不穩(wěn)定得現(xiàn)象，導(dǎo)致結(jié)果不正確。因此有必要進(jìn)行列主元技術(shù)，以最大可能的消除這種現(xiàn)象。這一技術(shù)要尋找行「，使得Ia（I）|二maxa（k-1）ik并將第r行和第k行的元素進(jìn)行交換以使得當(dāng)前的akk廠的數(shù)值比0要大的多。這種列主元的消去法的主要步驟如下:

.消元過程對k=1,2,…,n-1,進(jìn)行如下步驟。ik選主元，記Ia1=maxaikrk i>k'1a3很小，這說明方程的系數(shù)矩陣嚴(yán)重病態(tài)，給出警告，提示結(jié)果可能不對。2）交換增廣陣A的r,k兩行的元素。a-ajj o=k-n+1)3)計算消元jjakakj/,(i=k+1,...,n;j=k+1,……,n+1).回代過程對k=n,n-1,…,1,進(jìn)行如下計算kjjkjj=k-1jkk至此，完成了整個方程組的求解。(2)matlab源程序:%非壓縮,dat51.dat、dat53.datclear;clcfp=fopen('dat53.dat','rb');id=fread(fp,1,'int32');ver=fread(fp,1,'int32');N=fread(fp,1,'int32');q=fread(fp,1,'int32');p=fread(fp,1,'int32');fori=1:NA(i,:)=fread(fp,N,'float');endfori=1:Nd(i)=fread(fp,1,'float');end%正向消去過程fori=1:N-qfork=1:pll=A(i+k,i)/A(i,i);forj=i:i+qA(i+k,j)=A(i+k,j)-ll*A(i,j);endd(i+k)=d(i+k)-ll*d(i);endendfori=N-q+1:Nfork=1:N-ill=A(i+k,i)/A(i,i);forj=i:NA(i+k,j)=A(i+k,j)-ll*A(i,j);endd(i+k)=d(i+k)-ll*d(i);endend%回代過程x(N)=d(N)/A(N,N);fori=N-1:-1:1S=0;ifi+q>Ncv=N;%cv-criticalvalueelsecv=i+q;endforj=i+1:cvS=S+A(i,j)*x(j);endx(i)=(d(i)-S)/A(i,i);endx%壓縮,dat52.dat、dat54.datclear;clcfp=fopen('dat54.dat','rb');id=fread(fp,1,'int32');ver=fread(fp,1,'int32');N=fread(fp,1,'int32');q=fread(fp,1,'int32');p=fread(fp,1,'int32');fori=1:NA(i,:)=fread(fp,p+q+1,'float');endfori=1:Nd(i)=fread(fp,1,'float');end%正向消去過程fori=1:Nifi+p<Ncv=p;%cv-criticalvalueelsecv=N-i;endfork=1:cvll=A(i+k,p+1-k)/A(i,p+1);forj=p+1:p+q+1A(i+k,j-k)=A(i+k,j-k)-ll*A(i,j);endd(i+k)=d(i+k)-ll*d(i);endend%回代過程x(N)=d(N)/A(N,p+l);fori=N-l:-l:l;S=0;ii=i+l;jj=P+2；while(ii<=N)&&(jj<=p+q+l)S=S+A(i,jj)*x(ii);11=11+1;?? ??Tendx(i)=(d(i)-S)/A(i,p+l);endx⑶實(shí)驗結(jié)果：

,CcimmendAa.1400 1100 亂1印制 3,L-LOO2.1400 工180 3.110(1 2IQ兇 XL10Q3.1400 ，“忖 3,110(1 1 L*同Cabjina即而throueti20703.1300 1L4D02.WXi2.L400 2.140。 3.3300J.L4D0 皇1?0 3.L400W.1400 2.3100 2.liDOJ.1300 1L400Cal'Jina307]throuch2OS53.1400 11400 2.l*Xi 3.L400 3.14C-fl 2.3400 3.14噸 2.14u￡i 3.L40Q 第 3.3400 2.140Q 3.1400 3. L4Q0ColiJinz2036thrDUfiti2L0D3.1400 3.1400 2.l*Xi 2.L400 2.1400 3.]JQO 2.l40fl 3L1400 3.L400 2.14<i0 3.3100 2.140Q 3. ]JOO 2. L4D0Cabjin32101through2id33.1440 3. L4DD S.]-d?X' ，LM。 M140。 3.3400 ML40口 2.1400 3.L4DQ 3.1400 ，“口口 3.140口 3. 34QO 3, L4DQ匚乳皿匕2116>±irw曲2L?D3.1400 3. L4DD H.1JCC 3.L-LDO M140。 3.]400 3.L4DH Z.1400 3.L4D0 3.140D 3.]4iD0 3.1400 3. ]4Q0 3. L4D0C■口Ljoms213]>±irD-u￡h2L153m"QQ3,L4D03,1-SM 3,1400工3g3,1400 3,1400 3.L400 3,1400 3?“QQ 3.14003門卻0 3,L4dOCabjins2146throuih2L603”1切0 3,L4Q03.1.iWILMQ3,140。 1Ng3,L400 3.14W 3,L400 3,1400 1HE3,1400 1：NQO 3,L400>1非壓縮矩陣求解結(jié)果(部分)州nd日出'L5-21005.2LOO5?2100 5.2L0D 521005.2LOO52100 5.2100 5.2100 5.2M5.2100 5?2100 5.2100 5?2100Colimis$2141匚kruu曲力[前□-Z100 5.2LOO 5,2100 5.2100 52100 5.2100 521Q0 5.2100 5.21W 5,￡100 5.2100 a.ZIQO 5.2100 5,2100Colimis4131S6匚hrnu狗4J]705-2100 5.2LOO 5.2100 5.2100 52100 5.2LOO 52100 5.2100 5.2100 5,￡100 5.2LOO a-ElOO 5.2100 5.2100Ccdumai1217L匚krnugk門】的a-2100 5.2LOO 5,2100 5.2100 52100 5.2100 521Q0 5.2100 5.21W 5,￡100 5.2100 a.ZIQO 5.2100 5,2100Celuixia121￡I6匚hrnugh422005-2100 5.2LOO 5.2100 5.2100 52100 5.2LOO 52100 5.2100 5.2100 5,￡100 5.2LOO a-ElOO 5.2100 5.2100Ccduma1220L匚krnugk42215a-2100 5.2LOO 5,2100 5.2100 52100 5.2100 521Q0 5.2100 5.21W 5,￡100 5.2100 a.ZIQO 5.2100 5,2100roluma132L6uhrnugti4323D5.2100 5.210D 5u210O 5,2LO