数值分析实验.docx
《数值分析实验.docx》由会员分享,可在线阅读,更多相关《数值分析实验.docx(37页珍藏版)》请在冰点文库上搜索。
数值分析实验
《数值分析》
计算实习报告册
专业信息与计算科学
学号********
姓名张妃___
2012~2013年第一学期
实验一数值计算的工具Matlab
1.解释下MATLAB程序的输出结果
程序:
t=0.1
n=1:
10
e=n/10-n*t
答案:
t=0.1000小数点后输出四位
n=12345678910n取1到10
e=1.0e-015*
00-0.055500-0.1110-0.1110000
e等于-0.0555×
2.下面MATLAB程序的的功能是什么?
程序:
x=1;while1+x>1,x=x/2,pause(0.02),end
x=1;whilex+x>x,x=2*x,pause(0.02),end
x=1;whilex+x>x,x=x/2,pause(0.02),end
答案:
当x=1时,对x无限二分,到1+x>1结束,且循环之间间隔0.02秒
当x=1时,对
无限循环,到x+x>x结束,且循环之间间隔0.02秒
当x=1时,对x无限二分,到x+x>x结束,且循环之间间隔0.02秒
3.考虑下面二次代数方程的求解问题
公式
是熟知的,与之等价地有
,对于
,应当如何选择算法。
答案:
clc,clear
a=1;
b=100000000;
c=1;
x1=(-b+sqrt(b*b-4*a*c))/2*a
x2=(-b+sqrt(b*b-4*a*c))/2*a
x3=2*c/(-b+sqrt(b*b-4*a*c))
x4=2*c/(-b+sqrt(b*b-4*a*c))
d=x1*x1+100000000*x1+1
e=x2*x2+100000000*x2+1
f=x3*x3+100000000*x3+1
g=x4*x4+100000000*x4+1
x1=-7.4506e-009=-7.4506×
x2=-7.4506e-009
x3=-134217728
x4=-134217728
d=0.25494194030762
e=0.25494194030762
f=4.592625709481985e+015
g=4.592625709481985e+015
因为d,e比f,g更接近0所以公式
更好。
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
修改后程序
functions=powersin(x)
s=0;
x=pi/2;
t=x;
k=0;
n=1;
whiles+t~=s;
s=s+t;
t=-x^2/((n+1)*(n+2))*t;
n=n+2;
k=k+1;
end
K
结果
X=21*pi/2
k=60
ans=0.9999
X=11*pi/2
k=37
ans=-1.0000
X=pi/2
k=11
ans=1.0000
(a)解释上述程序的终止准则。
t=0时,终止
(b)对于
计算的进度是多少?
分别计算多少项?
ans=1.0000ans=-1.0000ans=0.999911,37,60
5.考虑调和级数
,它是微积分中的发散级数,在计算机上计算该级数的部分和,会得到怎么样的结果,为什么?
clc;clear
s=0;
n=20;
fori=n-1:
n+10000000000
t=1/i;
s=s+t;
end
S
结果
s=18.5697
fori=n-1:
n+100000
s=8.5952
fori=n-1:
n+10000000000000
s=18.5697
S会随着i的项数的增大而增大,到一定程度后不变。
因为在计算机中如果一个数很小,计算机会默认为0.
6.指数函数的级数展开
,如果对于
,用上述的级数近似计算指数函数的值,这样的算法结果是否会好,为什么?
1.clc,clear
s=1;
t=1;
x=-1;
fori=1:
100000
t=t*i;
tt=x^i/t;
s=s+tt;
end
s
X=1;s=2.7183x=-1;s=0.3679
2.functions=powersin(x)
s=0;
m=1;
t=1;
n=1;
whiles+t~=s;
s=s+t;
m=m*n;
t=-x^(n)/m;
n=n+2;
end
x=0,f=1;x=1,y=-0.4107;x=-1,f=2,4107
7.考虑数列
,它的统计平均定义为
,
标准差
数学上等价于
作为标准差的两种算法,你将如何评价他们的得与失。
1.clc,clear
a=input('请输入数据')
n=length(a);
sm=sum(a);
ave=sm/n;
b1=0;
fori=1:
n
b1=(a(i)-ave)^2+b1;
end
b1
b2=((1/(n-1))*b1)^(1/2)
结果
b1=133.9368
b2=11.5731
2.clc,clear
a=input('请输入数据')
n=length(a);
sm=sum(a);
ave=sm/n;
b3=0;
fori=1:
n
b3=(a(i))^2+b3;
end
b3
b1=((1/(n-1)*(b3-n*(ave^2))))^(1/2)
结果
b3=189545;
b1=11.8630
实验二插值法计算实习题
1.已知函数在下列个点的值为
0.2
0.4
0.6
0.8
1.0
0.98
0.92
0.81
0.64
0.38
试用4次插值多项式
及三次样条插值
(自然边界条件)对数据进行插值。
用图给出
,
及
。
一.对于多项式插值:
程序1:
(第一种编程)
(1)建立function函数
function[c,d]=poly1(x,y)
n=length(x);d=zeros(n,n);
d(:
1)=y';
forj=2:
n
fork=j:
n
d(k,j)=(d(k,j-1)-d(k-1,j-1))/(x(k)-x(k-j+1));
end
end
c=d(n,n);
fork=(n-1):
-1:
1
c=conv(c,poly(x(k)));
m=length(c);
c(m)=c(m)+d(k,k);
end
(2)建立一个主程序xunan1.m
clc
clear
poly1([0.2,0.4,0.6,0.8,1.0],[0.98,0.92,0.81,0.64,0.38])
结果:
ans=
-0.52080.8333-1.10420.19170.9800
得出所求的牛顿多项式为:
程序2:
(第二种编程)
x=[0.2,0.4,0.6,0.8,1];
y=[0.98,0.92,0.81,0.64,0.38];
a=polyfit(x,y,length(x)-1);%插值
poly2sym(a)%输出插值多项式
结果:
ans=
-1172812402960947/2251799813685248*x^4+1876499844737413/2251799813685248*x^3-4972724588554433/4503599627370496*x^2+6905519428633501/36028797018963968*x+8827055269646203/9007199254740992
画图:
程序1:
x=0.2;
y=-0.5208*x^4+0.8333*x^3-1.1042*x^2+0.1917*x+0.9800
x=0.2+0.08i;
y=-0.5208*x^4+0.8333*x^3-1.1042*x^2+0.1917*x+0.9800
x=0.2+0.08*11i;
y=-0.5208*x^4+0.8333*x^3-1.1042*x^2+0.1917*x+0.9800
x=0.2+0.08*10i;
y=-0.5208*x^4+0.8333*x^3-1.1042*x^2+0.1917*x+0.9800
结果:
y=0.9800
y=0.9847-0.0135i
y=1.2324-0.4306i
y=1.2334-0.3466i
程序2:
x=[0.20.2+0.08i0.2+0.08*11i0.2+0.08*10i];
y=[0.98000.9847-0.0135i1.2324-0.4306i1.2334-0.3466i]
plot(x,y,'-r+')
结果:
二.对于样条插值:
程序1:
x=0.2:
0.2:
1.0;
y=[0.98,0.92,0.81,0.64,0.38];
xx=linspace(0.2,1.0,21)
cs=spline(x,y,xx)
结果:
xx=
Columns1through12
0.20000.24000.28000.32000.36000.40000.44000.48000.52000.56000.60000.6400
Columns13through21
0.68000.72000.76000.80000.84000.88000.92000.96001.0000
cs=
Columns1through12
0.98000.97180.96170.94970.93580.92000.90220.88230.86030.83620.81000.7815
Columns13through21
0.75060.71680.68010.64000.59630.54880.49700.44090.380
clear
clc
x=0.2:
0.2:
1.0;
y=[0.98,0.92,0.81,0.64,0.38];
cs=spline(x,[0y0])
xx=linspace(0.2,1.0,21);
plot(x,y,'o',xx,ppval(cs,xx),'-')
画图:
程序:
clear
clc
x=0.2:
0.2:
1.0;
y=[0.98,0.92,0.81,0.64,0.38];
cs=spline(x,[0y0])
xx=[0.20.2+0.08i0.2+0.08*11i0.2+0.08*10i];
plot(x,y,'o',xx,ppval(cs,xx),'-')
2.在区间[-1,1]上分别取
用两组等距节点对龙格函数
作多项式插值及三次样条插值,对每个
值,分别画出插值函数及
的图形。
一.当
时的程序:
结果:
a=-1;b=1;n=10;
fori=1:
11
x(i)=a+(i-1)*(b-a)/n;
y(i)=1/(1+25*x(i)*x(i));
end
plot(x,y,'r-')
holdon
z=a:
(b-a)/(2*n):
b;
n=length(x);
forj=2:
n
fori=n:
-1:
j
y(i)=(y(i)-y(i-1))/(x(i)-x(i-j+1));
end
end
u=y(n);
m=length(z);
forj=1:
m
fori=n-1:
-1:
1
u=y(i)+u*(z(j)-x(i));
v(j)=u;
end
u=y(n);
end
plot(z,v,'b')
holdoff
二.当
时的程序:
a=-1;b=1;n=20
fori=1:
21
x(i)=a+(i-1)*(b-a)/n;
y(i)=1/(1+25*x(i)*x(i));
end
plot(x,y,'r-')
holdon
z=a:
(b-a)/(2*n):
b;
n=length(x);
forj=2:
n
fori=n:
-1:
j
y(i)=(y(i)-y(i-1))/(x(i)-x(i-j+1));
end
end
u=y(n);
m=length(z);
forj=1:
m
fori=n-1:
-1:
1
u=y(i)+u*(z(j)-x(i));
v(j)=u;
end
u=y(n);
end
plot(z,v,'b')
holdoff
结果:
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]上,两种插值哪个更精确?
一.
多项式插值
作图:
程序:
x=[01491625364964];
y=[0:
1:
8];
p=polyfit(x,y,8)
X=0:
1:
64;
Y=polyval(p,X);
plot(x,y,'r-',X,Y,'y-')
title('8次多项式')
xlabel('x')
ylabel('y')
1..建立function函数
functionf=nan(x,y,x0)
symst;
if(length(x)==length(y))
n=length(x);
else
disp('x!
=y')
return;
end
f=0;
for(i=1:
n)
l=y(i);
for(j=1:
i-1)
l=l*(t-x(j))/(x(i)-x(j));
end
f=f+l;
if(i==n)
if(nargin==3)
f=subs(f,'t',x0);
else
f=collect(f);
f=vpa(f,6);
end
end
end
2.建立一个主程序xunan1.m
clear
clc
x=[01491625364964];
y=[0:
1:
8];
f=nan(x,y)
结果:
f=
.764716e-12*t^8+.535301e-10*t^7+.161447e-7*t^6+.171358e-5*t^5+.136120e-3*t^4+.615224e-2*t^3+.132896*t^2+.860814*t
(2)
三次样条插值:
程序:
%得到三次样条拟合函数
%得到三次样条拟合函数
symsxs
x1=[01491625364964];
y1=[012345678];
x2=[0:
1:
64];
y2=spline(x1,y1,x2);
p=polyfit(x2,y2,3)
S=p
(1)+p
(2)*x+p(3)*x^2+p(4)*x^3
plot(x2,y2)
holdon
plot(x1,y1,'r-')
结果:
p=
0.0000-0.00420.25270.8876
S=
2288075067923491/73786976294838206464-2399112304472833/576460752303423488*x+4552380473376713/18014398509481984*x^2+999337332656867/1125899906842624*x^3
所以在[0,64]区间上样条插值更精确。
二.当[0,10]区间上
x=[01491625364964];
y=[0:
1:
8];
p=polyfit(x,y,8)
X=0:
1:
10;
Y=polyval(p,X);
plot(x,y,'r-',X,Y,'y-')
title('8次多项式')
xlabel('x')
ylabel('y')
结果:
程序:
symsxs
x1=[01491625364964];
y1=[012345678];
x2=[0:
1:
10];
y2=spline(x1,y1,x2);
p=polyfit(x2,y2,3)
S=p
(1)+p
(2)*x+p(3)*x^2+p(4)*x^3
plot(x2,y2)
holdon
plot(x1,y1,'r-')
title('三次样条')
结果:
由图可知在[0,10]区间上多项式插值更精确。
实验三函数逼近与快速傅立叶变换
1.对于给函数
在区间[-1,1]上去
,试求3次曲线拟合,试画出拟合曲线并打印方程,与实验二,第二章的计算实习题2的结果比较。
clc,clear
n=11;
fori=1:
n
x(i)=-1+0.2*(i-1);
y(i)=1/(1+25*x(i)*x(i));
end
A=polyfit(x,y,3)
z=polyval(A,x)
plot(x,y,'k+',x,z,'r')
结果:
A=
0.0000-0.5752-0.00000.4841
f(x)=-0.5752*X^2+0.4841
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次多项式的曲线拟合,再根据数据曲线图形,求一个另外函数的拟合曲线,用图示数据曲线及相应的三种拟合曲线。
clc,clear
x=[0.00.10.20.30.50.81.0]
y=[1.00.410.500.610.912.022.46]
A=polyfit(x,y,3)
z=polyval(A,x)
plot(x,y,'k+',x,z,'r')
结果:
A=
-6.622112.8147-4.65910.9266
f(x)=-6.6221*x^3+12.8147*X^2+-4.6591*x+0.9266
4次:
clc,clear
x=[0.00.10.20.30.50.81.0]
y=[1.00.410.500.610.912.022.46]
A=polyfit(x,y,4)
z=polyval(A,x)
plot(x,y,'k+',x,z,'r')
结果:
A=
2.8853-12.334816.2747-5.29870.9427
f(x)=2.8853*x^4-12.3348*x^3+16.2747*X^2+-5.2987*x+0.9427
根据数据曲线形状,我觉得它像是一个开口向上的二次函数的图像.
clc,clear
x=[0.00.10.20.30.50.81.0]
y=[1.00.410.500.610.912.022.46]
A=polyfit(x,y,2)
z=polyval(A,x)
plot(x,y,'k+',x,z,'r')
结果:
A=
3.1316-1.24000.7356
f(x)=3.1316*X^2+-1.2400*x+0.7356
实验五数值微分精度及步长的关系
实验目的:
数值计算中误差是不可避免的,希望同学能通过本实验初步认识数值分析中极端重要的两个概念:
截断误差和舍入误差,并认真体会误差对计算结果的影响.
问题提出:
设一元函数
,则
在
的导数定义为
实验内容:
计算在
的导数值可以用如下的差分给出其近似:
(1)
(2)
这也就给出了计算函数导数的两个简单算法.
请给出几个计算高阶导数的近似算法,并完成如下工作:
1.对同样的h,比较
(1)和
(2)的计算结果.
2.针对计算高阶导数的算法,比较h取不同值时的计算结果.
实验要求:
选择有代表性的函数
(最好选择多个),利用MATLAB提供的绘图工具画出该函数在某个区间的导数曲线
.再将数值计算的结果用MATLAB画出来,认真思考实验的结果.同学们还可以围绕这一问题设计其它的实验内容,特别认真的体会导数的结果。
对于问题1:
我们分别用算法
(1)、
(2)计算函数
在
上导数,并与
在
的值进行对比,画出图像(如下图1所示):
由图像可知:
其差距并不明显.于是我们便计算出分别用算法
(1)、
(2)计算函数
在
上导数的误差,并画出其图像,(如下图2所示,其中蓝色表示用算法
(1)计算得到的误差,红色表示用算法
(2)计算得到的误差):
图1
图2
由图像可知:
对于计算函数导数,
(2)是比
(1)更好的算法.
对于问题1实验结论的分析:
不论采用怎样的算法,计算结果通常都会有误差。
比如算法
(1),由Taylor公式,知:
所以有
利用上式来计算
,其截断误差为:
所以误差是存在的,并且当步长h越来越小时,
(1)的近似程度也越来越好。
类似地可以分析
(2)的截断误差为:
上述截断误差的分析表明
(2)是比
(1)更好的算法,因为对步长h(<<1),
(2)比
(1)更接近于
.
关于问题1的程序如下:
functiony=dsh1(fu,x,h)
y=(feval(fu,x+h)-feval(fu,x))/h;
y;
functiony=dsh2(fu,x,h)
y=(feval(fu,x+h)-feval(fu,x-h))/(2*h);
y;
clear,clc;
x=-pi:
0.1:
pi;
y=cos(x);
y1=dsh1('sin',x,0.01);
y2=dsh2('sin',x,0.01);
plot(x,y,'b*',x,y1,'ro',x,y2,'go')
pause
(1);
wucha1=abs(y1-y);
wucha2=abs(y2-y);
plot(x,wucha1,'b*',x,wucha2,'ro')
对于问题2:
我们用如下算法来计算二阶导数:
该算法程序如下:
functiony=dsh3(fu,x,h)
y=(feval(fu,x+2*h)-feval(fu,x)-feval(fu,x+h)+feval(fu,x-h))/(2*h*h);
y;
当h分别为0.