计算方法上机作业VIP

计算方法上机报告

姓名：

学号：

班级：

上课班级：

说明：

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

环境为MATLAB7.11.0,运行平台为Windows7。

1.对以下和式计算：

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

。，要求：

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

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

(1)算法思想

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

%修(4211\14

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

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

digits（位数）或vpa（变量，精度为数）

（2）算法结构

「凸-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(，请输入有效数字的位数m=');%输入有效数字的位数

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

break;

end

end;

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

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

（4）结果与分析

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

s=3.1415926536；

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

s=3.14159265358979323846264338328o

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

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

2.某通信公司在一次施工中，需要在水面宽度为20

米的河沟底部沿直线走向铺设一条沟底光缆。在铺设

光缆之前需要对沟底的地形进行初步探测，从而估计

所需光缆的长度，为工程预算提供依据。已探测到一

组等分点位置的深度数据（单位:米）如下表所示:

分点0123456

深度9.018.967.967.978.029.0510.13

分点78910111213

深度11.1812.2613.28133212.6111.2910.22

分点14151617181920

深度9.157.907.958.869.8110.8010.93

」请用合适的曲线拟合所测数据点;

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

曲线图；

(1)算法思想

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

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

所用数据点越多越好，可以采用分段插值方式，即用分段多项式代替单个多项

式作插值。分段多项式是由一些在相互连接的区间上的不同多项式连接而成的

一条连续曲线，其中三次样条插值方法是一种具有较好“光滑性”的分段插值

方法。

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

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

光缆从A点铺至B点，在某点处的深度为

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

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

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

L=I/(x)ds

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

fc=Ok(2)算法结构

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

，二〃,〃-1,…，丘1.1（"|一时1）/3-*心）="3Xi-Xon/l"

Fori=L2,・・・/l・14.1XJ+LX产力f+14.2

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

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

获取M的矩阵元素个数，存入m

8.ForA=2,3.…M&l。〃人-「豆2%/"尸"出3

乙-〃'乙-产丫人9.七/'内=”析10.Fork二6一1,小一2,…，110」

（3“此+1）/〃「Mk”获取X的元素个数存入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%产生从0到20含21个等分点的数组

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

n=length（x）;%等分点的数11

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);%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);%对三对角阵进行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);%追赶法求解样条参数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

dispf所需光缆长度为L=')；

disp(L);

figure

plot(x,y*,X,s,u)%绘制铺设河底光缆的曲线图

xlabe(位置；fontsize,,16);%标注坐标轴含义

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

title('铺设河底光缆的曲线图］fontsizH16);

grid;

(4)结果与分析

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

铺设河底光缆的曲线图

位置

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

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

3.假定某天的气温变化记录如下表所示，试用数据

拟合的方法找出这一天的气温变化的规律；试计算这

一天的平均气温，并试估计误差。.

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

方法做数据近似要用较高次的多项式，这不仅给计算带来困难，更主要的缺点

是误差很大。用“插值样条函数”做数据近似，虽然有很好的数值性质，且计

算量也不大，但存放参数Mi的量很大，且没有一个统一的数学公式来表示，也

带来了一些不便。另一方面，在有的实际问题中，用插值方法并不合适。当数

据点的数目很大时，要求P（x）通过所有数据点，可能会失去原数据所表示的规

律。如果数据点是由测量而来的，必然带有误差，插值法要求准确通过这些不

准确的数据点是不合适的。在这种情况下，不用插值标准而用其他近似标准更

加合理。通常情况下，是选取°|使£二最小，这就是最小二乘近似问题。

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

数、三次函数、四次函数以及指数型函数，计算相应的系数，估

算误差，并作图比较各种函数之间的区别。

（2）算法结构

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

作为其最后一列（第〃+1列）存放。结果在。中，P是误差2。

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

1.将“时刻值”存入数据点的个数存入川2.输入拟合多项式函数

"（*）的最高项次数人二〃-1,则拟合多项式函数为

P（x）=a%（x）+aM（x）+…+%"x）根据给定数据点确定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.[计算误差2]

m

。：用力]

，=E〃+i”、使用指数函数拟合时

现将指数函数进行变形：

将C=y,=x代入得：y=Qe-ba-c『对上式左右取对数得:

Iny=Ina-bc2^2bcx-bx2令

2

