数值分析实验.docx

数值分析实验.docx

《数值分析实验.docx》由会员分享，可在线阅读，更多相关《数值分析实验.docx（22页珍藏版）》请在冰豆网上搜索。

数值分析实验

数值分析实验（2014,9,16~10,28）

信计1201班，人数34人

数学系机房

数值分析

计算实习报告册

专业

学号

姓名

2014~2015年第一学期

实验一数值计算的工具Matlab

1.解释下MATLABS序的输出结果

程序：

t=0.1

n=1:

10

e=n/10-n*t

e的结果：

00-5.5511e-01700

-1.1102e-016-1.1102e-016000

2.下面MATLABS序的的功能是什么？

程序：

x=1;while1+x>1,x=x/2,pause（0.02）,end

用迭代法求出x=x/2,的最小值

x=1;whilex+x>x,x=2*x,pause（0.02）,end

用迭代法求出x=2*x,的值，使得2x>X

x=1;whilex+x>x,x=x/2,pause（0.02）,end

用迭代法求出x=x/2,的最小值，使得2x>X

3.考虑下面二次代数方程的求解问题

2

axbxc=0

公式x=电上4ac是熟知的，与之等价地有,对于

2a-b■b-4ac

a=1,b=100000000,c=1，应当如何选择算法。

应该用b~4ac计算，因为b与b2—4ac相近，两个相加减不宜

2au

做分母

357

4.函数sin（x）有幂级数展开sinx=x-x--■■

3!

5!

7!

利用幕级数计算sinx的MATLAB程序为functions=powersin（x）

s=0;

t=x;

n=1;

whiles+t~=s;

s=s+t；

t=-xA2/（（n+1）*（n+2））*t;

n=n+2；

end

t仁cputime;

pause（10）;

t2=cputime;

t0=t2-t1

（a）解释上述程序的终止准则。

当s+t=s，终止循环。

（b）对于x二二/2,1仁/2,21二/2计算的进度是多少？

分别计算多少项？

X=pi/2时，s=1.0000x=11pi/2时，s=-1.0000x=21pi/2时，s=0.9999Cputime分别是0.15630.04690.0156

5.考虑调和级数丄，它是微积分中的发散级数，在计算机上计算该级数的部

n三n

分和，会得到怎么样的结果，为什么？

functions=fun（n）

s=0;

t=1/n;

fori=1:

n

s=s+1/i;

end

当n=100时s=5.1874当n=80时s=4.9655当n=50时，s=4.4992

当n=10时，s=2.9290

23

6.

2!

3!

指数函数的级数展开ex=1x——...，如果对于x:

:

0，用上述的级数近

似计算指数函数的值，这样的算法结果是否会好，为什么?

function

s=powerexp（x）

s=1;n=1;t=1;

while

s+t〜=s;

t=（xAn）/factorial（n）;s=s+t;

n=n+1;

end

当x=-1时，s=0.3679当x=-2时，s=0.1353当x=-3时，s=0.0498

7.考虑数列xi,i=1,2,...n，它的统计平均定义为

11

一、"1n[21n2_2|2

标准差Q=|无（人-x）2数学上等价于坊=（Exi-nx）作为标准差

屮一li吕」n—1i二

的两种算法，你将如何评价他们的得与失。

clc,clear

x=randn（1,10000）;

n=length（x）;

a=sum（x）/n;

y1=sqrt（sum（（x-a）.A2）/（n-1））;

y2=sqrt（（sum（x.A2）-n*aA2）/（n-1））;

y1,y2后面的公式更好

改变m的值求出不同个数x标准差，没有好大差别

实验二插值法计算实习题

1.已知函数在下列个点的值为

Xi

0.2

0.4

0.6

0.8

1.0

f（Xi）

0.98

0.92

0.81

0.64

0.38

试用4次插值多项式p4（x）及三次样条插值S（x）（自然边界条件）对数据进行插

值。

用图给出｛（Xi，y）Xi=0.20.081,^0,1,11,10｝，p4（x）及S（x）。

functionf=lagfun（x）

a=[0.2,0.4,0.6,0.8,1.0];

b=[0.98,0.92,0.81,0.64,0.38];

fori=1:

5

L（i）=1;

forj=1:

5

ifj~=i

L（i）=L（i）*（x-a（j））/（a（i）-a（j））;

end

end

end

f=0;

fori=1:

5

f=f+L（i）*b（i）;

end

执行文件

