西安交大计算方法b大作业课件.docx
《西安交大计算方法b大作业课件.docx》由会员分享，可在线阅读，更多相关《西安交大计算方法b大作业课件.docx（21页珍藏版）》请在冰豆网上搜索。
西安交大计算方法b大作业课件
《计算方法B》上机
实验报告
学院：
机械工程学院
班级:
姓名：
学号：
2015年12月22日
1.计算以下和式：
，要求：
（1）若保留11个有效数字，给出计算结果，并评价计算的算法；
（2）若要保留30个有效数字，则又将如何进行计算。
实现思想：
以上问题出现了近似数相减的问题，为了减小误差，可分别求得减数之和以及被减数之和，最后将两者相减。
另外，减数与被减数求和均为同号计算，按照绝对值递增顺序相加可减小舍入误差。
此题中对有效数字有要求，因而计算时首先需要根据有效数字位数计算得出迭代次数，以保证计算值的精度。
源程序：
m=input（'输入有效数字个数m='）;
s0=1;s1=0;s2=0;n=0;
%判断迭代次数
whiles0>=0.5*10^-（m-1）
s0=4/（16^n*（8*n+1））-2/（16^n*（8*n+4））-1/（16^n*（8*n+5））-1/（16^n*（8*n+6））;
n=n+1;
end
%分别求解各项并求和
fork=n-1:
-1:
0
a1=4/（16^k*（8*k+1））;
a2=2/（16^k*（8*k+4））;
a3=1/（16^k*（8*k+5））;
a4=1/（16^k*（8*k+6））;
s1=a1+s1;
s2=a4+a3+a2+s2;
end
S=vpa（s1-s2,m）
实验结果：
11位有效数字计算结果如图1所示；30为有效数字计算结果如图2所示。
图1.11位有效数字计算结果图2.30为有效数字计算结果
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）估算所需光缆长度的近似值，并作出铺设河底光缆的曲线图；
算法思想：
由于题中所给点数为20，若采用高次多项式插值将产生很大的误差，所以拉格朗日或牛顿并不适用。
题中光缆为柔性，可光滑铺设于水底，鉴于此特性，采用三次样条插值插值法较为合适。
算法结构：
三次样条算法结构见《计算方法教程》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;
%%%带状矩阵求解（追赶法）
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
%%%绘图
plot（xx,s,'r-','linewidth',1.5）
disp（['¹光缆长度为ª',num2str（l（nn+1））,'Ã×']）
曲线图如图2-1所示，计算光缆长度如图2-2所示。
图2-1光缆插值曲线图
图2-1光缆计算长度显示
3.假定某天的气温变化记录如下表所示，试用数据拟合的方法找出这一天的气温变化的规律;试计算这一天的平均气温,并试估计误差。
时刻
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
实现思想：
此题中所给数据点数目较多，采用拉格朗日插值法或者牛顿插值法需要很高次的多项式，计算困难，误差大；采用样条插值计算量虽然不大，但是存放参数Mi的量很大，且没有一个统一的数学公式来表示，也不是很方便。
所以可考虑用最小二乘法进行拟合。
计算过程中，分别使用二次函数、三次函数以及四次函数，计算其相应的系数，估算误差并作图比较各个函数之间的区别。
算法结构：
（参考课本P123）
1.1[形成矩阵Qk]
1.2[变换Gk-1到Gk]
2.[求解三角方程]
3.[计算误差]
源代码：
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=length（x）;
n=input（'请输入函数的次数 '）;
plot（x,y,'k.',x,y,'-'）
grid;
hold on;
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
for j=1:
n
for i=1:
m
G（i,j）=x（1,i）^（j-1）;
end
end %%%建立矩阵Qk
for k=1:
n
for i=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来存放ω
for j=k+1:
m
w（j）=G（j,k）;
end
b（k）=c（k）*w（k）;
%%%变换矩阵Gk-1到Gk
G（k,k）=c（k）;
for j=k+1:
n+1
q=0;
for i=k:
m
q=w（i）*G（i,j）+q;
end
s=q/b（k）;
for i=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）;
for i=n-1:
（-1）:
1
for j=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
a
%%%回代过程
p=zeros（1,m）;
for j=1:
m
for i=1:
n
p（j）=p（j）+a（i）*x（j）^（i-1）;
end
end
plot（x,p,'r*',x,p,'-'）;
E2=0;%用E2来存放误差
%%%误差求解
for i=n+1:
m
E2=G（i,n+1）^2+E2;
end
E2=E2^0.5;
disp（'误差为'）;
disp（E2）;
t=0;
for i=1:
m
t=t+p（i）;
end
t=t/m;
%%%平均温度
disp（['平均温度为',num2str（t）,'℃']）
实验结果：
二次函数拟合，结果如下图所示
图3-1二次函数拟合结果
三次函数拟合，结果如下图所示
图3-2三次函数拟合结果
四次函数拟合，结果如下图所示
图3-3四次函数拟合结果
结果对比：
将二次函数、三次函数和四次函数拟合结果绘制在同一个坐标内，如图3-4所示。
其计算误差结果见表3-1所示。
图3-4拟合结果对比分析
4.设计算法，求出非线性方程
的所有实根，并使误差不超过
。
算法思想：
本题可采用牛顿法迭代求解,令
得带格式为
根据函数图像可以找出根的大致分布区间,带入不同的初值即可解出不同的根.
源代码:
functiony=f2（x）
y=6*x.^5-45*x.^2+20;%定义原函数
functiony=f3（x）
y=30*x^4-90*x;%定义原函数倒数
i=-5:
0.1:
5;
y=f2（i）;
plot（i,y）
holdon
plot（i,0,'-'）%画出原函数图像
%%Newton法求根
x1=input（'输入初值'）;
e=10^（-4）;%误差设定
Nmax=1000;%迭代最大次数限定
forn=1:
Nmax
f0=f2（x1）;
ifabs（f2（x1））fprintf（'输出的f（x）已经足够小'）;
x=x1;
break
else
F0=f3（x1）;
x=x1-f0/F0;
ifabs（x-x1）break
else
x1=x;
end
end
end
fprintf（'输出方程的根x=%2f',x）
计算结果:
函数图像如图4-1所示。
计算结果分别见图4-2所示。
图4-1函数图像
图4-2计算结果
根据带入不同的初值，可以求出不同的根，有图4-2可以看出，原函数的根大约有三个，分别是-0.654542、0.681174、1.870799。
5.线性方程组求解。
（1）编写程序实现大规模方程组的高斯消去法程序，并对所附的方程组进行求解。
所附方程组的类型为对角占优的带状方程组。
（2）针对本专业中所碰到的实际问题，提炼一个使用方程组进行求解的例子，并对求解过程进行分析、求解。
算法思想:
高斯消去法是利用现行方程组初等变换中的一种变换,将一个不为零的数乘到一个方程后加到另一个方程,使方程组变成同解的上三角方程组,然后再自下而上对上三角方程组求解。
算法结构：
1.读取二进制文件，存入计算矩阵
2.对矩阵进行初等变换，然后求解（详见计算方法教程第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’）
求解结果
对数据文件dat51求解，结果如下：
文件标识符为
id=F1E1D1A0
文件版本号为
ver=102
矩阵A的阶数
n=15
矩阵A的上带宽
q=3
矩阵A的下带宽
p=3
方程组的的解为：
1.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000
对数据文件dat52求解，结果如下：
文件标识符为
id=F1E1D1A0
文件版本号为
ver=202
矩阵A的阶数
n=20
矩阵A的上带宽
q=5
矩阵A的下带宽
p=5
方程组的的解为：
1.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000
对数据文件dat53求解，结果如下：
文件标识符为
id=F1E1D1A0
文件版本号为
ver=102
矩阵A的阶数
n=2160
矩阵A的上带宽
q=5
矩阵A的下带宽
p=5
方程组的解截图如图5-1所示（由于矩阵阶数较大，计算结果未能截取完整）。
图5-1数据文件dat53求解结果
对数据文件dat53求解，结果如下：
文件标识符为
id=F1E1D1A0
文件版本号为
ver=202
矩阵A的阶数
n=43240
矩阵A的上带宽
q=4
矩阵A的下带宽
p=4
方程组的解截图如图5-2所示（由于矩阵阶数较大，计算结果未能截取完整）。
图5-2数据文件dat54求解结果
升级会员 