計(jì)算方法上機(jī)作業(yè)

1、手及交通丈李XrAM JIAOTOFK UNIVERSITY計(jì)算方法上機(jī)報(bào)告姓 名:學(xué) 號(hào):班 級(jí):上課班級(jí)：.下載可編輯.說(shuō)明:MATLAB本次上機(jī)實(shí)驗(yàn)使用的編程語(yǔ)言是Matlab語(yǔ)言，編譯環(huán)境為7.11.0 ,運(yùn)行平臺(tái)為 Windows 7。1 4211一，、1. 對(duì)以下和式計(jì)算：0 凝8n18n48n58n6 , 要求： 若只需保留11個(gè)有效數(shù)字，該如何進(jìn)行計(jì)算；若要保留30個(gè)有效數(shù)字，則又將如何進(jìn)行計(jì)算；(1)算法思想1、根據(jù)精度要求估計(jì)所加的項(xiàng)數(shù)，可以使用后驗(yàn)誤差估計(jì)，通項(xiàng)為:1421114an -n -n-16 n 8n 18n 4 8n 5 8n 616n 8n 12、為了保證計(jì)

2、算結(jié)果的準(zhǔn)確性，寫程序時(shí)，從后向前計(jì)算;3、使用Matlab時(shí)，可以使用以下函數(shù)控制位數(shù)：digits(位數(shù))或vpa(變量，精度為數(shù))(2)算法結(jié)構(gòu)1. s 0;1421116n 8n 1 8n 4 8n 5 8n 6if2. for n 0,1 , 2, i10end;3. for n i,i 1,i 2, ,0s s t;.下載可編輯(3) Matlab源程序clear;clc;m=input('請(qǐng)輸入有效數(shù)字的位數(shù)m=');s=0;for n=0:50t=(1/16An)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6);if t<

3、=10A(-m) break ;end end;fprintf('需要將 n值加到 n=%dn' ,n-1);for i=n-1:-1:0t=(1/16Ai)*(4/(8*i+1)-2/(8*i+4)-1/(8*i+5)-1/(8*i+6);s=s+t;ends=vpa(s,m)(4)結(jié)果與分析當(dāng)保留11位有效數(shù)字時(shí)，需要將n值加到s =3.1415926536當(dāng)保留30位有效數(shù)字時(shí)，需要將n值加到s =3.14159265358979323846264338328通過(guò)上面的實(shí)驗(yàn)結(jié)果可以看出，通過(guò)從后往前計(jì)算， 果要求保留的準(zhǔn)確數(shù)字位數(shù)的要求。2.某通信公司在一次施工中，需要在

4、水面寬度為% 青除工作空間變量% 青除命令窗口命令% 俞入有效數(shù)字的位數(shù)喏U斷通項(xiàng)與精度的關(guān)系%!要將n值加到的數(shù)值%求和運(yùn)算%$制s的精度n=7,；n=22,o這種算法很好的保證了計(jì)算結(jié)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.15

5、7.907.958.869.8110.8010.93 請(qǐng)用合適的曲線擬合所測(cè)數(shù)據(jù)點(diǎn)； 預(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ì)光滑，緊貼地面，并且忽略水流對(duì)光纜

6、的沖擊。海底光纜線的長(zhǎng)度預(yù)測(cè)模型如下所示，光纜從A點(diǎn)鋪至B點(diǎn)，在某點(diǎn)處的深度為hi 0計(jì)算光纜長(zhǎng)度時(shí)，用如下公式:20L 0 f(x)ds_'2f (x) dx019 k 1f(x),1 f'(x)2dxy)2(2)算法結(jié)構(gòu)1.For0,1,2,2.3.4.5.6.7.8.9.1.1For2.1For4.14.24.3doVFor2.1.1x0hiMi1,2i n, n1, ,k(Mi Mi 1)/(xixi k) Mi1,2,1xi6Mi 1Modb。； n1,d1,n-1diMn；an；2Ci ;1 G a;2 bbn獲取M的矩陣元素個(gè)數(shù)，存入mFor8.18.28.3m/