x0=[0.2,0.4,0.6,0.8,1.0];

y0=[0.98,0.92,0.81,0.64,0.38];

plot（x0,y0,'o'）

holdon

gridon

fplot（'lagfun',[0,1]）;holdon

x=0:

0.1:

1;

plot（x,newton（x0,y0,x）,'r'）;

三次样条插值：

x=[0.2,0.4,0.6,0.8,1.0];

y=[0.98,0.92,0.81,0.64,0.38];

x1=0.2:

0.08:

1;

y仁interp1（x,y,x1,'spline'）

plot（x1,y1）

0.9“

0.8-

0.7f

0.6-

0.5-

2.在区间［-1,1］上分别取n=12,20用两组等距节点对龙格函数f（x）2

1+25x2多项式插值及三次样条插值，对每个n值，分别画出插值函数及f（x）的图形。

拉格朗日插值：

functiony=lagrl（xO,yO,x）

n=length（x0）;m=length（x）;

fori=1:

m

z=x（i）;

s=0.0;

fork=1:

n

p=1.0;

forj=1:

n

ifj〜=k

p=p*（z-x0（j））/（x0（k）-x0（j））;

end

ends=p*y0（k）+s;

end

y（i）=s;

end

再做12和20等距结点插值

forn=12:

8:

20x=-1:

0.01:

1;

y=1./（1+25.*x.A2）;

z=0*x;

x0=-1:

2/（n-1）:

1;y0=1./（1+25.*x0.A2）;

y1=lagrl（x0,y0,x）;plot（x,z,'r',x,y,'k:

',x,y1,'r'）gtext（['Lagr.',num2str（n）]）

holdon

end

title（'Lagrange'）legend（'Lagr插值','f（x）图象'）

图象

Lagrange

x

0

1

4

9

16

25

36

49

64

3.下列数据点的插值

y

0

1

2

3

4

5

6

7

8

可以得到平方根函数的近似值，在区间［0,64］上作图。

（1）用这九个点作8次多项式插值L8（x）.

（2）用三次样条（第一边界条件）程序求

从得到结果看在［0,64］上，哪个插值更精确；在区间［0,10］上，两种插值哪个更精确？

（1）多项式插值

x=［0,1,4,9,16,25,36,49,64］;

y=［0,1,2,3,4,5,6,7,8］;

a=polyfit（x,y,length（x）-1）

poly2sym（a）

答案-0.00000.0000-0.00000.0002-0.00500.0604

-0.38141.3257-0.0000

L（x）=-6345667567529957/19342813113834066795298816*xA8+507195785145098

3/75557863725914323419136*xA7-3204839575550849/590295810358705651712*xA6+8226197088139413/36893488147419103232*xA5-717795609662967/144115188075855872*xA4+2177199843684719/36028797018963968*xA3-3435436202510413/9007199254740992*xA2+5970618836686703/4503599627370496*x-7696702421972085/2475880078570760549798248448

画图x0=[0,1,4,9,16,25,36,49,64];

y0=sqrt（x）;

plot（x0,y0,'r-'）

holdony=[0,1,2,3,4,5,6,7,8];

a=polyfit（x,y,8）;xx=0:

0.1:

64;yy=polyval（a,xx）;plot（xx,yy,'b'）holdon

（3）三次样条

x=[0,1,4,9,16,25,36,49,64];

y=[0,1,2,3,4,5,6,7,8];

x1=0:

0.1:

64;

y1=interp1（x,y,x1,'spline'）

plot（x,sqrt（x）,'r-'）holdon

plot（x1,y1,'b'）holdon

实验三函数逼近与快速傅立叶变换

1

1.对于给函数f（x）2在区间［-1,1］上去x^-l•0.2i（i=1,2,川,10），试

1+25x

求3次曲线拟合，试画出拟合曲线并打印方程，与实验二，第二章的计算实

习题2的结果比较。

首先求出拟合曲线x=-1:

0.2:

1;

y=1./（25*x.A2+1）;

p=polyfit（x,y,3）

Y=-0.5752xA2+0.4814

再输入

x1=-1:

0.01:

1;y仁polyval（p,x1）;

2.由实验给出数据表

x

0

0.1

0.2

0.3

0.5

0.8

1

y

1.0

0.41

0.50

0.61

0.91

2.02

2.46

试求3次，4次多项式的曲线拟合，再根据数据曲线图形，求一个另外函数的拟

合曲线，用图示数据曲线及相应的三种拟合曲线。

