计算方法B上机报告.docx
- 文档编号:25131715
- 上传时间:2023-06-05
- 格式:DOCX
- 页数:28
- 大小:460.79KB
计算方法B上机报告.docx
《计算方法B上机报告.docx》由会员分享,可在线阅读,更多相关《计算方法B上机报告.docx(28页珍藏版)》请在冰豆网上搜索。
计算方法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℃。
函数形式为
四种函数图像的对比:
红色表示指数函数,绿色表示二次函数,青色表示三次函数,紫色表示四次函数,蓝色的折线表示数据点。
总结心得:
同上面的结果可以看出,随着次数的增加误差在逐渐减小,拟合曲线也更加接近数据点。
因此,在使用最小二乘进行数据拟合时应尽量使用较高次数的拟合函数。
通过对比图像,可以发现二次函数和指数型二次函数图像比较接近,那是因为指数型二次函数是对数据点求对数进行拟合的,因此二者比较接近。
给定的数据点在中间的位置时波动较大,造成了拟合曲线的误差也较大,但并没有影响到整体的规律。
在进行最小二次编程中,使用正交化方法中间变量较多容易造成混乱,计算量也较大,但是比使用法方程方法更稳定。
本题在编程过程中有一定的适用性,只需要对数据点进行更改就可以拟合出给定次数的多项式,为以后的数据拟合提供了方便。
4.设计算法,求出非线性方程
的所有根,并使误差不超过
。
问题分析:
本题采用牛顿迭代法计算方程根。
可令
,其迭代格式为
要使牛顿迭代法收敛,则必须寻找根的合理区间[a,b],使得在该区间内
即有根。
在选定的区间内函数
的一二阶导数
、
均不变号。
选定一个初值
,使
则牛顿迭代法收敛。
程序实现:
clear;clc;
%%%观察根的大致位置
a=-10:
0.1:
10;
fori=1:
201
y(i)=6*a(i)^5-45*a(i)^2+20;
end
plot(a,y);
grid;
k=1;
y(202)=0;
%%%使用二分法来确定迭代初值
fori=1:
201
ify(i)*y(i+1)<0
b(k)=a(i);
k=k+1;
continue
end
end
x=b;
n=length(x);
%%%牛顿迭代法计算方程根
dlt=1e-4;
fori=1:
n
forj=1:
100
X(i,j)=x(i)-(6*x(i)^5-45*x(i)^2+20)/(30*x(i)^4-90*x(i));
c=abs(X(i,j)-x(i));
ifc break end x(i)=X(i,j); end end disp('计算结果为: '); x 运行结果: 函数图像如下,可以观察根的大致位置。 计算结果: 总结心得: 在该题的计算中,首先画出了函数的图像可以初步观察出根的大致位置,有效缩减了求根区间。 为防止迭代进入死循环,还设置了最高迭代次数为100次,结果显示方程的三个根均在4次以内就可以迭代出,可以看出牛顿迭代法的求解速度相当快。 5.编写程序实现大规模方程组的列主元高斯消去法程序,并对所附的方程组进行求解。 针对本专业中所碰到的实际问题,提炼一个使用方程组进行求解的例子,并对求解过程进行分析、求解。 问题分析: 高斯消去法在整个消去过程中, 均不为零,但是若在某一步中,对角线上元的绝对值比其下方的某些元素 的绝对值小得多,即 那么,相应的数 将会是绝对值很大的数,这可能引起上溢而中断计算,即使不发生上溢,也会带来很大的舍入误差。 而列主元消去法的原则使 中,绝对值最大的一个换到(k,k)位置,这样就有效的避免了大数除以小数而造成的上溢现象。 列主元消去法的主要核心是查找最大的 并交换到 。 程序实现: (1)非压缩格式dat141.dat;dat142.dat clc;clear; %%%读取数据 fid=fopen('dat141.dat','r'); a=fread(fid,5,'long',0); n=a(3); a=fread(fid,[n,n+1],'float'); fori=1: n forj=1: n A(i,j)=a(j,i); end b(i)=a(i,n+1); end %%%查找主元并交换位置 fork=1: n form=k: n ifabs(A(m,k))~=max(abs(A(: k))) continue else r=A(k,: );A(k,: )=A(m,: );A(m,: )=r; r=b(k);b(k)=b(m);b(m)=r; continue end end fori=k+1: n m=A(i,k)/A(k,k); forj=1: n A(i,j)=A(i,j)-m*A(k,j); end b(i)=b(i)-m*b(k); end end %%%回代过程 x(n)=b(n)/A(n,n); b(n)=x(n); fork=n-1: (-1): 1 forj=k+1: n b(k)=b(k)-A(k,j)*x(j); end x(k)=b(k)/A(k,k); b(k)=x(k); end x=x' (2)压缩格式dat143.dat clear;clc; formatshort %%%读取数据 fid=fopen('dat143.dat','r'); d=fread(fid,5,'long',0); n=d(3);%矩阵维数 p=d(5);%上带宽 q=d(4);%下带宽 a=fread(fid,[p+q+1,n],'float'); b=fread(fid,[n,1],'float'); a=a'; d=b'; %%%消去过程 fork=1: n-1 fori=k+1: min(k+p,n) a(i,k-i+p+1)=a(i,k-i+p+1)/a(k,p+1); forj=k+1: min(k+q,n) a(i,j-i+p+1)=a(i,j-i+p+1)-a(i,k-i+p+1)*a(k,j-k+p+1); end b(i)=b(i)-b(k)*a(i,k-i+p+1); end end %%%回代过程 x(n)=b(n)/a(n,p+1); fork=n-1: -1: 1 sum=0; forj=k+1: min(k+q,n) sum=sum+a(k,j-k+p+1)*x(j); end x(k)=(b(k)-sum)/a(k,p+1); end x 运算结果: 1)非压缩格式dat141.dat;dat142.dat的运算结果如下 (2)压缩格式dat143.dat 总结心得: 所给数据的方程的根都相等,其中dat141.dat,为10阶方程,根均为1.0000;dat142.dat为2050阶,其根均为1.9600;dat143.dat为43512阶,其根均为3.1413。 在求解非压缩格式的方程根时,采用一般列主元高斯消去法,计算量大,耗时长。 在求解压缩格式的方程根时,由于阶数较大,运算时间太长,对程序进行了改进,只对非零元素进行计算,运算步数大量减少,速度得到了极大的提高。 本专业算例: 一底边绝热的正方形导热物体,内无热源,其余三边温度如下,求解该正方形内节点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。 运算结果和正确结果吻合。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 计算方法 上机 报告
![提示](https://static.bdocx.com/images/bang_tan.gif)