计算方法上机课后复习

计算方法上机报告

姓名:

学号:

班级:

上课班级:

说明：

本次上机实验使用的编程语言是Matlab语言，编译环境为MATLAB

7.11.0,运行平台为Windows7。

30],42]]、

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

①若只需保留11个有效数字，该如何进行计算；

②若要保留30个有效数字，则又将如何进行计算；

(1)算法思想

1、根据精度要求估计所加的项数，可以使用后验误差估计，通项为：

1(42114

a„=--------------------------------------------------<-------------<£；

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

2、为了保证计算结果的准确性，写程序时，从后向前计算；

3、使用Matlab时，可以使用以下函数控制位数：

digits(位数)或vpa(变量，精度为数)

(2)算法结构

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请输入有效数字的位数m=1）;%输入有效数字的位数

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)%判断通项与精度的关系

break;

end

end;

fprintf('需要将n值加到n=%d\n；n-1);%需要将n值加到的数值

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;%求和运算

end

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

（4）结果与分析

当保留11位有效数字时，需要将n值加到n=7,

s=3.1415926536；

当保留30位有效数字时，需要将n值加到n=22,

s=3.14159265358979323846264338328。

通过上面的实验结果可以看出，通过从后往前计算，这种算法很好的保证了计算结

果要求保留的准确数字位数的要求。

2.某通信公司在一次施工中，需要在水面宽度为20米的河沟底部沿直线走

向铺设一条沟底光缆。在铺设光缆之前需要对沟底的地形进行初步探测,

从而估计所需光缆的长度，为工程预算提供依据。己探测到一组等分点

位置的深度数据(单位:米)如下表所示：

分点0123456

深度9.018.967.967.978.029.0510.13

分点78910111213

深度11.1812.2613.2813.3212.6111.2910.22

分点14151617181920

深度9.157.907.958.869.8110.8C10.93

①请用合适的曲线拟合所测数据点；

②预测所需光缆长度的近似值，作出铺设河底光缆的曲线图；

(1)算法思想

如果使用多项式差值，则由于龙格现象，误差较大，因此，用相对较少的插值数据

点作插值，可以避免大的误差，但是如果又希望将所得数据点都用上，且所用数据点越

多越好，可以采用分段插值方式，即用分段多项式代替单个多项式作插值。分段多项式

是由一些在相互连接的区间上的不同多项式连接而成的一条连续曲线，其中三次样条插

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

在本题中，假设所铺设的光缆足够柔软，在铺设过程中光缆触地走势光滑，紧贴地

面，并且忽略水流对光缆的冲击。海底光缆线的长度预测模型如下所示，光缆从A点铺

至B点，在某点处的深度为

海底光缆线的长度预测模型

计算光缆长度时，用如下公式:

20

Jof(<x)ds

=x£+l

=Z，3)2+(Ay)2

(2)算法结构

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的矩阵元素个数，存入m

8.For%=2,3,…，机

8.1aj心二lk

8.2bk-lk-%n网

8.34-小九_]=九

9.7mMm=M,n

10.Fork=

10.1(九-q—)/4nM.

11.获取x的元素个数存入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;%产生从0到20含21个等分点的数组

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];%等分点位置的深度数据

n=length(x);%等分点的数目

N=length(X);

%%求三次样条插值函数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);%计算三对角阵系数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);%对三对角阵进行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);%追赶法求解样条参数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);%在各区间用三次样条插值函数计算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

%%计算所需光缆长度

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；-')%绘制铺设河底光缆的曲线图

xlabel('位置'「fontsize',16);%标注坐标轴含义

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

title（铺设河底光缆的曲线图；fontsize',16）;

grid;

（4）结果与分析

铺设海底光缆的曲线图如下图所示:

仿真结果表明，运用分段三次样条插值所得的拟合曲线能较准确地反映铺设光缆的

走势图，计算出所需光缆的长度为L=26.4844m。

3.假定某天的气温变化记录如下表所示，试用数据拟合的方法找出这一天

的气温变化的规律；试计算这一天的平均气温，并试估计误差。

时刻0123456789101112

平均气温15141414141516182020232528

时刻131415161718192021222324

平均气温313431292725242220181716

(1)算法思想

在本题中，数据点的数目较多。当数据点的数目很多时，用“多项式插值”方法做

数据近似要用较高次的多项式，这不仅给计算带来困难，更主要的缺点是误差很大。用

“插值样条函数”做数据近似，虽然有很好的数值性质，且计算量也不大，但存放参数

M的量很大，且没有一个统一的数学公式来表示，也带来了一些不便。另一方面，在有

的实际问题中，用插值方法并不合适。当数据点的数目很大时，要求P(x)通过所有数据

点，可能会失去原数据所表示的规律。如果数据点是由测量而来的，必然带有误差，插

值法要求准确通过这些不准确的数据点是不合适的。在这种情况下，不用插值标准而用

其他近似标准更加合理。通常情况下，是选取巴使E2最小，这就是最小二乘近似问题。

在本题中，采用“最小二乘法”找出这一天的气温变化的规律，使用二次函数、三

次函数、四次函数以及指数型函数c=四十.解，计算相应的系数，估算误差，并作图比

较各种函数之间的区别。

(2)算法结构

本算法用正交化方法求数据的最小二乘近似。假定数据以用来生成了G，并将

)，作为其最后一列(第〃+i列)存放。结果在々中，夕是误差七〉

