数值分析实验报告1.docx
- 文档编号:24646401
- 上传时间:2023-05-29
- 格式:DOCX
- 页数:21
- 大小:269.65KB
数值分析实验报告1.docx
《数值分析实验报告1.docx》由会员分享,可在线阅读,更多相关《数值分析实验报告1.docx(21页珍藏版)》请在冰豆网上搜索。
数值分析实验报告1
实验报告
实验项目名称插值法
实 验 室 数学实验室
所属课程名称数值逼近
实验类型算法设计
实验日期
班级
学号
姓名
成绩
实验概述:
【实验目的及要求】
本次实验的目的是熟练《数值分析》第二章“插值法”的相关内容,掌握三种插值方法:
牛顿多项式插值,三次样条插值,拉格朗日插值,并比较三种插值方法的优劣。
本次试验要求编写牛顿多项式插值,三次样条插值,拉格朗日插值的程序编码,并在MATLAB软件中去实现.
【实验原理】
《数值分析》第二章“插值法”的相关内容,包括:
牛顿多项式插值,三次样条插值,拉格朗日插值的相应算法和相关性质.
【实验环境】(使用的软硬件)
软件:
MATLAB2012a
硬件:
电脑型号:
联想Lenovo昭阳E46A笔记本电脑
操作系统:
Windows8专业版
处理器:
Intel(R)Core(TM)i3CPUM350@2.27GHz2。
27GHz
实验内容:
【实验方案设计】
第一步,将书上关于三种插值方法的内容转化成程序语言,用MATLAB实现;第二步,分别用牛顿多项式插值,三次样条插值,拉格朗日插值求解不同的问题.
【实验过程】(实验步骤、记录、数据、分析)
实验的主要步骤是:
首先分析问题,根据分析设计MATLAB程序,利用程序算出问题答案,分析所得答案结果,再得出最后结论。
实验一:
已知函数在下列各点的值为
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,yi),xi=0.2+0.08i,i=0,1,11,10},P4(x)及S(x)。
(1)首先我们先求牛顿插值多项式,此处要用4次牛顿插值多项式处理数据。
已知n次牛顿插值多项式如下:
Pn=f(x0)+f[x0,x1](x—x0)+f[x0,x1,x2](x-x0)(x—x1)+···+f[x0,x1,···xn](x—x0)···(x—xn—1)
我们要知道牛顿插值多项式的系数,即均差表中得部分均差.
在MATLAB的Editor中输入程序代码,计算牛顿插值中多项式系数的程序如下:
functionvarargout=newtonliu(varargin)
clear,clc
x=[0。
20.40。
60。
81。
0];
fx=[0.980.920。
810.640.38];
newtonchzh(x,fx);
functionnewtonchzh(x,fx)
%由此函数可得差分表
n=length(x);
fprintf(’*****************差分表*****************************\n');
FF=ones(n,n);
FF(:
1)=fx';
fori=2:
n
forj=i:
n
FF(j,i)=(FF(j,i-1)—FF(j-1,i-1))/(x(j)—x(j—i+1));
end
end
fori=1:
n
fprintf('%4。
2f’,x(i));
forj=1:
i
fprintf('%10.5f’,FF(i,j));
end
fprintf(’\n');
end
由MATLAB计算得:
xi
f(xi)
一阶差商
二阶差商
三阶差商
四阶差商
0。
20
0.980
0。
40
0。
920
-0.30000
0。
60
0。
810
-0。
55000
—0.62500
0。
80
0。
640
-0.85000
—0。
75000
—0.20833
1.00
0。
380
—1。
30000
—1。
12500
-0.62500
—0.52083
所以有四次插值牛顿多项式为:
P4(x)=0.98-0.3(x—0.2)—0。
62500(x-0.2)(x-0。
4)-0。
20833(x-0.2)(x—0.4)(x—0.6)—0。
52083(x—0.2)(x-0。
4)(x-0.6)(x—0.8)
(2)接下来我们求三次样条插值函数.
用三次样条插值函数由上题分析知,要求各点的M值:
三次样条插值函数计算的程序如下:
functiontgsanci(n,s,t)%n代表元素数,s,t代表端点的一阶导。
x=[0。
20.40.60.81.0];
y=[0.980.920.810.640.38];
n=5
forj=1:
1:
n-1
h(j)=x(j+1)—x(j);
end
forj=2:
1:
n-1
r(j)=h(j)/(h(j)+h(j-1));
end
forj=1:
1:
n-1
u(j)=1-r(j);
end
forj=1:
1:
n—1
f(j)=(y(j+1)—y(j))/h(j);
end
forj=2:
1:
n-1
d(j)=6*(f(j)—f(j-1))/(h(j-1)+h(j));
end
d
(1)=0
d(n)=0
a=zeros(n,n);
forj=1:
1:
n
a(j,j)=2;
end
r
(1)=0;
u(n)=0;
forj=1:
1:
n-1
a(j+1,j)=u(j+1);
a(j,j+1)=r(j);
end
b=inv(a)
m=b*d'
p=zeros(n—1,4);%p矩阵为S(x)函数的系数矩阵
forj=1:
1:
n-1
p(j,1)=m(j)/(6*h(j));
p(j,2)=m(j+1)/(6*h(j));
p(j,3)=(y(j)—m(j)*(h(j)^2/6))/h(j);
p(j,4)=(y(j+1)—m(j+1)*(h(j)^2/6))/h(j);
end
p
得到m=(0-1.6071—1。
0714—3.10710)T
即M0=0;M1=-1。
6071;M2=-1。
0714;M3=—3.1071;M4=0
则根据三次样条函数定义,可得:
S(x)=
接着,在CommandWindow里输入画图的程序代码,
下面是画牛顿插值以及三次样条插值图形的程序:
x=[0.20.40。
60。
81.0];
y=[0.980.920.810。
640。
38];
plot(x,y)
holdon
fori=1:
1:
5
y(i)=0.98-0。
3*(x(i)—0。
2)—0.62500*(x(i)—0.2)*(x(i)—0.4)—0.20833*(x(i)—0。
2)*(x(i)—0.4)*(x(i)-0.6)-0。
52083*(x(i)—0.2)*(x(i)-0.4)*(x(i)-0。
6)*(x(i)—0.8)
end
k=[011011]
x0=0。
2+0.08*k
fori=1:
1:
4
y0(i)=0.98—0。
3*(x(i)—0。
2)-0。
62500*(x(i)-0。
2)*(x(i)-0。
4)-0。
20833*(x(i)—0.2)*(x(i)—0.4)*(x(i)-0.6)-0。
52083*(x(i)—0。
2)*(x(i)—0.4)*(x(i)—0.6)*(x(i)-0。
8)
end
plot(x0,y0,'o',x0,y0)
holdon
y1=spline(x,y,x0)
plot(x0,y1,’o')
holdon
s=csape(x,y,’variational’)
fnplt(s,'r’)
holdon
gtext('三次样条自然边界’)
gtext('原图像')
gtext('4次牛顿插值')
运行上述程序可知:
给出的{(xi,yi),xi=0。
2+0.08i,i=0,1,11,10}点,S(x)及P4(x)图形如下所示:
实验二:
在区间[-1,1]上分别取
用两组等距节点对龙格函数
作多项式插值及三次样条插值,对每个
值,分别画出插值函数即
的图形。
我们先求多项式插值:
在MATLAB的Editor中建立一个多项式的M-file,输入如下的命令(如牛顿插值公式):
function[C,D]=newpoly(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
当n=10时,我们在CommandWindow中输入以下的命令:
clear,clf,holdon;
X=-1:
0。
2:
1;
Y=1./(1+25*X。
^2);
[C,D]=newpoly(X,Y);
x=—1:
0。
01:
1;
y=polyval(C,x);
plot(x,y,X,Y,’.’);
gridon;
xp=-1:
0。
2:
1;
z=1。
/(1+25*xp。
^2);
plot(xp,z,’r')
得到插值函数和f(x)图形:
当n=20时,我们在CommandWindow中输入以下的命令:
clear,clf,holdon;
X=-1:
0。
1:
1;
Y=1。
/(1+25*X.^2);
[C,D]=newpoly(X,Y);
x=—1:
0.01:
1;
y=polyval(C,x);
plot(x,y,X,Y,’。
’);
gridon;
xp=-1:
0。
1:
1;
z=1./(1+25*xp.^2);
plot(xp,z,’r')
得到插值函数和f(x)图形:
下面再求三次样条插值函数,在MATLAB的Editor中建立一个多项式的M—file,
输入下列程序代码:
functionS=csfit(X,Y,dx0,dxn)
N=length(X)-1;
H=diff(X);
D=diff(Y)./H;
A=H(2:
N-1);
B=2*(H(1:
N-1)+H(2:
N));
C=H(2:
N);
U=6*diff(D);
B
(1)=B
(1)-H
(1)/2;
U
(1)=U
(1)—3*(D
(1));
B(N-1)=B(N—1)-H(N)/2;
U(N—1)=U(N-1)—3*(-D(N));
fork=2:
N-1
temp=A(k—1)/B(k—1);
B(k)=B(k)—temp*C(k—1);
U(k)=U(k)—temp*U(k-1);
end
M(N)=U(N—1)/B(N—1);
fork=N-2:
—1:
1
M(k+1)=(U(k)-C(k)*M(k+2))/B(k);
end
M
(1)=3*(D
(1)-dx0)/H
(1)-M
(2)/2;
M(N+1)=3*(dxn—D(N))/H(N)—M(N)/2;
fork=0:
N-1
S(k+1,1)=(M(k+2)—M(k+1))/(6*H(k+1));
S(k+1,2)=M(k+1)/2;
S(k+1,3)=D(k+1)-H(k+1)*(2*M(k+1)+M(k+2))/6;
S(k+1,4)=Y(k+1);
end
当n=10时,我们在CommandWindow中输入以下的命令:
clear,clc
X=—1:
0.2:
1;
Y=1./(25*X。
^2+1);
dx0=0.0739644970414201;dxn=—0.0739644970414201;
S=csfit(X,Y,dx0,dxn)
x1=—1:
0。
01:
-0。
5;y1=polyval(S(1,:
),x1-X
(1));
x2=-0。
5:
0.01:
0;y2=polyval(S(2,:
),x2-X
(2));
x3=0:
0。
01:
0.5;y3=polyval(S(3,:
),x3-X(3));
x4=0.5:
0.01:
1;y4=polyval(S(4,:
),x4-X(4));
plot(x1,y1,x2,y2,x3,y3,x4,y4,X,Y,'.')
结果如图:
当n=20时,我们在CommandWindow中输入以下的命令:
clear,clc
X=—1:
0.1:
1;
Y=1./(25*X.^2+1);
dx0=0.0739644970414201;dxn=-0。
0739644970414201;
S=csfit(X,Y,dx0,dxn)
x1=-1:
0。
01:
—0。
5;y1=polyval(S(1,:
),x1-X
(1));
x2=-0.5:
0.01:
0;y2=polyval(S(2,:
),x2—X
(2));
x3=0:
0。
01:
0.5;y3=polyval(S(3,:
),x3-X(3));
x4=0。
5:
0.01:
1;y4=polyval(S(4,:
),x4-X(4));
plot(x1,y1,x2,y2,x3,y3,x4,y4,X,Y,’。
')
结果如图:
实验三:
下列数据点的插值
x
0
1
4
9
16
25
36
49
64
y
0
1
2
3
4
5
6
7
8
可以得到平方根函数的近似,在区间[0,64]上作图。
(1)用这9各点作8次多项式插值L8(x).
(2)用三次样条(自然边界条件)程序求S(x).从结果看在[0,64]上,那个插值更精确;在区间[0,1]上,两种哪个更精确?
L8(x)可由公式Ln(x)=
得出。
三次样条可以利用自然边界条件.写成矩阵:
其中
j=
i=
,dj=6f[xj-1,xj,xj+1],
n=
0=0d0=dn=0
l0(x)=
l1(x)=
l2(x)=
l3(x)=
l4(x)=
l5(x)=
l6(x)=
l7(x)=
l8(x)=
L8(x)=l1(x)+2l2(x)+3l3(x)+4l4(x)+5l5(x)+6l6(x)+7l7(x)+8l8(x)
求三次样条插值函数由MATLAB计算:
可得矩阵形式的线性方程组为:
在MATLAB中的Editor中输入程序代码,
以下是三次样条函数的程序代码:
functiontgsanci(n,s,t)%n代表元素数,s,t代表端点的一阶导。
y=[012345678];
x=[01491625364964];
n=9
forj=1:
1:
n—1
h(j)=x(j+1)-x(j);
end
forj=2:
1:
n-1
r(j)=h(j)/(h(j)+h(j—1));
end
forj=1:
1:
n-1
u(j)=1—r(j);
end
forj=1:
1:
n-1
f(j)=(y(j+1)—y(j))/h(j);
end
forj=2:
1:
n—1
d(j)=6*(f(j)—f(j—1))/(h(j—1)+h(j));
end
d
(1)=0
d(n)=0
a=zeros(n,n);
forj=1:
1:
n
a(j,j)=2;
end
r
(1)=0;
u(n)=0;
forj=1:
1:
n-1
a(j+1,j)=u(j+1);
a(j,j+1)=r(j);
end
b=inv(a)
m=b*d'
t=a
p=zeros(n—1,4);%p矩阵为S(x)函数的系数矩阵
forj=1:
1:
n—1
p(j,1)=m(j)/(6*h(j));
p(j,2)=m(j+1)/(6*h(j));
p(j,3)=(y(j)-m(j)*(h(j)^2/6))/h(j);
p(j,4)=(y(j+1)-m(j+1)*(h(j)^2/6))/h(j);
end
p
解得:
M0=0;M1=-0.5209;M2=0.0558;M3=—0.0261;M4=0。
0006;M5=-0.0029;M6=—0。
0008;M7=—-0。
0009;
M8=0,则三次样条函数:
S(x)=
下面进行画图,在CommandWindow中输入画图的程序代码:
%画图形比较那个插值更精确的函数:
x0=[01491625364964];
y0=[012345678];
x=0:
64;y=lagr1(x0,y0,x);
plot(x0,y0,’o’)
holdon
plot(x,y,’r’);
holdon;
pp=csape(x0,y0,'variational’)
fnplt(pp,’g');
holdon;
plot(x0,y0,':
b');holdon
%axis([0201]);%看[01]区间的图形时加上这条指令
gtext(’三次样条插值')
gtext(’原图像')
gtext('拉格朗日插值’)
%下面是求拉格朗日插值的函数
functiony=lagr1(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
拉格朗日插值函数与三次样条插值函数如图中所示,绿色实线条为三次样条插值曲线,蓝色虚线条为x=y2的曲线,另外一条红色线条为拉格朗日插值曲线。
图3-1为[01]的曲线,图3-2为[064]区间上的曲线。
图3-1
图3—2
由图3—1可以看出,红色的线条与蓝色虚线条几乎重合,所以可知拉格朗日插值函数的曲线更接近开平方根的函数曲线,在[0,1]朗格朗日插值更精确。
而在区间[0,64]上从图3—2中可以看出蓝色虚线条和绿色线条是几乎重合的,而红色线条在[30,70]之间有很大的振荡,所在在区间[0,64]三次样条插值更精确写.
【结论】(结果)
单个多项式高次插值效果并不理想,有龙格现象,偏差大,没有使用价值.而分段低次插值则精确度较高,拟合效果较好,而三次样条插值具有良好的收敛性与稳定性,与分段低次插值相比较光滑度更高,而且提供的信息也相对少一些.我们可以看到,在以上的三道实验题里,我们可以从图形中看出,三次样条的拟合程度是三种插值函数里最好的。
【小结】
通过此次实验,我对牛顿多项式插值,三次样条插值,拉格朗日插值有了更进一步的了解,知道了三次样条的拟合程度在高次的情况下更高,在理论上和应用上都有重要意义,在利用计算机编程软件进行高次插值的时候,我们可以多考虑利用三次样条进行插值.
指导教师评语及成绩:
成绩:
指导教师签名:
批阅日期:
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 分析 实验 报告