实验报告
实验项目名称实验室
实验项目名称
实验室
所属课程名称
实验类型
实验日期
数学实验室
数值逼近
算法设计
级
号
名
实验概述：
【实验目的及要求】
值方法：牛顿在
值方法：牛顿
在MATLAB件中
本次试验要求编写牛顿多项式插值，三次样条插值，拉格朗日插值的程序编码，并 去实现。
【实验原理】
《数值分析》第二章“插值法”的相关内容，包括：牛顿多项式插值，三次样条插值， 拉格朗日
插值的相应算法和相关性质。
【实验环境】（使用的软硬件） 软件：
MATLAB 2012a
硬件： 电脑型号：联想 Lenovo昭阳E46A笔记本电脑
操作系统： Win dows 8专业版
处理器：In tel ( R Core( TM i3 CPU M 350 @2.27GHz 2.27GHz
实验内容：
【实验方案设计】
第一步，将书上关于三种插值方法的内容转化成程序语言，用 MATLAB现;
第二步，分别用牛顿多项式插值，三次样条插值，拉格朗日插值求解不同的问题。
【实验过程】（实验步骤、记录、数据、分析）
实验的主要步骤是：首先分析问题，根据分析设计 MATLA程序，利用程序算出
问题答案，分析所得答案结果，再得出最后结论。
实验一：
已知函数在下列各点的值为
X
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）。
4次牛顿插值多项式处理数据。（1）
4次牛顿插值多项式处理数据。
已知n次牛顿插值多项式如下：
f=f (X 0) +f[X 0,X1] (X-X 0) + f[X O,X1,X2] (X-X 0)( X-X 1) + + f[x O,X1 ,为](X-X
0)我们要知道牛顿插值多项式的系数，即均差表中得部分均差。
在MATLAB勺Editor中输入程序代码，计算牛顿插值中多项式系数的程序如下:
fun cti on varargout=n ewt on liu(vararg in) clear,clc x=[0.2 0.4 0.6 0.8 1.0];
fx=[0.98 0.92 0.81 0.64 0.38];
n ewt on chzh(x,fx);
fun cti on n ewt on chzh(x,fx)
%由此函数可得差分表
n=len gth(x);************fprintfFF=on es( n,n);
n=len gth(x);
************
fprintf
FF=on es( n,n);
FF(:,1)=fx : for i=2: n
差分表
********************
*******
**\n'
);
for j=i:n
FF(j,i)=(FF(j,i-1)-FF(j-1,i-1))/(x(j)-x(j-i+1));
end
end for i=1: n
fprin tf('%4.2f,x(i));
for j=1:i
fprin tf('%10.5f,FF(i,j));
end fpri n tf('\n');
end
由MATLAB+算得:
X
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
所以有四次插值牛顿多项式为：
(x-0.2)(x-0.4)(x-0.6)(x-0.8 )
(2)接下来我们求三次样条插值函数。
用三次样条插值函数由上题分析知 ,要求各点的M值:
2
0
0
0
0
M
0
0.500
2
0.500
0
0
M
-3.7500
0
0.500
2
0.50
0
M
-4.5000
0
0
0.500
0 2
0.500
M
-6.7500
0
0
0
0
2
M
0
三次样条插值函数计算的程序如下
:
function tgsanci(n,s,t) %n 代表元素数，s,t 代表端点的一阶导。
x=[0.2 0.4 0.6 0.8 1.0];
y=[0.98 0.92 0.81 0.64 0.38];
n=5
for j=1:1: n-1
h(j)=x(j+1)-x(j);
end
for j=2:1: n-1
丽h(j)/(h(j)+ha-1));
end
for j=1:1: n-1
u(j)=1-r(j);
end
for j=1:1: n-1
f(j)=(y(j+1)-y(j))/h(j);
d(j)=6*(f(j)-f(j-1))/(h(j-1)+h(j));
end
d(1)=0 d(n )=0 a=zeros( n,n);
for j=1:1: n
a(i,i)=2;
end
TOC \o "1-5" \h \z 3 3
-1.33930.6 x )3 0.8929(x-0.4)—3 4.65360.6 x—)4.0857(x 0.4—) , x S(x)= [0406
3 3
-0.89290.8 X) -2.5893 X-0.6) 4.08570.8 x ) 3.3036x 0.6 ) , x ]
-2.58931.0 x ) 3-0(x 0.8 )33036(1.0 x ) 1.9(x 0.8) , x [0.8,1.0] [0.6,0.81
接着，在Comma nd Window里输入画图的程序代码，
下面是画牛顿插值以及三次样条插值图形的程序：
x=[0.2 0.4 0.6 0.8 1.0];
y=[0.98 0.92 0.81 0.64 0.38];
Plot(x,y)
hold on
for i=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=[0 1 10 11]
x0=0.2+0.08*k
for i=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( xO,yO,'o',xO,yO )
hold on
y1=spli n e(x,y,xO)
plot(x0,y1,'o')
hold on
s=csape(x,y,'variatio nal')
fnplt(sy)
hold on
gtext('三次样条自然边界')
gtext('原图像')
gtext('4 次牛顿插值')
运行上述程序可知：给出的｛ (X, yi), xi=0.2+0.08i , i=0 , 1,11,10｝点，S (x)及 P4 (x)图
形如下所示：
实验
在区间［-1,1 :上分别取n 10,20用两组等距节点对龙格函数
f(x)
1
1 25x2
作多项式插值及
得到插值函数和f
得到插值函数和f (x)图形:
三次样条插值，对每个 n值，分别画出插值函数即 f(x)的图形。
我们先求多项式插值：
在MATLAB的Editor中建立一个多项式的 M-file,输入如下的命令(如牛顿插值公式)
fun ctio n [ C,D]=n ewpoly(X,Y)
n=le n gth(X);
D=zeros (n,n)
D(:,1)=Y'
for j=2: n
for k=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);
for k=( n-1):-1:1
C=co n v(C,poly(X(k)))
m=le n gth(C);
C(m)= C(m)+D(k,k);
end
当n=10时，我们在 Comma nd Win dow中输入以下的命令： clear,clf,hold on;
X=-1:021;
Y =1./(1+25*X. A2);
[C,D]=newpoly(X, Y);
x=-1:0.01:1;
y=polyval(C,x);
plot(x,y,X, Y,'.');
grid on;
xp=-1:0.2:1;
z=1./(1+25*xp.A2);
plot(xp, z,'r')
当n=20时，我们在 Comma nd Win dow中输入以下的命令：clear,clf,hold on;
X=-1:0.1:1;
Y =1丿(1+25*X. A2);
[C,D]=newpoly(X, Y);
x=-1:0.01:1;
y=polyval(C,x);
plot(x,y,X, Y,'.');
grid on;
xp=-1:0.1:1;
z=1./(1+25*xp.A2);
plot(xp, z,'r')
得到插值函数和f (x)图形:
F面再求三次样条插值函数，在MATLAB的
F面再求三次样条插值函数，在
MATLAB的Editor中建立一个多项式的
M- file,
输入下列程序代码:
fun ctio n S=csfit(X,Y,dxO,dx n)
N=le n gth(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));
for k=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);
for k=N-2:-1:1
1
1
1
1
M(k+1)=(U(k)-C(k)*M(k+2))/B(k);
end
M(1)=3*(D(1)-dxO)/H(1)-M(2)/2; M(N+1)=3*(dxn-D(N))/H(N)-M(N)/2;
for k=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
clear,clc当 n=10 时，我们在 Comma nd Win dow 中输入以下的命令 : X=-1:0.2:1;
clear,clc
A
Y=1./(25*X. A2+1);
dx0= 0.0739644970414201;dx n= -0.0739644970414201; S=csfit(X,Y,dx0,dx n)
x1=-1:0.01:-0.5;y 仁 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,'.') 结果如图：
clear,clc
X=-1:0.1:1;
Y =1./(25*X92+1);
dxO= 0.0739644970414201;dx n= -0.0739644970414201;
S=csfit(X，丫 ,dxO,dx n)
x1= -1:0.01:-0.5;y 仁 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,'.')
结果如图：
1 0 8 0 6 O 4 0 ? 0 O 2 0 A G 6 O 8
实验三:
F列数据点的插值
X
0J
1
4
9
16
25
36
49
64
y
0
1
2
3
4
5
6
7
8
诃以得到平方根函数的近似，*区间 [0,64]上作图。
(1)用这9务点作8次多项式扌rfi伯 L 8(x),
[0,1:上，两(2)用三次样条(自然边界条件)程序求S (x)。从结果看在[0,64
[0,1:上，两
n
L8（X）可由公式 L n（x）= y klk（xj ）得出
k 0
三次样条可以利用自然边界条件。写成矩阵:
2
0
0 0 0
M
d°
1
2
1 0 0
M
d1
0
2
2 2 0
M
d2
0
0
3 2 3
M
d3
0
0
0 4 2
M
d4
h j-1
hj
其中j =
hj-1 hi
，i=
h
i -1 h
，dj=6f[x
n= 0=0 d 0=dn=0
O 匕-1)(Xx -0)Q4)( X 19Xx 4)(X6)(x 9)(￡5)( x 25)3欽 x36)(49499( x 娜
TOC \o "1-5" \h \z = 但16(0 0)(1@)(0 1)(9员0 4)(16>(0 9)(165)(25)(36)(0 3649)60~~494(16 64)
(x -0 )(xx-0)(xl)( x19((x 4)佑)(x 9)(x>5)(1X)(x 36)36)( x 49)(x x 64)64)
h(xx) (251 0)0)(25 1)(25 9)(4)(256)(1 9)(25)(1 16)3250 36)(25 649)(25 64)
= (x-0嬉-1)发(% )(x 4)(x?)(x 9)()25)(1X)(x36)( 25)()499()(x 64) 64)
12( x) (36 00)(36 1)(卒)(39)(44)(366)(4 9)(鉤(4 16)3360 25)936—49)36 64)"
=6(X) (x -0)(x x-0)1x(x1)(x 4)(x 41x)(x 9)(25)(1X)(x36)( 25)( x4360xx 6464)
= (90)(90)(49 1)(91)(44)(94)(496)(9 9)(25)(9 16)349(9 25949 34)(49 64)
| 3( X) (X - -1)(x 4)( x 9)( x 16)(x 25)( X 36)( X 49)
=7(X) (64 0)(x)(64 1)(64 4)(64 9)(64 16)(64 25)(64 36)(64 49)"
|8(X)
L8(
L8( x)= 1 1 (X)+2 I 2( x)+3 I 3( x)+4 I
4(x)+5 I
5( x)+6 I 6(x)+7 I 7(x)+8 I
8(X)
2.0000
0
0
0
0
0
0
0
0
Mo
0
0.2500
2.0000
0.7500
0
0
0
0
0
0
M !
1.0
0
0.3750
2.0000
0.6250
0
0
0
0
0
0.1
0
0
0.4167
2.0000
0.5833
0
0
0
0
0.0286
0
0
0
0.4375
2.0000
0.5625
0
0
0.0119
0
0
0
0
0.4500
2.0000
0.5500
0
0.0061
0
0
0
0
0
0.4583
2.0000
0.5417
0.0035
0
0
0
0
0
0
0.4634
2.0000
0.5357 M
7
0.0022
0
0
0
0
0
0
0
0
2.0000
0
求三次样条插值函数由 MATLA时算:
可得矩阵形式的线性方程组为
%n代表元素数，
%n代表元素数， s,t代表端点的一阶导
fun cti on tgsa nci(n, s,t)
y=[0 1 2 3 4 5 6 7 8];
x=[0 1 4 9 16 25 36 49 64];
n=9
for j=1:1: n-1
h(j)=x(j+1)-x(j); |
end
for j=2:1: n-1
r(j)=h(j)/(h(j)+h(j-1)); |
end
for j=1:1: n-1
u(j)=1-r(j);
end
for j=1:1: n-1
f(j)=(y(j+1)-y(j))/h(j); end
for j=2:1: n-1 d(j)=6*(f(j)-f(j-1))/(h(j-1)+ha)) end
d(1)=0
d(n )=0
a=zeros (n,n);
for j=1:1: n
a(j,j)=2;
end
r(1)=0; u( n )=0;
for j=1:1: n-1
a(j+1,j)=u(j+1);
a(j,j+1)=r(j);
end
b=i n v(a)
for j=1:1: n-1 p(j,1)=m(j)/(6*h(j)); p(j,2)=ma+1)/(6*h(j));
P(j , 3)=(y(j)-m(j)*(h(j)A2/6))/h(j); P (j ,
A
4)=(y(j+1)-m(j+1)*(h(j) 2/6))/h(j); end
P
解得:
MO=O;M i =-0.5209;M 2=0.0558;M 3=-0.0261;M 4=0.0006;M 5=-0.0029;M 6=-0.0008;M 7=--0.0009;
-0), x [0,1]0(1 x )3 0.0868(x 0) 3 0(1 x
-0), x [0,1]
0.0019(9 x )3
0.0009(x -4 )3
x [4,9]
-0.0006(16 x )3
0(x 9)3 0.4590
(16
x) - 0.5708(x-9)
，x [9,16]
0(25
x)3 - 0.0001
(x
16)3-0.4436 (
25
x) 0.5600(x -16)，
x [16,25]
0(36
x)3- 0( x
25)3
0.4599 (36
x)
0.5470(x -25) ，x
[25,36]
0(48
x)3- 0( x
36)3
0.4633 (48
x)
0.5404(x -36) ，x
[36,48]
0(64
x)3 0(x
48)3
0.4689 ( 64
x)
0.5333(x -48) ，x
[48,64]
-0.0289(4 x )3 - 0.003 ( x -1)3
0.5938( 4 x ) 0.6388 (x 1 )
x [1,4]
M=0,则三次样条函数：
%画图形比较那个插值更精确的函数：
x0=[0 1 4 9 16 25 36 49 64];
y0=[0 1 2 3 4 5 6 7 8]; 0.3535( 9 x ) 0.6271 ( x 4 )
S(x)=x=0：64;y=lagr1(x0,y0,x);
plot(x0,y0,'o')
hold on
plot(x,y,'r');
hold on;
F面进行画图，在 Comma nd Win dow 中输入画图的程序代码：
pp=csape(x0,y0,'variatio nal')
fnplt(pp,'g');
hold on;
plot(x0,y0,':b');hold on
%axis([0 2 0 1]); %看[0 1]区间的图形时加上这条指令
gtext ('三次样条插值')
gtext ('原图像')
gtext ('拉格朗日插值')
%F面是求拉格朗日插值的函数
0
0.2
0.4
0.6
0.8
1
100
"T"
L
~r~
L
80
—
60
拉格
朗日插值
40
■
20
T fW
0
■
g 橡
壬人 vr
1 IH
on
I
c
I
0 10 20 30 40 50 60 70
图3 — 2
由图3-1可以看出，红色的线条与蓝色虚线条几乎重合， 所以可知拉格朗日插值函数的曲线
更接近开平方根的函数曲线，在 ［0, 1:朗格朗日插值更精确。而在区间 ［0, 64］上从图3-2中可以
看出蓝色虚线条和绿色线条是几乎重合的，而红色线条在 ［30, 70］之间有很大的振荡，所在在区
间［0, 64］三次样条插值更精确写。
【结论】（结果）
分段低次插值则单个多项式高次插值效果并不理想，有龙格现象，偏差大，没有使用价值。而
分段低次插值则
稳定性，与分段低次插值相比较精确度较高，拟合效果较好，而三次样条插值具有良好的收敛性与 光滑度更高，而且提供的信息也相对少一些。我们 可以看到，在以上的三道实验题里，我们可以从 图形中看出，三次样条的拟合程度 是三种插值函数里最好的。
稳定性，与分段低次插值相比较
小结】
通过此次实验，我对牛顿多项式插值，三次样条插值，拉格朗日插值有了更进 一步的了解，知 道了三次样条的拟合程度在高次的情况下更高，在理论上和应用上 都有重要意义，在利用计算机编 程软件进行高次插值的时候，我们可以多考虑利用 三次样条进行插值。
指导教师评语及成绩
成绩 : 指导教师签名 : 批阅日期：
- 下载文档
- 收藏
- 0
推荐访问:实验报告 数值 实验 报告 最新数值分析实验报告176453