計(jì)算方法上機(jī)報(bào)告

姓名：

學(xué)號(hào)：

班級(jí)：

上課班級(jí)：

說(shuō)明：

本次上機(jī)實(shí)驗(yàn)使用的編程語(yǔ)言是Matlab語(yǔ)言，編譯

環(huán)境為MATLAB7.11.0,運(yùn)行平臺(tái)為Windows7。

1.對(duì)以下和式計(jì)算：

T-8n+4-8n+5-8n+6).、

。，要求：

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

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

(1)算法思想

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

%修(4211\14

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

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

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

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

「凸-1_____?______1______1_\

n

Ls=0;16\8n+l8n+48nf58n+6/;

2.forn=0,12?-JiftW10end;

3.for"二ij—1J-2,???，0§=5+￡；(3)Matlab源程序

clear;%清除工作空間變量

cic;%清除命令窗I」命令

m=input(，請(qǐng)輸入有效數(shù)字的位數(shù)m=');%輸入有效數(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<=10N-m)%判斷通項(xiàng)與精度的關(guān)系

break;

end

end;

fprintf。需要將n值加到n=%d\n',n-l);%需要將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)算

s=vpa（s,m）%控制s的精度

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

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

s=3.1415926536；

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

s=3.14159265358979323846264338328o

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

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

2.某通信公司在一次施工中，需要在水面寬度為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.28133212.6111.2910.22

分點(diǎn)14151617181920

深度9.157.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ì)光纜的沖擊。海底光纜線的長(zhǎng)度預(yù)測(cè)模型如下所示，

光纜從A點(diǎn)鋪至B點(diǎn)，在某點(diǎn)處的深度為

海底光纜線的長(zhǎng)度預(yù)測(cè)模型

計(jì)算光纜長(zhǎng)度時(shí)，用如下公式:

f20=f2°/(x)Vl+/(x)2dx

L=I/(x)ds

=ErO(xWl+f(X)2dx=^V(Ax)2+(Ay)2

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

[Fori=0，12…,〃L1y尸M'2.Fork=l，22.lFor

，二〃,〃-1,…，丘1.1（"|一時(shí)1）/3-*心）="3Xi-Xon/l"

Fori=L2,???/l?14.1XJ+LX產(chǎn)力f+14.2

力Hd/g>6H4）="l-C尸四；2nb436%A%

donM^anM=Co2=>bo甲/%;aw。％."『從/尸入7.

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

8.ForA=2,3.…M&l?！ㄈ?「豆2%/"尸"出3

乙-〃'乙-產(chǎn)丫人9.七/'內(nèi)=”析10.Fork二6一1,小一2,…，110」

（3“此+1）/〃「Mk”獲取X的元素個(gè)數(shù)存入s

一

12.1=A13.Fori=L2,…，5-113.1if*4~theni=k；break

elsei”=k14,Xk-Xjnh；X—nX；XT1=

_3^3

22

XX/>'—，

[M—T+M4仇石）》+以一如可）打/…

（3

Matlab源程序

clear;

x=0:l:201%產(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ù)11

N=length（X）;

%%求三次樣條插值函數(shù)s(x)

M=y;

fork=2:3;%計(jì)算二階差商并存放在M中

fori=n:-l:k;

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

end

end

計(jì)算三對(duì)角陣系數(shù)及右端向量

h(l)=x(2)-x(l);%azb,cd

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)=0;

d(l)=0;

d(n)=0;

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

yi(D=d(i);

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

在各區(qū)間月三次樣條插值函數(shù)計(jì)算點(diǎn)處的值

H=x(k+l)-x(k);%X

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

%%計(jì)算所需光纜長(zhǎng)度

L=0;%計(jì)算所需光纜長(zhǎng)度

fori=2:N

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

end

dispf所需光纜長(zhǎng)度為L(zhǎng)=')；

disp(L);

figure

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

xlabe(位置；fontsize,,16);%標(biāo)注坐標(biāo)軸含義

ylabel('深度/m'「fontsizH16);

title('鋪設(shè)河底光纜的曲線圖］fontsizH16);

grid;

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

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

鋪設(shè)河底光纜的曲線圖

位置

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

纜的走勢(shì)圖，計(jì)算出所需光纜的長(zhǎng)度為L(zhǎng)=26.4844m0

3.假定某天的氣溫變化記錄如下表所示，試用數(shù)據(jù)

擬合的方法找出這一天的氣溫變化的規(guī)律；試計(jì)算這

一天的平均氣溫，并試估計(jì)誤差。.

在本題中，數(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ù)Mi的量很大，且沒有一個(gè)統(tǒng)一的數(shù)學(xué)公式來(lái)表示，也

帶來(lái)了一些不便。另一方面，在有的實(shí)際問(wèn)題中，用插值方法并不合適。當(dāng)數(shù)

據(jù)點(diǎn)的數(shù)目很大時(shí)，要求P（x）通過(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)更

加合理。通常情況下，是選取°|使￡二最小，這就是最小二乘近似問(wèn)題。

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

數(shù)、三次函數(shù)、四次函數(shù)以及指數(shù)型函數(shù)，計(jì)算相應(yīng)的系數(shù)，估

算誤差，并作圖比較各種函數(shù)之間的區(qū)別。

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

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

作為其最后一列（第〃+1列）存放。結(jié)果在。中，P是誤差2。

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

1.將“時(shí)刻值”存入數(shù)據(jù)點(diǎn)的個(gè)數(shù)存入川2.輸入擬合多項(xiàng)式函數(shù)

"（*）的最高項(xiàng)次數(shù)人二〃-1,則擬合多項(xiàng)式函數(shù)為

P（x）=a%（x）+aM（x）+…+%"x）根據(jù)給定數(shù)據(jù)點(diǎn)確定G

For/=0,l,2,--Ji-lFor…，印2.1%=氏42.2"尸氏"13.

Fork=l,2,???M3.1［形成矩陣。打

m

l/2

-sgn(gkk)(^g^)=>o

q,,一on3k

3.1.13.1.2ukkA3.1.3For

/=k+l,k+2,…，刖3.L3.1。小叼3.1.403ks3.2[變換Gi到

G&]

m

03由

3.2.1"ng**For/=+…,+I3.2.21=k

3.2.3For曰，丘1,???,m3.2.3.1。/'"尸況/4.[解三角方程

R。一八勺

4.19n.n+1^9nn^0n4.2For/=H-1,H-2,???,I4.2.1

n

⑸,eNg/R/g，尸5E2

六注15.[計(jì)算誤差2]

m

。：用力]

，=E〃+i”、使用指數(shù)函數(shù)擬合時(shí)

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

將C=y,=x代入得：y=Qe-ba-c『對(duì)上式左右取對(duì)數(shù)得:

Iny=Ina-bc2^2bcx-bx2令

2

z=lny,a0=lna-bcfal=2bc,a2=-b則可得多項(xiàng)式：

2

Z=aQlalX\a2X⑶Matlab源程序

clear;%清除工作空間變量

cic;%清除命令窗口命令

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ù)存入m

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

fprint"一天的平均氣溫為T=%An-J);%求一天的平均氣溫

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

h=input(，請(qǐng)輸入擬合多項(xiàng)式的最高項(xiàng)次數(shù)=);%根據(jù)給定數(shù)據(jù)點(diǎn)生成矩陣G

n=h+l;

G=n；

forj=O:(n-l)

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

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

end

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

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

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(kzk)-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);%計(jì)算誤差e

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

disp(A);

fprintfC使用％d次函數(shù)擬合的誤差是％f：

t=0:0.05:24;

A=fliplr(A);%將系數(shù)數(shù)組左右翻轉(zhuǎn)

Y=poly2sym(A);%將系數(shù)數(shù)組轉(zhuǎn)化為多項(xiàng)式

subsfX'x'A);

Y=double(ans);

figure(l)

plot(XM'k”YWb%繪制擬合多項(xiàng)式函數(shù)圖形

xlabe(時(shí)刻，);%標(biāo)注坐標(biāo)軸含義

ylabel。平均氣溫，)；

title([num2str(n?l)J次函數(shù)的最小二乘曲線訕

grid;

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

yy=log(y);

n=3;

G=[]；

GG=n；

forj=O:(n-l)

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

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

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

GG=vertcat(GG,gg);%g亞宜連接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)==O;

sgn=O;

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(i,j)=G(i,j)+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(M+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);%計(jì)算誤差e

指數(shù)函數(shù)的系數(shù)是：%輸出系數(shù)及誤差

fprintf('\na=%f,b=%f,c=%f\a/b/c);e

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」);%繪制擬合指數(shù)函數(shù)圖形

xlabel。時(shí)亥山);%標(biāo)注坐標(biāo)軸含義

ylabel。平均氣溫1)；