z=lny,a0=lna-bcfal=2bc,a2=-b则可得多项式：

2

Z=aQlalX\a2X⑶Matlab源程序

clear;%清除工作空间变量

cic;%清除命令窗口命令

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

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

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

fprint"一天的平均气温为T=%An-J);%求一天的平均气温

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

h=input(，请输入拟合多项式的最高项次数=);%根据给定数据点生成矩阵G

n=h+l;

G=n；

forj=O:(n-l)

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

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

end

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

fori=l:m%将数据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);%解三角方程求系数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);

fprintfC使用％d次函数拟合的误差是％f：

t=0:0.05:24;

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

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

subsfX'x'A);

Y=double(ans);

figure(l)

plot(XM'k”YWb%绘制拟合多项式函数图形

xlabe(时刻，);%标注坐标轴含义

ylabel。平均气温，)；

title([num2str(n・l)J次函数的最小二乘曲线讪

grid;

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

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';%转置得到矩阵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)==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);%解三角方程求系数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);%计算误差e

指数函数的系数是：%输出系数及误差

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

fprintf('\n使用指数函数拟合的误差是:%f,,e);

t=0:0.05:24;

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

figure(2)

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

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

ylabel。平均气温1)；

title(「指数函数的最小二乘曲线']);

grid;

(4)结果与分析

a、二次函数：

一天的平均气温为：21.2000

2次函数的系数：8.30632.6064-0.0938

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

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

b、三次函数：

一天的平均气温为：21.2000

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

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

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

3次的额的晟小二乘由线

35

0510152025

时刻

C、四次函数：

一天的平均气温为：21.2000

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

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

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

4次困数的最小二乘曲线

35

30

黑

r

骡

正

25

20

d、指数函数：

一天的平均气温为：21.2000

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

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

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

指数由数的最小二系曲线

35

0510152025

时刻

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

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

最小，说明结果最准确。

4.设计算法，求出非线性方程6d-45x2+20=0的所有根，

并使误差不超过

(1)算法思想

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

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

二分法是产生一串区间，使新区间/(“”是旧区间/("的一个子区间，其长

度是八"的一半，且有一个端点是八”的一个端点。由区间⑸，让+%

确定区间的方法是计算区间的中点

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

/（什也口⑹,xW+2）］,否则取/仅7=［小+2r（k+叫，重复这一过程即可。

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

（2）算法结构

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;%非线性方程组的表达式

g=U;

fori=100:1: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,xO)*subs(f,x,(xO+>:l)/2)>eps

xO=(xO+xl)/2;

else

xl=(x0+xl)/2;

end

end

root=xO%输出方程的根

end

(4)结果与分析

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

差要求。

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

程序，并对所附的方程组进行求解。针对本专业中所

碰到的实际问题，提炼一个使用方程组进行求解的例

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

（1）算法思想

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

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

而上对上三角方程组求解。

列主元消去法是当高斯消元到第k步时，从k列的””以下（包括0豚）的

各元素中选出绝对值最大的，然后通过行交换将其交换到的位置上。交换系

数矩阵中的两行（包括常数项），只相当于两个方程的位置交换了，因此，列选

主元不影响求解的结果。

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

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

消去过程中，适当交换方程的顺序对保证消去过程能顺利进行及计算解的精确

度都是有必要的，交换方程的原则是使嘘…中，绝对值

最大的一个换到（k,k）位置而成为第k步消去的主元，这就是列主元Gauss消

去法。

（2）算法结构

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

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

（1）文件头部分，其结构：

typedefstruct

（

longintid;

longintver;

longintn;

)

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

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

版本号说明

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

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

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

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

n:表示方程的阶数

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

typedefstruct

(

longintq;//为上带宽

longintp;//为下带宽

)

(3)系数矩阵

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

为n*n,类型为float型；

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

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

(4)右端系数

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

3、二进制文件的读取：

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

a=fread(f,3,Juinf)%读取头文件，3-读取前3个，若读取压缩格式的，头

文件为5个

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

id=dec2hex(a(l));

ver=dec2hex(a(2));%这两句是进行进制转换，读取id与ver

1.A的阶数="2.Fork=123,・小一12.1找满足

|a，*|二max|。伏|的下标〃*

/>*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〉产4

n

B/an/(4N人产j)/an=0

Xn

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

Matlab源程序

dear;%清除工作空间变量

cic;%清除命令窗口命令

%%读取系数矩阵

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

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

name=strcat(p,f);

file=fopen(name,'r');

,读取二进制头文件

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

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

fprint"文件标识符为。；

id

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

fprint"文件版本号为，);

ver

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

fprint"矩阵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(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('错误!,)

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

计算并将其放到中

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);%回代过程

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

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

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)%开始消去过程

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

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);%回代过程

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

