计算方法上机作业之令狐文艳创作
计算方法上机报告
令狐文艳
姓名：
学号：
班级：
上课班级：
说明：
本次上机实验使用的编程语言是Matlab 语言，编译环境为MATLAB 7.11.0，运行平台为Windows 7。 1. 对以下和式计算：
∑
∞
??
?
??+-+-+-+=0681581482184161n n n n S n ，要求：
①若只需保留11个有效数字，该如何进行计算； ②若要保留30个有效数字，则又将如何进行计算； （1） 算法思想
1、根据精度要求估计所加的项数，可以使用后验误差估计，通项为：
1421114
16818485861681n n n a n n n n n ε??=
---<< ?
+++++??；
2、为了保证计算结果的准确性，写程序时，从后向前计算；
3、使用Matlab 时，可以使用以下函数控制位数： digits(位数)或vpa(变量，精度为数)
（2）算法结构
1.;0=s
??? ??+-+-+-+=681581482184161n n n n t n
；
2.
for 0,1,2,,n i =???
if 10m
t -≤
end;
3.
for ,1,2,,0n i i i =--???
;s s t =+
（3）Matlab 源程序
clear; %清除工作空间变量 clc; %清除命令窗口命令
m=input('请输入有效数字的位数m='); %输入有效数字的位数 s=0; for n=0:50
t=(1/16^n)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6));
if t<=10^(-m) %判断通项与精度的关系
break;
end
end;
fprintf('需要将n值加到n=%d\n',n-1); %需要将n值加到的数值
for i=n-1:-1:0
t=(1/16^i)*(4/(8*i+1)-2/(8*i+4)-1/(8*i+5)-1/(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米的河沟底部
沿直线走向铺设一条沟底光缆。在铺设光缆之前需要对沟底的地形进行初步探测，从而估计所需光缆的长度，为工程预算提供依据。已探测到一组等分点位置的深度数据(单位:米)如下表所示：
①
②预测所需光缆长度的近似值，作出铺设河底光缆的曲线图；（1）算法思想
如果使用多项式差值，则由于龙格现象，误差较大，因此，用相对较少的插
值数据点作插值，可以避免大的误差，但是如果又希望将所得数据点都用上，且所用数据点越多越好，可以采用分段插值方式，即用分段多项式代替单个多项式作插值。分段多项式是由一些在相互连接的区间上的不同多项式连接而成的一条连续曲线，其中三次样条插值方法是一种具有较好“光滑性”的分段插值方法。
在本题中，假设所铺设的光缆足够柔软，在铺设过程中光缆触地走势光滑，紧贴地面，并且忽略水流对光缆的冲击。海底光缆线的长度预测模型如下所示，光缆从A 点铺至B 点，在某点处的深度为i h 。
海底光缆线的长度预测模型
计算光缆长度时，用如下公式：
20
()L f x ds
=
?
20
'20
()1()f x f x dx =+?
19
1
'20()1()k k
k f x f x dx
+==+∑?
22
()()x y =?+?（2）算法结构
1.
For n i ,,2,1,0???= 1.1 i
i M y ?
2.
For 2,1=k
2.1 For k n n i ,,1, -=
2.1.1
i
k i i i i M x x M M ?----)/()(1
3. 1
01h x x ?-
4.
For 1-,,2,1n i = 4.1 1
1++?-i i i h x x
4.2 b
a c c h h h i i i i i i ??-?+++2;1;)/(11
4.3
i
i d M ?+16
5.
000;;c M d M d n n ???λ
n
n n b a b ???2;;20μ
6. 1111,γμ??d b
7. 获取M 的矩阵元素个数，存入m 8.
For m k ,,3,2 = 8.1 k
k k l a ?-1/μ 8.2 k
k k k c l b μ??-1- 8.3
k
k k k l d γγ??-1-
9. m
m m M ?μγ/
10.
For 1,,2,1 --=m m k 10.1
k
k k k k M M c ??-+μγ/)(1
11. 获取x 的元素个数存入s
12. k ?1
13.
For 1,,2,1-=s i
13.1 if i
x x ≤~ then k i ?；break
else k i ?+1
14.
x
x x x x x h x x k k k k ?~;~;
11?-?-?---
y
h x h M y x h M y x M x M k k k k k k ~/]?)6()6(6?6[2
211331
?-+-++---
（3）Matlab 源程序
clear;
clc;
x=0:1:20; %产生从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); %等分点的数目
N=length(X);
%% 求三次样条插值函数s(x)
M=y;
for k=2:3; %计算二阶差商并存放在M中
for i=n:-1:k;
M(i)=(M(i)-M(i-1))/(x(i)-x(i-k+1));
end
end
h(1)=x(2)-x(1); %计算三对角阵系数a,b,c及右端向量d
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);
end
M(1)=0; %选择自然边界条件
M(n)=0;
b(1)=2;
b(n)=2;
c(1)=0;
a(n)=0;
d(1)=0;
d(n)=0;
u(1)=b(1); %对三对角阵进行LU分解
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);
end
M(n)=y1(n)/u(n);%追赶法求解样条参数M(i)
for k=n-1:-1:1;
M(k)=(y1(k)-c(k)*M(k+1))/u(k);
end
s=zeros(1,N);
for m=1:N;
k=1;
for i=2:n-1
if X(m)<=x(i);
k=i-1;
break;
else
k=i;
end
end
H=x(k+1)-x(k); %在各区间用三次样条插值函数计算X点处的值 x1=x(k+1)-X(m);
x2=X(m)-x(k); s(m)=(M(k)*(x1^3)/6+M(k+1)*(x2^3)/6+(y(k)-
(M(k)*(H^2)/6))*x1+(y(k+1)-(M(k+1)*(H^2)/6))*x2)/H;
end
%% 计算所需光缆长度
L=0; %计算所需光缆长度
for i=2:N
L=L+sqrt((X(i)-X(i-1))^2+(s(i)-s(i-1))^2);
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.假定某天的气温变化记录如下表所示，试用数据拟合的方法找出这一天的气温变化的规律；试计算这一天的平均气温，并试估计误差。
时刻
1
2
3
4
5
6
7
8
9 10 11 12
平均气温 15 14 14 14 14 15 16 18 20 20 23 25 28 时刻
13 14 15 16 17 18 19 20 21 22 23 24
平均气温 31 34 31 29 27 25 24 22 20 18 17 16
（1在本题中，数据点的数目较多。当数据点的数目很多时，用“多项式插值”方法做数据近似要用较高次的多项式，这不仅给计算带来困难，更主要的缺点是误差很大。用“插值样条函数”做数据近似，虽然有很好的数值性质，且计算量也不大，但存放参数Mi 的量很大，且没有一个统一的数学公式来表示，也带来了一些不便。另一方面，在有的实际问题中，用插值方法并不合适。当数据点的数目很大时，要求)(x p 通过所有数据点，可能会失去原数据所表示的规律。如果数据点是由测量而来的，必然带有误差，插值法要求准确通过这些不准确的数据点是不合适的。在这种情况下，不用插值标准而用其他近似标准更加合理。通常情况下，是选取i 使2E 最小，这就是最小二
乘近似问题。
在本题中，采用“最小二乘法”找出这一天的气温变化的规律，使用二次函数、三次函数、四次函数以及指数型函数2
)(c t b ae C --=，计算相应的系
数，估算误差，并作图比较各种函数之间的区别。
（2）算法结构
本算法用正交化方法求数据的最小二乘近似。假定数据以用来生成了G ，并
将y 作为其最后一列（第1+n 列）存放。结果在a 中，ρ是误差2
2E 。
I 、使用二次函数、三次函数、四次函数拟合时
1. 将“时刻值”存入
i
x ，数据点的个数存入m
2.
输入拟合多项式函数)(x p 的最高项次数1-=n h ，则拟合多项式函数为
)()()()(2211x g x g x g x p n n ααα+???++=,根据给定数据点确定G
For
1,,2,1,0-???=n j
For m i ,,2,1???=
2.1 1
,+?j i j
i g x
2.2
1
,+?n i i g y
3.For n k ,,2,1???=
3.1 [形成矩阵
k
Q ]
3.1.1 σ
?-∑=m
k
i ik kk g g 2
/12))(sgn(
3.1.2
k
kk g ωσ?-
3.1.3 For m k k j ,,2,1???++=
3.1.3.1
j
jk g ω?
3.1.4
β
σω?k 3.2 [变换1
-k G 到
k
G ]
3.2.1
kk g ?σ
For 1,,,2,1+???++=n n k k j
3.2.2
t
g m
k
i ij i ?∑=βω/)(
3.2.3 For m k k i ,,1,???+=
3.2.3.1
ij
i ij g t g ?+ω
4. [解三角方程1h Ra =] 4.1
n
nn n n a g g ?+/1.
4.2For 1,,2,1???--=n n i 4.2.1 i
ii n
i j j ij
n i a g x g
g ?-
∑+=+/][1
1.
5.
[计算误差22E ]
ρ
?∑+=+m
n i n i g
1
2
1
,
II 、使用指数函数拟合时
现将指数函数进行变形： 将y C =,x t
=代入2
)(c t b ae
C --=得:2
)(c x b ae
y --=
对上式左右取对数得:2
22ln ln bx bcx bc a y -+-=
令
b
bc bc a y z -==-==21202ln ln ααα，，，
则可得多项式:
2
210x x z ααα++=
（3）Matlab 源程序
clear; %清除工作空间变量 clc; %清除命令窗口命令 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(1:m))/m;
fprintf('一天的平均气温为T=%f\n',T); %求一天的平均气温 %% 二次、三次、四次函数的最小二乘近似
h=input('请输入拟合多项式的最高项次数='); %根据给定数据点生成矩阵G n=h+1; G=[];
for j=0:(n-1)
g=x.^j; %g(x)按列排列
G=vertcat(G,g); %g垂直连接G
end
G=G'; %转置得到矩阵G
for i=1:m %将数据y作为G的最后一列（n+1列）
G(i,n+1)=y(i);
end
G;
for k=1:n %形成矩阵Q(k)
if G(k,k)>0;
sgn=1;
elseif G(k,k)==0;
sgn=0;
else sgn=-1;
end
sgm=-sgn*sqrt(sum(G(k:m,k).^2));
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; %变换Gk-1到Gk
for j=k+1:n+1
t=sum(w(k:m)*G(k:m,j))/bt;
for i=k:m;
G(i,j)=G(i,j)+t*w(i);
end
end
end
A (n)=G(n,n+1)/G(n,n); %解三角方程求系数A for i=n-1:-1:1
A (i)=(G(i,n+1)-sum(G(i,i+1:n).*A (i+1:n)))/G(i,i);
end
e=sum(G(n+1:m,n+1).^2); %计算误差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(1)
plot(x,y,'k*',t,Y,'r-'); %绘制拟合多项式函数图形
xlabel('时刻'); %标注坐标轴含义ylabel('平均气温');
title([num2str(n-1),'次函数的最小二乘曲线']);
grid;
%% 指数函数的最小二乘近似
yy=log(y);
n=3;
G=[];
GG=[];
for j=0:(n-1)
g=x.^j; %g(x)按列排列
G=vertcat(G,g); %g垂直连接G
gg=t.^j; %g(x)按列排列
GG=vertcat(GG,gg); %g垂直连接G
end
G=G'; %转置得到矩阵G
for i=1:m %将数据y作为G的最后一列（n+1列）
G(i,n+1)=yy(i);
end
G;
for k=1:n %形成矩阵Q(k)
if G(k,k)>0;
sgn=1;
elseif G(k,k)==0;
sgn=0;
else sgn=-1;
end
sgm=-sgn*sqrt(sum(G(k:m,k).^2));
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; %变换Gk-1到Gk
for j=k+1:n+1
t=sum(w(k:m)*G(k:m,j))/bt;
for i=k:m;
G(i,j)=G(i,j)+t*w(i);
end
end
end
A(n)=G(n,n+1)/G(n,n); %解三角方程求系数A
for i=n-1:-1:1
A (i)=(G(i,n+1)-sum(G(i,i+1:n).*A (i+1:n)))/G(i,i);
end
b=-A(3);
c=A(2)/(2*b);
a=exp(A(1)+b*(c^2));
G(n+1:m,n+1)=exp(sum(G(n+1:m,n+1).^2));
e=sum(G(n+1:m,n+1).^2); %计算误差e
fprintf('\n指数函数的系数是：a=%f,b=%f,c=%f',a,b,c);%输出系数及误差e fprintf('\n使用指数函数拟合的误差是：%f',e);
t=0:0.05:24;
YY=a.*exp(-b.*(t-c).^2);
figure(2)
plot(x,y,'k*',t,YY,'r-'); %绘制拟合指数函数图形xlabel('时刻'); %标注坐标轴含义
ylabel('平均气温');
title(['指数函数的最小二乘曲线']);
grid;
（4）结果与分析
a、二次函数：
一天的平均气温为： 21.2000
2次函数的系数: 8.3063 2.6064 -0.0938
使用2次函数拟合的误差是：280.339547
二次函数的最小二乘曲线如下图所示：
b、三次函数：
一天的平均气温为： 21.2000
3次函数的系数: 13.3880 -0.2273 0.2075 -0.0084 使用3次函数拟合的误差是： 131.061822
三次函数的最小二乘曲线如下图所示：
c、四次函数：
一天的平均气温为： 21.2000
4次函数的系数: 16.7939 -3.7050 0.8909 -0.0532 0.0009 使用4次函数拟合的误差是：59.04118
四次函数的最小二乘曲线如下图所示：
d、指数函数：
一天的平均气温为： 21.2000
指数函数的系数是：a=26.160286，b=0.004442，c=14.081900 使用指数函数拟合的误差是：57.034644
指数函数的最小二乘曲线如下图所示：
通过上述几种拟合可以发现，多项式的次数越高，计算拟合的效果越好，误差越小，说明结果越准确；同时，指数多项式拟合的次数虽然不高，但误差最小，说明结果最准确。
4.设计算法，求出非线性方程
52645200x x -+=的所有根，并使误差不超过4
10-。
（1）算法思想
首先，研究函数的形态，确定根的范围；通过剖分区间的方法确定根的位置，然后利用二分法的基本原理进行求解，找到满足精度要求的解。
二分法是产生一串区间，使新区间)
1(+k I
是旧区间)
(k I
的一个子区间，其
长度是)(k I 的一半，且有一个端点是)
(k I 的一个端点。由区间]
,[)1()()
(+=k k k x x I
确定区间)
1(+k I
的方法是计算区间)
(k I
的中点
)
(2
1)1()()
2(+++=
k k k x x x
若0)()()2()(<+k k x f x f ，则取],[)2()()1(++=k k k x x I ，否则取
],[)1()2()1(+++=k k k x x I ，重复这一过程即可。显然，每次迭代使区间长度减小一半，故二分法总是收敛的。
（2）算法结构
1.
)0()(f x f ?；
1)1()(f x f ? 2.If 0
10>f f then stop
3.If
2
0||ε 4.If 21||ε 1(x 作为根； stop 5.x x x ?+)(21)1() 0( 6.If ||||) 1(1)1(x x x ε<-then 输出x 作为根； stop 7.f x f ?)( 8.If 2||ε 9.If 01 9.1) 0(x x ?； f f ? else 9.2 ) 1(x x ?；1f f ? 10.go to 5 （3）Matlab 源程序 x=-100:100; y=6*(x.^5)-45*(x.^2)+20; %非线性方程组的表达式 g=[]; for i=-100:1:100 %确定根所在的区间 k=i+1; if (y(x==i).*y(x==k) f=6*x^5-45*x^2+20; n=length(g); %确定根的个数 for j=1:n x0=g(j); %求根区间左端点 x1=g(j)+1; %求根区间右端点 while (x1-x0)>=10^(-4) if subs(f,x,x0)*subs(f,x,(x0+x1)/2)>eps x0=(x0+x1)/2; else x1=(x0+x1)/2; end end root=x0 %输出方程的根 end （4）结果与分析 该非线性方程组有三个实根，分别为1.8708，0.6812，-0.6545，且满足误差要求。 5.编写程序实现大规模方程组的列主元高斯消去法程序，并对所附的方程组进行求解。针对本专业中所碰到的实际问题，提炼一个使用方程组进行求解的例子，并对求解过程进行分析、求解。 （1）算法思想 高斯消去法是利用现行方程组初等变换中的一种变换，即用一个不为零的数乘一个方程后加只另一个方程，使方程组变成同解的上三角方程组，然后再自下而上对上三角方程组求解。 列主元消去法是当高斯消元到第k步时，从k列的kk a以下（包括kk a）的各元素中选出绝对值最大的，然后通过行交换将其交换到kk a的位置上。交换系数矩阵中的两行（包括常数项），只相当于两个方程的位置交换了，因此，列选主元不影响求解的结果。 程序的核心就是高斯列主元消去法。根据教材提供的算法，编写列主元消去法的子函数与适应于超大规模超出系统内存的方程组的改编程序。同时，在Gauss消去过程中，适当交换方程的顺序对保证消去过程能顺利进行及计算解 的精确度都是有必要的，交换方程的原则是使 ) , ,1 , ()1 (n k k i k ik + = - α 中，绝对 值最大的一个换到（k,k)位置而成为第k步消去的主元，这就是列主元Gauss 消去法。 （2）算法结构 1、数据文件的文件名为：文件名+.dat 2、数据文件中的数据为二进制记录结构，分为以下四个部分： (1)文件头部分，其结构： typedef struct { long int id; long int ver; long int n; } 其中：id:为该数据文件的标识，值为0xF1E1D1A0，即为：十六进制的F1E1D1A0 ver:为数据文件的版本号，值为16进制数据, n: (2)文件头2:此部分说明为条状矩阵的上下带宽,结构： typedef struct { long int q; // 为上带宽 long int p; // 为下带宽 } (3)系数矩阵 a.如存贮格式非为压缩方式，则按行方式存贮系数矩阵中的每一个元素，个数为n*n，类型为float型； b.如果存贮格式是压缩方式,则按行方式存贮,每行中只存放上下带宽内的非零元素,即，每行中存贮的最多元素为p+q+1个。 (4)右端系数 按顺序存贮右端系数的每个元素，个数为n个，类型为float型 3、二进制文件的读取： f=fopen('fun003.dat','r'); %打开文件，.dat文件放在 m 文件同一目录下， a=fread(f,3,'uint') %读取头文件，3-读取前 3 个，若读取压缩格式的，头文件 为 5 个 b=fread(f,inf,'float'); %读取剩下的文件，float 型 id=dec2hex(a(1)); ver=dec2hex(a(2)); %这两句是进行进制转换，读取 id 与ver 1. A 的阶数n ? 2. For 1,,3,2,1-=n k 2.1找满足 k ik k i k k μααμ的下标≥=max 2.2For n j ,,2,1 = 2.2.1 j kj k μ αα? 2.3 k k μ ββ? 2.4 For n k k i ,,2,1 ++= 2.4.1 ik kk ik ααα? 2.4.2 For n k k j ,,2,1 ++= 2.4.2.1 ij kj ik ij αααα?- 2.4.3 i k ik i ββαβ?- n nn n x ?αβ/ For 1,,2,1 --=n n k ∑+=?- n k j k kk j kj k x x 1 /)(αα β （3）Matlab 源程序 clear; %清除工作空间变量 clc; %清除命令窗口命令 %% 读取系数矩阵 [f,p]=uigetfile('*.dat','选择数据文件'); %读取数据文件 num=5; %输入系数矩阵文件头的个数 name=strcat(p,f); file=fopen(name,'r'); head=fread(file,num,'uint'); %读取二进制头文件 id=dec2hex(head(1)); %读取标识符 fprintf('文件标识符为');