实验一---函数插值方法报告.doc

实验一 函数插值方法报告

一、问题提出
对于给定的一元函数的n+1个节点值。试用Lagrange公式求其插值多项式或分段二次Lagrange插值多项式。
数据如下：
（1）

0.4
0.55
0.65
0.80
0.95
1.05

0.41075
0.57815
0.69675
0.90
1.00
1.25382

求五次Lagrange多项式，和分段三次插值多项式，计算,
的值。（提示：结果为, ）
（2）

1
2
3
4
5
6
7

0.368
0.135
0.050
0.018
0.007
0.002
0.001
试构造Lagrange多项式，计算的,值。（提示：结果为, ）
二、要求
1、 利用Lagrange插值公式
编写出插值多项式程序；
2、 给出插值多项式或分段三次插值多项式的表达式；
3、 根据节点选取原则，对问题（2）用三点插值或二点插值，其结果如何；
4、 对此插值问题用Newton插值多项式其结果如何。Newton插值多项式如下：

其中：
三、目的和意义
1、 学会常用的插值方法，求函数的近似表达式，以解决其它实际问题；
2、 明确插值多项式和分段插值多项式各自的优缺点；
3、 熟悉插值方法的程序编制；
4、 如果绘出插值函数的曲线，观察其光滑性。
四、实验学时：2学时
五、实验步骤：
1．进入C或matlab开发环境；
2．根据实验内容和要求编写程序；
3．调试程序；
4．运行程序；
5．撰写报告,讨论分析实验结果.

解：
一、编写插值函数结构程序
Lagrange插值多项式M文件：lagrange1.m
function [A1,LN,L1,B1]=lagrange1(X,Y)
m=length(X); LN=ones(m,m);
for k=1: m
x1=1;
for i=1:m
if k~=i
x1=conv(x1,poly(X(i)))/(X(k)-X(i));
end
end
L1(k,:)=x1; B1(k,:)=poly2sym (x1)
end
A1=Y*L1;LN=Y*B1

分段三次艾尔米特插值多项式的M文件：Hermite3.m
function [f,ff] = Hermite3(x,y,y1)
syms t;
f = 0.0;
if(length(x) == length(y))
if(length(y) == length(y1))
n = length(x);
else
disp('y和y的导数的维数不相等');
return;
end
else
disp('x和y的维数不相等！ ');
return;
end
for i=1:n
h = 1.0;
a = 0.0;
for j=1:n
if( j ~= i)
h = h*(t-x(j))^2/((x(i)-x(j))^2);
a = a + 1/(x(i)-x(j));
end
end
f = f + h*((x(i)-t)*(2*a*y(i)-y1(i))+y(i));
end
ff = subs(f,'t');
（1）、求五次Lagrange多项式，和分段三次插值多项式。
在主显示区，输入五次Lagrange多项式程序：
>> X=[0.4 0.55 0.65 0.80 0.95 1.05];
>> Y=[0.41075 0.57815 0.69675 0.90 1.00 1.25382];
>> [A1,LN,L1,B1]=lagrange1(X,Y)
>> plot(X,A1);
>> F=poly2sym(A1)

运行后，输出五次Lagrange多项式的结果：
A1 =

121.6264 -422.7503 572.5667 -377.2549 121.9718 -15.0845

F =

(2139673480305281*x^5)/17592186044416 - (1859275536318005*x^4)/4398046511104 + (9836621836743*x^3)/17179869184 - (414796119737013*x^2)/1099511627776 + (2145751274873259*x)/17592186044416 - 1061478972867847/70368744177664
拉格朗日插值多项式的图如下：

在主显示区，输入分段三次艾尔米特插值多项式的程序：
>> x=[0.4 0.55 0.65 0.80 0.95 1.05];
>> y=[0.41075 0.57815 0.69675 0.90 1.00 1.25382];
>> y1=[2.3440 0.9032 1.4329 0.9903 0.9170 5.1439];
>> [f,ff] = Hermite3(x,y,y1);
>> ff
运行后，分段三次艾尔米特插值多项式的输出结果：

ff =

