计算方法B上机报告.docx
《计算方法B上机报告.docx》由会员分享，可在线阅读，更多相关《计算方法B上机报告.docx（24页珍藏版）》请在冰豆网上搜索。
计算方法B上机报告
计算方法B上机报告
第1题
某通信公司在一次施工中，需要在水面宽度为20米的河沟底部沿直线走向铺设一条沟底光缆。
在铺设光缆之前需要对沟底的地形进行初步探测，从而估计所需光缆的长度，为工程预算提供依据。
已探测到一组等分点位置的深度数据（单位:
米）如下表所示：
分点
0
1
2
3
4
5
6
深度
9.01
8.96
7.96
7.97
8.02
9.05
10.13
分点
7
8
9
10
11
12
13
深度
11.18
12.26
13.28
13.32
12.61
11.29
10.22
分点
14
15
16
17
18
19
20
深度
9.15
7.90
7.95
8.86
9.81
10.80
10.93
（1）请用合适的曲线拟合所测数据点；
（2）估算所需光缆长度的近似值，并作出铺设河底光缆的曲线图；
问题分析和算法思想：
本题的主要目的是对21个测量数据进行拟合，同时对拟合曲线进行线积分即可得到河底光缆长度的近似值，可以用的插值方法很多：
多项式插值、Lagrange插值、Newton插值、三次样条插值等。
由于数值点较多时，采用高次多项式插值将产生很大的误差，用拉格朗日插值多项式会出现龙格现象。
故为了将所有的数据点都用上，且题中光缆为柔性，可光滑铺设于水底，鉴于此特性，采用三次样条插值的方法较为合适。
计算光缆长度近似值，只需将每两点之间的距离算出，然后依次相加，所得的折线长度，即为光缆长度的近似值。
光缆长度计算公式：
算法结构：
三次样条算法结构见《计算方法教程》P110。
源程序：
clear;clc;
x=0:
20;
y=[9.018.967.967.978.029.0510.1311.1812.2613.2813.3212.6111.2910.229.157.907.958.869.8110.8010.93];
d=y;
plot（x,y,'k.','markersize',15）
holdon
%%%计算二阶差商
fork=1:
2
fori=21:
-1:
（k+1）
d（i）=（d（i）-d（i-1））/（x（i）-x（i-k））;
end
end
%%%假定d的边界条件，采用自然三次样条
fori=2:
20
d（i）=6*d（i+1）;
end
d
（1）=0;
d（21）=0;
%%%追赶法求解带状矩阵的m值
a=0.5*ones（1,21）;
b=2*ones（1,21）;
c=0.5*ones（1,21）;
a
（1）=0;c（21）=0;
u=ones（1,21）;
u
（1）=b
（1）;
r=c;
yy
（1）=d
（1）;
%%%追的过程
fork=2:
21
l（k）=a（k）/u（k-1）;
u（k）=b（k）-l（k）*r（k-1）;
yy（k）=d（k）-l（k）*yy（k-1）;
end
%%%赶的过程
m（21）=yy（21）/u（21）;
fork=20:
-1:
1
m（k）=（yy（k）-r（k）*m（k+1））/u（k）;
end
%%%利用插值点画出拟合曲线
k=1;
nn=100;
xx=linspace（0,20,nn）;
l=0;
forj=1:
nn
fori=2:
20
ifxx（j）<=x（i）
k=i;
break;
else
k=i+1;
end
end
h=1;
xbar=x（k）-xx（j）;
xmao=xx（j）-x（k-1）;
s（j）=（m（k-1）*xbar^3/6+m（k）*xmao^3/6+（y（k-1）-m（k-1）*h^2/6）*xbar+（y（k）-m（k）*h^2/6）*xmao）/h;
sp（j）=-m（k-1）*（x（k）-xx（j））^2/（2*h）+m（k）*（xx（j）-x（k-1））^2/（2*h）+（y（k）-y（k-1））/h-（m（k）-m（k-1））*h/6;
l（j+1）=（1+sp（j）^2）^0.5*（20/nn）+l（j）;%利用第一类线积分求河底光缆的长度
end
%%%绘图
title（'光缆插值曲线'）
xlabel（'分点'）
ylabel（'深度'）
plot（xx,s,'r-','linewidth',1.5）
grid
disp（['所需光缆长度为',num2str（l（nn+1））,'米']）
运行结果：
第2题
假定某天的气温变化记录如下表所示，试用数据拟合的方法找出这一天的气温变化的规律;试计算这一天的平均气温,并试估计误差。
时刻
0
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
问题分析及算法思路：
由于题中所给数据点有25个，此时采用多项式插值的方法做数据近似需要用较高次的多项式，这不仅给计算带来困难，而且有余项带来的误差很大；若采用样条插值计算量虽然不大，但是参数的个数很多，而且没有一个统一的数学公式来表示，对计算带来了很大的麻烦。
所以可考虑使用最小二乘法进行拟合。
在本题中，采用最小二乘法找出这一天的气温变化的规律，依次使用二次函数、三次函数以及四次函数进行拟合，计算其相应的系数，估算误差，并作图比较三个函数之间的区别。
算法结构：
最小二乘法算法结构见《计算方法教程》P123。
源程序：
x=0:
24;
y=[15141414141516182020232528313431292725242220181716];
m=length（x）;
n=input（'请输入函数的次数'）;
plot（x,y,'k.',x,y,'-'）
grid;
holdon;
n=n+1;
G=zeros（m,n+1）;
G（:
n+1）=y';
c=zeros（1,n）;%建立c来存放σ
q=0;
f=0;
b=zeros（1,m）;%建立b用来存放β
%%%形成矩阵G
forj=1:
n
fori=1:
m
G（i,j）=x（1,i）^（j-1）;
end
end
%%%建立矩阵Qk
fork=1:
n
fori=k:
m
c（k）=G（i,k）^2+c（k）;
end
c（k）=-sign（G（k,k））*（c（k）^0.5）;
w（k）=G（k,k）-c（k）;%建立w来存放ω
forj=k+1:
m
w（j）=G（j,k）;
end
b（k）=c（k）*w（k）;
%%%变换矩阵Gk-1到Gk
G（k,k）=c（k）;
forj=k+1:
n+1
q=0;
fori=k:
m
q=w（i）*G（i,j）+q;
end
s=q/b（k）;
fori=k:
m
G（i,j）=s*w（i）+G（i,j）;
end
end
end
%%%求解三角方程Rx=h1
a（n）=G（n,n+1）/G（n,n）;
fori=n-1:
（-1）:
1
forj=i+1:
n
f=G（i,j）*a（j）+f;
end
a（i）=（G（i,n+1）-f）/G（i,i）;%%%a（i）存放各级系数
f=0;
end
%%%拟合后的回代过程
p=zeros（1,m）;
forj=1:
m
fori=1:
n
p（j）=p（j）+a（i）*x（j）^（i-1）;
end
end
plot（x,p,'r*',x,p,'-'）;
E2=0;%用E2来存放误差
%%%误差求解
fori=n+1:
m
E2=G（i,n+1）^2+E2;
end
E2=E2^0.5;
disp（'误差为'）;
disp（E2）;
t=0
fori=1:
m
t=t+p（i）;
end
t=t/m;%%%平均温度
disp（['平均温度为',num2str（t）,'℃']）
运行结果：
二次函数时，系数a0=8.3063，a1=2.6064，a2=-0.0938，误差为16.7433，平均温度为21.2℃。
函数形式为：
三次函数时，系数a0=13.3880，a1=-0.2273，a2=0.2075，a3=-0.0084，误差为11.4482，平均温度为21.2℃。
函数形式：
四次函数时，系数a0=16.7939，a1=-3.7050，a2=0.8909，a3=-0.0532，a4=0.0009，误差为7.6838，平均温度为21.2℃。
函数形式为：
第3题
已知非线性方程
在[2,3]中有根。
请设计算法，求出该根，并使求出的根的误差不超过
。
问题分析以及算法思路：
本题采用复化梯形求积公式计算方程根，可令
。
算法结构：
牛顿迭代法算法结构见《计算方法教程》P196。
源程序：
clear;
clc;
x0=2;%设置求解区间
x1=3;
n=20;%复化梯形求积公式的步数
e=0.0004;
f0=T3Sub（x0,n）;
f1=T3Sub（x1,n）;
if（f0*f1）>0%判断在区间是否有
error（'Noanswer,becausef0*f1>0'）;
end
fori=1:
201
x=（x0+x1）/2;
ifabs（x1-x）disp（['方程的解',num2str（x）]）;
disp（['迭代步数',num2str（i）]）;
disp（['误差',num2str（abs（x1-x））]）;
break;
end
f=T3Sub（x,n）;
if（f*f0）<0
x1=x;
f1=f;
else
x0=x;
f0=f;
end
end
function[f]=T3Sub（x,n）
%复化梯形求积公式
h=pi/n;
f=0;
fori=1:
n
u=i*h;
f=f+cos（x*sin（u-h））+cos（x*sin（u））;
end
f=（h*f）/（2*pi）;
end
运行结果：
第4题
线性方程组求解
（1）编写程序实现大规模方程组的高斯消去法程序，并对所附的方程组进行求解。
所附方程组的类型为对角占优的带状方程组。
（2）针对本专业中所碰到的实际问题，提炼一个使用方程组进行求解的例子，并对求解过程进行分析、求解。
算法思想:
高斯消去法是利用现行方程组初等变换中的一种变换,将一个不为零的数乘到一个方程后加到另一个方程,使方程组变成同解的上三角方程组,然后再自下而上对上三角方程组求解。
算法结构：
1.读取二进制文件，存入计算矩阵
2.对矩阵进行初等变换，然后求解（计算方法教程高斯消去法及列主元高斯消去法算法）
源代码：
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（'文件标识符为：
'）;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=1:
n,
forj=1:
n,
a（i,j）=A（（i-1）*n+j）;%求系数矩阵a（i,j）
end
end
b=zeros（n,1）;
fori=1:
n,
b（i）=A（n*n+i）;
end
fork=1:
n-1,%列主元高斯消去法
m=k;
fori=k+1:
n,%寻找主元
ifabs（a（m,k））m=i;
end
end
ifa（m,k）==0%遇到条件终止
disp（'错误！
'）
return
end
forj=1:
n,%交换元素位置得主元
t=a（k,j）;
a（k,j）=a（m,j）;
a（m,j）=t;
t=b（k）;
b（k）=b（m）;
b（m）=t;
end
fori=k+1:
n,%计算l（i,k）并将其放到a（i,k）中
a（i,k）=a（i,k）/a（k,k）;
forj=k+1:
n
a（i,j）=a（i,j）-a（i,k）*a（k,j）;
end
b（i）=b（i）-a（i,k）*b（k）;
end
end
x=zeros（n,1）;%回代过程
x（n）=b（n）/a（n,n）;
fork=n-1:
-1:
1,
x（k）=（b（k）-sum（a（k,k+1:
n）*x（k+1:
n）））/a（k,k）;
end
end
%%对压缩带状矩阵进行求解
ifver=='202',%高斯消去法
m=p+q+1;
a=zeros（n,m）;
fori=1:
1:
n
forj=1:
1:
m
a（i,j）=A（（i-1）*m+j）;
end
end
b=zeros（n,1）;
fori=1:
1:
n
b（i）=A（n*m+i）;%求b（i）
end
fork=1:
1:
（n-1）%开始消去
ifa（k,（p+1））==0
disp（'错误！
'）;
break;
end
st1=n;
if（k+p）st1=k+p;
end
fori=（k+1）:
1:
st1
a（i,（k+p-i+1））=a（i,（k+p-i+1））/a（k,（p+1））;
forj=（k+1）:
1:
（k+q）
a（i,j+p-i+1）=a（i,j+p-i+1）-a（i,k+p-i+1）*a（k,j+p-k+1）;
end
b（i）=b（i）-a（i,k+p-i+1）*b（k）;
end
end
x=zeros（n,1）;%回代
x（n）=b（n）/a（n,p+1）;
sum=0;
fork=（n-1）:
-1:
1
sum=b（k）;
st2=n;
if（k+q）st2=k+q;
end
forj=（k+1）:
1:
st2
sum=sum-a（k,j+p-k+1）*x（j）;
end
x（k）=sum/a（k,p+1）;
sum=0;
end
end
disp（'方程组的的解为：
'）%输出解
disp（x）
运行结果：
首先是对测试文件dat20171/dat20172的求解结果，具体如下：
经过测试矩阵的初步测试，利用程序对接下来的两个数据文件进行处理，由于dat20173/dat20174解的个数较多，在这里只截取不分解作为示意，具体求解如下：
本专业算例
一底边绝热的形导热物体，无热源，其余三边温度如下，求解该形节点1、2、3、4的温度。
解：
导热物体无热源，物性为常数的控制方程可以写成
考虑到空间步长相同时，离散方程可以写成：
式中
从而有：
所有的节点方程可整理为
用列主元高斯消去法对上述矩阵进行求解，其程序如下：
clear;clc;
A=[-4110-70;1-401-50;10-31-15;011-3-10];
[m,n]=size（A）;
fori=1:
（m-1）
max=i;
fork=i+1:
m
ifabs（A（k,i））>abs（A（max,i））
max=k;
end
end
temp=A（i,i:
m）;A（i,i:
m）=A（max,i:
m）;A（max,i:
m）=temp;
forj=（i+1）:
m
A（j,:
）=A（j,:
）-A（i,:
）*A（j,i）/A（i,i）;
end
end
disp（'回代求解'）
x（m）=A（m,n）/A（m,m）;
fork=（m-1）:
-1:
1
x（k）=（A（k,n）-A（k,k+1:
m）*x（k+1:
m）'）/A（k,k）;
end
x
运算结果：
T1=28.7368，T2=24.2632，T3=20.6842，T4=18.3158。
运算结果和正确结果吻合。
升级会员 