计算方法综合实践.docx
《计算方法综合实践.docx》由会员分享,可在线阅读,更多相关《计算方法综合实践.docx(18页珍藏版)》请在冰豆网上搜索。
计算方法综合实践
1.应用自己熟悉的算法语言编写程序,使之尽可能具有通用性。
2.上机前充分准备,复习有关算法,写出计算步骤,反复检查,调试程序。
(注:
在练习本上写,不上交)
3.完成计算后写出实验报告,容包括:
算法步骤叙述,变量说明,程序清单,输出计算结果,结构分析和小结等。
(注:
具体题目具体分析,并不是所有的题目的实验报告都包含上述容!
)
4.独立完成,如有雷同,一律判为零分!
5.上机期间不允许做其他任何与课程设计无关的事情,否则被发现一次扣10分,被发现三次判为不及格!
非特殊情况,不能请假。
旷课3个半天及以上者,直接判为不及格。
一、基本技能训练
1、误差分析
实验1.2误差传播与算法稳定性
实验目的:
体会稳定性在选择算法中的地位。
误差扩的算法是不稳定的,是我们所不期望的;误差衰减的算法是稳定的,是我们努力寻求的,这是贯穿本课程的目标。
问题提出:
考虑一个简单的由积分定义的序列
显然
。
当
时,
而对于
时,利用分部积分易得
另一方面,我们有
实验容:
由以上递推关系,我们可得到计算序列
的两种方法。
(I)
(II)
symsnIn5In6In7;
In5=vpa((exp(-1)),5);
In6=vpa((exp(-1)),6);
In7=vpa((exp(-1)),7);
fprintf('%.5f%.6f%.7f\n',eval(In5),eval(In6),eval(In7));
forn=2:
10
In5=vpa((1-n*In5),5);
In6=vpa((1-n*In6),6);
In7=vpa((1-n*In7),7);
fprintf('%.5f%.6f%.7f\n',eval(In5),eval(In6),eval(In7));
end
五位六位七位
0.367880.3678790.3678794
0.264240.2642410.2642411
0.207280.2072770.2072766
0.170890.1708930.1708934
0.145530.1455330.1455329
0.126800.1268020.1268024
0.112380.1123840.1123835
0.100930.1009320.1009320
0.091610.0916120.0916123
0.083880.0838770.0838771
symsnEn;
En5=vpa(0,5);
En6=vpa(0,6);
En7=vpa(0,7);
fprintf('%.5f%.6f%.7f\n',eval(En5),eval(En6),eval(En7));
forn=10:
-1:
2
En5=vpa(((1-En5)/n),5);
En6=vpa(((1-En6)/n),6);
En7=vpa(((1-En7)/n),7);
fprintf('%.5f%.6f%.7f\n',eval(En5),eval(En6),eval(En7));
end
五位六位七位
0.000000.0000000.0000000
0.100000.1000000.1000000
0.100000.1000000.1000000
0.112500.1125000.1125000
0.126790.1267860.1267857
0.145540.1455360.1455357
0.170890.1708930.1708929
0.207280.2072770.2072768
0.264240.2642410.2642411
0.367880.3678790.3678795
(1)由所得数据可以了解,两种算法随着n的增大,误差越来越大,而第二种算法随着n的减小,数据越来越精确,且三种有效数字结果大致相同,所以第二种更精确。
(2)算法一E1的误差e1,由于En=1-n*En-1,n=2,3,4…..,通过推算可得误差
=
,n越大,其误差
就越大,那么最后
算出的结果是e1的n!
倍;
算法二
的误差
,由于
,推倒可得误差
=
,
比
缩小了(N-n)!
倍,则n减少,所以
就越小。
所以算法一越往后算,一步步的误差会使误差越大,而算法二由于是从后面递推回来,其误差会被缩小。
所以算法二更优。
(3)算法一中,n增大,误差增大,算法二中,n减小,误差减小。
所以当某一步发生误差
后随着n变大,算法一的误差越来越大,而算法二由于是往回推,所以误差变小。
(4)通过以上可一推出,算法一随着n变大,容易使误差越来越大,而算法二会使误差变小,所以算法二更加稳定。
2、求解非线性方程
3、插值
实验目的:
掌握Lagrange插值法和Newton插值法
问题提出:
,已知
的函数值表如下
x
0
0.1
0.2
0.3
0.4
0.5000
0.5398
0.5793
0.6179
0.7554
用插值法求
和
的近似值。
实验容:
(1)分别用Lagrange插值法和Newton插值法编程求解;
(2)求出插值多项式系数,对比计算结果。
拉格朗日差值
functionyh=lagrange(x,y,xh)
tic
n=length(x);
m=length(xh);
x=x(:
);
y=y(:
);
xh=xh(:
);
yh=zeros(m,1);
c1=ones(1,n-1);
c2=ones(m,1);
fori=1:
n,
xp=x([1:
i-1i+1:
n]);
yh=yh+y(i)*prod((xh*c1-c2*xp')./(c2*(x(i)*c1-xp')),2);
end
toc
x=[00.10.20.30.4];
y=[0.50000.53980.57930.61790.7554];
xh=[0.130.36];
lagrange(x,y,xh)
时间已过0.026678秒。
ans=
0.5537
0.6780
多项式为
((5*xh-1)*((5*xh)/2-1)*(10*xh-1)*((10*xh)/3-1))/2+(5793*xh*(5*xh-2)*(10*xh-1)*(10*xh-3))/2000-(6179*xh*(5*xh-1/2)*(10*xh-2)*(10*xh-4))/3000+(3777*xh*(5*xh-1)*(10*xh-3)*((10*xh)/3-1/3))/2000-(2699*xh*(5*xh-3/2)*(10*xh-2)*((10*xh)/3-4/3))/500
作图
x0=[00.10.20.30.4];
y0=[0.50000.53980.57930.61790.7554];
xh=[0:
0.01:
0.4];
y=lagrange(x0,y0,xh);
时间已过0.000158秒。
>>plot(xh,y);holdon;
牛顿差值
functionf=Newton(x,y,x0,x1)
tic
symst;
if(length(x)==length(y))
n=length(x);
c(1:
n)=0.0;
else
disp('x和y的维数不相等!
');
return;
end
f=y
(1);
y1=0;
l=1;
for(i=1:
n-1)
for(j=i+1:
n)
y1(j)=(y(j)-y(i))/(x(j)-x(i));
end
c(i)=y1(i+1);
l=l*(t-x(i));
f=f+c(i)*l;
y=y1;
end
f=simplify(f);
g=subs(f,'t',x0)
g1=subs(f,'t',x1)
A=zeros(n,n-1);
A=[y',A];
forj=2:
n
fori=j:
n
A(i,j)=(A(i,j-1)-A(i-1,j-1))/(x(i)-x(i+1-j));
end
end
disp('差商表为');
disp(A);
toc
x0=[00.10.20.30.4];
y0=[0.50000.53980.57930.61790.7554];
Newton(x0,y0,0.13,0.36)
g=
9228481/0000000
g1=
91741/00000
差商表为
1.0e+04*
00000
0.00000.0004000
-0.0000-0.0004-0.004100
-0.0000-0.00010.00160.01900
0.00420.04190.21010.69481.6896
时间已过0.368255秒。
ans=
(251*t^4)/6-(97*t^3)/+(501*t^2)/00+(7*t)/00+1/2
作图
>>t=[0:
0.01:
0.4];
>>f=[];
>>f=(251.*t.^4)./6.-(97.*t.^3)./.+(501.*t.^2)./00.+(7.*t)./00.+1/2;
>>plot(t,f);holdon;
两种方法作图曲线为
拉格朗日插值法时间复杂度为O(n),牛顿插值法时间复杂度为O(n^2),实际上拉格朗日算法计算时间为0.000158,牛顿算法计算时间0.368255。
说明拉格朗日算法时间复杂度低。
4、数值积分
实验目的:
掌握Romberg积分法
实验容:
用Romberg积分法计算下列定积分:
function[R,y]=Romberg(a,b,n)
k=1;
while1
T=zeros(k+1,k+1);
T(1,1)=1/2*(b-a)*(f(a)+f(b));
fori=1:
k
h=(b-a)/2^i;
s=0;
forj=1:
2^(i-1)
s=s+f(a+(2*j-1)*h);
end
T(i+1,1)=T(i,1)/2+h*s;
end
forj=1:
k
c=1/(4^j-1);
form=j:
k
T(m+1,j+1)=T(m+1,j)+c*(T(m+1,j)-T(m,j));
end
end
ifabs(T(k,k)-T(k+1,k+1))break;
end
k=k+1;
end
y=T(k,k);
R=T;
输入积分公式
(1)
functiony=f(x)
y=1/(1+x^4);
[R,y]=Romberg(0,1,10)
R=
0.75000
0.84560.8775
y=
0.7500
(2)
functiony=f(x)
y=cos(sin(x));
[R,y]=Romberg(0,pi/4,10)
R=
0.69120
0.70990.7161
y=
0.6912
(3)
functiony=f(x)
y=exp(x)/x;
[R,y]=Romberg(1,2,10)
R=
3.20640
3.09713.0607
y=
3.2064
(4)
functiony=f(x)
y=x.^x*log(3+x);
[R,y]=Romberg(2,3,10)
R=
27.40760
22.127120.3669
y=
27.4076
龙贝格计算量较小
二、提高技能训练
1、
(2)实验目的与要求
(1)了解Kepler方程;
(2)能够用牛顿迭代方法计算Kepler方程;
(3)通过调节轨道偏心率e查看运算收敛情况。
实验容及数据来源
编写解kepler方程的方法,给定偏心率为0.01、平近点角为32度时计算出偏近点角
。
clear;
x0=32;
symsE;
f=inline('E-0.01*sin(E)-32');
g=inline('1-0.01*cos(E)');
x1=x0-f(x0)/g(x0);
y=abs(x1-x0);
x0=x1;i=i+1;
whiley>10^(-5)
x1=x0-f(x0)/g(x0);
y=abs(x1-x0);
x0=x1;
i=i+1;
end
E=vpa(x1)
结果:
E=32.83
(3)通过调节轨道偏心率e查看运算收敛情况;
e0.00111.82.210
E32.0005532.9999133.5355930.8838131.35097
导0.999161.013181.93911-0.89582-8.97892
可以看出,e在0~1之间,收敛;在1~1.8之间,发散;在1.8~2.2左右,收敛;2.2之后,发散.
2、
(日照时间分布)已知某地区在不同月份的平均日照时间的观测数据如下表所示
月份
1
2
3
4
5
6
7
8
9
10
11
12
日照/(h/月)
80.9
67.2
67.1
50.5
32.0
33.6
35.6
46.8
52.3
62.0
64.1
71.2
试分析日照时间的变化规律。
x=1:
1:
12;
y=[80.967.267.150.532.033.636.646.852.362.064.171.2];
x1=1:
0.01:
12;
y1=interp1(x,y,x1);
y2=interp1(x,y,x1,'spline');
plot(x,y,'go',x1,y1,'r-',x1,y2,'b')
其中红线是分段插值得到的,因此为折线,绿色是三次样条曲线,‘o‘处是插值节点。
从图中可以看出,日照时间从一月开始减少至5月,之后又开始增加至12月。
(水道测量问题)水深数据如下表所示,设船的吃水深度为1.5m,问在
的区域,哪些地方船只避免进入,画船只避免进入的区域。
x
129.0
140.0
108.5
88.0
185.5
195.5
105.5
y
7.5
141.5
28.0
147.0
22.5
137.5
85.5
z
1.22
2.44
1.83
2.44
1.83
2.44
2.44
x
157.5
107.5
77.0
81.0
162.0
162.0
117.5
y
-6.5
-81.5
3.0
56.5
-66.5
84
-38.5
z
2.74
2.74
2.44
2.44
2.74
1.22
2.74
x=[129.0140.0108.588.0185.5195.5105.5157.5107.577.081.0162.0162.0117.5];
y=[7.5141.528.0147.022.5137.585.5-6.5-81.53.056.5-66.584-38.5];
z=-[1.222.441.832.441.832.442.442.742.742.442.442.741.222.74];
x1=linspace(min(x),max(x),40);
y1=linspace(min(y),max(y),40);
[Xi,Yi]=meshgrid(x1,y1);
[Xi,Yi,Zi]=griddata(x,y,z,Xi,Yi,'cubic');
mesh(Xi,Yi,Zi);
xlabel('X');
ylabel('Y');
zlabel('Z');
figure;
[c,h]=contour(Xi,Yi,Zi,[-1.2:
-0.1:
-2.8]);
clabel(c,h);
xlabel('X');
ylabel('Y');
三、本课程设计的心得体会(500字左右)
通过这次计算方法的课程设计,我比过去得到了更多的收获。
不仅仅是知识的学习,更明白了巩固旧知识的重要性,明白了什么是温故而知新。
也懂得了人际关系的重要性,保持好的人际关系,遇到问题,别人才乐于伸出援手。
在试验的进行中,我们使用matlab对问题进行求解。
刚开始的时候,使用matlab解决此类问题有一些难度,包括如何操作,具体算法编排都不是很了解,但是通过互相探讨,查找资料,渐渐的掌握了一些技巧,能够解决一些基础问题。
同时我们也明白了maatlab的强大,它对数学数据和图形的处理能力是c比不上的。
在本次实验中,我的计算方法知识得到了很好的巩固与更新。
在不断的出错,不断的查找,不断的求索与修改之中,过去知识中的漏洞与迷茫渐渐不复存在,也明白了自己掌握的东西还是太少,太浅薄,需要努力寻找新的东西。
而寻找东西的能力恰恰也在本次实验中得到了培养。
由于一次次的求解失误,只能努力的寻找解决方法,书本,网络,同学,这些都在本次试验中运用到。
课程设计已经结束,但是我收获了软件操作知识,数学知识,以及查找资料等能力,这些都是我对专业,对课程有了一个更加深入的了解。