I、使用二次函数、三次函数、四次函数拟合时

1.将“时刻值”存入数据点的个数存入加

2.输入拟合多项式函数P&)的最高项次数〃=1,则拟合多项式函数为

p(x)=(x)+a2g2(x)+•••+angn(x),根据给定数据点确定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履严

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.「解三角方程网=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使用指数函数拟合时

现将指数函数进行变形：

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

对上式左右取对数得:Iny=In6(-fee2+Ibex-bx1

2

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

2

贝！可得多项式：z=4+a”+a2x

(3)Matlab源程序

clear;%清除工作空间变量

cic;%清除命令窗口命令

x=0:24;%将时刻值存入数组

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);%将数据点的个数存入m

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

fprintf。/的平均气温为T=%f\n',T);%求一天的平均气温

%%二次、三次、四次函数的最小二乘近似

h二input（'请输入拟合多项式的最高项次数二’）；%根据给定数据点生成矩阵G

n=h+l；

G=[]；

forj=0:(n-l)

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

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

end

G=G';%转置得到矩阵G

fori=l:m%将数据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);%解三角方程求系数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次函数的系数是:；h);%输出系数a及误差e

disp(A);

fprintf(使用％d次函数拟合的误差是%f:\h,e);

t=0:0.05:24;

A=fliplr(A);%将系数数组左右翻转

Y=poly2sym(A);%将系数数组转化为多项式