6400000000*(t - 4/5)^2*(t - 11/20)^2*(t - 13/20)^2*(t - 19/20)^2*(t - 21/20)^2*((2240245151070481*t)/140737488355328 - 52393133567890089/8796093022208000))/184041 - (16000000*((6348013345609171*t)/140737488355328 - 85523418631741336287/1759218604441600000)*(t - 2/5)^2*(t - 4/5)^2*(t - 11/20)^2*(t - 13/20)^2*(t - 19/20)^2)/169 + (16000000*((4105617466549689*t)/281474976710656 - 5238387122042657959/703687441776640000)*(t - 2/5)^2*(t - 4/5)^2*(t - 13/20)^2*(t - 19/20)^2*(t - 21/20)^2)/9 - (256000000*((35097*t)/10000 - 46347/12500)*(t - 2/5)^2*(t - 11/20)^2*(t - 13/20)^2*(t - 19/20)^2*(t - 21/20)^2)/81 - (400000000*((13147*t)/20000 - 449611/400000)*(t - 2/5)^2*(t - 4/5)^2*(t - 11/20)^2*(t - 19/20)^2*(t - 21/20)^2)/81 - (10000000000*((84913*t)/11000 - 1833347/220000)*(t - 2/5)^2*(t - 4/5)^2*(t - 11/20)^2*(t - 13/20)^2*(t - 21/20)^2)/9801
分段三次艾尔米特插值多项式的图如下：

(2)、试构造Lagrange多项式。
结果为,
在主显示区，输入程序：
>> X=[1 2 3 4 5 6 7];
>> Y=[0.368 0.135 0.050 0.018 0.007 0.002 0.001];
>> [A1,LN,L1,B1]=lagrange1(X,Y)
>> plot(X,A1);
>> F=poly2sym(A1)
运行后，输出结果的Lagrange多项式的结果：
A1 =

0.0001 -0.0016 0.0186 -0.1175 0.4419 -0.9683 0.9950

F =

(4304240283865561*x^6)/73786976294838206464 - (7417128346304051*x^5)/4611686018427387904 + (223*x^4)/12000 - (2821*x^3)/24000 + (994976512675275*x^2)/2251799813685248 - (19367*x)/20000 + 199/200
Lagrange多项式的图如下：

二、计算函数值

计算函数值的主程序：lagrangezhi.m
function [y,R]=lagrangezhi(X,Y,x,M)
n=length(X); m=length(x);
for i=1:m
z=x(i);s=0.0;
for k=1:n
p=1.0; q1=1.0; c1=1.0;
for j=1:n
if j~=k
p=p*(z-X(j))/(X(k)-X(j));
end
q1=abs(q1*(z-X(j)));c1=c1*j;
end
s=p*Y(k)+s;
end
y(i)=s;
end
R=M*q1/c1;
（1）、计算、 的值。
在主显示区，输入程序：
>> x=0.596; M=1; X=[0.4,0.55,0.65,0.80,0.95,1.05];
>> Y=[0.41075,0.57815,0.69675,0.90,1.00,1.25382];
>> [y,R]=lagrangezhi(X,Y,x,M)
运行结果：
y =
0.6257
R =
2.2170e-008
在主显示区，输入程序：
>> x=0.99; M=1; X=[0.4,0.55,0.65,0.80,0.95,1.05];
>> Y=[0.41075,0.57815,0.69675,0.90,1.00,1.25382];
>> [y,R]=lagrangezhi(X,Y,x,M)
运行结果：
y =
1.0542
R =
5.5901e-008
（2）、计算、的值
在主显示区，输入程序：
>> x=1.8; M=1; X=[1,2,3,4,5,6,7];
>> Y=[0.368,0.135,0.050,0.018,0.007,0.002,0.001];
>> [y,R]=lagrangezhi(X,Y,x,M)
运行结果：
y =
0.1648
R =
0.0059
在主显示区，输入程序：
>> x=6.15; M=1; X=[1,2,3,4,5,6,7];
>> Y=[0.368,0.135,0.050,0.018,0.007,0.002,0.001];
>> [y,R]=lagrangezhi(X,Y,x,M)
运行结果：
y =
0.0013
R =
0.0042
三、Newton插值多项式
Newton插值多项式主程序M文件：Newton.m
function [A,C,L,wcgs,Cw]= Newton(X,Y)
n=length(X); A=zeros(n,n); A(:,1)=Y';
s=0.0; p=1.0; q=1.0; c1=1.0;
for j=2:n
for i=j:n
A(i,j)=(A(i,j-1)- A(i-1,j-1))/(X(i)-X(i-j+1));
end
b=poly(X(j-1));q1=conv(q,b); c1=c1*j; q=q1;
end
C=A(n,n); b=poly(X(n)); q1=conv(q1,b);
for k=(n-1):-1:1
C=conv(C,poly(X(k))); d=length(C); C(d)=C(d)+A(k,k);
end
L(k,:)=poly2sym(C); Q=poly2sym(q1);
syms M
wcgs=M*Q/c1; Cw=q1/c1;
在主显示区，输入的程序:
>> x=[0.4 0.55 0.65 0.80 0.95 1.05];
>> y=[041075 0.57815 0.69675 0.90 1.00 1.25382];
>> [A,C,L,wcgs,Cw]= Newton(x,y)
>> syms x;
>> ezplot(L,[0 1.1]);
运行结果如下，得到
A =