先输入代码

x=[0.00.10.20.30.50.81.0];

p3=polyfit（x,y,3）

p4=polyfit（x,y,4）

z=[0.0:

0.1:

1.0];

y3=polyval（p3,z）;

y4=polyval（p4,z）;

plot（x,y,'b',z,y3,'r',z,y4,'g'）

结果-6.6221e+0001.2815e+001-4.6591e+0009.2659e-001

2.8853e+000-1.2335e+0011.6275e+001-5.2987e+0009.4272e-001

实验四数值微分与数值积分

1.用不同数值分析方法计算积分ixlnxdx二…4,

‘09

（1）取不同步长h。

分别用复合梯形及复合辛普森求积计算积分，给出

误差中关于h的函数，并与积分精确值比较两个公式的精度，是否存在一个最小的h，使得精度不能再被改善？

（2）用龙贝格求积计算完成问题

（1）

（3）用自适应辛普森积分，使得精度达到10J。

复合梯形公式

建立m文件clc,clear;

x=0.0000001:

0.1:

1;

y=sqrt（x）.*log（x）;

Fx=trapz（x,y）

error=Fx+4/9

结果Fx=-0.4123

error=0.0321

h=0.01，Fx=-0.4431，error=-0.0014；h=0.001，Fx=-0.4444，

error=-5.4875e-005；h=0.0001，Fx=-0.4444，error=-2.0338e-006；h=0.00001，

Fx=-0.4444，error=-6.4496e-008.

没有一个最小的h使精度不能再改变

复合辛普森求积

Fx=quad（inline（'sqrt（x）*log（x）'）,0.000001,1,0.0001）

error=Fx+4/9

结果Fx=-0.4440

error=4.3273e-004

h=0.00001，Fx=-0.4444error=-6.3537e-005；

h=0.000001，Fx=-0.4444，error=-2.9938e-006；

h=0.0000001，Fx=-0.4444，error=-3.0657e-007；

h=0.0000000000001，Fx=-0.4444error=-9.6548e-009;

h=0.00000000000001，Fx=-0.4444error=-9.6548e-009

存在一个最小的h使精度不能在改善，且h在0.0000000000001与0.00000000000001之间。

龙贝格求积

function[q,n]=shiyan42（f,a,b）

M=1;

abs0=10;k=0;

T=zeros（1,1）;h=b-a;

T（1,1）=（h/2）*（subs（f,a）+subs（f,b））;

whileabs0>0.0001

k=k+1;

h=h/2;p=0;

fori=1:

M

x=a+h*（2*i-1）;

p=p+subs（f,x）;

end

T（k+1,1）=T（k,1）/2+h*p;

M=2*M;

forj=1:

k

T（k+1,j+1）=（（4Aj）*T（k+1,j）-T（k,j））/（4Aj-1）;

end

abs0=abs（T（k+1,j+1）-T（k,j））;

end

q=T（k+1,k+1）;n=k;

在命令窗口输入[Fx,n]=shiyan42（'sqrt（x）*log（x）',10A（-8）,1）自适应辛普森积分

Fx=quad（inline（'sqrt（x）.*log（x）'）,0.000001,1）

结果Fx=-0.4444

实验五数值微分精度及步长的关系

实验目的：

数值计算中误差是不可避免的，希望同学能通过本实验初步认识数值分析中极端重要的两个概念：

截断误差和舍入误差，并认真体会误差对计算结果的影响

问题提出：

设一元函数f：

R>R，贝Uf在xo的导数定义为

实验内容:

计算在x0的导数值可以用如下的差分给出其近似:

这也就给出了计算函数导数的两个简单算法.

请给出几个计算高阶导数的近似算法，并完成如下工作:

1.对同样的h,比较

（1）和

（2）的计算结果.

2.针对计算高阶导数的算法，比较h取不同值时的计算结果

实验要求:

选择有代表性的函数f（x）（最好选择多个），利用MATLAB提供的绘图工具画出该函数在某个区间的导数曲线f⑸（x）.再将数值计算的结果用MATLAB画出来，认真思考实验的结果.同学们还可以围绕这一问题设计其它的实验内容，特别认真的体会导数的结果。

这是n阶导数的程序function[C,L,dyk,k]=ndaolag（X,Y,n）

m=length（X）;n1=m;L=ones（m,m）;

fork=1:

m

v=1;

fori=1:

m

ifk~=i

