计算物理学（数值分析）上机实验答案2、插值法

实验二、插值法

​ 插值法是函数逼近的一种重要方法，它是数值积分、微分方程数值解等数值计算的基础与工具，其中多项式插值是最常用和最基本的方法。拉格朗日插值多项式的优点是表达式简单明确，形式对称，便于记忆，它的缺点是如果想要增加插值节点，公式必须整个改变，这就增加了计算工作量。而牛顿插值多项式对此做了改进，当增加一个节点时只需在原牛顿插值多项式基础上增加一项，此时原有的项无需改变，从而达到节省计算次数、节约存储单元、应用较少节点达到应有精度的目的。

一、实验目的

1、理解插值的基本概念，掌握各种插值方法，包括拉格朗日插值和牛顿插值等， 注意其不同特点；

2、通过实验进一步理解并掌握各种插值的基本算法。

二、算法实例

1. 与插值有关的 MATLAB 函数

(一) POLY2SYM 函数

调用格式一：poly2sym ©

调用格式二：f1=poly2sym(C,‘V’) 或 f2=poly2sym(C, sym (‘V’) ),

(二) POLYVAL 函数

调用格式：Y = polyval(P,X)

(三) POLY 函数

调用格式：Y = poly (V)

(四) CONV 函数

调用格式：C =conv (A, B)

(五) DECONV 函数

调用格式：[Q,R] =deconv (B,A)

(六) roots(poly(1:n))命令

调用格式：roots(poly(1:n))

(七) det(a*eye(size (A)) - A)命令

调用格式：b=det(a*eye(size (A)) - A)

2.线性插值及其 MATLAB 程序

例 2.1

​ 已知函数 f (x) 在[1,3] 上具有二阶连续导数，$|f^{\prime \prime }(x)|\le 5$，且满足条件 f (1) = 1, f (3) = 2 。求线性插值多项式和函数值 f (1.5) ，并估计其误差。

X=[1,3];Y=[1,2];
l01= poly(X(2))/( X(1)- X(2));
l11= poly(X(1))/( X(2)- X(1));
l0=poly2sym(l01),l1=poly2sym(l11)
P=l01*Y(1)+l11*Y(2),L=poly2sym (P)
x=1.5;Y = polyval(P,x)
%运行后输出基函数 l0 和 l1 及其插值多项式的系数向量 P（略）、插值多项式L和插值Y、为
%L = x/2 + 1/2		Y = 1.2500
M=5;R1=M*abs((x-X(1))* (x-X(2)))/2
%运行后输出误差限为R1 = 1.8750

例 2.2

​ 求函数$f(x)=e^{-x}$在[0,1] 上线性插值多项式，并估计其误差。

 X=[0,1];Y=exp(-X)
l01= poly(X(2))/( X(1)- X(2));
l11= poly(X(1))/( X(2)- X(1));
l0=poly2sym (l01),l1=poly2sym (l11),
P = l01* Y(1)+l11* Y(2),L=poly2sym (P)
%运行后输出基函数 l0 和 l1 及其插值多项式的系数向量 P 和插值多项式 L为
% l0 =1 - x		l1 =x
% P =-0.6321    1.0000
% L =1 - (1423408956596761*x)/2251799813685248
M=1;x=0:0.001:1; R1=M*max(abs((x-X(1)).*(x-X(2))))./2
%运行后输出误差限为R1 = 0.1250

3.抛物线插值及其 MATLAB 程序

例 2.3

​ 求将区间 [0, π/2] 分成n 等份(n = 1,2) ，用 y = f (x) = cosx 产生n +1个节点，然后分别作线性插值函数 P1 (x) 和抛物线插值函数 P2 (x) .用它们分别计算 cos (π/6) (取四位有效数字)，并估计其误差。

