数值分析实验.docx
- 文档编号:29782940
- 上传时间:2023-07-26
- 格式:DOCX
- 页数:24
- 大小:153.80KB
数值分析实验.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
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 分析 实验
![提示](https://static.bdocx.com/images/bang_tan.gif)