v=conv（v,poly（X（i）））/（X（k）-X（i））;

end

end

L1（k,:

）=v;l（k,:

）=poly2sym（v）;

end

C=Y*L1;L=Y*I;

symsxdyk

fork=1:

n

k；

dyk=diff（L,x,k）

end在命令窗□输

X=[pi/6,pi/4,pi/3,5*pi/12,pi/2];Y=[0.05,0.7071,0.8660,0.9659,1];

[C,L,dyk,k]=ndaolag（X,Y,4）

Fori=1:

4

symsix,dyi=diff（sin（x）,x,i）

end

实验六数值积分误差传播与算法稳定性

实验目的：

体会算法稳定性在选择算法中的地位.

问题提出：

考虑一个简单的用积分定义的序列

1

En二0xnexdx,n二0,1,2,|1|,

求解En,n=0,1,2,3,HI

实验内容：

由递推关系En=1-nEn4,n=1,2,3川I，可以得到计算巳=[xnexdx,n=0,1,2,111,积分序列洽〉的两种算法.

算法一

初始值为E°=1/e

（1）

递推公式：

En=1-En4,n-1,2,IH

算法

初始值为En=0（3）

递推公式：

EnJ-,n二N,N_1,N_2,||（,3,2,1（4）

n

实验要求：

分别用算法一，算法二，并分别在计算中采用5位、6位和7位有效数字,

请判断哪种算法能给出更精确的结果。

算法一：

a=input（'inputyouxiaoweishua:

'）;

symsnin

in=vpa（（exp（-1））,a）

forn=2:

10

in=vpa（（1-n*in）,a）

end

结果a=5

.36788.26424.20728.17088.14560.12640.11520.7840e-1

.29440-1.9440

算法二：

functionin=noib

b=input（'inputyouxiaoweishua:

'）;

symsnin

in=vpa（0,b）

forn=10:

-1:

2

in=vpa（（（1-in）/n）,b）

end

结果a=5

0..10000.10000.11250.12679.14554.17089.20728

.26424.36788

实验七多项式求根病态问题

实验目的:

希望同学们通过本次实验对问题的病态性有一个初步的体会.

问题提出：

考虑一个高次的代数多项式

20

P（x）=（x-1）（x-2）||]（x-20）~H（x-k）

（1）

kA

显然该多项式的全部根为1,2,,,，20,共计20个，且每个根都是单重的（也

称为简单的），现在考虑该多项式的一个扰动

（2）

P（x）亠7.x19=0

其中；是一个非常小的数，这相当于考虑

（1）中的x19的系数有一个小的扰动,我们比较

（1），

（2）根的情况

实验内容：

为了实现方便，我们先介绍两个MATLAB函数：

“roots”和“poly”.对

u=roots（a）

aF,a?

川am）

则u的各分量是多项式方程

qxna2XnJt川a.xa.1=0

的全部根.而函数

b=poly（v）

若其中v是一个n维向量，该函数的输出b是一个n1维向量，它是以

v的各分量为根的多项式的系数，可见“roots”和“poly”是两个互逆的运算.

ve二zeros（1,21）;

ve

（2）二ess;

roots（poly（1:

20）ve）

上述简单的MATLAB程序便得到

（2）的全部根，程序中“es€即使

（2）中的；

实验要求:

选择充分小的ess，重复上述实验，记录结果的变化，还可以选择

（2）中扰动项改为；x18或者其它形式，观察会有怎样的现象出现。

functiont_charp1_1

%数值实验1.1变态问题

%输入：

［020］之间的扰动项及小的扰动常数%输出：

加扰动都得到的全

部根clc

result=inputdlg（{'请输入扰动项：

在［020］之间的整数：

'},'charpt1-1',1,{'19'}）;

Numb=str2num（char（result））;

if（（Numb>20）|（Numb<0））errordlg（'请输入正确的扰动项：

［020］之间

的整数!

'）;return;end

result=inputdlg（{'请输入（01）间的扰动常数：

'},'charpt1_1',1,{'0.00001'}）;

ess=str2num（char（result））;ve=zeros（1,21）;ve（21-Numb）=ess;root=roots（poly（1:

20）+ve）;

disp（［'对扰动项’,num2str（Numb）,'加扰动’,num2str（ess）,'得到的全

部根为：

’］）；

disp（num2str（root））;

结果分析：

请输入扰动项：

在［020］之间的整数：

19请输入（01）间的

扰动常数：

0.00001