计算方法综合实践讲解.docx

计算方法综合实践讲解.docx

《计算方法综合实践讲解.docx》由会员分享，可在线阅读，更多相关《计算方法综合实践讲解.docx（16页珍藏版）》请在冰豆网上搜索。

计算方法综合实践讲解

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=

249341921721956909228481/450359962737049600000000

g1=

1192694200291039391741/175********41600000000

差商表为

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-（2269814212194733997*t^3）/90071992547409920+（12474970967816329501*t^2）/2702159776422297600+（659777345409770407*t）/4503599627370496000+1/2

作图

>>t=[0:

0.01:

0.4];

>>f=[];

>>f=（251.*t.^4）./6.-（2269814212194733997.*t.^3）./90071992547409920.+（12474970967816329501.*t.^2）./2702159776422297600.+（659777345409770407.*t）./4503599627370496000.+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.005560568887858607922680675983

（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'）;

3、本课程设计的心得体会（500字左右）

通过这次计算方法的课程设计，我比过去得到了更多的收获。

不仅仅是知识的学习，更明白了巩固旧知识的重要性，明白了什么是温故而知新。

也懂得了人际关系的重要性，保持好的人际关系，遇到问题，别人才乐于伸出援手。

在试验的进行中，我们使用matlab对问题进行求解。

刚开始的时候，使用matlab解决此类问题有一些难度，包括如何操作，具体算法编排都不是很了解，但是通过互相探讨，查找资料，渐渐的掌握了一些技巧，能够解决一些基础问题。

同时我们也明白了maatlab的强大，它对数学数据和图形的处理能力是c比不上的。

在本次实验中，我的计算方法知识得到了很好的巩固与更新。

在不断的出错，不断的查找，不断的求索与修改之中，过去知识中的漏洞与迷茫渐渐不复存在，也明白了自己掌握的东西还是太少，太浅薄，需要努力寻找新的东西。

而寻找东西的能力恰恰也在本次实验中得到了培养。

由于一次次的求解失误，只能努力的寻找解决方法，书本，网络，同学，这些都在本次试验中运用到。

课程设计已经结束，但是我收获了软件操作知识，数学知识，以及查找资料等能力，这些都是我对专业，对课程有了一个更加深入的了解。