实验二、插值法
插值法是函数逼近的一种重要方法,它是数值积分、微分方程数值解等数值计算的基础与工具,其中多项式插值是最常用和最基本的方法。拉格朗日插值多项式的优点是表达式简单明确,形式对称,便于记忆,它的缺点是如果想要增加插值节点,公式必须整个改变,这就增加了计算工作量。而牛顿插值多项式对此做了改进,当增加一个节点时只需在原牛顿插值多项式基础上增加一项,此时原有的项无需改变,从而达到节省计算次数、节约存储单元、应用较少节点达到应有精度的目的。
一、实验目的
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 ′ ′ ( x ) ∣ ≤ 5 | f ^ { \prime \prime } ( x ) | \leq 5 ∣f′′(x)∣≤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 f ( x ) = e ^ { - x } 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
即 R 3 ≈ 0.04 f ( 4 ) ( ξ ) , ξ ∈ ( − 2.00 , 2.00 ) R 3 \approx 0.04 f ^ { ( 4 ) } ( \xi ) , \quad \xi \in ( - 2.00,2.00 ) R3≈0.04f(4)(