方程组一：

文件标识符为id=F1E1D1AO

文件版本号为ver=102

矩阵A的阶数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

方程组二：

文件标识符为id=F1E1D1AO

文件版本号为ver=102

矩阵A的阶数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

方程组三：

文件标识符为id=F1E1D1AO

文件版本号为ver=202

矩阵A的阶数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

实际应用:

以求得三阶系统的第二阶固有频率，现通过特征方程来求解其主振型，其中

A=[7227;113011;7227],b=[000],根据振动理论归一化理论,取x(2)=L

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

1]

收获与体会

首先，非常感谢老师百忙之中给我们讲授《计算方法》这门课程，使我对数

值计算有了一个全新的认识，在课堂的学习中，我对数值计算方法有了一个基

本的了解，但是这门课程要经过上机练习才能很好的掌握。

根据老师给定的题目，我运用Matlab语言对题目进行了编程，并对计算结

果进行了详细的分析。由于这是我第一次深入接触Matlab语言，在编程的过程

中也遇到了不少困难，于是就去图书馆找资料，或者从网上查询每一条语句的

用法，一句一句的编写程序，最终编写完了所有的程序，我的自信一下子提高

了，享受到了劳动成果的滋味。

这次编程实践，是我学完理论课程之后对自己该方面的能力的一次很好的检

验，从开始的算法思路，到运行调试，以及另人兴奋的可用程序，都是一个很

好的学习和锻炼的过程。使我巩固了原有的理论知识，培养了我灵活运用和组

合集成所学过知识及技能来分析解决实际问题的能力。使我体会到自身知识和

能力能在实际中的应用和发挥。不但可以激发创新意识，还可以开发创造能力、

培养沟通能力。

报告中存在的不当之处，还请老师批评指正。

计算方法上机报告

姓名：

学号：

班级：

上课班级：

说明:

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

环境为MATLAB7.1L0,运行平台为Windows7。

1.对以下和式计算：

211\

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

0，要求：.

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

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

（1）算法思想

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

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

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

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

digits（位数）或vpa（变量，精度为数）

（2）算法结构

…。尸击（高-晟-壶-麻）.

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

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

clear;%清除工作空间变量

de;%清除命令窗口命令

m=input(，请输入有效数字的位数m=');%输入有效数字的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)%判断通项与精度的关系

break;

end

end;

fprintf。需要将n值加到n=%d\rf,n-l);%需要将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.14159265358979323846264338328c

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

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

2.某通信公司在一次施工中，需要在水面宽度为20

米的河沟底部沿直线走向铺设一条沟底光缆。在铺设

光缆之前需要对沟底的地形进行初步探测，从而估计

所需光缆的长度，为工程预算提供依据。己探测到一

组等分点位置的深度数据(单位:米)如下表所不:

分点0123456

深度9.018.967.967.978.029.0510.13

分点78910111213

深度11.1812.2613.28133212.6111.2910.22

分点14151617181920

深度9.157.907.958.869.8110.8010.93

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

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

曲线图；

(1)算法思想

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

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

所用数据点越多越好，可以采用分段插值方式，即用分段多项式代替单个多项

式作插值。分段多项式是由一些在相互连接的区间上的不同多项式连接而成的

一条连续曲线，其中三次样条插值方法是一种具有较好“光滑性”的分段插值

方法。

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

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

光缆从A点铺至B点，在某点处的深度为

，，

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

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

「2。=f27(x)VwWdx

L=If(x)ds

Jo