0.4108 0 0 0 0 0
0.5782 1.1160 0 0 0 0
0.6967 1.1860 0.2800 0 0 0
0.9000 1.3550 0.6760 0.9900 0 0
1.0000 0.6667 -2.2944 -7.4261 -15.3020 0
1.2538 2.5382 7.4861 24.4514 63.7551 121.6264

C =

121.6264 -422.7503 572.5667 -377.2549 121.9718 -15.0845

L =

(8558693921221117*x^5)/70368744177664 - (3718551072636019*x^4)/8796093022208 + (5036350380412441*x^3)/8796093022208 - (3318368957896111*x^2)/8796093022208 + (536437818718315*x)/4398046511104 - 8491831782942691/562949953421312

wcgs =

(M*(x^6 - (22*x^5)/5 + (1583*x^4)/200 - (3721*x^3)/500 + (542206127247039*x^2)/140737488355328 - (4682696525551953*x)/4503599627370496 + 4111390143022055/36028797018963968))/720

Cw =

0.0014 -0.0061 0.0110 -0.0103 0.0054 -0.0014 0.0002

牛顿插值多项式的图如下：

在主显示区，输入的程序:
>> x=[1 2 3 4 5 6 7];
>> y=[0.368 0.135 0.050 0.018 0.007 0.002 0.001];
>> [A,C,L,wcgs,Cw]= Newton(x,y)
>> syms x;
>> ezplot(L,[0 8]);
运行结果如下，得到：
A =

0.3680 0 0 0 0 0 0
0.1350 -0.2330 0 0 0 0 0
0.0500 -0.0850 0.0740 0 0 0 0
0.0180 -0.0320 0.0265 -0.0158 0 0 0
0.0070 -0.0110 0.0105 -0.0053 0.0026 0 0
0.0020 -0.0050 0.0030 -0.0025 0.0007 -0.0004 0
0.0010 -0.0010 0.0020 -0.0003 0.0005 -0.0000 0.0001

C =

0.0001 -0.0016 0.0186 -0.1175 0.4419 -0.9683 0.9950

L =

(8608480567731121*x^6)/147573952589676412928 - (7417128346304047*x^5)/4611686018427387904 + (223*x^4)/12000 - (2821*x^3)/24000 + (7959812101402191*x^2)/18014398509481984 - (19367*x)/20000 + 199/200

wcgs =

(M*(x^7 - 28*x^6 + 322*x^5 - 1960*x^4 + 6769*x^3 - 13132*x^2 + 13068*x - 5040))/5040

Cw =

0.0002 -0.0056 0.0639 -0.3889 1.3431 -2.6056 2.5929 -1.0000
牛顿插值多项式的图如下

四、结果讨论和分析
通过上机实验，我对主要函数常见的插值方法有了更深的认识。在编辑拉格朗日插值、牛顿插值以及分段插值MATLAB程序过程中，深入了解了它们各自的算法，为我以后在工程中应用插值函数奠定了基础。
根据得出的数据知道，与拉格朗日插值相比牛顿插值的计算量要小，程序设计方便，并且牛顿插值曲线相对要光滑些。分段三次埃尔米特插值与分段低次插值函数都有一致收敛性但是低次插值要求给出节点上的导数值，其光滑性较差，改进这种插值以克服其缺点要用三次样条插值。三次样条插值具有良好的收敛性和稳定性，又有二阶光滑度。