数值分析实验Word文档下载推荐.docx
《数值分析实验Word文档下载推荐.docx》由会员分享，可在线阅读，更多相关《数值分析实验Word文档下载推荐.docx（37页珍藏版）》请在冰点文库上搜索。
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
修改后程序
x=pi/2;
k=0;
s=s+t;
t=-x^2/（（n+1）*（n+2））*t;
n=n+2;
k=k+1;
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
n=20;
fori=n-1:
n+10000000000
t=1/i;
S
s=18.5697
n+100000
s=8.5952
n+10000000000000
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;
s
X=1;
s=2.7183x=-1;
s=0.3679
2.functions=powersin（x）
m=1;
m=m*n;
t=-x^（n）/m;
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;
n
b1=（a（i）-ave）^2+b1;
b1
b2=（（1/（n-1））*b1）^（1/2）
b1=133.9368
b2=11.5731
2.clc,clear
b3=0;
b3=（a（i））^2+b3;
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:
d（k,j）=（d（k,j-1）-d（k-1,j-1））/（x（k）-x（k-j+1））;
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）;
（2）建立一个主程序xunan1.m
clc
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）%输出插值多项式
-1172812402960947/2251799813685248*x^4+1876499844737413/2251799813685248*x^3-4972724588554433/4503599627370496*x^2+6905519428633501/36028797018963968*x+8827055269646203/9007199254740992
画图：
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;
x=0.2+0.08*11i;
x=0.2+0.08*10i;
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+'
二．对于样条插值：
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=
0.98000.97180.96170.94970.93580.92000.90220.88230.86030.83620.81000.7815
Columns13through21
0.75060.71680.68010.64000.59630.54880.49700.44090.380
cs=spline（x,[0y0]）
xx=linspace（0.2,1.0,21）;
o'
xx,ppval（cs,xx）,'
-'
xx=[0.20.2+0.08i0.2+0.08*11i0.2+0.08*10i];
2.在区间[-1,1]上分别取
用两组等距节点对龙格函数
作多项式插值及三次样条插值，对每个
值，分别画出插值函数及
的图形。
一．当
时的程序：
a=-1;
b=1;
n=10;
11
x（i）=a+（i-1）*（b-a）/n;
y（i）=1/（1+25*x（i）*x（i））;
r-'
holdon
z=a:
（b-a）/（2*n）:
b;
n=length（x）;
forj=2:
fori=n:
-1:
j
y（i）=（y（i）-y（i-1））/（x（i）-x（i-j+1））;
u=y（n）;
m=length（z）;
forj=1:
m
fori=n-1:
1
u=y（i）+u*（z（j）-x（i））;
v（j）=u;
u=y（n）;
plot（z,v,'
b'
holdoff
二．当
n=20
21
3.下列数据点的插值
x
4
9
16
25
36
49
64
y
2
3
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:
64;
Y=polyval（p,X）;
X,Y,'
y-'
title（'
8次多项式'
xlabel（'
x'
ylabel（'
y'
1..建立function函数
functionf=nan（x,y,x0）
symst;
if（length（x）==length（y））
else
disp（'
x!
=y'
return;
f=0;
for（i=1:
n）
l=y（i）;
for（j=1:
i-1）
l=l*（t-x（j））/（x（i）-x（j））;
f=f+l;
if（i==n）
if（nargin==3）
f=subs（f,'
t'
x0）;
else
f=collect（f）;
f=vpa（f,6）;
2.建立一个主程序xunan1.m
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:
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）
plot（x1,y1,'
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]区间上
10;
10];
三次样条'
由图可知在[0,10]区间上多项式插值更精确。
实验三函数逼近与快速傅立叶变换
1.对于给函数
在区间[-1,1]上去
，试求3次曲线拟合，试画出拟合曲线并打印方程，与实验二，第二章的计算实习题2的结果比较。
n=11;
x（i）=-1+0.2*（i-1）;
y（i）=1/（1+25*x（i）*x（i））;
A=polyfit（x,y,3）
z=polyval（A,x）
k+'
x,z,'
r'
A=
0.0000-0.5752-0.00000.4841
f（x）=-0.5752*X^2+0.4841
2.由实验给出数据表
0.1
0.3
0.5
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.500.610.912.022.46]
-6.622112.8147-4.65910.9266
f（x）=-6.6221*x^3+12.8147*X^2+-4.6591*x+0.9266
4次:
A=polyfit（x,y,4）
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
根据数据曲线形状，我觉得它像是一个开口向上的二次函数的图像.
A=polyfit（x,y,2）
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）;
clear,clc;
x=-pi:
0.1:
pi;
y=cos（x）;
y1=dsh1（'
sin'
x,0.01）;
y2=dsh2（'
b*'
x,y1,'
ro'
x,y2,'
go'
pause
（1）;
wucha1=abs（y1-y）;
wucha2=abs（y2-y）;
plot（x,wucha1,'
x,wucha2,'
对于问题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）;
当h分别为0.
升级会员 