X=[0,pi/2];Y=cos(X);
l01= poly(X(2))/( X(1)- X(2));
l11= poly(X(1))/( X(2)- X(1));
l0=poly2sym (l01),l1=poly2sym (l11)
P = l01* Y(1)+ l11* Y(2), L=poly2sym (P)
x=pi/6;Y = polyval(P,x)
M=1;x=pi/6; R1=M*abs((x-X(1))*(x-X(2)))/2
%输出结果为
% l0 =1 - (5734161139222659*x)/9007199254740992 
% l1 =(5734161139222659*x)/9007199254740992
% P =-0.6366    1.0000 
% L =1 - (5734161139222659*x)/9007199254740992
% Y =0.6667
% R1 =0.2742

X=0:pi/4:pi/2;Y=cos(X);
l01= conv (poly(X(2)),poly(X(3)))/(( X(1)- X(2))* ( X(1)- X(3)));
l11= conv (poly(X(1)),poly(X(3)))/(( X(2)- X(1))* ( X(2)- X(3)));
l21= conv (poly(X(1)),poly(X(2)))/(( X(3)- X(1))* ( X(3)- X(2)));
l0=poly2sym (l01),l1=poly2sym (l11),l2=poly2sym (l21)
P = l01* Y(1)+ l11* Y(2) + l21* Y(3), L=poly2sym (P)
x=pi/6;Y=polyval(P,x)
M=1;x=pi/6; R2=M*abs((x-X(1))*(x-X(2)) *(x-X(3)))/6
%输出结果为
% l0 =(228155022448185*x^2)/281474976710656 - (2150310427208497*x)/1125899906842624 + 1
% l1 =(5734161139222659*x)/2251799813685248 - (228155022448185*x^2)/140737488355328
% l2 =(228155022448185*x^2)/281474976710656 - (5734161139222659*x)/9007199254740992
% P =-0.3357   -0.1092    1.0000
% L =1 - (6048313895780875*x^2)/18014398509481984 - (7870612110600739*x)/72057594037927936
% Y =0.8508
% R2 =0.0239

4.n 次拉格朗日（Lagrange）插值及其 MATLAB 程序

例 2.4

​ 给出节点数据 f (−2.00) = 17.00 ， f (0.00) = 1.00 ， f (1.00) = 2.00 ，f (2.00) = 17.00，作三次拉格朗日插值多项式计算 f (0.6)，并估计其误差。

X=[-2,0,1,2]; Y =[17,1,2,17];
p1=poly(X(1)); p2=poly(X(2));
p3=poly(X(3)); p4=poly(X(4));
l01= conv ( conv (p2, p3), p4)/(( X(1)- X(2))* ( X(1)- X(3)) * ( X(1)- X(4)));
l11= conv ( conv (p1, p3), p4)/(( X(2)- X(1))* ( X(2)- X(3)) * ( X(2)- X(4)));
l21= conv ( conv (p1, p2), p4)/(( X(3)- X(1))* ( X(3)- X(2)) * ( X(3)- X(4)));
l31= conv ( conv (p1, p2), p3)/(( X(4)- X(1))* ( X(4)- X(2)) * ( X(4)- X(3)));
l0=poly2sym (l01),l1=poly2sym (l11),l2=poly2sym (l21), l3=poly2sym (l31)
P = l01* Y(1)+ l11* Y(2) + l21* Y(3) + l31* Y(4),L=poly2sym (P)
x=0.6; Y = polyval(P,x);
syms M; x=0.6;
R3=M*abs((x-X(1))*(x-X(2)) *(x-X(3)) *(x-X(4)))/24
%运行后输出为
% l0 =- x^3/24 + x^2/8 - x/12 
% l1 =x^3/4 - x^2/4 - x + 1 
% l2 =- x^3/3 + (4*x)/3 
% l3 =x^3/8 + x^2/8 - x/4
% P =1     4    -4     1 
% L =x^3 + 4*x^2 - 4*x + 1
% Y =0.2560 
% R3 =(91*M)/2500

即$R3\approx 0.04f^{(4)}($