10 多项式插值

1、1第一节第一节 多项式插值多项式插值2 Approximation theory is concerned with how functions can best be approximated with simpler functions, and with quantitatively characterizing the errors introduced thereby. 逼近理论研究如何将函数利用一组简单函数近似表征，并逼近理论研究如何将函数利用一组简单函数近似表征，并定量分析逼近过程中产生的误差。定量分析逼近过程中产生的误差。 数值逼近包括两大类：插值和拟合数值逼近包括两大类：插值

2、和拟合 应用：应用： 问题问题1：已知美国：已知美国1900年年2000年人口数，分析美国人口变化规律，年人口数，分析美国人口变化规律，预测美国未来人口数目；预测美国未来人口数目； 问题问题2：已知某元件在电压为：已知某元件在电压为v1vn时，其电流为时，其电流为i1in，求该元件，求该元件的伏安特性曲线。的伏安特性曲线。3已经测得在某处海洋不同深度处的水温如下：深度（M） 466 741 950 1422 1634水温（oC）7.04 4.28 3.40 2.54 2.13根据这些数据，希望合理地估计出其它深度（如500米，600米，1000米）处的水温举例这就是本章要讨论的“插值问题”4

3、当极少量的高纯氢气在高真空玻璃管中，加入高电压使之放电，管中发出光束，使这种光经过分光作用，在可见光区得到四条颜色不同的谱线。5 实验总结： 1885年，瑞士一位中学数学教师J.J.Balmar(巴尔默)指出，上述谱线的频率符合下列公式： 由此公式可算出： 当n3时，是H的频率 当n4时，是H的频率 当n5时，是H的频率 当n6时，是H的频率 1522113.289 10 ()2vHzn6 插值（插值（interpolate） 已知函数在已知函数在xi处的值为处的值为 yi ，求，求 f (x)，使之满足：，使之满足： yi = f (xi) 其中其中， f (x)为为插值函数插值函数， xi

4、处为处为插值节点插值节点，插值节点的区间称为，插值节点的区间称为插值区间插值区间， yi = f (xi)为为插值条件插值条件。 拟合（拟合（fit） 已知函数在已知函数在xi处的值为处的值为 yi ，求，求 f (x)，使之满足：，使之满足： e =yi - f (xi) 在给定的在给定的准则准则下最小。下最小。 差异：差异： 插值函数必插值函数必须经过插值点。须经过插值点。拟合函数不必拟合函数不必经过拟合点。经过拟合点。7 在已知在已知 的前提下，有多少函数满足的前提下，有多少函数满足yi = f (xi) ？ 给定任意一组给定任意一组，存在无穷多函数满足，存在无穷多函数满足yi = f

5、(xi) ，因，因此，在解决插值问题前，必须首先明确所采用的插值函数。此，在解决插值问题前，必须首先明确所采用的插值函数。 常用插值函数：常用插值函数： 多项式函数；多项式函数；xn sinc函数等；函数等；sin（x）/ x8 问题描述：问题描述： 给定插值点给定插值点，构造多项式函数，构造多项式函数 Pn(x) = a0 + a1x + a2x2+ anxn，使之满足：，使之满足： Pn(xi) = yi (i = 0,1,2,n)。 如何计算：如何计算： 多项式多项式Pn(x) 由其多项式系数由其多项式系数a0 , a1, a2, , an决定，只需要求解多项决定，只需要求解多项式系数，

6、即可获得该插值多项式。式系数，即可获得该插值多项式。 将将Pn(xi) = yi写为矩阵形式可得：写为矩阵形式可得：1001 000001111 1011111010.1.1.1nnnnnnnnnnnnnnnnaya xa xayxxya xa xayxxaaya xa xayxx求解该线性方程组即可得到多项式的系数求解该线性方程组即可得到多项式的系数范德蒙矩阵（范德蒙矩阵（ Vandermonde）9 该线性方程组有解吗，解唯一吗？该线性方程组有解吗，解唯一吗？ 唯一性定理：通过唯一性定理：通过n+1个节点的个节点的n阶插值多项式存在且唯一。阶插值多项式存在且唯一。0011.1.1det()

