計算方法上機(jī)課后復(fù)習(xí)

計算方法上機(jī)報告

姓名:

學(xué)號:

班級:

上課班級:

說明：

本次上機(jī)實(shí)驗(yàn)使用的編程語言是Matlab語言，編譯環(huán)境為MATLAB

7.11.0,運(yùn)行平臺為Windows7。

30],42]]、

1.對以下和式計算：s=?16718〃+】―8〃+4_8〃+5_8〃+6J，要求：

①若只需保留11個有效數(shù)字，該如何進(jìn)行計算；

②若要保留30個有效數(shù)字，則又將如何進(jìn)行計算；

(1)算法思想

1、根據(jù)精度要求估計所加的項(xiàng)數(shù)，可以使用后驗(yàn)誤差估計，通項(xiàng)為：

1(42114

a?=--------------------------------------------------<-------------<￡；

〃8+14n4851+8J6”川毋8

2、為了保證計算結(jié)果的準(zhǔn)確性，寫程序時，從后向前計算；

3、使用Matlab時，可以使用以下函數(shù)控制位數(shù)：

digits(位數(shù))或vpa(變量，精度為數(shù))

(2)算法結(jié)構(gòu)

1.5=0；

1(42___111

"一方18>+1-8〃+4-8〃-5-8〃+6『

2.for〃=0,1,2,??，

ift<W-m

end;

3.for〃=

s=s+/;

(3)Matlab源程序

clear;%清除工作空間變量

cic;%清除命令窗口命令

m=inputC請輸入有效數(shù)字的位數(shù)m=1）;%輸入有效數(shù)字的位數(shù)

s=0;

forn=0:50

t=(l/16An)*(4/(8*n+l)-2/(8*n+4)-l/(8*n+5)-l/(8*n+6));

ift<=10A(-m)%判斷通項(xiàng)與精度的關(guān)系

break;

end

end;

fprintf('需要將n值加到n=%d\n；n-1);%需要將n值加到的數(shù)值

fori=n-l:-l:0

t=(l/16Ai)*(4/(8*i+l)-2/(8*i+4)-l/(8*i+5)-l/(8*i+6));

s=s+t;%求和運(yùn)算

end

s=vpa(s,m)%控制s的精度

（4）結(jié)果與分析

當(dāng)保留11位有效數(shù)字時，需要將n值加到n=7,

s=3.1415926536；

當(dāng)保留30位有效數(shù)字時，需要將n值加到n=22,

s=3.14159265358979323846264338328。

通過上面的實(shí)驗(yàn)結(jié)果可以看出，通過從后往前計算，這種算法很好的保證了計算結(jié)

果要求保留的準(zhǔn)確數(shù)字位數(shù)的要求。

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.8C10.93

①請用合適的曲線擬合所測數(shù)據(jù)點(diǎn)；

②預(yù)測所需光纜長度的近似值，作出鋪設(shè)河底光纜的曲線圖；

(1)算法思想

如果使用多項(xiàng)式差值，則由于龍格現(xiàn)象，誤差較大，因此，用相對較少的插值數(shù)據(jù)

點(diǎn)作插值，可以避免大的誤差，但是如果又希望將所得數(shù)據(jù)點(diǎn)都用上，且所用數(shù)據(jù)點(diǎn)越

多越好，可以采用分段插值方式，即用分段多項(xiàng)式代替單個多項(xiàng)式作插值。分段多項(xiàng)式

是由一些在相互連接的區(qū)間上的不同多項(xiàng)式連接而成的一條連續(xù)曲線，其中三次樣條插

值方法是一種具有較好“光滑性”的分段插值方法。

在本題中，假設(shè)所鋪設(shè)的光纜足夠柔軟，在鋪設(shè)過程中光纜觸地走勢光滑，緊貼地

面，并且忽略水流對光纜的沖擊。海底光纜線的長度預(yù)測模型如下所示，光纜從A點(diǎn)鋪

至B點(diǎn)，在某點(diǎn)處的深度為

海底光纜線的長度預(yù)測模型

計算光纜長度時，用如下公式:

20

Jof(<x)ds

=x￡+l

=Z，3)2+(Ay)2

(2)算法結(jié)構(gòu)

1.Fori=0,1,2,-??,n

i.iy=K

2.Fork=1,2

2.1Fori=n,n-\,--,k

2.1.1(%-〃-)/-=M

3.=4

4.Fori=1,2,--,n-l

4.1%+i-玉n%

4.2hj+l/(//,.+%)=>q;1-q=>%2nb

4.36A/.+1=>&

5.dQ=>Af0;dn=>Mn;c0

2nbo;//〃=>cin;2=>blt

6.b、=>M，4=>7]

7.獲取M的矩陣元素個數(shù)，存入m

8.For%=2,3,…，機(jī)

8.1aj心二lk

8.2bk-lk-%n網(wǎng)

8.34-小九_]=九

9.7mMm=M,n

10.Fork=

10.1(九-q—)/4nM.

11.獲取x的元素個數(shù)存入s

12.In%

13.ForZ=

13.1ifx<Xjthenink；break

elsei+1nZ

14.演一4_[=>//；xk-x=>x\x-x^=>x

r3X3h2h2

[MATT+%T+(”T-M-)x+(y-M—)x]/h^y

666k6k

(3)Matlab源程序

clear;

cic;

x=0:l:20;%產(chǎn)生從0到20含21個等分點(diǎn)的數(shù)組

X=0:0.2:20;

y=[9.01r8.96,7.96,7.97,8.02,9.05,10.13,11.18,12.26,13.28,1332,12.61.11.29,10.22,9.157.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;%計算二階差商并存放在M中

fori=n:-l:k;

M(i)=(M(i)-M(i-l))/(x(i)-x(i-k+l));

end

end

h(l)=x(2)-x(l);%計算三對角陣系數(shù)a,b,c及右端向量d

fori=2:n-l;

h(i)=x(i+l)-x(i);

c(i)=h(i)/(h(i)+h(i-l));

a(i)=l-c(i);

b(i)=2;

d(i)=6*M(i+l);

end

M(l)=0;%選擇自然邊界條件

M(n)=0;

b(l)=2;

b(n)=2;

c(l)=O;

a(n)=O;

d(l)=O;

d(n)=O;

u(l)=b(l);%對三對角陣進(jìn)行LU分解

yl(D=d(l)；

fork=2:n;

l(k)=a(k)/u(k-l);

u(k)=b(k)-l(k)*c(k-l);

yl(k)=d(k)-l(k)*yl(k-l);

end

M(n)=yl(n)/u(n);%追趕法求解樣條參數(shù)M(i)

fork=n-l:-l:l;

M(k)=(yl(k)-c(k)*M(k+l))/u(k);

end

s=zeros(l,N);

form=l:N;

k=l;

fori=2:n-l

ifX(m)<=x(i);

k=i-l;

break;

else

k=i;

end

end

H=x(k+l)-x(k);%在各區(qū)間用三次樣條插值函數(shù)計算X點(diǎn)處的值

xl=x(k+l)-X(m);

x2=X(m)-x(k);

s(m)=(M(k)*(xlA3)/6+M(k+l)*(x2A3)/6+(y(k)-(M(k)*(HA2)/6))*xl+(y(k+l)-(M(k+l)*(HA2)/

6))*x2)/H;

end

%%計算所需光纜長度

L=0;%計算所需光纜長度

fori=2:N

L=L+sqrt((X(i)-X(i-l))A2+(s(i)-s(i-l))A2);

end

disp('所需光纜長度為L=');

disp(L);

figure

plot(x,y,,*,,X,s；-')%繪制鋪設(shè)河底光纜的曲線圖

xlabel('位置'「fontsize',16);%標(biāo)注坐標(biāo)軸含義

ylabel（'深度/m','fontsize',16）;

title（鋪設(shè)河底光纜的曲線圖；fontsize',16）;

grid;

（4）結(jié)果與分析

鋪設(shè)海底光纜的曲線圖如下圖所示:

仿真結(jié)果表明，運(yùn)用分段三次樣條插值所得的擬合曲線能較準(zhǔn)確地反映鋪設(shè)光纜的

走勢圖，計算出所需光纜的長度為L=26.4844m。

3.假定某天的氣溫變化記錄如下表所示，試用數(shù)據(jù)擬合的方法找出這一天

的氣溫變化的規(guī)律；試計算這一天的平均氣溫，并試估計誤差。

時刻0123456789101112

平均氣溫15141414141516182020232528

時刻131415161718192021222324

平均氣溫313431292725242220181716

(1)算法思想

在本題中，數(shù)據(jù)點(diǎn)的數(shù)目較多。當(dāng)數(shù)據(jù)點(diǎn)的數(shù)目很多時，用“多項(xiàng)式插值”方法做

數(shù)據(jù)近似要用較高次的多項(xiàng)式，這不僅給計算帶來困難，更主要的缺點(diǎn)是誤差很大。用

“插值樣條函數(shù)”做數(shù)據(jù)近似，雖然有很好的數(shù)值性質(zhì)，且計算量也不大，但存放參數(shù)

M的量很大，且沒有一個統(tǒng)一的數(shù)學(xué)公式來表示，也帶來了一些不便。另一方面，在有

的實(shí)際問題中，用插值方法并不合適。當(dāng)數(shù)據(jù)點(diǎn)的數(shù)目很大時，要求P(x)通過所有數(shù)據(jù)

點(diǎn)，可能會失去原數(shù)據(jù)所表示的規(guī)律。如果數(shù)據(jù)點(diǎn)是由測量而來的，必然帶有誤差，插

值法要求準(zhǔn)確通過這些不準(zhǔn)確的數(shù)據(jù)點(diǎn)是不合適的。在這種情況下，不用插值標(biāo)準(zhǔn)而用

其他近似標(biāo)準(zhǔn)更加合理。通常情況下，是選取巴使E2最小，這就是最小二乘近似問題。

在本題中，采用“最小二乘法”找出這一天的氣溫變化的規(guī)律，使用二次函數(shù)、三

次函數(shù)、四次函數(shù)以及指數(shù)型函數(shù)c=四十.解，計算相應(yīng)的系數(shù)，估算誤差，并作圖比

較各種函數(shù)之間的區(qū)別。

(2)算法結(jié)構(gòu)

本算法用正交化方法求數(shù)據(jù)的最小二乘近似。假定數(shù)據(jù)以用來生成了G，并將

)，作為其最后一列(第〃+i列)存放。結(jié)果在々中，夕是誤差七〉

I、使用二次函數(shù)、三次函數(shù)、四次函數(shù)擬合時

1.將“時刻值”存入數(shù)據(jù)點(diǎn)的個數(shù)存入加

2.輸入擬合多項(xiàng)式函數(shù)P&)的最高項(xiàng)次數(shù)〃=1,則擬合多項(xiàng)式函數(shù)為

p(x)=(x)+a2g2(x)+???+angn(x),根據(jù)給定數(shù)據(jù)點(diǎn)確定G

Forj=0,1,2,-?H-1

Fori=1,2,-sm

2.1x/ngjj+i

2.2y,=*gj,〃+[

3.Fork=1,2,-??,n

3.1[形成矩陣Q]

3.1.1-sgn(gQ(Z履嚴(yán)

3.1.2gyon①￡

3.1.3Forj=k-\-\,k+2,--,in

3.1.3.1gjk=>coj

3.1.4oG)k=>p

3.2［變換Gj到GJ

321ongkk

ForJ=k+l,k+2,???,〃,〃+l

“I

3.2.20①8)/=i

i=k

3.2.3Fori=k,k+l,…,m

3.2.3.1gij+igngij

4.「解三角方程網(wǎng)=41

4］gtt.n+Jgnn=

4.2Fori=n—1,H—2,--,l

1［4.”+i-Zg/j1/g"=>ai

，乙?1J=i+1

5.［計算誤差E；］

>：gJ〃+I=>P

IK使用指數(shù)函數(shù)擬合時

現(xiàn)將指數(shù)函數(shù)進(jìn)行變形：

符C=yJ=%代入C=四加與得：y=ae-b(x-c)

對上式左右取對數(shù)得:Iny=In6(-fee2+Ibex-bx1

2

令z=lny,a()=]na-bc,a、=2bc,a2=-b

2

貝！可得多項(xiàng)式：z=4+a”+a2x

(3)Matlab源程序

clear;%清除工作空間變量

cic;%清除命令窗口命令

x=0:24;%將時刻值存入數(shù)組

y=[15,14,14,14,14,15,16,18,20,20,23,25,28,3134,31,29,27,25,24,22,20,18,17,16];

[-,m]=size(x);%將數(shù)據(jù)點(diǎn)的個數(shù)存入m

T=sum(y(l:m))/m;

fprintf。/的平均氣溫為T=%f\n',T);%求一天的平均氣溫

%%二次、三次、四次函數(shù)的最小二乘近似

h二input（'請輸入擬合多項(xiàng)式的最高項(xiàng)次數(shù)二’）；%根據(jù)給定數(shù)據(jù)點(diǎn)生成矩陣G

n=h+l；

G=[]；

forj=0:(n-l)

g=x.Aj；%g（x）按列排列

G=vertcat(G,g);%9垂直連接G

end

G=G';%轉(zhuǎn)置得到矩陣G

fori=l:m%將數(shù)據(jù)y作為G的最后一列（n+1列）

G(i,n+l)=y(i);

end

G;

fork=l:n%形成矩陣Q(k)

ifG(k,k)>0;

sgn=l;

elseifG(k,k)==0;

sgn=0;

elsesgn=-l;

end

sgm=-sgn*sqrt(sum(G(k:m,k).A2));

w=zeros(l,n);

w(k)=G(k,k)-sgm;

forj=k+l:m

w(j)=G(j,k);

end

bt=sgm*w(k);

G(k,k)=sgm;%變換Gk-1到Gk

forj=k+l:n+l

t=sum(w(k:m)*G(k:mj))/bt;

fori=k:m;

G(ij)=G(ij)+t*w(i);

end

end

end

A(n)=G(n,n+l)/G(n,n);%解三角方程求系數(shù)A

fori=n-l:-l:l

A(i)=(G(i,n+l)-sum(G(i,i+l:n).*A(i+l:n)))/G(i,i);

end

e=sum(G(n+l:m,n+l).A2);%計算誤差e

fprintf('%d次函數(shù)的系數(shù)是:；h);%輸出系數(shù)a及誤差e

disp(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(l)

plot(x,y,'k*',t,Y,%繪制擬合多項(xiàng)式函數(shù)圖形

xlabel(時刻');%標(biāo)注坐標(biāo)軸含義

ylabel('平均氣溫')；

title([num2str(n-l)二次函數(shù)的最小二乘曲線']);

grid;

%%指數(shù)函數(shù)的最小二乘近似

yy=log(y);

n=3;

G=[];

GG=[];

forj=O:(n-l)

g=x.Aj；%g(x)按列排列

G=vertcat(G,g);%9垂直連接G

gg=t.Aj；%g(x)按列排列

GG=vertcat(GG,gg);%9垂直連接G

end

G=G';%轉(zhuǎn)置得到矩陣G

fori=l:m%將數(shù)據(jù)y作為G的最后一列(n+1列)

G(i,n+l)=yy(i);

end

G;

fork=l:n%形成矩陣Q(k)

ifG(k,k)>0;

sgn=l;

elseifG(k,k)==0;

sgn=0;

elsesgn=-l;

end

sgm=-sgn*sqrt(sum(G(k:m,k).A2));

w=zeros(l,n);

w(k)=G(k,k)-sgm;

forj=k+l:m

w(j)=G(j,k);

end

bt=sgm*w(k);

G(k,k)=sgm;%變換Gk-1到Gk

forj=k+l:n+l

t=sum(w(k:m)*G(k:mj))/bt;

fori=k:m;

G(ij)=G(i,j)+t*w(i)；

end

end

end

%解三角方程求系數(shù)

A(n)=G(n,n+l)/G(n,n);A

fori=n-l:-l:l

A(i)=(G(i/n+l)-sum(G(ij+l:n).*A(i+l:n)))/G(i,i);

end

b=-A⑶；

c=A(2)/(2*b);

a=exp(A(l)+b*(cA2))；

A

G(n+l:m/n+l)=exp(sum(G(n+l:m,n+l).2));

e=sum(G(n+l:m,n+l).A2);%計算誤差e

fprintf('\n指數(shù)函數(shù)的系數(shù)是:a=%f,b=%f,c=%f1,afb,c);%輸出系數(shù)及誤差e

fprintf（'\n使用指數(shù)函數(shù)擬合的誤差是：%f1,e）;

t=0:0.05:24;

YY=a.*exp(-b.*(t-c).A2);

figure⑵

plot(x,y/k*',t,YY,'r-');%繪制擬合指數(shù)函數(shù)圖形

xlabel。時亥!T);%標(biāo)注坐標(biāo)軸含義

ylabel('平均氣溫');

title（「指數(shù)函數(shù)的最小二乘曲線】）;

grid;

（4）結(jié)果與分析

a、二次函數(shù):

一天的平均氣溫為:21.2000

2次函數(shù)的系數(shù):8.30632.6064-0.0938

使用2次函數(shù)擬合的誤差是：280.339547

二次函數(shù)的最小二乘曲線如下圖所示:

2次函數(shù)的最小二乘曲線

35

30

25

r喟

處20

H-

10

10152025

時刻

b、三次函數(shù)：

一天的平均氣溫為:21.2000

3次函數(shù)的系數(shù)：13.3880-0.22730.2075-0.0084

使用3次函數(shù)擬合的誤差是：131.061822

三次函數(shù)的最小二乘曲線如下圖所示：

3次困數(shù)的最小二乘曲線

時刻

C、四次函數(shù)：

一天的平均氣溫為：21.2000

4次函數(shù)的系數(shù)：16.7939-3.70500.8909-0.05320.0009

使用4次函數(shù)擬合的誤差是：59.04118

四次函數(shù)的最小二乘曲線如下圖所示：

4次函數(shù)的最小二乘曲線

35

30

25

20

15-

1°0

10152025

時刻

d、指數(shù)函數(shù)：

一天的平均氣溫為:21.2000

指數(shù)函數(shù)的系數(shù)是:a=26.160286,b=0.004442,c=14.081900

使用指數(shù)函數(shù)擬合的誤差是：57.034644

指數(shù)函數(shù)的最小二乘曲線如下圖所示:

指數(shù)困數(shù)的最小二乘曲線

35

0510152025

時刻

通過上述幾種擬合可以發(fā)現(xiàn)，多項(xiàng)式的次數(shù)越高，計算擬合的效果越好，誤差越

小，說明結(jié)果越準(zhǔn)確；同時，指數(shù)多項(xiàng)式擬合的次數(shù)雖然不高，但誤差最小，說明結(jié)果

最準(zhǔn)確。

4.設(shè)計算法，求出非線性方程6爐-45/+20=°的所有根，并使誤差不超過

KF4。

（1）算法思想

首先，研究函數(shù)的形態(tài)，確定根的范圍；通過剖分區(qū)間的方法確定根的位置，然后

利用二分法的基本原理進(jìn)行求解，找到滿足精度要求的解。

二分法是產(chǎn)生一串區(qū)間，使新區(qū)間產(chǎn)+D是舊區(qū)間/（幻的一個子區(qū)間，其長度是/⑷

的一半，且有一個端點(diǎn)是/㈤的一個端點(diǎn)。由區(qū)間/陰二次叱膽+"］確定區(qū)間/的的方法

是計算區(qū)間八幻的中點(diǎn)

“2）,（”）+—+D）

2

若/（產(chǎn)）。（產(chǎn)坳）＜0,則取/伏劃=口叫公『，否則取產(chǎn)；儼+2"叫，重復(fù)這一過

程即口」。顯然，每次迭代使區(qū)間長度減小一半，故二分法總是收斂的。

（2）算法結(jié)構(gòu)

L"))"；⑴)n力

2.Iffof\>0thenstop

3.31yolv邑then輸出大(。)作為根;stop

4.If|工|<&then輸出產(chǎn)作為根;stop

l(x(0)+x(,))=>x

,2

6.IfIx(])x\<1x⑴Ithen愉出X作為根；stop

7./(%)=/

8.Ifthen輸出X作為根；

9.If/,/<0then

9.1xni/n/o

else

9.2x=>x(iy.f=>f\

10.goto5

(3)Matlab源程序

x=-100:100;

y=6*(x.A5)-45*(x.A2)+20;%非線性方程組的表達(dá)式

g=[];

fori=-100:l:100%確定根所在的區(qū)間

k=i+l;

if(y(x==i).*y(x==k)<eps)%區(qū)間長度為1

g=[gil;

end

end

symsx;

f=6*xA5-45*xA2+20;

n=length(g);%確定根的個數(shù)

forj=l:n

xO=g(j);%求根區(qū)間左端點(diǎn)

xl=g(j)+l;%求根區(qū)間右端點(diǎn)

while(xl-x0)>=10A(-4)

ifsubs(f,x,x0)*subs(f,x,(x0+xl)/2)>eps

x0=(x0+xl)/2;

else

xl=(x0+xl)/2;

end

end

root=xO%輸出方程的根

end

(4)結(jié)果與分析

該非線性方程組有三個實(shí)根，分別為1.8708,0.6812,-0.6545,且滿足誤差要求。

5.編寫程序?qū)崿F(xiàn)大規(guī)模方程組的列主元高斯消去法程序，并對所附的方程

組進(jìn)行求解。針對本專業(yè)中所碰到的實(shí)際問題，提煉一個使用方程組進(jìn)行

求解的例子，并對求解過程進(jìn)行分析、求解。

(1)算法思想

高斯消去法是利用現(xiàn)行方程組初等變換中的一種變換，即用一個不為零的數(shù)乘一個

方程后加只另一個方程，使方程組變成同解的上三角方程組，然后再自下而上對上三角

方程組求解。

列主元消去法是當(dāng)高斯消元到第欠步時，從女列的&以下〔包括火火)的各元素中

選出絕對值最大的，然后通過行交換將其交換到外.的位置上。交換系數(shù)矩陣中的兩行

(包括常數(shù)項(xiàng))，只相當(dāng)于兩個方程的位置交換了，因此，列選主元不影響求解的結(jié)果。

程序的核心就是高斯列主元消去法。根據(jù)教材提供的算法，編寫列主元消去法的子

函數(shù)與適應(yīng)于超大規(guī)模超出系統(tǒng)內(nèi)存的方程組的改編程序。同時，在Gauss消去過程中，

適當(dāng)交換方程的順序?qū)ΡＷC消去過程能順利進(jìn)行及計算解的精確度都是有必要的，交換

方程的原則是使a十&+1,…中，絕對值最大的一個換到(k,k)位置而成為第k

步消去的主元，這就是列主元Gauss消去法。

(2)算法結(jié)構(gòu)

1、數(shù)據(jù)文件的文件名為：文件名+.dat

2、數(shù)據(jù)文件中的數(shù)據(jù)為二進(jìn)制記錄結(jié)構(gòu)，分為以下四個部分：

(1)文件頭部分，其結(jié)構(gòu)：

typedefstruct

(

longintid;

longintver;

longintn;

)

其中：id:為該數(shù)據(jù)文件的標(biāo)識，值為OxFIElDlAO,即為：十六進(jìn)制的F1E1D1A0

ver:為數(shù)據(jù)文件的版本號，值為16進(jìn)制數(shù)據(jù)，

版本號說明

0x101系數(shù)矩陣為非壓縮格式稀疏矩陣

0x102系數(shù)矩陣為非壓縮格式帶狀對角陣

0x201系數(shù)矩陣為壓縮格式稀疏矩陣

0x202系數(shù)矩陣為壓縮格式帶狀對角陣

n:表示方程的階數(shù)

(2)文件頭2:此部分說明為條狀矩陣的上下帶寬，結(jié)構(gòu)：

typedefstruct

(

longintq;//為上帶寬

longintp;//為下帶寬

J

(3)系數(shù)矩陣

a.如存貯格式非為壓縮方式，則按行方式存貯系數(shù)矩陣中的每一個元素，個數(shù)為

n*n,類型為float型；

b.如果存貯格式是壓縮方式，則按行方式存貯，每行中只存放上下帶寬內(nèi)的非零元素,

即，每行中存貯的最多元素為p+q+1個。

⑷右端系數(shù)

按順序存貯右端系數(shù)的每個元素，個數(shù)為n個，類型為float型

3、二進(jìn)制文件的讀?。?/p>

f=fopen('fun003.dat,,'r,);%打開文件，.dat文件放在m文件同一目錄下，

a=fread(f3；uinf)%讀取頭文件，3-讀取前3個，若讀取壓縮格式的，頭文件為5個

b=fread(f,inf,'float");%讀取剩下的文件，float型

id=dec2hex(a(l));

ver=dec2hex(a(2));%這兩句是進(jìn)行進(jìn)制轉(zhuǎn)換，讀取id與ver

1.A的階數(shù)

2.Fork=1,2,3,-

2.1找滿足以J=i暝的下標(biāo)4

2.2For/=1,2,…,〃

2.2.1akj<=>aftJ

23瓦o然

2.4Fori=A+1,k+2,…,〃

2.4.1—―na法

%

2.4.2Forj=k+l,k+2,…,n

2.4.2.1%-。晨)?=ai)

2.4.3

Fork=〃一一2,…J

n

ax

(4—Zkjj)/。灶nxk

y-i+l

(3)Matlab源程序

clear;%清除工作空間變量

cic;%清除命令窗口命令

%%讀取系數(shù)矩陣

[f,p]=uigetfile('*.dat；選擇數(shù)據(jù)文件')；%讀取數(shù)據(jù)文件

num=5;%輸入系數(shù)矩陣文件頭的個數(shù)

name=strcat(p,f);

file=fopen(name/r');

head=fread(file,num；uint,);%讀取二進(jìn)制頭文件

id=dec2hex(head(l));%讀取標(biāo)識符

fprintf('文件標(biāo)識符為')；

id

ver=dec2hex(head(2));%讀取版本號

fprintf('文件版本號為');

ver

n=head(3);%讀取階數(shù)

fprintf(矩陣A的階數(shù)’)；

n

q=head(4);%上帶寬

fprintf(矩陣A的上帶寬)；

q

p=head(5);%下帶寬

fprintf(矩陣A的下帶寬)；

P

dist=4*num;

fseek(file,dist,'bof');%把句柄值轉(zhuǎn)向第六個元素開頭處

%讀取二進(jìn)制文件，獲取系數(shù)矩陣

[A/count]=fread(file/inf,'float');

fclose(file);%關(guān)閉二進(jìn)制頭文件

%%對非壓縮帶狀矩陣進(jìn)行求解

ifver=='102\

a=zeros(n,n);

fori=l:n,

forj=l:n,

%求系數(shù)矩陣

a(ij)=A((i-l)*n+j);a(i,j)

end

end

b=zeros(n,l);

fori=l:n,

b(i)=A(n*n+i);

end

fork=l:n-l,%列主元高斯消去法

m=k;

fori=k+l:n,%尋找主元

ifabs(a(m,k))<abs(a(i,k))

m=i;

end

end

ifa(m,k)==0%遇到條件終止

dispel!')

return

end

forj=l:n,%交換元素位置得主元

t=a(k,j);

a(kj)=a(m,j);

a(mj)=t;

t=b(k);

b(k)=b(m);

b(m)=t;

end

fori=k+l:n,%計算l(i,k)并將其放至i」a(i,k)中

a(i,k)=a(i,k)/a(k,k);

forj=k+l:n

a(ij)=a(ij)-a(i,k)*a(k,j);

end

b(i)=b(i)-a(i,k)*b(k);

end

end

x=zeros(n,l);%回代過程

x(n)=b(n)/a(n,n);

fork=n-l:-l:l,

x(k)=(b(k)-sum(a(k,k+l:n)*x(k+l：n)))/a(k,k);

end

end

%%對壓縮帶狀矩陣進(jìn)行求解

ifver=='2O2；%高斯消去法

m=p+q+l;

a=zeros(n,m);

fori=l:l:n

forj=l:l:m

a(ij)=A((i-l)*m+j);

end

end

b=zeros(n,l);

fori=l:l:n

b(i)=A(n*m+i);%求可）

end

fork=l:l:(n-l)%開始消去過程

ifa(k,(p+l))==0

disp(㈱吳！1)；

break;

end

stl=n;

if(k+p)<n

stl=k+p;

end

fori=(k+l):l:stl

a(i,(k+p-i+l))=a(i,(k+p-i+l))/a(k,(p+l))；