数值分析实验讲解.docx
《数值分析实验讲解.docx》由会员分享，可在线阅读，更多相关《数值分析实验讲解.docx（24页珍藏版）》请在冰点文库上搜索。
数值分析实验讲解
数值分析实验（2014,9,16~10,28）
信计1201班，人数34人
数学系机房
数值分析
计算实习报告册
专业
学号
姓名
2014~2015年第一学期
实验一数值计算的工具Matlab
1.解释下MATLAB程序的输出结果
程序：
t=0.1
n=1:
10
e=n/10-n*t
e的结果：
00-5.5511e-01700-1.1102e-016-1.1102e-016000
2.下面MATLAB程序的的功能是什么？
程序：
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.考虑下面二次代数方程的求解问题
公式
是熟知的，与之等价地有
，对于
，应当如何选择算法。
应该用
计算，因为b与
相近，两个相加减不宜做分母
4.函数
有幂级数展开
利用幂级数计算
的MATLAB程序为
functions=powersin（x）
s=0;
t=x;
n=1;
whiles+t~=s;
s=s+t；
t=-x^2/（（n+1）*（n+2））*t；
n=n+2；
end
t1=cputime;
pause（10）;
t2=cputime;
t0=t2-t1
（a）解释上述程序的终止准则。
当s+t=s，终止循环。
（b）对于
计算的进度是多少？
分别计算多少项？
X=pi/2时，s=1.0000x=11pi/2时，s=-1.0000x=21pi/2时，s=0.9999
Cputime分别是0.15630.04690.0156
5.考虑调和级数
，它是微积分中的发散级数，在计算机上计算该级数的部分和，会得到怎么样的结果，为什么？
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
6.指数函数的级数展开
，如果对于
，用上述的级数近似计算指数函数的值，这样的算法结果是否会好，为什么？
functions=powerexp（x）
s=1;
n=1;
t=1;
whiles+t~=s;
t=（x^n）/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.考虑数列
，它的统计平均定义为
，
标准差
数学上等价于
作为标准差的两种算法，你将如何评价他们的得与失。
clc,clear
x=randn（1,10000）;
n=length（x）;
a=sum（x）/n;
y1=sqrt（sum（（x-a）.^2）/（n-1））;
y2=sqrt（（sum（x.^2）-n*a^2）/（n-1））;
y1,y2后面的公式更好
改变m的值求出不同个数x标准差，没有好大差别
实验二插值法计算实习题
1.已知函数在下列个点的值为
0.2
0.4
0.6
0.8
1.0
0.98
0.92
0.81
0.64
0.38
试用4次插值多项式
及三次样条插值
（自然边界条件）对数据进行插值。
用图给出
，
及
。
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;
y1=interp1（x,y,x1,'spline'）
plot（x1,y1）
2.在区间[-1,1]上分别取
用两组等距节点对龙格函数
作多项式插值及三次样条插值，对每个
值，分别画出插值函数及
的图形。
拉格朗日插值：
functiony=lagrl（x0,y0,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
end
s=p*y0（k）+s;
end
y（i）=s;
end
再做12和20等距结点插值
forn=12:
8:
20
x=-1:
0.01:
1;
y=1./（1+25.*x.^2）;
z=0*x;
x0=-1:
2/（n-1）:
1;
y0=1./（1+25.*x0.^2）;
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）图象'）
图象
3.下列数据点的插值
x
0
1
4
9
16
25
36
49
64
y
0
1
2
3
4
5
6
7
8
可以得到平方根函数的近似值，在区间[0,64]上作图。
（1）用这九个点作8次多项式插值
.
（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*x^8+5071957851450983/75557863725914323419136*x^7-3204839575550849/590295810358705651712*x^6+8226197088139413/36893488147419103232*x^5-717795609662967/144115188075855872*x^4+2177199843684719/36028797018963968*x^3-3435436202510413/9007199254740992*x^2+5970618836686703/4503599627370496*x-7696702421972085/2475880078570760549798248448
画图x0=[0,1,4,9,16,25,36,49,64];
y0=sqrt（x）;
plot（x0,y0,'r-'）
holdon
x=[0,1,4,9,16,25,36,49,64];
y=[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,1]上去
，试求3次曲线拟合，试画出拟合曲线并打印方程，与实验二，第二章的计算实习题2的结果比较。
首先求出拟合曲线x=-1:
0.2:
1;
y=1./（25*x.^2+1）;
p=polyfit（x,y,3）
得出-1.6723e-016-5.7518e-0011.4919e-0174.8412e-001
Y=-0.5752x^2+0.4814
再输入
x1=-1:
0.01:
1;y1=polyval（p,x1）;
plot（x,y,'--*',x1,y1）,grid
2.由实验给出数据表
0
0.1
0.2
0.3
0.5
0.8
1
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];
y=[1.00.410.50.610.912.022.46];
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.用不同数值分析方法计算积分
，
（1）取不同步长h。
分别用复合梯形及复合辛普森求积计算积分，给出误差中关于h的函数，并与积分精确值比较两个公式的精度，是否存在一个最小的h，使得精度不能再被改善？
（2）用龙贝格求积计算完成问题
（1）。
（3）用自适应辛普森积分，使得精度达到
。
复合梯形公式
建立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）=（（4^j）*T（k+1,j）-T（k,j））/（4^j-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）',10^（-8）,1）
自适应辛普森积分
Fx=quad（inline（'sqrt（x）.*log（x）'）,0.000001,1）
结果Fx=-0.4444
实验五数值微分精度及步长的关系
实验目的：
数值计算中误差是不可避免的，希望同学能通过本实验初步认识数值分析中极端重要的两个概念：
截断误差和舍入误差，并认真体会误差对计算结果的影响.
问题提出：
设一元函数
，则
在
的导数定义为
实验内容：
计算在
的导数值可以用如下的差分给出其近似：
（1）
（2）
这也就给出了计算函数导数的两个简单算法.
请给出几个计算高阶导数的近似算法，并完成如下工作：
1.对同样的h，比较
（1）和
（2）的计算结果.
2.针对计算高阶导数的算法，比较h取不同值时的计算结果.
实验要求：
选择有代表性的函数
（最好选择多个），利用MATLAB提供的绘图工具画出该函数在某个区间的导数曲线
.再将数值计算的结果用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*l;
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）
递推公式：
（2）
算法二
初始值为
（3）
递推公式：
（4）
实验要求：
分别用算法一，算法二，并分别在计算中采用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
实验七多项式求根病态问题
实验目的：
希望同学们通过本次实验对问题的病态性有一个初步的体会.
问题提出：
考虑一个高次的代数多项式
（1）
显然该多项式的全部根为1,2，……，20,共计20个，且每个根都是单重的（也称为简单的），现在考虑该多项式的一个扰动
（2）
其中
是一个非常小的数，这相当于考虑
（1）中的
的系数有一个小的扰动，我们比较
（1），
（2）根的情况
实验内容：
为了实现方便，我们先介绍两个MATLAB函数：
“roots”和“poly”.对
u=roots（a）
若其中
为
维的向量，则该函数的输出
为一个
维的向量.设
则
的各分量是多项式方程
的全部根.而函数
若其中
是一个
维向量，该函数的输出
是一个
维向量，它是以
的各分量为根的多项式的系数，可见“roots”和“poly”是两个互逆的运算.
上述简单的MATLAB程序便得到
（2）的全部根，程序中“ess”即使
（2）中的
实验要求：
选择充分小的ess，重复上述实验，记录结果的变化，还可以选择
（2）中扰动项改为
或者其它形式，观察会有怎样的现象出现。
function t_charp1_1
%数值实验1.1变态问题
%输入：
[0 20]之间的扰动项及小的扰动常数 %输出：
加扰动都得到的全部根 clc
result=inputdlg（{'请输入扰动项：
在[0 20]之间的整数：
'},'charpt 1-1',1,{'19'}）;
Numb=str2num（char（result））;
if（（Numb>20）|（Numb<0）） errordlg（'请输入正确的扰动项:
[0 20]之间的整数!
'）;return;end
result=inputdlg（{'请输入（0 1）间的扰动常数：
'},'charpt 1_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））;
结果分析：
请输入扰动项：
在[0 20]之间的整数：
19 请输入（0 1）间的扰动常数：
0.00001
升级会员 