title(「指數(shù)函數(shù)的最小二乘曲線']);

grid;

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

a、二次函數(shù)：

一天的平均氣溫為：21.2000

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

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

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

b、三次函數(shù)：

一天的平均氣溫為：21.2000

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

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

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

3次的額的晟小二乘由線

35

0510152025

時(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

黑

r

騾

正

25

20

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

時(shí)刻

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

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

最小，說(shuō)明結(jié)果最準(zhǔn)確。

4.設(shè)計(jì)算法，求出非線性方程6d-45x2+20=0的所有根，

并使誤差不超過(guò)

(1)算法思想

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

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

二分法是產(chǎn)生一串區(qū)間，使新區(qū)間/(“”是舊區(qū)間/("的一個(gè)子區(qū)間，其長(zhǎng)

度是八"的一半，且有一個(gè)端點(diǎn)是八”的一個(gè)端點(diǎn)。由區(qū)間⑸，讓+%

確定區(qū)間的方法是計(jì)算區(qū)間的中點(diǎn)

/f（x⑸…D）若/（/⑼/（/2））<0,則取

/（什也口⑹,xW+2）］,否則取/僅7=［小+2r（k+叫，重復(fù)這一過(guò)程即可。

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

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

1./■（晨°））=/。」（晨1））=/12.iJQOthenstop

3.Ifl/d<￡2then輸出X（°）作為根；stop

4.Ifl/J<￡2then輸出X。）作為根；stop

5.If1*°）-*|<￡山叫then輸出*作為根；stop

7./（*）=/

8.IfVl<￡2then輸出X作為根；

9.If/』<°then

9.1X=>X（°）.f^/oelse

9.2X=X（0；/n/110.goto5

（3）Matlab源程序

x=-100:100;

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

g=U;

fori=100:1:100%確定根所在的【乂問(wèn)

k=i+l;

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

g=[gil；

end

end

symsx;

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

n=length(g);%確定根的個(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,xO)*subs(f,x,(xO+>:l)/2)>eps

xO=(xO+xl)/2;

else

xl=(x0+xl)/2;

end

end

root=xO%輸出方程的根

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)行分析、求解。

（1）算法思想

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

一個(gè)方程后加只另一個(gè)方程，使方程組變成同解的上三角方程組，然后再自下

而上對(duì)上三角方程組求解。

列主元消去法是當(dāng)高斯消元到第k步時(shí)，從k列的””以下（包括0豚）的

各元素中選出絕對(duì)值最大的，然后通過(guò)行交換將其交換到的位置上。交換系

數(shù)矩陣中的兩行（包括常數(shù)項(xiàng)），只相當(dāng)于兩個(gè)方程的位置交換了，因此，列選

主元不影響求解的結(jié)果。

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

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

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

度都是有必要的，交換方程的原則是使噓…中，絕對(duì)值

最大的一個(gè)換到（k,k）位置而成為第k步消去的主元，這就是列主元Gauss消

去法。

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

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

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

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

typedefstruct

（

longintid;

longintver;

longintn;

)

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

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

版本號(hào)說(shuō)明

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

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

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

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

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

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

typedefstruct

(

longintq;//為上帶寬

longintp;//為下帶寬

)

(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、二進(jìn)制文件的讀?。?/p>

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

a=fread(f,3,Juinf)%讀取頭文件，3-讀取前3個(gè)，若讀取壓縮格式的，頭

文件為5個(gè)

b=fread(f,inf,ffloat*);%讀取剩下的文件,float型

id=dec2hex(a(l));

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

1.A的階數(shù)="2.Fork=123,?小一12.1找滿足

|a，*|二max|。伏|的下標(biāo)〃*

/>*2,2For/=L2,…，〃2.2.1

ik

-很

八2.4Fori=k+l,k+2,???，〃2.4.1kk2.4.2For

j=k+l,k+2/2.4.2.1ail~aikakl^atl2.4.3%0〉產(chǎn)4

n

B/an/(4N人產(chǎn)j)/an=0

Xn

Dn/°nnFork=H-1,H-2,-,1尸k+1(3)

Matlab源程序

dear;%清除工作空間變量

cic;%清除命令窗口命令

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

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

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

name=strcat(p,f);

file=fopen(name,'r');

,讀取二進(jìn)制頭文件

head=fread(file/num/uint');%

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

fprint"文件標(biāo)識(shí)符為。；

id

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

fprint"文件版本號(hào)為，);

ver

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

fprint"矩陣A的階數(shù))；

n

q=head(4);%匕帶寬

fprintf(知陣A的上帶寬)；

q

p=head(5);%下帶寬

fprintf。矩陣A的下帶寬);

P

dist=4*num;

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

(A,count]=fread(file,inf/float");%讀取二進(jìn)制文件，獲取系數(shù)矩陣

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

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

ifver==,102；

a=zeros(n,n);

fori=l:n,

forj=l:n,

a(ij)=A((i-l)*n+j);%求系數(shù)矩陣a(ij)

end

end

b=zeros(n,l)；

for1=1:0,

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

end

for%列主元高斯消去法

m=k;

%尋找主元

fori=k+l:nz

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

m=i;

end

end

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

disp('錯(cuò)誤!,)

return

end

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

t=a(k,j);

a(k,j)=a(mj);

a(mj)=t;

t=b(k);

b(k)=b(m);

b(m)=t;

end

計(jì)算并將其放到中

forl=k+l:nz%1(1,k)a(l,k)

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

forj=k+l:n

a(iJ)=a(ij)-a(izk)*a(kj);

end

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

end

end

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

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

fork=n-l:-l:lz

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

end

end

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

ifver==，202—%高斯消去法

m=p+q+l;

a=zeros(n,m);

fori=l:l:n

forj=l:l:m

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

end

end

b=zeros(nzl);

forl=l:l:n

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

end

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

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

disp(，錯(cuò)誤!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));

forj=(k+l):l:(k+q)

a(ij+p-i+l)=a(i/j+p-i+l)-a(i/k+p-i+l)*a(kj+p-k+l);

end

b(i)=b(i)-a(i,k+p-i+l)*b(k);

end

end

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

x(n)=b(n)/a(n,p+l);

sum=0;

fork=(n-l):-l:l

sum=b(k);

st2=n;

if(k+q)<n

st2=k+q;

end

forj=(k+l):l:st2

sum=sum-a(kzj+p-k+l)*x(j);

end

x(k)=sum/a(k,p+l);

sum=0;

end

end

dispf方程組的的解為：)％輸出方程組的解

disp(x)

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

方程組一：

文件標(biāo)識(shí)符為id=F1E1D1AO

文件版本號(hào)為ver=102

矩陣A的階數(shù)n=15

矩陣A的上帶寬q=6

矩陣A的下帶寬p=6

方程組的的解為：

1.0000

1.0000

1.0000

1.0000

1.0000

1.0000

1.0000

1.0000

1.0000

1.0000

1.0000

1.0000

1.0000

1.0000

1.0000

方程組二：

文件標(biāo)識(shí)符為id=F1E1D1AO

文件版本號(hào)為ver=102

矩陣A的階數(shù)n=2050

矩陣A的上帶寬q=6

矩陣A的下帶寬p=5

方程組的的解為：

1.9600

1.9600

1.9600

1.9600

1.9600

1.9600

1.9600

1.9600

1.9600

1.9600

方程組三：

文件標(biāo)識(shí)符為id=F1E1D1AO

文件版本號(hào)為ver=202

矩陣A的階數(shù)n=43512

矩陣A的上帶寬q-4

矩陣A的下帶寬p=3

方程組的的解為：

3.1413

3.1413

3.1413

3.1413

3.1413

3.1413

3.1413

3.1413

3.1413

3.1413

3.1413

實(shí)際應(yīng)用:

以求得三階系統(tǒng)的第二階固有頻率，現(xiàn)通過(guò)特征方程來(lái)求解其主振型，其中

A=[7227;113011;7227],b=[000],根據(jù)振動(dòng)理論歸一化理論,取x(2)=L

計(jì)算出x(l)、x(2)o即可得到第二階的主振型。利用Gauss法求解可得x=[T0

1]

收獲與體會(huì)

首先，非常感謝老師百忙之中給我們講授《計(jì)算方法》這門課程，使我對(duì)數(shù)

值計(jì)算有了一個(gè)全新的認(rèn)識(shí)，在課堂的學(xué)習(xí)中，我對(duì)數(shù)值計(jì)算方法有了一個(gè)基

本的了解，但是這門課程要經(jīng)過(guò)上機(jī)練習(xí)才能很好的掌握。

根據(jù)老師給定的題目，我運(yùn)用Matlab語(yǔ)言對(duì)題目進(jìn)行了編程，并對(duì)計(jì)算結(jié)

果進(jìn)行了詳細(xì)的分析。由于這是我第一次深入接觸Matlab語(yǔ)言，在編程的過(guò)程

中也遇到了不少困難，于是就去圖書館找資料，或者從網(wǎng)上查詢每一條語(yǔ)句的

用法，一句一句的編寫程序，最終編寫完了所有的程序，我的自信一下子提高

了，享受到了勞動(dòng)成果的滋味。

這次編程實(shí)踐，是我學(xué)完理論課程之后對(duì)自己該方面的能力的一次很好的檢

驗(yàn)，從開始的算法思路，到運(yùn)行調(diào)試，以及另人興奮的可用程序，都是一個(gè)很

好的學(xué)習(xí)和鍛煉的過(guò)程。使我鞏固了原有的理論知識(shí)，培養(yǎng)了我靈活運(yùn)用和組

合集成所學(xué)過(guò)知識(shí)及技能來(lái)分析解決實(shí)際問(wèn)題的能力。使我體會(huì)到自身知識(shí)和

能力能在實(shí)際中的應(yīng)用和發(fā)揮。不但可以激發(fā)創(chuàng)新意識(shí)，還可以開發(fā)創(chuàng)造能力、

培養(yǎng)溝通能力。

報(bào)告中存在的不當(dāng)之處，還請(qǐng)老師批評(píng)指正。

計(jì)算方法上機(jī)報(bào)告

姓名：

學(xué)號(hào)：

班級(jí)：

上課班級(jí)：

說(shuō)明:

本次上機(jī)實(shí)驗(yàn)使用的編程語(yǔ)言是Matlab語(yǔ)言，編譯

環(huán)境為MATLAB7.1L0,運(yùn)行平臺(tái)為Windows7。

1.對(duì)以下和式計(jì)算：

211\

8/>+4+58n+6/—、

0，要求：.

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

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

（1）算法思想

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

%=育18〃+「8〃+4-8"+5-8〃+6尸16〃8“+14\

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

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

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

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

…。尸擊（高-晟-壺-麻）.

2.for門=0,12…Jift&l°end;

3.for-2,???,0s=S+t；(3)Matlab源程序

clear;%清除工作空間變量

de;%清除命令窗口命令

m=input(，請(qǐng)輸入有效數(shù)字的位數(shù)m=');%輸入有效數(shù)字的I

5=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\rf,n-l);%需要將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ù)字時(shí)，需要將n值加到n=7,

s=3.1415926536；

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

s=3.14159265358979323846264338328c

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

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

2.某通信公司在一次施工中，需要在水面寬度為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.28133212.6111.2910.22

分點(diǎn)14151617181920

深度9.157.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ì)光纜的沖擊。海底光纜線的長(zhǎng)度預(yù)測(cè)模型如下所示，

光纜從A點(diǎn)鋪至B點(diǎn)，在某點(diǎn)處的深度為

，，

海底光纜線的長(zhǎng)度預(yù)測(cè)模型

計(jì)算光纜長(zhǎng)度時(shí)，用如下公式:

「2。=f27(x)VwWdx

L=If(x)ds

Jo

=￡廣/。){1+/3(^=￡43)2+(")2

女=0*(2)算法結(jié)構(gòu)

1.Fori=0,12…,〃i.l"'MbForhlZ.lFor

i=nn-l-k,i.i-x.J=>M

ttt2r/3XL-X0=>/I14

Fori=L2,???/l?14.1XJ+LX產(chǎn)力f+14.2

人+1/(力什%+1)=C,；1-C產(chǎn)a.2nb436MH4nd巧

.產(chǎn)明汨產(chǎn)乂〃40=>。02=>冊(cè);〃產(chǎn)%;2=36力=41必=>h7

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

8.ForA=2,3,…川8.1%/人-尸豆？%〃“產(chǎn)乙8.3

九?〃九_(tái)戶49.丫析/("""?］。Fork二6一1,小一2,…，110」

(七一入'+1)/"尸刈勺］.獲取x的元素個(gè)數(shù)存入s

XX

12.1=A13.Forf=l,2r-,5-113.1if-*then.break

elsei”=k14.Xkfk?Lh；Xk-X^X;XT

_3

VV萬(wàn)/—b4A，

必一后+咐石+石)》+以一如石)打/2⑴

Matlab源程序

clear;

de;

x=0:l:20;%產(chǎn)生從0至l]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:-l:k;

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

end

h(l)=x(2)-x(l);%計(jì)算三對(duì)角陣系數(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(1)=O;%選擇自然邊界條件

M(n)=O;

b(l)=2;

b(n)=2;

c(l)=O;

a(n)=O;

d(l)=O;

d(n)=O;

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

yl(l)=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ù)計(jì)算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

%%計(jì)算所需光纜改度

L=0;%計(jì)算所需光纜長(zhǎng)度

fori=2:N

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

end

disp(，所需光纜長(zhǎng)度為L(zhǎng)=');

disp(L);

figure

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

xlabelf位置；fontsizH16）;%標(biāo)注坐標(biāo)軸含義

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

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

grid;

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

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

鋪設(shè)河底光纜的曲線圖

位置

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

纜的走勢(shì)圖，計(jì)算出所需光纜的長(zhǎng)度為L(zhǎng)=26.4844mo

3.假定某天的氣溫變化記錄如下表所示，試用數(shù)據(jù)

擬合的方法找出這一天的氣溫變化的規(guī)律；試計(jì)算這

一天的平均氣溫，并試估計(jì)誤差。

（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ù)近似，雖然有很好的數(shù)值性質(zhì)，且計(jì)

算量也不大，但存放參數(shù)Mi的量很大，且沒有一個(gè)統(tǒng)一的數(shù)學(xué)公式來(lái)表示，也

帶來(lái)了一些不便。另一方面，在有的實(shí)際問(wèn)題中，用插值方法并不合適。當(dāng)數(shù)

據(jù)點(diǎn)的數(shù)目很大時(shí)，要求P（x）通過(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)更

加合理。通常情況下，是選取使2二最小，這就是最小二乘近似問(wèn)題。

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

數(shù)、三次函數(shù)、四次函數(shù)以及指數(shù)型函數(shù)計(jì)算相應(yīng)的系數(shù)，估

算誤差，并作圖比較各種函數(shù)之間的區(qū)別。

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

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

作為其最后一列（第〃+1列）存放。結(jié)果在°中，P是誤差2。

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

1.將“時(shí)刻值”存入?，數(shù)據(jù)點(diǎn)的個(gè)數(shù)存入川2,輸入擬合多項(xiàng)式函數(shù)

"（*）的最高項(xiàng)次數(shù)力二〃-1,則擬合多項(xiàng)式函數(shù)為

P（x）=a/（x）+"孫⑴+…"叫⑶根據(jù)給定數(shù)據(jù)點(diǎn)確定G

For/=0,1,2,???,/!-Ipor…，印2.1%=氏"12.2"尸九"13.

ForA=l,2,???/】3.1［形成矩陣。打

一sg〃(gn)(，g3"2no

3.1.103.1.29"一0n3A3.L3For

/=k+l,A+2,…M3.1.3.1%=叼3.1.4[變換到

GA]

m

03由3/(Jnt

32.1°n%AFor/二k+l，k+2,…L+13.2.20

3.2.3Fori=鼠k+1,…,m3.2.3.1。/“尸氏/4.[解三角方程

R。二61]

4.19小〃+1"〃=%4.2For/=n-1,n-2,--?,14.2.1

n

|g，/「E。力】/g〃=%2

尸—15.[計(jì)算誤差2]

m

E%一”

，=葉1"、使用指數(shù)函數(shù)擬合時(shí)

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

將C=y'=x代入得：y=oe-b(—婚對(duì)上式左右取對(duì)數(shù)得:

Iny=\na-bc2^2bcx-bx2令

2

z=lny,a0=lna-bctal=2bc,a2=-b則可得多項(xiàng)式：

2

Z=aQlalXIa2X(3)Matlab源程序

clear;%清除工作空間變顯

cic;%清除命令窗口命令

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ù)存入m

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

fprintfC一大的平均氣溫為T=%f\N,T);%求一天的平均氣溫

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

h=input(，請(qǐng)臨入擬合多項(xiàng)式的最高項(xiàng)次數(shù)力;％根據(jù)給定數(shù)據(jù)點(diǎn)生成矩陣G

n=h+l;

G=n；

forj=0:(n-l)

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

G=vertcat(G,g);%g垂克連接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伙)

ifG(k,k)>0;

sgn=l;

elseifG(k,k)==0;

sgn=0;

elsesgn="l;

end

A

sgm=-sgn*sqrt(sum(G(k:m/k).2));

w=zeros(l,n);

w(k)=G(kzk)-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:m,j))/bt;

fori=k:m;

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

end

end

end

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

A(n)=G(n,n+l)/G(nzn);%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

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

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

disp(A);

fphntf。使用％d次函數(shù)擬合的i關(guān)差是％f：

t=0:0.05:24;

A=fliplr(A);%將系數(shù)數(shù)組左右翻轉(zhuǎn)

Y=poly2sym(A);%將系數(shù)數(shù)組轉(zhuǎn)化為多項(xiàng)式

subs(X'x'zt);

Y=double(ans);

figure(l)

plot(x,y/k*，MK);%繪制擬合多項(xiàng)式函數(shù)圖形

xlabel。時(shí)刻1%標(biāo)注坐標(biāo)軸含義

ylabel('平均氣溫,)；

title([nuE2str(nr),'次函數(shù)的最小二乘曲線,]);

grid;

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

yy=log(y);

n=3;

G=n；

GG=[];

forj=O:(n-l)

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

G=vertcat(G,g);%g偃門.連接G

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

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

end

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

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

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

end

G;

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

ifG(k,k)>0;

sgn=l;

elseifG(kzk)==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(kzk)=sgm;%Gk-1Gk

forj=k+l:n+l

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

fori=k:m;

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

end

end

end

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

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

fori=n-l:-l:l

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

end

b=-A(3);

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);%l\算誤差e

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

fprintf('\n使用指數(shù)函數(shù)擬合的誤差是：％f,e);

t=0:0.05:24;

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

figure。)

plot(x,y「『,t,YY=r」);%繪制擬合指數(shù)函數(shù)圖形

xlabel(時(shí)亥『);%標(biāo)注坐標(biāo)軸含義

ylabelC平均氣溫匕

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

grid;

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

a、二次函數(shù)：

一天的平均氣溫為：21.2000

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

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

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

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

35

25

明

r2O

史

時(shí)

b、三次函數(shù)：

一天的平均氣溫為：21.2000

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

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

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

3次的額的晟小二乘由線

35

0510152025

時(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

黑

r

騾

正

25

20

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