subs(Y/x',t);

Y=double(ans);

figure(l)

plot(x,y,'k*',t,Y,%绘制拟合多项式函数图形

xlabel(时刻');%标注坐标轴含义

ylabel('平均气温')；

title([num2str(n-l)二次函数的最小二乘曲线']);

grid;

%%指数函数的最小二乘近似

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';%转置得到矩阵G

fori=l:m%将数据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

%解三角方程求系数

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指数函数的系数是:a=%f,b=%f,c=%f1,afb,c);%输出系数及误差e

fprintf（'\n使用指数函数拟合的误差是：%f1,e）;

t=0:0.05:24;

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

figure⑵

plot(x,y/k*',t,YY,'r-');%绘制拟合指数函数图形

xlabel。时亥!T);%标注坐标轴含义

ylabel('平均气温');

title（「指数函数的最小二乘曲线】）;

grid;

（4）结果与分析

a、二次函数:

一天的平均气温为:21.2000

2次函数的系数:8.30632.6064-0.0938

使用2次函数拟合的误差是：280.339547

二次函数的最小二乘曲线如下图所示:

2次函数的最小二乘曲线

35

30

25

r喟

处20

H-

10

10152025

时刻

b、三次函数：

一天的平均气温为:21.2000

3次函数的系数：13.3880-0.22730.2075-0.0084

使用3次函数拟合的误差是：131.061822

三次函数的最小二乘曲线如下图所示：

3次困数的最小二乘曲线

时刻

C、四次函数：

一天的平均气温为：21.2000

4次函数的系数：16.7939-3.70500.8909-0.05320.0009

使用4次函数拟合的误差是：59.04118

四次函数的最小二乘曲线如下图所示：

4次函数的最小二乘曲线

35

30

25

20

15-

1°0

10152025

时刻

d、指数函数：

一天的平均气温为:21.2000

指数函数的系数是:a=26.160286,b=0.004442,c=14.081900

使用指数函数拟合的误差是：57.034644

指数函数的最小二乘曲线如下图所示:

指数困数的最小二乘曲线

35

0510152025

时刻

通过上述几种拟合可以发现，多项式的次数越高，计算拟合的效果越好，误差越

小，说明结果越准确；同时，指数多项式拟合的次数虽然不高，但误差最小，说明结果

最准确。

4.设计算法，求出非线性方程6炉-45/+20=°的所有根，并使误差不超过

KF4。

（1）算法思想

首先，研究函数的形态，确定根的范围；通过剖分区间的方法确定根的位置，然后

利用二分法的基本原理进行求解，找到满足精度要求的解。

二分法是产生一串区间，使新区间产+D是旧区间/（幻的一个子区间，其长度是/⑷

的一半，且有一个端点是/㈤的一个端点。由区间/阴二次叱胆+"］确定区间/的的方法

是计算区间八幻的中点

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

2

若/（产）。（产坳）＜0,则取/伏划=口叫公『，否则取产；俨+2"叫，重复这一过

程即口」。显然，每次迭代使区间长度减小一半，故二分法总是收敛的。

（2）算法结构

L"))"；⑴)n力

2.Iffof\>0thenstop

3.31yolv邑then输出大(。)作为根;stop

4.If|工|<&then输出产作为根;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;%非线性方程组的表达式

g=[];

fori=-100:l:100%确定根所在的区间

k=i+l;

if(y(x==i).*y(x==k)<eps)%区间长度为1

g=[gil;

end

end

symsx;

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

n=length(g);%确定根的个数

forj=l:n

xO=g(j);%求根区间左端点

xl=g(j)+l;%求根区间右端点

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)结果与分析

该非线性方程组有三个实根，分别为1.8708,0.6812,-0.6545,且满足误差要求。

5.编写程序实现大规模方程组的列主元高斯消去法程序，并对所附的方程

组进行求解。针对本专业中所碰到的实际问题，提炼一个使用方程组进行

求解的例子，并对求解过程进行分析、求解。

(1)算法思想

高斯消去法是利用现行方程组初等变换中的一种变换，即用一个不为零的数乘一个

方程后加只另一个方程，使方程组变成同解的上三角方程组，然后再自下而上对上三角

方程组求解。

列主元消去法是当高斯消元到第欠步时，从女列的&以下〔包括火火)的各元素中

选出绝对值最大的，然后通过行交换将其交换到外.的位置上。交换系数矩阵中的两行

(包括常数项)，只相当于两个方程的位置交换了，因此，列选主元不影响求解的结果。

程序的核心就是高斯列主元消去法。根据教材提供的算法，编写列主元消去法的子

函数与适应于超大规模超出系统内存的方程组的改编程序。同时，在Gauss消去过程中，

适当交换方程的顺序对保证消去过程能顺利进行及计算解的精确度都是有必要的，交换

方程的原则是使a十&+1,…中，绝对值最大的一个换到(k,k)位置而成为第k

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

(2)算法结构

1、数据文件的文件名为：文件名+.dat

2、数据文件中的数据为二进制记录结构，分为以下四个部分：

(1)文件头部分，其结构：

typedefstruct

(

longintid;

longintver;

longintn;

)

其中：id:为该数据文件的标识，值为OxFIElDlAO,即为：十六进制的F1E1D1A0

ver:为数据文件的版本号，值为16进制数据，

版本号说明

0x101系数矩阵为非压缩格式稀疏矩阵

0x102系数矩阵为非压缩格式带状对角阵

0x201系数矩阵为压缩格式稀疏矩阵

0x202系数矩阵为压缩格式带状对角阵

n:表示方程的阶数

(2)文件头2:此部分说明为条状矩阵的上下带宽，结构：

typedefstruct

(

longintq;//为上带宽

longintp;//为下带宽

J

(3)系数矩阵

a.如存贮格式非为压缩方式，则按行方式存贮系数矩阵中的每一个元素，个数为

n*n,类型为float型；

b.如果存贮格式是压缩方式，则按行方式存贮，每行中只存放上下带宽内的非零元素,

即，每行中存贮的最多元素为p+q+1个。

⑷右端系数

按顺序存贮右端系数的每个元素，个数为n个，类型为float型

3、二进制文件的读取：

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));%这两句是进行进制转换，读取id与ver

1.A的阶数

2.Fork=1,2,3,-

2.1找满足以J=i暝的下标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;%清除命令窗口命令

%%读取系数矩阵

[f,p]=uigetfile('*.dat；选择数据文件')；%读取数据文件

num=5;%输入系数矩阵文件头的个数

name=strcat(p,f);

file=fopen(name/r');

head=fread(file,num；uint,);%读取二进制头文件

id=dec2hex(head(l));%读取标识符

fprintf('文件标识符为')；

id

ver=dec2hex(head(2));%读取版本号

fprintf('文件版本号为');

ver

n=head(3);%读取阶数

fprintf(矩阵A的阶数’)；

n

q=head(4);%上带宽

fprintf(矩阵A的上带宽)；

q

p=head(5);%下带宽

fprintf(矩阵A的下带宽)；

P

dist=4*num;

fseek(file,dist,'bof');%把句柄值转向第六个元素开头处

%读取二进制文件，获取系数矩阵

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

fclose(file);%关闭二进制头文件

%%对非压缩带状矩阵进行求解

ifver=='102\

a=zeros(n,n);

fori=l:n,

forj=l:n,

%求系数矩阵

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

%%对压缩带状矩阵进行求解

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))；