7、k 2,3, ,mak / k 1 lkbk -lk Ck 1dk - lk k 1m Mm10. For k m 1,m 2, ,110.1 ( k ck Mk i)/ k Mk11. 獲取x的元素個(gè)數(shù)存入s12. 1 k13. For i 1,2,s 113.1 if Xi then i k; breakelse i 1 kh一八;Xkxx;XXk1?X3x3h2h2Mk* Mk否(yk1也1百丈(yk川石闋加(3) Matlab源程序clear;clc;x=0:1:20;爐生從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

8、,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);%求三次樣條插值函數(shù)s(x)M=y;for k=2:3;for i=n:-1:k;M(i)=(M(i)-M(i-1)/(x(i)-x(i-k+1);endendh(1)=x(2)-x(1);for 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;

9、M(n)=0;嗡分點(diǎn)位置的深度數(shù)據(jù)%等分點(diǎn)的數(shù)目%十算二階差商并存放在W%十算三對(duì)角陣系數(shù)a,b,c及右端向量d處擇自然邊界條件b(1)=2；b(n)=2；c(1)=0；a(n)=0;d(1)=0;d(n)=0;u(1)=b(1)； y1(1)=d(1);for k=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)； 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

10、i=2:n-1if X(m)<=x(i);k=i-1;break ;else%寸三對(duì)角陣進(jìn)行LU分解%追趕法求解樣條參數(shù)M(i)k=i; endendH=x(k+1)-x(k);%B各區(qū)間用三次樣條插值函數(shù)計(jì)算X自處的值x1=x(k+1)-X(m);x2=X(m)-x(k);s(m)=(M(k)*(x1A3)/6+M(k+1)*(x2A3)/6+(y(k)-(M(k)*(HA2)/6)*x1+(y(k+1)-(M(k+1)*(HA2)/6)*x2)/H;end%計(jì)算所需光纜長(zhǎng)度L=0;for i=2:N% 十算所需光纜長(zhǎng)度L=L+sqrt(X(i)-X(i-1)A2+(s(i)-s(i-

11、1)A2); end disp('所需光纜長(zhǎng)度為L(zhǎng)=');disp(L);figure plot(x,y, '*' ,X,s,'-') 線圖%繪制鋪設(shè)河底光纜的曲癡注坐標(biāo)軸含義xlabel('位置，'fontsize' ,16);ylabel('深度/m','fontsize' ,16);title('鋪設(shè)河底光纜的曲線圖'，'fontsize' ,16);grid；(4)結(jié)果與分析鋪設(shè)海底光纜的曲線圖如下圖所示：仿真結(jié)果表明，運(yùn)用分段三次樣條插值所得的擬合曲

12、線能較準(zhǔn)確地反映鋪設(shè)光纜的走勢(shì)圖，計(jì)算出所需光纜的長(zhǎng)度為L(zhǎng)=26.4844m。3.假定某天的氣溫變化記錄如下表所示，試用數(shù)據(jù)擬合的方法找出這一天的氣溫變化的規(guī)律；試計(jì)算這一天的平均氣溫，并試估計(jì)誤差。時(shí)刻0123456789101112平均氣溫15141414141516182020232528時(shí)刻131415161718192021222324平均氣溫313431292725242220181716(1)算法思想在本題中，數(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ù)近似，雖然有

13、很好的數(shù)值性質(zhì)，且計(jì)算量也不大，但存放參數(shù)Mi的量很大，且沒(méi)有一個(gè)統(tǒng)一的數(shù)學(xué)公式來(lái)表示，也帶來(lái)了一些不便。另一方面， 在 有的實(shí)際問(wèn)題中，用插值方法并不合適。當(dāng)數(shù)據(jù)點(diǎn)的數(shù)目很大時(shí)，要求P通過(guò)所有數(shù)據(jù)點(diǎn)，可能會(huì)失去原數(shù)據(jù)所表示的規(guī)律。 如果數(shù)據(jù)點(diǎn)是由測(cè)量而來(lái)的，必然帶有誤差， 插值法要求準(zhǔn)確通過(guò)這些不準(zhǔn)確的數(shù)據(jù)點(diǎn)是不合適的。在這種情況下，不用插值標(biāo)準(zhǔn)而用其他近似標(biāo)準(zhǔn)更加合理。通常情況下，是選取i使E2最小，這就是最小二乘近似問(wèn)題。在本題中，采用“最小二乘法”找出這一天的氣溫變化的規(guī)律，使用二次函數(shù)、三、2次函數(shù)、四次函數(shù)以及指數(shù)型函數(shù) C aeb(tc),計(jì)算相應(yīng)的系數(shù)，估算誤差，并作圖比 較各

14、種函數(shù)之間的區(qū)別。(2)算法結(jié)構(gòu)本算法用正交化方法求數(shù)據(jù)的最小二乘近似。假定數(shù)據(jù)以用來(lái)生成了G ,并將y作為其最后一列(第n 1列)存放。結(jié)果在a中， 是誤差E；。I、使用二次函數(shù)、三次函數(shù)、四次函數(shù)擬合時(shí)1 .將“時(shí)刻值”存入Xi ,數(shù)據(jù)點(diǎn)的個(gè)數(shù)存入 m2 .輸入擬合多項(xiàng)式函數(shù)p(x)的最高項(xiàng)次數(shù)h n 1，則擬合多項(xiàng)式函數(shù)為 p(x) igl(x)2g2(x)ngn(x),根據(jù)給定數(shù)據(jù)點(diǎn)確定 GFor j 0,1,2, ,n 1For i 1,2, m2.1 xgi,j12.2 Vgi,n13 .Fork1,2,n3.1 形成矩陣Qk3.1.1 sgn(gkk)(gik)i k3.1.2

15、gkkk3.1.3 For j k 1,k3.1.3.1 gjk3.1.4 k3.2 變換 Gk 1 到Gk3.2.1 gkkFor j k 1,k 2,m3.2.2 ( igj)/ ti k3.2.3 For i k, k 1,3.2.3.1 gjt i4 .解三角方程Rah14 1 gn.n 1 / g nnan5 .2 For i n 1,n 2, ,1n6 2 1gi.n 1 g ij x j / g ii5 .計(jì)算誤差E2m2g i,n 11 n 1II、使用指數(shù)函數(shù)擬合時(shí)現(xiàn)將指數(shù)函數(shù)進(jìn)行變形：2將C y,t X代入C aeb( )得:對(duì)上式左右取對(duì)數(shù)得：2, ,m,n 1,mgja

16、ib(x c)2y ae22ln y lna bc 2bcx bx令 z lny,0 Ina bc2, 1 2bc,則可得多項(xiàng)式:2z 01x2X(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);T=sum(y(1:m)/m;fprintf(' 一天的平土!氣溫為 T=%fn' ,T);%二次、三次、四次函數(shù)的最小二乘近似 h=input

17、('請(qǐng)輸入擬合多項(xiàng)式的最高項(xiàng)次數(shù)=');n=h+1;G=;for j=0:(n-1)g=x，j;G=vertcat(G,g);endG=G'for i=1:mG(i,n+1)=y(i);endG;for k=1:nif G(k,k)>0;sgn=1;elseif G(k,k)=0;sgn=0;else sgn=-1; endsgm=-sgn*sqrt(sum(G(k:m,k).A2);w=zeros(1,n);w(k)=G(k,k)-sgm;for j=k+1:m w(j)=G(j,k);endbt=sgm*w(k);G(k,k尸sgm;for j=k+1:n+1

18、t=sum(w(k:m)*G(k:m,j)/bt;for i=k:m;%各數(shù)據(jù)點(diǎn)的個(gè)數(shù)存入m%求一天的平均氣溫% 艮據(jù)給定數(shù)據(jù)點(diǎn)生成矩陣%g(x)按列排列%亞直連接G%專置得到矩陣G%各數(shù)據(jù)y作為GW最后一列嘛成矢I陣Q(k)射換Gk-1到Gk(n+1 歹U)G(i,j尸G(i,j)+t*w(i);endendendA (n)=G(n,n+1)/G(n,n);for i=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).A2);%由三角方程求系數(shù)Afprintf( disp(A); f

19、printf('%d次函數(shù)的系數(shù)是：，h);%計(jì)算誤差e喻出系數(shù)a及誤差e使用歆函數(shù)擬合的誤差是f: ' ,h,e);t=0:0.05:24;A=fliplr(A);Y=poly2sym(A);subs(Y, 'x' ,t);Y=double(ans);figure。)plot(x,y, 'k*' ,t,Y, 'r-');xlabel('時(shí)刻);ylabel('平均氣溫');title(num2str(n-1),'次函數(shù)的最小二乘曲線');grid;%指數(shù)函數(shù)的最小二乘近似yy=log(y)

20、;n=3;G=；GG=;for j=0:(n-1)g=x，j;G=vertcat(G,g);gg=t.Aj；GG=vertcat(GG,gg);endG=G'for i=1:mG(i,n+1)=yy(i);%各系數(shù)數(shù)組左右翻轉(zhuǎn)%各系數(shù)數(shù)組轉(zhuǎn)化為多項(xiàng)式儂制擬合多項(xiàng)式函數(shù)圖形 兩注坐標(biāo)軸含義%g(x)按列排列%亞直連接G%g(x)按列排列%我直連接G%專置得到矩陣G%各數(shù)據(jù)y作為GW最后一列(n+1歹U)endG;fork=1:n%形成矢I陣Q(k)if G(k,k)>0;sgn=1；elseif G(k,k)=0;sgn=0;else sgn=-1; endsgm=-sgn*sqr

21、t(sum(G(k:m,k2);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;for j=k+1:n+1t=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);for i=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*(cA2)

22、;G(n+1:m,n+1)=exp(sum(G(n+1:m,n+1)，2);e=sum(G(n+1:m,n+1).A2);fprintf( 'n 指數(shù)函數(shù)的系數(shù)是：a=%f,b=%f,c=%f fprintf( 'n使用指數(shù)函數(shù)擬合的誤差是：f ,e);t=0:0.05:24;YY=a.*exp(-b.*(t-c).A2);figure(2)plot(x,y, 'k*' ,t,YY, 'r-');xlabel('時(shí)刻);ylabel('平均氣溫');title('指數(shù)函數(shù)的最小二乘曲線');grid;儂換G

23、k-1到Gk%由三角方程求系數(shù)A%十算誤差e,a,b,c);%俞出系數(shù)及誤差e%繪制擬合指數(shù)函數(shù)圖形%標(biāo)注坐標(biāo)軸含義(4)結(jié)果與分析a、二次函數(shù)：一天的平均氣溫為：21.20002 次函數(shù)的系數(shù)：8.3063 2.6064 -0.0938使用2次函數(shù)擬合的誤差是：280.339547.下載可編輯.二次函數(shù)的最小二乘曲線如下圖所示:2次曲船為品子二暴曲孤20353O2SM佑to 明lr£前b、三次函數(shù):一天的平均氣溫為：21.20003 次函數(shù)的系數(shù)：13.3880 -0.2273 0.2075 -0.0084使用3次函數(shù)擬合的誤差是：131.061822三次函數(shù)的最小二乘曲線如下圖所

24、示：03次甌故的昂小二桀曲第iac、四次函數(shù)：一天的平均氣溫為：21.20004 次函數(shù)的系數(shù)：16.7939 -3.7050 0.8909 -0.0532 0.0009使用4次函數(shù)擬合的誤差是：59.04118四次函數(shù)的最小二乘曲線如下圖所示：.下載可編輯351生 id05ID152025時(shí)荊4次打毀的我匚、乘曲綠5 P2 Zd、指數(shù)函數(shù)：一天的平均氣溫為：21.2000指數(shù)函數(shù)的系數(shù)是：a=26.160286, b=0.004442, c=14.081900使用指數(shù)函數(shù)擬合的誤差是：57.034644指數(shù)函數(shù)的最小二乘曲線如下圖所示：指卻由荊劉易曾二乘曲姓L+*一*，、/7 :/+*LAi

25、*4 y * */F1011110610IS2025時(shí)邦通過(guò)上述幾種擬合可以發(fā)現(xiàn)，多項(xiàng)式的次數(shù)越高，計(jì)算擬合的效果越好，誤差越 小，說(shuō)明結(jié)果越準(zhǔn)確；同時(shí)，指數(shù)多項(xiàng)式擬合的次數(shù)雖然不高，但誤差最小，說(shuō)明結(jié)果 最準(zhǔn)確。524.設(shè)計(jì)算法，求出非線性方程6x 45x 20。的所有根，并使誤差不超過(guò)10(1)算法思想首先，研究函數(shù)的形態(tài)，確定根的范圍；通過(guò)剖分區(qū)間的方法確定根的位置，然后 利用二分法的基本原理進(jìn)行求解，找到滿足精度要求的解。x,x(k °確定區(qū)間I(k1)的二分法是產(chǎn)生一審區(qū)間，使新區(qū)間是舊區(qū)間I (k)的一個(gè)子區(qū)間，其長(zhǎng)度是I的一半，且有一個(gè)端點(diǎn)是I(心的一個(gè)端點(diǎn)。由區(qū)間I(

26、k)方法是計(jì)算區(qū)間I(k)的中點(diǎn)(k 2) x-(x(k) x(k1)2(k) (k 2)(k 1)(k) (k 2)(k 1) (k 2) (k 1)-右f(x )f(x ) 0,則取I x ,x ，畬則取I x ,x ，生旻這也程即可。顯然，每次迭代使區(qū)間長(zhǎng)度減小一半，故二分法總是收斂的。(2)算法結(jié)構(gòu)1. f(x(0)f0； f(x(1)fi2. If f0 f10 then stop3. If|f0|2 then輸出x(0)作為根；stop4. If|f1|2 then輸出x作為根；stop5. 2(x(0) x(1) x6. If |x(1) x |1 | x(1) | then 輸

27、出 x 作為根；stop7. f (x) f8. If |f | 2 then輸出x作為根；9. If f1 f 0 then9.1 xx(0)； ff。else(1)9.2 x x . ff110. go to 5(3) Matlab源程序線性方程組的表達(dá)式%確定根所在的區(qū)間%區(qū)間長(zhǎng)度為1嘛定根的個(gè)數(shù)x=-100:100;y=6*(x.A5)-45*(x.A2)+20;g=；for i=-100:1:100k=i+1;if (y(x=i).*y(x=k)<eps)g=g i； endend syms x;f=6*xA5-45*xA2+20;n=length(g);for j=1:nx0

28、=g(j);%求根區(qū)間左端點(diǎn)x1=g(j)+1；9家根區(qū)間右端點(diǎn)while (x1-x0)>=10A(-4) if subs(f,x,x0)*subs(f,x,(x0+x1)/2)>epsx0=(x0+x1)/2; else x1=(x0+x1)/2; end endroot=x0喻出方程的根end(4)結(jié)果與分析該非線性方程組有三個(gè)實(shí)根，分別為 1.8708, 0.6812 , -0.6545 ,且滿足誤差要求。5.編寫程序?qū)崿F(xiàn)大規(guī)模方程組的列主元高斯消去法程序，并對(duì)所附的方程組進(jìn)行求解。針對(duì)本專業(yè)中所碰到的實(shí)際問(wèn)題，提煉一個(gè)使用方程組進(jìn)行求解的例子，并對(duì)求解過(guò)程進(jìn)行分析、求解。

29、(1)算法思想高斯消去法是利用現(xiàn)行方程組初等變換中的一種變換，即用一個(gè)不為零的數(shù)乘一個(gè)方程后加只另一個(gè)方程，使方程組變成同解的上三角方程組，然后再自下而上對(duì)上 三角方程組求解。列主元消去法是當(dāng)高斯消元到第 k步時(shí)，從k列的akk以下(包括akk)的各元素中選出絕對(duì)值最大的，然后通過(guò)行交換將其交換到akk的位置上。交換系數(shù)矩陣中的兩行(包括常數(shù)項(xiàng))，只相當(dāng)于兩個(gè)方程的位置交換了，因此，列選主元不影響求解的結(jié)果。程序的核心就是高斯列主元消去法。根據(jù)教材提供的算法，編寫列主元消去法的子函數(shù)與適應(yīng)于超大規(guī)模超出系統(tǒng)內(nèi)存的方程組的改編程序。同時(shí)，在Gaus0肖去過(guò)程中，適當(dāng)交換方程的順序?qū)ΡＷC消去過(guò)程能

30、順利進(jìn)行及計(jì)算解的精確度都是有必要的，交換方程的原則是使ikk1)(ik,k 1, ,n)中，絕對(duì)值最大的一個(gè)換到(k,k)位置而 成為第k步消去的主元，這就是列主元 Gauss消去法。(2)算法結(jié)構(gòu)1、數(shù)據(jù)文件的文件名為：文件名+.dat2、數(shù)據(jù)文件中的數(shù)據(jù)為二進(jìn)制記錄結(jié)構(gòu)，分為以下四個(gè)部分：(1)文件頭部分，其結(jié)構(gòu)：typedef structlong int id;long int ver;long int n;其中：id:為該數(shù)據(jù)文件的標(biāo)識(shí)，值為 0xF1E1D1A0即為：十六進(jìn)制的F1E1D1A0ver:為數(shù)據(jù)文件的版本號(hào)，值為16進(jìn)制數(shù)據(jù)，版本號(hào)說(shuō)明0x101系數(shù)矩陣為非壓縮格式稀

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

32、進(jìn)制文件的讀?。篺=fopen('fun003.dat','r'); %打開(kāi)文件，.dat文件放在m文件同一目錄下,a=fread(f,3,'uint') %讀取頭文件，3-讀取前3個(gè)，若讀取壓縮格式的，頭文件為5個(gè)b=fread(f,inf,'float'); %讀取剩下白文件，float 型id=dec2hex(a(1);ver=dec2hex(a(2); %這兩句是進(jìn)行進(jìn)制轉(zhuǎn)換，讀取 id 與ver1. A的階數(shù) n2. For k 1,2,3, ,n 12.1找滿足kkmaxik的下標(biāo)k2.2 For j 1,2, ,n2

33、.2.1kjk j2.3 k2.4 Fori k 1,k 2, , n2.4.1ikikkk2.4.2 For j k 1,k 2, , n2.4.2.1ij ik kjij2.4.3i - ik k in/nnxnFor k n 1, n 2, ,1n( kkj xj )/ kk xk(3) Matlab源程序clear;clc;%讀取系數(shù)矩陣f,p=uigetfile( '*.dat','選擇數(shù)據(jù)文件);num=5;name=strcat(p,f);file=fopen(name, 'r');head=fread(file,num, 'uin

34、t');id=dec2hex(head(1);fprintf('文件標(biāo)識(shí)符為)； idver=dec2hex(head(2);fprintf('文件版本號(hào)為)；vern=head(3);fprintf('矩陣A勺階數(shù))； nq=head(4);fprintf('矩陣A勺上帶寬)；qp=head(5);fprintf('矩陣A勺下帶寬)；pdist=4*num;fseek(file,dist,'bof );頭處A,count=fread(file,inf, 'float');矩陣fclose(file);%對(duì)非壓縮帶狀矩陣進(jìn)

35、行求解if ver= '102', a=zeros(n,n);for i=1:n,for j=1:n,a(i,j尸A(i-1)*n+j);endendb=zeros(n,1);for i=1:n, b(i)=A(n*n+i);end% 青除工作空間變量% 青除命令窗口命令% 賣取數(shù)據(jù)文件% 俞入系數(shù)矩陣文件頭的個(gè)數(shù)取二進(jìn)制頭文件%取標(biāo)識(shí)符% 賣取版本號(hào)%賣取階數(shù)%上帶寬%r帶寬叫巴句柄值轉(zhuǎn)向第六個(gè)元素開(kāi)取二進(jìn)制文件，獲取系數(shù)%關(guān)閉二進(jìn)制頭文件%求系數(shù)矩陣a(i,j)%主元高斯消去法for k=1:n-1,m=k;for i=k+1:n, if abs(a(m,k)<abs

36、(a(i,k) m=i;endendif a(m,k)=0disp( '錯(cuò)誤！) returnendfor 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,中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);endendx=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)

37、*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:nfor j=1:1:ma(i,j)=A(i-1)*m+j); end end b=zeros(n,1);for i=1:1:n%尋找主元%遇到條件終止眇換元素位置得主元%十算l(i,k)并將其放到a(i,k)刎代過(guò)程喻斯消去法麻 a(i,j).下載可編輯麻 b(i)%開(kāi)始消去過(guò)程b(i)=A(n*m+i);endfor k=1:1:(n-1)if a(k,(p+1)=0disp('錯(cuò)誤！');break ;endst1=n;if (k+p)<n st1=k+p;endfor i=(k+1):1:st1a(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); endb(i)=b(i)-a(i,k+p-i+1)*b(k);endend%回代過(guò)程x=zeros(n,1);x