7、0.1nnijijnnnxxxxxxxx范德蒙矩阵的行列式的值为范德蒙矩阵的行列式的值为xi- xj 的连乘积，当的连乘积，当 xixj时，该时，该行列式的值不为零，即线性方程组有解，且存在唯一解行列式的值不为零，即线性方程组有解，且存在唯一解证明：证明：10 例：已知 x = 0 , 0.2, 0.4, 0.6, 0.8, 1.0，y = 1.0, 1.17, 1.70, 2.59, 3.93, 6.00，求5阶多项式插值。00.20.40.60.810123456711clear allclose allclcx = 0 : 0.2 : 1% display(y = 3 * x.5 - 4

8、 * x.4 + 2 * x.3 + 4 * x.2 + 1) y = 3 * x.5 - 4 * x.4 + 2 * x.3 + 4 * x.2 + 1 A = vander(x)display(strcat(cond(A) = , num2str(cond(A) display( ) display(P = A y. )P = A y. ;P = P. plot(x, y, o) hold onxx = 0 : 0.01 : 1;yy = P(1) * xx.5 + P(2) * xx.4 + P(3) * xx.3 + P(4) * xx.2 + P(5) * xx.1 + P(6);p

9、lot(xx, yy)12 通过解方程组求得插值多项式通过解方程组求得插值多项式pn(x)的方法并不可的方法并不可取取.这是因为当这是因为当n较大时解方程组的计算量较大，而较大时解方程组的计算量较大，而且方程组系数矩阵的条件数一般较大（可能且方程组系数矩阵的条件数一般较大（可能是病态是病态方程组）。方程组）。为此我们必须从其它为此我们必须从其它途径来求途径来求Pn(x)：不通过求解方程组而不通过求解方程组而获得插值多项式获得插值多项式13基本思想：基本思想：在n次多项式空间P Pn n中找一组合适的基函数中找一组合适的基函数 0 0( (x),x), 1 1(x), (x), 3 3（x),x

10、),使使pn(x)=a0 0 0( (x)x) +a1 1 1(x)(x) +an 3 3（x)x)不同的基函数的选取导致不同的不同的基函数的选取导致不同的插值方法插值方法Lagrange插值插值Newton插值插值14n = 1使得使得可见可见 P1(x) 是过是过 ( x0 , y0 ) 和和 ( x1, y1 ) 两点的直线。两点的直线。)()(0010101xxxxyyyxP 101xxxx 010 xxxx = y0 + y1l0(x)l1(x) 10)(iiiyxlLagrangeLagrange插值插值求求 n 次多项式次多项式 使得使得nnnxaxaaxP 10)( ),0,1

11、,niiP xy in已知已知 x0 , x1 ; y0 , y1 ，求求101( )P xaa x100111(),()P xy P xy使得15 多项式插值的拉格朗日多项式表示：多项式插值的拉格朗日多项式表示： 给定插值点给定插值点，其插值多项式可表示为：，其插值多项式可表示为：0 01 10( )( )( ) .( )( )nnn nk kkP xy l xyl xy l xy l x0001111()().()( )=().()()()()()()iiiiinjnijiij knkjixxxxxxl xxxxxxxxxxxxxxx，其中101101111101 ( )()().() (

12、 )().()().()( ) ( )()( ) ( )()()nnniiiiiiinnininnnkkknkixxxxxxxxxxxxxxxxxl xxxxxP xyxxx：令则16 性质1：拉格朗日多项式在本插值点值为1，其它插值点值为0。1()0kjkjl xkj-0.2-0.100.10.2-0.6-0.4-0.200.20.40.60.81-0.1-0.0500.050.100.20.40.60.81二阶拉格朗日差值函数四阶拉格朗日差值函数01010110021201010210120122021( ) ( )()()()()( ) ( ) ()()()()()() ( )()()x

