西安交通大學(xué)計(jì)算方法B上機(jī)報(bào)告

計(jì)算方法上機(jī)報(bào)告姓名：學(xué)號(hào)：班級(jí)：能動(dòng)上課班級(jí)：題目及求解：一、對(duì)下列和式計(jì)算：，規(guī)定：

=1\*GB3①若只需保存11個(gè)有效數(shù)字，該如何進(jìn)行計(jì)算；

=2\*GB3②若要保存30個(gè)有效數(shù)字，則又將如何進(jìn)行計(jì)算；1算法思想（1）根據(jù)精度規(guī)定預(yù)計(jì)所加的項(xiàng)數(shù)，能夠使用后驗(yàn)誤差預(yù)計(jì)，通項(xiàng)為：；（2）為了確保計(jì)算成果的精確性，寫(xiě)程序時(shí)，從后向前計(jì)算；（3）使用Matlab時(shí)，能夠使用下列函數(shù)控制位數(shù)：digits(位數(shù))或vpa(變量，精度為數(shù))2算法構(gòu)造；forifend;for3Matlab源程序clear;%去除工作空間變量clc;%去除命令窗口命令m=input('請(qǐng)輸入有效數(shù)字的位數(shù)m=');%輸入有效數(shù)字的位數(shù)s=0;forn=0:50t=(1/16^n)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6));ift<=10^(-m)%判斷通項(xiàng)與精度的關(guān)系break;endend;fprintf('需要將n值加到n=%d\n',n-1);%需要將n值加到的數(shù)值fori=n-1:-1:0t=(1/16^i)*(4/(8*i+1)-2/(8*i+4)-1/(8*i+5)-1/(8*i+6));s=s+t;%求和運(yùn)算ends=vpa(s,m)%控制s的精度4成果與分析若保存11位有效數(shù)字，則n=7，此時(shí)求解得：s=3.；若保存30位有效數(shù)字時(shí)，則n=22,此時(shí)求解得：s=3.338328。通過(guò)上面的實(shí)驗(yàn)成果能夠看出，通過(guò)從后往前計(jì)算，這種算法較好的確保了計(jì)算成果規(guī)定保存的精確數(shù)字位數(shù)的規(guī)定。二、某通信公司在一次施工中，需要在水面寬度為20米的河溝底部沿直線走向鋪設(shè)一條溝底光纜。在鋪設(shè)光纜之前需要對(duì)溝底的地形進(jìn)行初步探測(cè)，從而預(yù)計(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\*GB3①請(qǐng)用適宜的曲線擬合所測(cè)數(shù)據(jù)點(diǎn)；=2\*GB3②預(yù)測(cè)所需光纜長(zhǎng)度的近似值，作出鋪設(shè)河底光纜的曲線圖；1算法思想如果使用多項(xiàng)式差值，則由于龍格現(xiàn)象，誤差較大，因此，用相對(duì)較少的插值數(shù)據(jù)點(diǎn)作插值，能夠避免大的誤差，但是如果又但愿將所得數(shù)據(jù)點(diǎn)都用上，且所用數(shù)據(jù)點(diǎn)越多越好，能夠采用分段插值方式，即用分段多項(xiàng)式替代單個(gè)多項(xiàng)式作插值。分段多項(xiàng)式是由某些在互相連接的區(qū)間上的不同多項(xiàng)式連接而成的一條持續(xù)曲線，其中三次樣條插值辦法是一種含有較好“光滑性”的分段插值辦法。在本題中，假設(shè)所鋪設(shè)的光纜足夠柔軟，在鋪設(shè)過(guò)程中光纜觸地走勢(shì)光滑，緊貼地面，并且無(wú)視水流對(duì)光纜的沖擊。海底光纜線的長(zhǎng)度預(yù)測(cè)模型如圖2-1所示，光纜從A點(diǎn)鋪至B點(diǎn)，在某點(diǎn)處的深度為。圖2-1海底光纜線的長(zhǎng)度預(yù)測(cè)模型計(jì)算光纜長(zhǎng)度時(shí)，用以下公式：2算法構(gòu)造1）For1.12）For2.1For2.1.13）4）For4.14.24.35）6）7）獲取M的矩陣元素個(gè)數(shù)，存入m8）For8.18.28.39）10）For10.111）獲取x的元素個(gè)數(shù)存入s12）13）For13.1ifthen；breakelse14）3Matlab源程序clear;%去除工作空間變量clc;%去除命令窗口命令x=0:1:20;%產(chǎn)生從0到20含21個(gè)等分點(diǎn)的數(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];%等分點(diǎn)位置的深度數(shù)據(jù)n=length(x);%等分點(diǎn)的數(shù)目N=length(X);%%求三次樣條插值函數(shù)s(x)M=y;fork=2:3;%計(jì)算二階差商并寄存在M中fori=n:-1:k;M(i)=(M(i)-M(i-1))/(x(i)-x(i-k+1));endendh(1)=x(2)-x(1);%計(jì)算三對(duì)角陣系數(shù)a,b,c及右端向量dfori=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);%對(duì)三對(duì)角陣進(jìn)行LU分解y1(1)=d(1);fork=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)fork=n-1:-1:1;M(k)=(y1(k)-c(k)*M(k+1))/u(k);ends=zeros(1,N);form=1:N;k=1;fori=2:n-1ifX(m)<=x(i);k=i-1;break;elsek=i;endendH=x(k+1)-x(k);%在各區(qū)間用三次樣條插值函數(shù)計(jì)算X點(diǎn)處的值x1=x(k+1)-X(m);x2=X(m)-x(k);s(m)=(M(k)*(x1^3)/6+M(k+1)*(x2^3)/6+(y(k)-(M(k)*(H^2)/6))*x1+(y(k+1)-(M(k+1)*(H^2)/6))*x2)/H;end%%計(jì)算所需光纜長(zhǎng)度L=0;%計(jì)算所需光纜長(zhǎng)度f(wàn)ori=2:NL=L+sqrt((X(i)-X(i-1))^2+(s(i)-s(i-1))^2);enddisp('所需光纜長(zhǎng)度為L(zhǎng)=');disp(L);figureplot(x,y,'*',X,s,'-')%繪制鋪設(shè)河底光纜的曲線圖xlabel('位置','fontsize',16);%標(biāo)注坐標(biāo)軸含義ylabel('深度/m','fontsize',16);title('鋪設(shè)河底光纜的曲線圖','fontsize',16);grid;4成果與分析鋪設(shè)海底光纜的曲線圖如圖2-2所示：圖2-2鋪設(shè)河底光纜曲線圖仿真成果表明，運(yùn)用分段三次樣條插值所得的擬合曲線能較精確地反映鋪設(shè)光纜的走勢(shì)圖，計(jì)算出所需光纜的長(zhǎng)度為L(zhǎng)=26.4844m。三、假定某天的氣溫變化統(tǒng)計(jì)以下表所示，試用數(shù)據(jù)擬合的辦法找出這一天的氣溫變化的規(guī)律；試計(jì)算這一天的平均氣溫，并試預(yù)計(jì)誤差。時(shí)刻0123456789101112平均氣溫15141414141516182020232528時(shí)刻131415161718192021222324平均氣溫3134312927252422201817161算法思想在本題中，數(shù)據(jù)點(diǎn)的數(shù)目較多。當(dāng)數(shù)據(jù)點(diǎn)的數(shù)目諸多時(shí)，用“多項(xiàng)式插值”辦法做數(shù)據(jù)近似要用較高次的多項(xiàng)式，這不僅給計(jì)算帶來(lái)困難，更重要的缺點(diǎn)是誤差很大。用“插值樣條函數(shù)”做數(shù)據(jù)近似，即使有較好的數(shù)值性質(zhì)，且計(jì)算量也不大，但寄存參數(shù)的量很大，且沒(méi)有一種統(tǒng)一的數(shù)學(xué)公式來(lái)表達(dá)，也帶來(lái)了某些不便。另首先，在有的實(shí)際問(wèn)題中，用插值辦法并不適宜。當(dāng)數(shù)據(jù)點(diǎn)的數(shù)目很大時(shí)，規(guī)定通過(guò)全部數(shù)據(jù)點(diǎn)，可能會(huì)失去原數(shù)據(jù)所示的規(guī)律。如果數(shù)據(jù)點(diǎn)是由測(cè)量而來(lái)的，必然帶有誤差，插值法規(guī)定精確通過(guò)這些不精確的數(shù)據(jù)點(diǎn)是不適宜的。在這種狀況下，不用插值原則而用其它近似原則更加合理。普通狀況下，是選用使最小，這就是最小二乘近似問(wèn)題。在本題中，采用“最小二乘法”找出這一天的氣溫變化的規(guī)律，使用二次函數(shù)、三次函數(shù)、四次函數(shù)以及指數(shù)型函數(shù)，計(jì)算對(duì)應(yīng)的系數(shù)，估算誤差，并作圖比較多個(gè)函數(shù)之間的區(qū)別。2算法構(gòu)造本算法用正交化辦法求數(shù)據(jù)的最小二乘近似。假定數(shù)據(jù)以用來(lái)生成了，并將作為其最后一列（第列）寄存。成果在中，是誤差。=1\*ROMANI、使用二次函數(shù)、三次函數(shù)、四次函數(shù)擬合時(shí)1. 將“時(shí)刻值”存入，數(shù)據(jù)點(diǎn)的個(gè)數(shù)存入2. 輸入擬合多項(xiàng)式函數(shù)的最高項(xiàng)次數(shù)，則擬合多項(xiàng)式函數(shù)為,根據(jù)給定數(shù)據(jù)點(diǎn)擬定ForFor2.12.23. For3.1[形成矩陣]3.1.13.1.23.1.3For3.1.3.13.1.43.2[變換到]3.2.1For3.2.23.2.3For3.2.3.14.[解三角方程]4.14.2 For4.2.15.[計(jì)算誤差]=2\*ROMANII、使用指數(shù)函數(shù)擬合時(shí)現(xiàn)將指數(shù)函數(shù)進(jìn)行變形：將,代入得:對(duì)上式左右取對(duì)數(shù)得:令則可得多項(xiàng)式:（3）Matlab源程序clear;%去除工作空間變量clc;%去除命令窗口命令x=0:24;%將時(shí)刻值存入數(shù)組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);%將數(shù)據(jù)點(diǎn)的個(gè)數(shù)存入mT=sum(y(1:m))/m;fprintf('一天的平均氣溫為T(mén)=%f\n',T);%求一天的平均氣溫%%二次、三次、四次函數(shù)的最小二乘近似h=input('請(qǐng)輸入擬合多項(xiàng)式的最高項(xiàng)次數(shù)=');%根據(jù)給定數(shù)據(jù)點(diǎn)生成矩陣Gn=h+1;G=[];forj=0:(n-1)g=x.^j;%g(x)按列排列G=vertcat(G,g);%g垂直連接GendG=G';%轉(zhuǎn)置得到矩陣Gfori=1:m%將數(shù)據(jù)y作為G的最后一列（n+1列）G(i,n+1)=y(i);endG;fork=1:n%形成矩陣Q(k)ifG(k,k)>0;sgn=1;elseifG(k,k)==0;sgn=0;elsesgn=-1;endsgm=-sgn*sqrt(sum(G(k:m,k).^2));w=zeros(1,n);w(k)=G(k,k)-sgm;forj=k+1:mw(j)=G(j,k);endbt=sgm*w(k);G(k,k)=sgm;%變換Gk-1到Gkforj=k+1:n+1t=sum(w(k:m)*G(k:m,j))/bt;fori=k:m;G(i,j)=G(i,j)+t*w(i);endendendA(n)=G(n,n+1)/G(n,n);%解三角方程求系數(shù)Afori=n-1:-1:1A(i)=(G(i,n+1)-sum(G(i,i+1:n).*A(i+1:n)))/G(i,i);ende=sum(G(n+1:m,n+1).^2);%計(jì)算誤差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(A);%將系數(shù)數(shù)組轉(zhuǎn)化為多項(xiàng)式subs(Y,'x',t);Y=double(ans);figure(1)plot(x,y,'k*',t,Y,'r-');%繪制擬合多項(xiàng)式函數(shù)圖形xlabel('時(shí)刻');%標(biāo)注坐標(biāo)軸含義ylabel('平均氣溫');title([num2str(n-1),'次函數(shù)的最小二乘曲線']);grid;%%指數(shù)函數(shù)的最小二乘近似yy=log(y);n=3;G=[];GG=[];forj=0:(n-1)g=x.^j;%g(x)按列排列G=vertcat(G,g);%g垂直連接Ggg=t.^j;%g(x)按列排列GG=vertcat(GG,gg);%g垂直連接GendG=G';%轉(zhuǎn)置得到矩陣Gfori=1:m%將數(shù)據(jù)y作為G的最后一列（n+1列）G(i,n+1)=yy(i);endG;fork=1:n%形成矩陣Q(k)ifG(k,k)>0;sgn=1;elseifG(k,k)==0;sgn=0;elsesgn=-1;endsgm=-sgn*sqrt(sum(G(k:m,k).^2));w=zeros(1,n);w(k)=G(k,k)-sgm;forj=k+1:mw(j)=G(j,k);endbt=sgm*w(k);G(k,k)=sgm;%變換Gk-1到Gkforj=k+1:n+1t=sum(w(k:m)*G(k:m,j))/bt;fori=k:m;G(i,j)=G(i,j)+t*w(i);endendendA(n)=G(n,n+1)/G(n,n);%解三角方程求系數(shù)Afori=n-1:-1:1A(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*(c^2));G(n+1:m,n+1)=exp(sum(G(n+1:m,n+1).^2));e=sum(G(n+1:m,n+1).^2);%計(jì)算誤差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('時(shí)刻');%標(biāo)注坐標(biāo)軸含義ylabel('平均氣溫');title(['指數(shù)函數(shù)的最小二乘曲線']);grid;4成果與分析a二次函數(shù)：一天的平均氣溫為：21.二次函數(shù)的系數(shù):8.30632.6064-0.0938使用二次函數(shù)擬合的誤差是：280.339547二次函數(shù)的最小二乘曲線以下圖3-1所示： 圖3-1二次函數(shù)的最小二乘曲線b三次函數(shù)：一天的平均氣溫為：21.三次函數(shù)的系數(shù):13.3880-0.22730.2075-0.0084使用三次函數(shù)擬合的誤差是：131.061822三次函數(shù)的最小二乘曲線如圖3-2所示： 圖3-2三次函數(shù)的最小二乘曲線c四次函數(shù)：一天的平均氣溫為：21.四次函數(shù)的系數(shù):16.7939-3.70500.8909-0.05320.0009使用四次函數(shù)擬合的誤差是：59.04118四次函數(shù)的最小二乘曲線如圖3-3所示：圖3-3四次函數(shù)的最小二乘曲線d指數(shù)函數(shù)：一天的平均氣溫為：21.指數(shù)函數(shù)的系數(shù)是：a=26.160286，b=0.004442，c=14.081900使用指數(shù)函數(shù)擬合的誤差是：57.034644指數(shù)函數(shù)的最小二乘曲線以下圖所示：圖3-4指數(shù)函數(shù)的最小二乘曲線通過(guò)上述幾個(gè)擬合能夠發(fā)現(xiàn)，多項(xiàng)式的次數(shù)越高，計(jì)算擬合的效果越好，誤差越小，闡明成果越精確；同時(shí)，指數(shù)多項(xiàng)式擬合的次數(shù)即使不高，但誤差最小，闡明成果最精確。四、設(shè)計(jì)算法，求出非線性方程的全部根，并使誤差不超出。1算法思想首先，研究函數(shù)的形態(tài)，擬定根的范疇；通過(guò)剖分區(qū)間的辦法擬定根的位置，然后運(yùn)用二分法的基本原理進(jìn)行求解，找到滿足精度規(guī)定的解。二分法是產(chǎn)生一串區(qū)間，使新區(qū)間是舊區(qū)間的一種子區(qū)間，其長(zhǎng)度是的二分之一，且有一種端點(diǎn)是的一種端點(diǎn)。由區(qū)間擬定區(qū)間的辦法是計(jì)算區(qū)間的中點(diǎn)若，則取，否則取，重復(fù)這一過(guò)程即可。顯然，每次迭代使區(qū)間長(zhǎng)度減小二分之一，故二分法總是收斂的。2算法構(gòu)造1） ；2） Ifthenstop3） Ifthen輸出作為根；stop4） Ifthen輸出作為根；stop5） 6） Ifthen輸出作為根；stop7） 8） Ifthen輸出作為根；9） Ifthen9.1；else9.2；10）goto53Matlab源程序clear;%去除工作空間變量clc;%去除命令窗口命令x=-100:100;y=6*(x.^5)-45*(x.^2)+20;%非線性方程組的體現(xiàn)式g=[];fori=-100:1:100%擬定根所在的區(qū)間k=i+1;if(y(x==i).*y(x==k)<eps)%區(qū)間長(zhǎng)度為1g=[gi];endendsymsx;f=6*x^5-45*x^2+20;n=length(g);%擬定根的個(gè)數(shù)forj=1:nx0=g(j);%求根區(qū)間左端點(diǎn)x1=g(j)+1;%求根區(qū)間右端點(diǎn)while(x1-x0)>=10^(-4)ifsubs(f,x,x0)*subs(f,x,(x0+x1)/2)>epsx0=(x0+x1)/2;elsex1=(x0+x1)/2;endendroot=x0%輸出方程的根end4成果與分析該非線性方程組有三個(gè)實(shí)根，分別為1.8708，0.6812，-0.6545，且滿足誤差規(guī)定。五、編寫(xiě)程序?qū)崿F(xiàn)大規(guī)模方程組的列主元高斯消去法程序，并對(duì)所附的方程組進(jìn)行求解。針對(duì)本專業(yè)中所碰到的實(shí)際問(wèn)題，提煉一種使用方程組進(jìn)行求解的例子，并對(duì)求解過(guò)程進(jìn)行分析、求解。1算法思想高斯消去法是運(yùn)用現(xiàn)行方程組初等變換中的一種變換，即用一種不為零的數(shù)乘一種方程后加只另一種方程，使方程組變成同解的上三角方程組，然后再自下而上對(duì)上三角方程組求解。列主元消去法是當(dāng)高斯消元到第步時(shí)，從列的下列（涉及）的各元素中選出絕對(duì)值最大的，然后通過(guò)行交換將其交換到的位置上。交換系數(shù)矩陣中的兩行（涉及常數(shù)項(xiàng)），只相稱于兩個(gè)方程的位置交換了，因此，列選主元不影響求解的成果。程序的核心就是高斯列主元消去法。根據(jù)教材提供的算法，編寫(xiě)列主元消去法的子函數(shù)與適應(yīng)于超大規(guī)模超出系統(tǒng)內(nèi)存的方程組的改編程序。同時(shí)，在Gauss消去過(guò)程中，適宜交換方程的次序?qū)Υ_保消去過(guò)程能順利進(jìn)行及計(jì)算解的精確度都是有必要的，交換方程的原則是使中，絕對(duì)值最大的一種換到（k,k)位置而成為第k步消去的主元，這就是列主元Gauss消去法。2算法構(gòu)造1、數(shù)據(jù)文獻(xiàn)的文獻(xiàn)名為：文獻(xiàn)名+.dat2、數(shù)據(jù)文獻(xiàn)中的數(shù)據(jù)為二進(jìn)制統(tǒng)計(jì)構(gòu)造，分為下列四個(gè)部分：(1)文獻(xiàn)頭部分，其構(gòu)造： typedefstruct{longintid;longintver;longintn;}其中：id:為該數(shù)據(jù)文獻(xiàn)的標(biāo)記，值為0xF1E1D1A0，即為：十六進(jìn)制的F1E1D1A0ver:為數(shù)據(jù)文獻(xiàn)的版本號(hào)，值為16進(jìn)制數(shù)據(jù),版本號(hào)闡明0x101系數(shù)矩陣為非壓縮格式稀疏矩陣0x102系數(shù)矩陣為非壓縮格式帶狀對(duì)角陣0x201系數(shù)矩陣為壓縮格式稀疏矩陣0x202系數(shù)矩陣為壓縮格式帶狀對(duì)角陣n:表達(dá)方程的階數(shù)(2)文獻(xiàn)頭2:此部分闡明為條狀矩陣的上下帶寬,構(gòu)造：typedefstruct{longintq;//為上帶寬longintp;//為下帶寬}(3)系數(shù)矩陣a.如存貯格式非為壓縮方式，則按行方式存貯系數(shù)矩陣中的每一種元素，個(gè)數(shù)為n*n，類型為float型；b.如果存貯格式是壓縮方式,則按行方式存貯,每行中只寄存上下帶寬內(nèi)的非零元素,即，每行中存貯的最多元素為p+q+1個(gè)。(4)右端系數(shù)按次序存貯右端系數(shù)的每個(gè)元素，個(gè)數(shù)為n個(gè)，類型為float型3、二進(jìn)制文獻(xiàn)的讀?。篺=fopen('fun003.dat','r');%打開(kāi)文獻(xiàn)，.dat文獻(xiàn)放在m文獻(xiàn)同一目錄下，a=fread(f,3,'uint')%讀取頭文獻(xiàn)，3-讀取前3個(gè)，若讀取壓縮格式的，頭文獻(xiàn)為5個(gè)b=fread(f,inf,'float');%讀取剩余的文獻(xiàn)，float型id=dec2hex(a(1));ver=dec2hex(a(2));%這兩句是進(jìn)行進(jìn)制轉(zhuǎn)換，讀取id與ver1.A的階數(shù)2.For2.1找滿足2.2 For2.2.12.32.4For2.4.12.4.2For2.4.2.12.4.3For（3）Matlab源程序clear;%去除工作空間變量clc;%去除命令窗口命令%%讀取系數(shù)矩陣[f,p]=uigetfile('*.dat','選擇數(shù)據(jù)文獻(xiàn)');%讀取數(shù)據(jù)文獻(xiàn)num=5;%輸入系數(shù)矩陣文獻(xiàn)頭的個(gè)數(shù)name=strcat(p,f);file=fopen(name,'r');head=fread(file,num,'uint');%讀取二進(jìn)制頭文獻(xiàn)id=dec2hex(head(1));%讀取標(biāo)記符fprintf('文獻(xiàn)標(biāo)記符為');idver=dec2hex(head(2));%讀取版本號(hào)fprintf('文獻(xiàn)版本號(hào)為');vern=head(3);%讀取階數(shù)fprintf('矩陣A的階數(shù)');nq=head(4);%上帶寬fprintf('矩陣A的上帶寬');qp=head(5);%下帶寬fprintf('矩陣A的下帶寬');pdist=4*num;fseek(file,dist,'bof');%把句柄值轉(zhuǎn)向第六個(gè)元素開(kāi)頭處[A,count]=fread(file,inf,'float');%讀取二進(jìn)制文獻(xiàn)，獲取系數(shù)矩陣fclose(file);%關(guān)閉二進(jìn)制頭文獻(xiàn)%%對(duì)非壓縮帶狀矩陣進(jìn)行求解ifver=='102',a=zeros(n,n);fori=1:n,forj=1:n,a(i,j)=A((i-1)*n+j);%求系數(shù)矩陣a(i,j)endendb=zeros(n,1);fori=1:n,b(i)=A(n*n+i);endfork=1:n-1,%列主元高斯消去法m=k;fori=k+1:n,%尋找主元ifabs(a(m,k))<abs(a(i,k))m=i;endendifa(m,k)==0%碰到條件終止disp('錯(cuò)誤！')returnendforj=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;endfori=k+1:n,%計(jì)算l(i,k)并將其放到a(i,k)中a(i,k)=a(i,k)/a(k,k);forj=k+1:na(i,j)=a(i,j)-a(i,k)*a(k,j);endb(i)=b(i)-a(i,k)*b(k);endendx=zeros(n,1);%回代過(guò)程x(n)=b(n)/a(n,n);fork=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)行求解ifver=='202',%高斯消去法m=p+q+1;a=zeros(n,m);fori=1:1:nforj=1:1:ma(i,j)=A((i-1)*m+j);%求a(i,j)endendb=zeros(n,1);fori=1:1:nb(i)=A(n*m+i);%求b(i)endfork=1:1:(n-1)%開(kāi)始消去過(guò)程ifa(k,(p+1))==0disp('錯(cuò)誤！');break;endst1=n;if(k+p)<nst1=k+p;endfori=(k+1):1:st1a(i,(k+p-i+1))=a(i,(k+p-i+1))/a(k,(p+1));forj=(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);endb(i)=b(i)-a(i,k+p-i+1)*b(k);endendx=zeros(n,1);