数值分析第五章

1、拉格朗日插值公式拉格朗日插值公式均差与牛顿插值公式均差与牛顿插值公式埃尔米特插值函数埃尔米特插值函数样条插值方法样条插值方法数值分析 Ch5引例引例.函数函数 sin x 的逼近的逼近(1) 线性函数逼近线性函数逼近 y0 = x(2)泰勒级数逼近泰勒级数逼近 y1(x)= x x3/3! + x5/5! (3)抛物线逼近抛物线逼近(error=0.0559) y2=4x( x)/2(4)帕特逼近帕特逼近(error=0.0036)0123400.511.50123400.511.50123400.51)536411088(1555122260166320)(4253xxxxxxP 00.20

2、.40.60.811.21.41.61.82-0.4-0.200.20.40.6(1)复杂函数的计算复杂函数的计算;(2)函数表中非表格点计算函数表中非表格点计算(3)光滑曲线的绘制光滑曲线的绘制;(4)提高照片分辩率算法提高照片分辩率算法(5)定积分的离散化处理定积分的离散化处理;(6)微分方程的离散化处理微分方程的离散化处理;(7)积分方程的离散化处理积分方程的离散化处理;插值法的应用背景插值法的应用背景:510155101524682468已知已知f(x)在点在点xi上的函数值上的函数值 yi=f(xi), (i=0,1,2,n)则称则称 P(x) 为为 f(x) 的的插值函数插值函数.

3、 称称 x0, x1, , xn为为 插值结点插值结点; 称称 f(x) 为为被插值函数被插值函数. 如果如果 P(x)=a0 + a1x + anxn满足满足: P(xk)= yk (k = 0,1,n)设设 f(x)C a , b, 取点取点 a x0 x1xnb代数插值问题代数插值问题插值条件插值条件点点,则满足插值条件则满足插值条件 P(xk)= yk (k = 0,1,n)的的 n 次插值多项式次插值多项式 P(x)=a0 + a1x + anxn存在而且是唯一的存在而且是唯一的。 nnnnnnnnnyxaxaayxaxaayxaxaa101111000010证明证明: 由插值条件由

4、插值条件P(x0)= y0P(x1)=y1P(xn)=yn定理定理5.1 若插值结点若插值结点x0, x1, xn 是是 (n+1)个互异个互异nnnnnxxxxxxA111|1100 方程组系数矩阵取行列式方程组系数矩阵取行列式故方程组有唯一解故方程组有唯一解.从而插值多项式从而插值多项式 P(x) 存在而且是唯一的存在而且是唯一的.例例5.1 误差函数表可构造误差函数表可构造6次插值函数次插值函数0)(0 jijinxxx 0 0.5000 1.0000 1.5000 2.0000 2.5000 3.0000y 0 0.5205 0.8427 0.9661 0.9953 0.9996 1.

5、0000)()(001010 xxxxyyyxL 过两点直线方程过两点直线方程求满足求满足: L(x0)=y0 , L(x1)=y1的线性函数的线性函数 L(x)已知函数表已知函数表 x x0 x1 f(x) y0 y1例例 求求 的近似值的近似值(函数值函数值: 10.7238)11510.7143)100115(100121101110115 01010110)(,)(xxxxxlxxxxxl 记记当当 x0 x x1 时时 0l0(x)1, 0l1(x)1x x0 x1l0(x) 1 0l1(x) 0 11100)()()(yxlyxlxL y0 y1 = 1 0y0 + 0 1y100

6、111010)(yxxxxyxxxxxL 对称形式对称形式二次插值问题二次插值问题 x x0 x1 x2 f(x) y0 y1 y2已知函数表已知函数表求函数求函数 L(x)=a0 + a1x + a2 x2 满足满足: L(x0)=y0 , L(x1)=y1, L(x2)=y2y0 y1 y2 = 1 0 0y0 + 0 1 0y1+ 0 0 1y2L(x)=l0(x)y0+l1(x)y1+l2(x)y2，)()()(2010210 xxxxxxxxxl 二次插值函数二次插值函数: L(x)=l0(x)y0+l1(x)y1+l2(x)y2，xx0 x1x2l0(x) 1 0 0l0(x) 1

7、00l1(x) 010l2(x) 0 0 1L(x) y0y1y2 xx0 x1 x2)()()(2101201xxxxxxxxxl )()()(1202102xxxxxxxxxl 二次插值基函数图形二次插值基函数图形00.51-0.500.5100.51-0.500.5100.5100.51取取 x0 =0, x1 = 0.5, x2 = 100.51-0.500.51l0(x)=2(x 0.5)(x 1);l1(x)= 4 x(x 1);l2(x) = 2(x 0.5)x拉格朗日插值公式拉格朗日插值公式插值条件插值条件:L(xk)= yk (k = 0,1,n)nnnyxlyxlyxlxL

8、)()()()(1100 )()()()()()()(110110nkkkkkknkkkxxxxxxxxxxxxxxxxxl 其中其中,第第k (k=0,1,，n)个插值基函数个插值基函数nkjjjkjkxxxxxl 0)()()(或或:Runge反例反例: , (-5x5)211)(xxf -5-4-3-2-1012345-0.500.511.52L10(t) f(t) f(x)取取xk= 5+k 计算计算: f(xk) (k=0,1,10) 构造构造L10(x).取取:tk= 5+0.05k (k=0,1,200),计算计算: L10(tk)x=-5:5;y=1./(1+x.2);t=-5

9、:0.05:5;y1=1./(1+t.2);n=length(t);for i=1:n z=t(i);s=0; for k=1:11 Lk=1;u=x(k); for j=1:11 if j=k,Lk=Lk*(z-x(j)/(u-x(j); end end s=s+Lk*y(k); end y2(i)=s;endplot(x,y,ko,t,y1,t,y2,r)00111010)(yxxxxyxxxxxL 两点线性插值两点线性插值插值余项插值余项(误差误差): R(x) = f(x) L(x) 由插值条件由插值条件,知知 R(x)=C(x) (x x0)(x x1)即即 f(x) L(x) =

10、C(x) (x x0)(x x1) C(x) = ? ax0 x1xnb则对任何则对任何xa , b, 满足满足 Ln(xk) = f(xk) 的的 n 次次插值多项式插值多项式Ln(x) 的误差的误差)()!1()()()()(1)1(xnfxLxfxRnnnnn )()()(101nnxxxxxxx 其中其中,),(ban 且与且与x有关有关定理定理5.2 设设 f(x)Ca, b, 且且 f (x) 在在(a, b)内具有内具有n+1阶导数阶导数, 取插值结点取插值结点证明证明: 记记 n+1(x) =(x x0)(x x1)(x xn)f(x) Ln(x)= C(x) n+1(x)取定

11、取定 x(a, b), 设设 t( a, b ). 构造函数构造函数 )()()()(1tCtLtftFnn 显然显然, F(x) = 0, F(xj) = 0, ( j = 0,1,n ) 由插值条件由插值条件Ln(xk) = f(xk) (k = 0,1,n)知存在知存在C(x)使得使得 F(t) 有有(n+2)个相异零点个相异零点. 根据根据Rolle定理定理, F(t)在区间在区间(a, b)内至少有内至少有 (n +1)个相异零点个相异零点. 0)!1()()1( nCfn )!1()()1( nfCn )()!1()()()(1)1(xnfxLxfnnn 依此类推依此类推,F (

12、n+ 1 )(t) 在区间在区间 ( a, b ) 内至少有内至少有一个零点一个零点。故存在故存在 (a, b), 使使F(n+1)( )=0)()()()()1(1)1()1()1(tCtLtftFnnnnnn 3232752)(xxxxP 例例 分析三次多项式插值误差分析三次多项式插值误差 x 0 1 2 3P3(x) -2 -2 -4 40)()4(3 xP)3)(2)(1()(! 41)()4(33 xxxxPxR 0)(3 xR32332752)()(xxxxPxL 例例 取被插值函数为正弦函数取被插值函数为正弦函数 f(x) = sin x，取三点做取三点做二次插值。二次插值。 2

13、/ x 0 Sin x 0 1 06/ | )(2/(|cos|2 xxxR22/ )(4)( xxxL 486142|322 R思考题思考题: : 在区间在区间 内增加插值结点是否会内增加插值结点是否会导致导致RungeRunge现象现象, 0 例例5.3 设设 y = f(x) 在区间在区间 a, b上有连续上有连续,且且 f (x) 在在 (a, b)内具有内具有2阶导数阶导数,已知已知f (x)在区间端点处在区间端点处的值的值.如果当如果当x (a, b)时时,有有|f (x)|M. 试证明试证明21)(8| )(|abMxR 证明证明 由由Lagrange插值误差定理插值误差定理)(

14、2)()()()(11bxaxfxLxfxR 令令h(x) = |( x a )( x b )| 4)()2()(max2abbahxhbxa 21)(8| )(|abMxR 应用应用: 考虑制做考虑制做 sin x 在在0, , 上等距结点的函数表上等距结点的函数表, ,要要求用线性插值计算非表格点数据时求用线性插值计算非表格点数据时, ,能准确到小数后两能准确到小数后两位位, ,问函数表中自变量数据的步长问函数表中自变量数据的步长h应取多少为好应取多少为好？ 解解：设应取的步长为设应取的步长为h , 则则 xj = jh ( j = 0,1,n). 当当 x(xj , xj+1)时时88)

15、(| )(|max| )(|2211hxxxfxRjjxxxjj sin)(sin)(1sin11jjjjxxxxxxhx 2210218 h h 0.2 只须只须取取x0, x1, x2,求二次函数求二次函数 P(x)=a0 + a1(x x0) + a2 (x x0)(x x1)满足条件满足条件 P(x0)=f(x0), P(x1)=f(x1), P(x2)=f(x2) )()()()()()(21202202101011000 xfxxxxaxxaaxfxxaaxfa 插值条件引出关于插值条件引出关于a0, a1, a2方程方程牛顿插值问题牛顿插值问题010110)()(,xxxfxfx

16、xf 121221)()(,xxxfxfxxf 021021210,xxxxfxxfxxxf 解下三角方程组过程中引入符号解下三角方程组过程中引入符号a0 = f(x0), a1 = fx0, x1, a2 = fx0, x1, x2P(x)= f(x0) + fx1, x2(x x0) + fx0, x1, x2(x x0)(x x1)牛顿插值公式牛顿插值公式:定义定义5.3 若已知函数若已知函数 f(x) 在点在点 x0,x1,xn 处的处的值值 f(x0), f(x1), , f(xn).如果如果 i j ,则则( j = 0,1,n-1 ) 一阶均差一阶均差n阶均差阶均差010110,

17、xxxxfxxfxxxfnnnn jjjjjjxxxfxfxxf 111)()(,二阶均差二阶均差jjjjjjjjjxxxxfxxfxxxf 212121,( j = 0,1,n-2 ) x- 2-1013y-56-16-2-24例例 由函数表由函数表求各阶均差求各阶均差xf(x) 一阶均差一阶均差二阶均差二阶均差三阶均差三阶均差-2-56-1-16 400-2 14 -131-2 0 -7 23 4 3 1 2解解:按公式计算一阶均差按公式计算一阶均差、二阶均差二阶均差、三阶均差三阶均差MATLAB程序计算程序计算x=-2 -1 0 1 3;y=-56 -16 -2 -2 4;data=x,

18、y;dy=y;n=length(x);for k=1:n-1 dx=x(k+1:n)-x(1:n-k) dy=diff(dy)./dx; f=zeros(n,1);f(k+1:n)=dy; data=data,f;enddata -2 -56 -1 -16 40 0 -2 14 -13 1 -2 0 -7 2 3 4 3 1 2 0函数值的计算函数值的计算:N3(x) = 56 + (x + 2) 40 +(x + 1) 13 +2 xN3(x) = 56 + 40(x + 2) 13(x + 2)(x + 1) + 2(x + 2)(x + 1) x 22680560 a5439401 a7

19、6132 a23 a3232752)(xxxxP -2-10123-60-50-40-30-20-10010根据代数插值存在唯一性定理根据代数插值存在唯一性定理, n 次牛顿插值公式恒次牛顿插值公式恒等于等于n次拉格朗日插值公式次拉格朗日插值公式,误差余项也相等，即误差余项也相等，即)()!1()()(1)1(xnfxRnnn 牛顿多项式算法牛顿多项式算法: 记插值节点为记插值节点为x0,x1,xn, f(x)的各阶差商为的各阶差商为 f0, f1, f2, , fn sfn 计算计算 s fk+s*(x-xk) (k=n-1,n-2, , 0)(1)(3) N(x)=s)()()(101nn

20、xxxxxxx 取插值结点取插值结点: ax0 x1xnb满足满足Ln(xk)=f(xk)的的 n 次多项式插值余项次多项式插值余项)()!1()()()()(1)1(xnfxLxfxRnnnnn )()()(101nnxxxxxxx 其中其中,选取选取: x0, x1 , xn , 使使min| )(|max1 xnbxa 结论结论: 切比雪夫多项式切比雪夫多项式Tn+1(x)的全部零点。的全部零点。拉格朗日插值余项拉格朗日插值余项 n+1阶切比雪夫多项式阶切比雪夫多项式: Tn+1=cos(n+1) cos = x 代入得代入得 Tn+1( x ) = cos(n+1) arccos x

21、)1(2)12(cos( nkxk 即即2)12(arccos)1( kxn( k=0,1,n )取取 f(x)C1, 1, 令令 x = cos , 则有则有 1, 1 0, 将将g( ) = f(cos )展开成余弦级数展开成余弦级数 10cos21)(nnnaag 切比雪夫结点切比雪夫结点211)(xxf 例例1. 函数函数取等距插值结点取等距插值结点: -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5 x-5, 5 11(x)=(x+5)(x+4)(x+3)(x+2)(x+1)x(x-1)(x-2)(x-3)(x-4)(x-5)(!11)()()(11)11(1

22、0 xfxLxfn 11(x) -4.9491 -4.5482 -3.7787 -2.7032 -1.4087 0.0000 1.4087 2.7032 3.7787 4.5482 4.9491在在-5, 5区间上区间上,取取11个切比雪夫结点个切比雪夫结点)22)12(cos(5 kxk( k=10, 9, 8, , 1, 0 ) 11(x)=(x x0)(x x1)(x x2)(x x10) 11(x) 00)(yxf 11)(yxf 00)(mxf 11)(mxf 插值函数插值函数 H(x)= a0 + a1x + a2x2 + a3x3插值条件插值条件三次三次Hermite插值问题插值

23、问题xx0 x1H(x)y0y1H(x)m0m1设设,)(104xxCxf 已知已知-0.4-0.200.20.40.60.811.21.41.6-0.6-0.4-0.200.20.40.60.81例例2. 插值条件插值条件:求求3次插值函数次插值函数.)1()(1xxCxH 解解: 设设积分积分,得得利用函数插值条件利用函数插值条件, 得得 C2 = 0 , C1 = 6 H(x) = 3x2 2x3 = (3 2x)x2 x 0 1 H(x) 0 1 H(x) 0 02321)23(61)(CxxCxH 利用基函数表示利用基函数表示Hermite插值插值)()()()()(11001100

24、 xmxmxyxyxH x0 x1 10 00 01 0 0)(0 x )(1x )(0 x )(1x x)(1x )(0 x x0 x1 00 10 0 0 0 1)(1x )(0 x x令令010 xxxx 则则0111xxxx 20110100)(21()(xxxxxxxxx 20100111)(21()(xxxxxxxxx )1()(1 CH21)23()( 1)23()(20 设设2216/ )23()(CCH 则则2)1(21 22)1(2)1( )21)(1(2 )21()1(2 )(1 )(0 0 1 10 00 0 1 0 0)(1 )(0 201011)()(xxxxxxx

25、 210)()(xxxxCxH 设设1)()(2010 xxCxH201100)()(xxxxxxx )(1x )(0 x x0 x1 00 10 0 0 0 1)(1x )(0 x x1)()(2011 xxCxH)()()(120 xxxxCxH 设设210)4(3)(! 4)()()()(xxxxfxHxfxR 两点两点Hermite插值的误差估计式插值的误差估计式证明证明: 由插值条件知由插值条件知构造辅助函数构造辅助函数2120)()()()()(xtxtxCtHtftF 利用利用 f(x) H(x)=C(x)(x x0)2(x x1)22120)()()(xxxxxCxR 取取 x

26、 异于异于 x0 和和 x1, 有有0)()(00 xRxR0)()(11 xRxR2120)()()()()(xtxtxCtHtftF 0) ! 4)()()()4()4( xCfF ! 4)()()4( fxC 210)4()(! 4)()(xxxxfxR 显然显然,F(t) 有三个零点有三个零点 x0, x, x1, 由由Roll定理知定理知,存在存在 两个零点两个零点 t0, t1. )(tF )(tF 故故 有四个相异零点有四个相异零点1100 xtxtx 反复应用反复应用 Roll 定理定理, 得得 F(4)(t) 知一个零点设为知一个零点设为 分段线性插值分段线性插值插值节点满足

27、插值节点满足: x0 x1xn 已知已知yj=f (xj) ( j= 0,1,2,n)1111)( jjjjjjjjhyxxxxyxxxxxL( j= 0,1,n-1)xxj，xj+1时时, 线性插值函数线性插值函数分段三次分段三次Hermite插值插值 jjjjjjjhyxxxxxxxxxH2111)(21()( ( j= 0,1,2,n-1)已知函数值和导数值已知函数值和导数值 )(),(jjjjxfmxfy xxj，xj+112111)(21( jjjjjjjyxxxxxxxxjjjjjmxxxxxx211)( 1211)( jjjjjmxxxxxx( j= 0,1,2,n)一元函数线性

28、插值基一元函数线性插值基01 2/ )1(1 N2/ )1(2 N二元函数双线性插值基二元函数双线性插值基1 4/ )1)(1(1 4/ )1)(1(2 4/ )1)(1(3 4/ )1)(1(4 E )1, 1(1 P)1 , 1(3P)1 , 1(4 P)1, 1(2 P1 , 1 N1 1, 0N2 0, 1 01000110)()(111 NN )()(223 NN )()(214 NN )()(111 NN )()(122 NN 10001010 00101001 00010101 例例1. 飞机机翼剖面图飞机机翼剖面图 1.1.数据采集数据采集2. 数据插值处理数据插值处理X 0

29、-0.4552 -0.6913 -0.8640 -0.9689 -0.9996Y 0 0.3285 0.3467 0.2716 0.1408 -0.0160T=1:6;t=1:.2:6;x=spline(T,X,t);y=spline(T,Y,t);例例2:龙格函数的插值逼近龙格函数的插值逼近 -505-0.500.51211)(xxf -505-0.500.517结点等距插值结点等距插值7结点切比雪夫插值结点切比雪夫插值7结点样条插值结点样条插值7结点埃尔米特插值结点埃尔米特插值-50500.51-50500.51利用龙格函数的数据表做样条插值第一步利用龙格函数的数据表做样条插值第一步 x-

30、5.0 -3.33 -1.66 0 1.66 3.33 5.0 y0.038 0.082 0.264 1.0 0.264 0.082 0.038m 0.014 -0.0054 0.4142 0. -0.4142 0.0054-0.014 y 0.0140.045 0.233 0 -0.233 -0.045 -0.014估算结点处导数值估算结点处导数值 mk, 由三对角方程组求解得出由三对角方程组求解得出 54321543214114114114114FFFFFmmmmmhyyFkkk)(311 (k=1，2，3，4，5，6)， h=10/6 定义定义 5.4 给定区间给定区间a , b上的一个

31、分划上的一个分划: a = x0 x1 xn = b已知已知 f(xj) = yj (j = 0,1,n), 如果如果 ,),(,),(,),()(1212101nnnxxxxSxxxxSxxxxSxS满足满足: (1) S(x)在在 xj，xj+ 1上为三次多项式上为三次多项式; (2) S”(x)在区间在区间a，b上连续上连续; (3) S(xj) = yj ( j = 0,1,n). 则称则称 S(x)为三次样条插值函数为三次样条插值函数.当当xxj , xj+ 1 ( j= 0,1,n-1 )时时 Sj(x)= aj + bj x + cj x2 + dj x3插值条件插值条件: S(

32、xj) = yj ( j = 0,1,n)连续性条件连续性条件: S(xj+0) =S(xj-0) ( j = 1,n-1) S(xj+0) =S(xj-0) ( j = 1,n-1) S”(xj+0) =S”(xj-0) ( j = 1,n-1)由样条定义由样条定义,可建立方程可建立方程(4n-2)个！个！ n个三次多项式个三次多项式, 待定系数共待定系数共4n个个！方程数少于未知数个数方程数少于未知数个数 ？(1)自然边界条件自然边界条件: S”(x0)=0, S”(xn)=0例例 5.7 已知已知f(1) = 1, f(0) = 0, f(1) = 1.求求1，1 上的三次自然样条上的三

33、次自然样条(满足自然边界条件满足自然边界条件). 解解 设设 ,)(1001222232112131xdxcxbxaxdxcxbxaxS则有则有: a1+b1c1+d1=1,d1=0, a2+b2+c2+d2=1 d1=d2, c1=c2, b1=b2 (2)周期边界条件周期边界条件: S(x0)=S(xn), S”(x0)=S”(xn)(3)固定边界条件固定边界条件: S(x0)=f (x0), S(xn)=f (xn)由自然边界条件由自然边界条件: 6a1+2b1=0, 6a2+2b2=0 解方程组解方程组,得得 a1=-a2=1/2, b1=b2=3/2, c1=c2=d1=d2=0,)(1023210123212323xxxxxxxS问题的解问题的解 -1-0.8-0.6-0.4-0.200.20.40.60.8100.20.40.60.