=£广/。){1+/3(^=£43)2+(")2

女=0*(2)算法结构

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产力f+14.2

人+1/(力什%+1)=C,；1-C产a.2nb436MH4nd巧

.产明汨产乂〃40=>。02=>册;〃产%;2=36力=41必=>h7

获取M的矩阵元素个数，存入m

8.ForA=2,3,…川8.1%/人-尸豆？%〃“产乙8.3

九・〃九_户49.丫析/("""?］。Fork二6一1,小一2,…，110」

(七一入'+1)/"尸刈勺］.获取x的元素个数存入s

XX

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

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

_3

VV万/—b4A，

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

Matlab源程序

clear;

de;

x=0:l:20;%产生从0至l]20含21个等分点的数组

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

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

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(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);%对三对角阵进行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);%追赶法求解样条参数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))^2);

end

disp(，所需光缆长度为L=');

disp(L);

figure

plot(x,y,*,X,s,'-')%绘制铺设河底光缆的曲线图

xlabelf位置；fontsizH16）;%标注坐标轴含义

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

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

grid;

（4）结果与分析

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

铺设河底光缆的曲线图

位置

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

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

3.假定某天的气温变化记录如下表所示，试用数据

拟合的方法找出这一天的气温变化的规律；试计算这

一天的平均气温，并试估计误差。

（1）算法思想

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

方法做数据近似要用较高次的多项式，这不仅给计算带来困难，更主要的缺点

是误差很大。用“插值样条函数”做数据近似，虽然有很好的数值性质，且计

算量也不大，但存放参数Mi的量很大，且没有一个统一的数学公式来表示，也

带来了一些不便。另一方面，在有的实际问题中，用插值方法并不合适。当数

据点的数目很大时，要求P（x）通过所有数据点，可能会失去原数据所表示的规

律。如果数据点是由测量而来的，必然带有误差，插值法要求准确通过这些不

准确的数据点是不合适的。在这种情况下，不用插值标准而用其他近似标准更

加合理。通常情况下，是选取使2二最小，这就是最小二乘近似问题。

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

数、三次函数、四次函数以及指数型函数计算相应的系数，估

算误差，并作图比较各种函数之间的区别。

（2）算法结构

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

作为其最后一列（第〃+1列）存放。结果在°中，P是误差2。

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

1.将“时刻值”存入〜，数据点的个数存入川2,输入拟合多项式函数

"（*）的最高项次数力二〃-1,则拟合多项式函数为

P（x）=a/（x）+"孙⑴+…"叫⑶根据给定数据点确定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.[计算误差2]

m

E%一”

，=叶1"、使用指数函数拟合时

现将指数函数进行变形：

将C=y'=x代入得：y=oe-b(—婚对上式左右取对数得:

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

2

z=lny,a0=lna-bctal=2bc,a2=-b则可得多项式：

2

Z=aQlalXIa2X(3)Matlab源程序

clear;%清除工作空间变显

cic;%清除命令窗口命令

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

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

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

fprintfC一大的平均气温为T=%f\N,T);%求一天的平均气温

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

h=input(，请临入拟合多项式的最高项次数力;％根据给定数据点生成矩阵G

n=h+l;

G=n；

forj=0:(n-l)

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

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

end

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

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

解三角方程求系数

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);%计算误差e

fprintf('%d次函数的系数是：)h);%输出系数a及误差e

disp(A);

fphntf。使用％d次函数拟合的i关差是％f：

t=0:0.05:24;

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

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

subs(X'x'zt);

Y=double(ans);

figure(l)

plot(x,y/k*，MK);%绘制拟合多项式函数图形

xlabel。时刻1%标注坐标轴含义

ylabel('平均气温,)；

title([nuE2str(nr),'次函数的最小二乘曲线,]);

grid;

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

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

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

解三角方程求系数

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

fprintf('\n使用指数函数拟合的误差是：％f,e);

t=0:0.05:24;

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

figure。)

plot(x,y「『,t,YY=r」);%绘制拟合指数函数图形

xlabel(时亥『);%标注坐标轴含义

ylabelC平均气温匕

title(「指数函数的最小二乘曲线1);

grid;

(4)结果与分析

a、二次函数：

一天的平均气温为：21.2000

2次函数的系数：8.30632.6064-0.0938

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

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

2次函数的最小二乘由线

35

25

明

r2O

史

时

b、三次函数：

一天的平均气温为：21.2000

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

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

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

3次的额的晟小二乘由线

35

0510152025

时刻

C、四次函数：

一天的平均气温为：21.2000

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

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

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

4次困数的最小二乘曲线

35

30

黑

r

骡

正

25

20

d、指数函数：

一天的平均气温为：21.2000

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

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

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

指数由数的最小二系曲线

35

0510152025