计算方法B上机报告.docx
《计算方法B上机报告.docx》由会员分享，可在线阅读，更多相关《计算方法B上机报告.docx（24页珍藏版）》请在冰豆网上搜索。
计算方法B上机报告
计算方法B上机实习报告
计算方法B上机实习报告
1.对以下和式计算：
，要求：
（1）若只需保留11个有效数字，该如何进行计算；
（2）若要保留30个有效数字，则又将如何进行计算；
问题分析：
在该题中S的每一项存在两个相近的数相减的问题，因此为了避免有效数字损失，最好是改变运算顺序，分别将正数和负数相加，然后再将其和相加。
另外，sn中有多个负数相加，可以按照绝对值递增的顺序求和，以减少舍入误差的影响。
同时，为了避免大数吃小数的问题，本题先计算出保留目标有效数字所需
要的迭代次数，然后采用倒序相加的方法实现。
程序实现：
clear;clc;
m=input（'请输入要保留的有效数字位数：
'）;
s1=0;
s2=0;
k=0;
s=1;
%%%%判断多需要的迭代次数
whiles>=0.5*10^-（m-1）s=4/（16^k*（8*k+1））-（2/（16^k*（8*k+4））+1/（16^k*（8*k+5））+1/（16^k*（8*k+6）））;
k=k+1;
end
%%%%正负数分别按照绝对值递增的顺序倒序相加
forn=（k-1）:
-1:
0
a1=4/（16^n*（8*n+1））;
a2=2/（16^n*（8*n+4））;
a3=1/（16^n*（8*n+5））;
a4=1/（16^n*（8*n+6））;
s1=a1+s1;
s2=a4+a3+a2+s2;
end
S=s1-s2;
S=vpa（S,m）
运算结果：
总结心得：
在计算求和问题中，应特别注意相近数相减的问题，这样会造成有效数字灾难性的损失。
另外在两个数量级相差较大的数字相加减时，较小数的有效数字会被丧失，因此要按照从小到大的顺序相加。
在上题计算中分别对正负相采用倒序相加，这样就有效的避免了“大数吃小数”的问题。
2.某通信公司在一次施工中，需要在水面宽度为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）预测所需光缆长度的近似值，并作出铺设河底光缆的曲线图；
问题分析：
本题的主要目的是对测量数据进行拟合，同时对拟合曲线进行线积分即可得到河底光缆长度的近似值。
由于数值点较多时，使用拉格朗日差值多项式会出现龙格现象。
为了将所有的数据点都用上，采用分段差插法，本题使用三次样条插值。
算法思想：
样条函数在每个子区间上是三次多项式，它的二阶导数必是一次多项式。
若用记在处的二阶导数。
则在区间上
式中
（1）
对上式进行两次积分得
（2）
它的一阶导数为
（3）
满足连续性条件，即
上式和（3）式得
（4）
用差商记号，并记
（4）式可以写成
方程组可以写成如下形式
自然样条插值条件为
在估计河底光缆长度时使用第一类线积分
程序实现：
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
%%%绘图
plot（xx,s,'r-','linewidth',1.5）
grid
disp（['所需光缆长度为',num2str（l（nn+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
问题分析：
对一天的气温进行数据拟合，可以考虑使用最小二乘的二次函数、三次函数、四次函数以及指数函数。
问题的难度是求解各种拟合函数的系数。
利用法方程求解最小二乘系数时，方程的解不够稳定，本题利用正交化方法。
程序实现：
clear;clc;
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
a
%%%拟合后的回代过程
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）,'℃']）
运行结果：
二次函数时，系数a1=8.3063,a2=2.6064,a3=-0.0938。
误差为16.7433，平均温度为21.2℃。
函数形式为
三次函数时，系数a1=13.3880,a2=-0.2273,a3=0.2075,a4=-0.0084,误差为11.4482，平均温度为21.2℃。
函数形式
四次函数时，系数a1=16.7939,a2=-3.7050,a3=0.8909,a4=-0.0532,a5=0.0009,误差为7.6838，平均温度为21.2℃。
函数形式为
当为指数函数时，假定指数函数的形式为。
只需将y值求对数，其主体部分不变，求得的系数为a1=2.3835,a2=0.1251,a3=-0.0044,误差为14.6867，平均温度为20.9399℃。
函数形式为
四种函数图像的对比：
红色表示指数函数，绿色表示二次函数，青色表示三次函数，紫色表示四次函数，蓝色的折线表示数据点。
总结心得：
同上
升级会员 