13、xxxlxl xxxxxxxxxxxxxlxl xxxxxxxxxxxxxlxxxxx：一阶插值公式二阶插值公式17 例：已知lg10=1,lg15=1.1761,lg20=1.3010,利用一次、二次多项式插值计算 lg12的近似值。0110 01 11020201101( )(20)( )(10)10 201020 101011.3010 ( )( )( )(20)(10)1010 xxxl xxl xxP xy l xyl xxx一阶插值，选择 和 则插值基函数为：0121020(15) (20)1( )(15)(20)(10 15) (1020)50(10) (20)1( )(10)(

14、20)(15 10) (1520)25(10) (15)1( )(10)(15)(20 10) (20 15)50 xxxlxxxxxl xxxxxlxxx二阶插值，选择 和 则插值基函数为：1820 01 12 212 ( )( )( )( )1 (20)(15)501.1761 (10)(20)251.3010 (10)(15)50 (12)1.0602 (12)1.0766 log(12)1.0792P xy lxy l xy lxxxxxxxPP最后得到: 10121416182011.051.11.151.21.251.3 log10P1(x)P2(X)X19clear allclo

15、se allclc% 插值x = 10, 15, 20y = 1,1.1761,1.3010 xx = 10 : 1 : 20;% 一阶插值yy1 = y(1) * (xx - x(3) / (x(1) - x(3) + y(3) * (xx - x(1) / (x(3) - x(1);% 一阶插值l1 = (xx - x(2) .* (xx - x(3) / (x(1) - x(2) / (x(1) - x(3);l2 = (xx - x(1) .* (xx - x(3) / (x(2) - x(1) / (x(2) - x(3);l3 = (xx - x(1) .* (xx - x(2)

16、/ (x(3) - x(1) / (x(3) - x(2);yy2 = y(1) * l1 + y(2) * l2 + y(3) * l3;figurehold onplot(xx, log10(xx), r);plot(xx, yy1, g);plot(xx, yy2, black);plot(x, y, b)20 定理：( )(1)01(1)(1)10( ) ( ) , ( ),.,( ) , ,( , ), ( )( )( )( )( )( )(1)!(1)!nnnnnnnnnnjjyf xnfxa bfxa baxxxbPxnxa ba bffR xf xPxxxxnn：()设函数的

17、阶导数在上连续，在（）存在，节点是 次拉格朗日插值多项式，则对任意的必存一点使插值余项1012012( )1( )()()2( )2( )()()()6fnR xxxxxfnR xxxxxxx时时(1)(1)110( ) ( ), ( , ),|(1)!nnnnniifxfxMMxa bxxn ，通常 未知，一般取的上界则有误差小于21 例：113sin, sin, sin624322已知 分别利用 sin x 的1次、2次 Lagrange 插值计算 sin 50 并估计误差。1111/6/4/41/615 ( ), ()0.77614/6/42/4/6182/4/3/31/435 ( ),

18、 () 0.76008/4/3/3/42182 exexininxxxPxPxxxPxP一阶插值，选择 和则插值基函数为：一阶插值，选择 和则插值基函数为：二阶插值，插值基函数为：225 ( ) , () 0.76543185 sin()=0.76604, 0.01001 0.005960.0006118PxP：，：，：外推误差 内插误差二阶误外推误差内插差 误差二阶误差22clear allclose allclcformat long% 插值x = 1/6, 1/4, 1/3 * pi;y = 1 / 2,1 / sqrt(2), sqrt(3) / 2;xx = 1/6 : 1/60 :

19、 1/3;xx = xx * pi;xx50 = 50 / 180 * pi;% 一阶插值 - exyy1ex = y(1) * (xx - x(2) / (x(1) - x(2) + y(2) * (xx - x(1) / (x(2) - x(1);yy150ex = y(1) * (xx50 - x(2) / (x(1) - x(2) + y(2) * (xx50 - x(1) / (x(2) - x(1)% 一阶插值 - inyy1in = y(2) * (xx - x(3) / (x(2) - x(3) + y(3) * (xx - x(2) / (x(3) - x(2);yy150i

20、n = y(2) * (xx50 - x(3) / (x(2) - x(3) + y(3) * (xx50 - x(2) / (x(3) - x(2)% % 一阶插值% l1 = (xx - x(2) .* (xx - x(3) / (x(1) - x(2) / (x(1) - x(3);% l2 = (xx - x(1) .* (xx - x(3) / (x(2) - x(1) / (x(2) - x(3);% l3 = (xx - x(1) .* (xx - x(2) / (x(3) - x(1) / (x(3) - x(2);% % yy2 = y(1) * l1 + y(2) * l2

21、 + y(3) * l3;% % l1 = (xx50 - x(2) .* (xx50 - x(3) / (x(1) - x(2) / (x(1) - x(3);% l2 = (xx50 - x(1) .* (xx50 - x(3) / (x(2) - x(1) / (x(2) - x(3);% l3 = (xx50 - x(1) .* (xx50 - x(2) / (x(3) - x(1) / (x(3) - x(2);% yy250 = y(1) * l1 + y(2) * l2 + y(3) * l3figurehold onplot(xx * 180 / pi, sin(xx), r)

22、;plot(xx * 180 / pi, yy1ex, g);plot(xx * 180 / pi, yy1in, black);% plot(xx * 180 / pi, yy2, black);plot(x * 180 / pi, y, b)303540455055600.50.550.60.650.70.750.80.85 sin外推内插节点23是否插值的节点越多，是否插值的节点越多，多项式插值越精确？多项式插值越精确？是否多项式的阶数越高，是否多项式的阶数越高，多项式插值越精确？多项式插值越精确？2421( ),-5,5,2,4,8,10( )1f xxnf xnLagrangex设取

23、作的 次插值多项式25-5-4-3-2-1012345-1.5-1-0.500.511.52-5-4-3-2-1012345-1.5-1-0.500.511.52-5-4-3-2-1012345-1.5-1-0.500.511.52-5-4-3-2-1012345-1.5-1-0.500.511.52-5-4-3-2-1012345-1.5-1-0.500.511.52-5-4-3-2-1012345-1.5-1-0.500.511.52n=2n=4n=6n=8n=10f(x)=1/(1+x2)-5-4-3-2-1012345-1.5-1-0.500.511.52n=2n=4n=6n=8n=1

24、0f(x)=1/(1+x2)不同次数的不同次数的Lagrange插值多项式的比较图插值多项式的比较图Runge现象现象26-50500.20.40.60.81-505-0.4-0.200.20.40.60.811.2-505-1.5-1-0.500.511.5-505-0.500.511.522阶4阶8阶10阶注意： 上述插值结果是正确的多项式插值结果，而且是唯一的多项式插值结果。Runge现象是多项式插值本身的缺陷而非误差。27%lagrangen.mfunction y=lagrangen(x0,y0,x)n=length(x0);m=length(x);for i=1:m z=x(i);

25、s=0; for k=1:n L=1; for j=1:n if j=k L=L*(z-x0(j)/(x0(k)-x0(j); end end s=s+L*y0(k); end y(i)=s;endy; Lagrange插值多项式求插值的Matlab程序.28%Chazhibijiao.mx=-5:0.1:5;z=0*x;y=1./(1+x.2);plot(x,z,k,x,y,r)axis(-5 5 -1.5 2);pause,hold onfor n=2:2:10 x0=linspace(-5,5,n+1); y0=1./(1+x0.2); x=-5:0.1:5; y1=lagrangen(

26、x0,y0,x); plot(x,y1), pauseendy2=1./(1+x0.2);y=interp1(x0,y2,x);plot (x,y,k),hold offgtext(n=2),gtext(n=4),gtext(n=6)gtext(n=8),gtext(n=10)gtext(f(x)=1/(1+x2)比较不同的插值多项式次数对插值的影响比较不同的插值多项式次数对插值的影响29x=-5:5;y=1./(1+x.2);t=-5:.05:5;y0=1./(1+t.2);p=polyfit(x,y,10);y1=polyval(p,t);plot(t,y0,x,y,o,t,y1,.)-5

27、-4-3-2-1012345-0.500.511.52Range现象演示（利用现象演示（利用matlab函数）函数）30 在不少实际问题中，对插值不但要求在节点上函数值相等函数值相等而且还要求它的导数值也相等导数值也相等。01.,( ),(0,1,2,., )( )( )( ),(0,1,2,., )( )( )niiiiiiiiaxxxb yf xinH xHmfxHxf xinxm数学描述：设在节点要求插值多项式满足：10210( ) , 1. , ( )( ),( )( ) 21 ( ) ( )( )( )( )1 ( )1 2()niiiinniiiiiiikikf xC a bnxx

28、a bH xf xH xfxnHermiteHxf xxfxxxxxxx：构造定理：给定和个不同的节点 ，,则满足条件：的最小阶多项式为阶插值多项式其中，220( )( )() ( )niiiik ilxxxx lx310-110-11(-).( -)( -).( -)( ) ()0,()0,(),( )1()(-).(-)(-).(-)()0,()0()(iiniiiiiiiniijijiijiiiixxx xx xx xl xxxxxxxxxl xijxxxxijaxxbx、目证明：步骤一由于拉格朗日基函数 满足条件标：构造： 构造满足条件: ，22()( )( )21 0, ( )()(

29、 )1( )( )( )2() ( )012()0,()( )1( )0( )0-2 (ijiiijiiiiiiiiiiiiiiiiiiiiiixxijlxxnxaxb lxxl xal xaxb lxaxbaxalxlxx，则为次多项式，且满足条件：为满，足条件：则整理得，解，得：22)12( )( )(-2 ( )12( ) ( )1-2( -) ( )( )i iiiiii iiiiiiibxlxxlx xxlxlxx x lxlx ，320-110-11200(-).( -)( -).( -) ( )(-).(-)1-( ) (-).(-) ( )1 ( )1 2 ()( ) iini

30、iiiiiiniiniiikikknkikk iiixxx xx xx xl xxxxxxxxxlxxxxlxxxxxx：、目标：构造证明（续）：将求导，得 即：满足条件:步骤二22( )() ( )( )21 ( ()0,()0,(),( )0( )1()0()0 ()1( )( )1 ( =0 )ijijiiiiiiiiijijiiiiiiiixxc xxxijxxxxijxlxxnxxclx ，构造则为次多项式，且满足条件：为满足条件：即 ，： 2( )() ( )iixxx lx3322102 -121220021221012121 ( ). ( )2.0.1.1 .2. 10.2.

31、10nnnnnnnnnnnnnnH xa xa xaH xna xaaxxaxxnxanx设阶多项式为：则：构造矩阵可得：000.nnyyyay34*2121*2121*2121( )( ) ( )( )-( ) ()()-()0,(0,1,2,., ) ()()-()0,(0,1,2,., )nnknknkknknkHermiteH xHxHermitexHxHxxHxHxknxHxHxkn。0n唯一性定理：插值存在且唯一证：假设和均满足插值条件，于是 在节点x ,.,x b,有： ( )22( )21kxxnxn 在每个节点 上均有二重根，即有个根。 与是次多项式矛盾。思路二： 证明矩阵形

32、式中节点矩阵满秩。35(22)21101(22)21( )22( ), , , ( )( )-( ) ( )( -)( -).( )( )(22 !-)()nnnnnnHermiteyf xnfxa bxa bHermiteR xf xHxxx xx xx xfxn插值余项定理： 函数的 阶导数在（）存在，对任意的插值余项其中：Hermite插值也存在插值也存在Runge现象。现象。361231231231.3, 1.6, 1.90.6200, 0.4554, 0.2818-0.5220, -0.5699, -0.5811xxxyyyyyyHermite已知： 求插值公式。2120001022

33、021110120122021()()50175152100175( ) ( )()()99999()()100320247200320( ) ( )()()99999()()5( )()()xxxxlxxxlxxxxxxxxxxl xxxlxxxxxxxxxxlxxxxx步骤1、写出拉格朗日多项式和其导数：2220145104100145 ( )99999xxlxx3722200000222111112222222202( )( )50175152( )1 2() () ( )(1012)()999100320247( )1 2() () ( )(1)()99950145104( )1 2(

34、) () ( )10(2)()999 (HermiteH xH xHxxx lxlxxxxHxxx lxlxxxHxxx lxlxxxxH步骤 、写出多项式和：2220022211122222250175152) ( )(1.3)()999100320247 ( ) ( )(1.6)()99950145104 ( ) ( )(1.9)()999xxx lxxxxHxxx lxxxxHxxx lxxxx3850011220011223( )( )( )( ) ( ) ( ) ( )HermiteHxy Hxy Hxy HxyHxyHxyHx步骤 、写出多项式插值结果：1.31.41.51.61.

35、71.81.92-0.58-0.57-0.56-0.55-0.54-0.53-0.52 f(x)y1.31.41.51.61.71.81.920.250.30.350.40.450.50.550.60.65 f(x)yHermiteHermite插值结果，保留了插值结果，保留了f f(x)(x)的特征的特征391.31.41.51.61.71.81.920.250.30.350.40.450.50.550.60.65 f(x)y1.31.41.51.61.71.81.92-0.6-0.59-0.58-0.57-0.56-0.55-0.54-0.53-0.52 f(x)y多项式插值结果，没有保留

36、了多项式插值结果，没有保留了f f(x)(x)的特征的特征40clear allclose allclc t = 1.3, 1.6, 1.9;y1 = 0.6200, 0.4554, 0.2818;dy = -0.5220, -0.5699, -0.5812; A = zeros(6); for iii = 1 : 3 for jjj = 1 : 6 A(iii, jjj) = t(iii) (6 - jjj); endend for iii = 4 : 6 for jjj = 1 : 5 A(iii, jjj) = (6 - jjj) * t(iii - 3) (5 - jjj); ende

37、nd y = y1, dy p = A y.;x = 1.3 :0.001 : 1.9;yy = polyval(p,x); plot(x, yy)hold onplot(t, y1, o) dyy = diff(yy) * 1000;figureplot(x(1 : 600), dyy)hold onplot(t, dy, o) figureyy2 = polyval(polyfit(t, y1, 2), x)plot(x, yy2, r)hold onplot(t, y1, o) dyy2 = diff(yy2) * 1000;figureplot(x(1 : 600), dyy2)hol

38、d onplot(t, dy, o)415432543210432543211152433311222133035( )( )5432()()() ()()()hermitH xa xa xa xa xa xaH xa xa xa xa xaH xyaH xyaH xyaH xyaH xyaH xya直接通过求解线性方程组计算。考虑 个插值节点，则插值为 阶：建立方程组：2.533.544.555.56-1-0.500.51X: 5.655Y: -0.5326 hermit节点二阶多项式sin（x）sin函数的三阶多项式插值和Hermite插值42clear allclose allclc t

39、 = 1 : 3;t = t / 3 * pi / 0.5; y1 = sin(t)dy = cos(t) A = zeros(6); for iii = 1 : 3 for jjj = 1 : 6 A(iii, jjj) = t(iii) (6 - jjj); endend for iii = 4 : 6 for jjj = 1 : 5 A(iii, jjj) = (6 - jjj) * t(iii - 3) (5 - jjj); endend y = y1, dy p = A y.x = 1 : 30;x = x / 30 * pi / 0.5;yy = polyval(p,x) plot(x, yy)hold onplot(t, y1, o) yy2 = polyval(polyfit(t, y1, 2), x)plot(x, yy2, r) plo