1/602022-11-7趣例趣例1:图像放大图像放大2/602022-11-73/602022-11-7Non-damagedDamaged趣例趣例2:图像修复图像修复4/602022-11-7趣例趣例3:数据可视化数据可视化http:/ Paul de Casteljau 和和Pierre Bzier,随后美国通用汽随后美国通用汽车的其它人一起推动了现在称为三次样条和车的其它人一起推动了现在称为三次样条和Bzier 样条的建立。样条是通过样条的建立。样条是通过很少很少的控制点就能够生成复杂平滑曲线的控制点就能够生成复杂平滑曲线的方法。的方法。6/602022-11-7趣例趣例5:游戏与电影游戏与电影Ref:http:/ 如果一个函数如果一个函数P(x)满足满足P(xi)=yi (i=0,n),那么那么函数函数P(x)插值了一系列数据点插值了一系列数据点(x0,y0),(xn,yn)，其中其中P(x)称为称为插值函数插值函数,点点x0,xn称为称为插值节点插值节点。x0 x1x2x3x4xP(x)函数是描述自然界客观规律的重要工具。函数是描述自然界客观规律的重要工具。8/602022-11-7 选择多项式函数的理由选择多项式函数的理由:计算方面计算方面多项式函数是计算机最基本的函多项式函数是计算机最基本的函数数,计算多项式函数的值只需用加和乘运算计算多项式函数的值只需用加和乘运算,且且积分和微分均非常方便。积分和微分均非常方便。理论方面理论方面多项式函数简单明了的数学性质。多项式函数简单明了的数学性质。有一个简单的原理可以说明什么时候存在给定有一个简单的原理可以说明什么时候存在给定次数的插值多项式。次数的插值多项式。插值函数类的选择插值函数类的选择:9/602022-11-7100010()()yyP xyxxxx 过两点直线方程过两点直线方程已知函数表求满足已知函数表求满足:P(x0)=y0和和P(x1)=y1的线性函数的线性函数 P(x)x x0 x1 y y0 y1引例引例 求求 的近似值的近似值11510.7143)100115(100121101110115 真实值真实值:10.723810线性插值函数线性插值函数x0 x1(x0，y0)(x1，y1)P(x)可见可见 是过是过 和和 两点的直线。两点的直线。11抛物插值函数抛物插值函数x0 x1x2因过三点的二次曲线为抛物线因过三点的二次曲线为抛物线,故称为抛物插值。故称为抛物插值。12/602022-11-7则称则称 P(x)为为 插值多项式插值多项式,称称 x0,x1,xn为为 插插值节点。值节点。如果如果 P(x)=a0+a1x+anxn满足满足 P(xk)=yk (k=0,1,n)考虑区间考虑区间a,b上上(n+1)个点个点a x0 x1xnb。插值条件插值条件 nnnnnnnnnyxaxaayxaxaayxaxaa101111000010 由插值条件由插值条件P(x0)=y0P(x1)=y1P(xn)=yn13/602022-11-7范德蒙范德蒙(Vandermonde)矩阵矩阵001111 1nnnnnxxxxAxx 14/602022-11-7100010()()yyP xyxxxx 过两点直线方程过两点直线方程已知函数表求满足已知函数表求满足:P(x0)=y0 和和 P(x1)=y1的线性函数的线性函数 P(x)。x x0 x1 y y0 y101010110()xxxxP xyyxxxx 15/602022-11-7The voyage of discovery is not in seeking new landscapes but in having new eyes.-Marcel Proust16/602022-11-710010101,()()xxxxxxxlxxlx 记记x x0 x1l0(x)1 0l1(x)0 10011()()()P xlx ylx y 01010110()xxxxP xyyxxxx 00111010()0,()011P xyyP xyyyy 17/602022-11-7I=imread(yao.png);J=imread(li.png);for alpha=1:-0.01:0 K=alpha*I+(1-alpha)*J;pause(0.3),imshow(K,)end0011()()()P xlx ylx y 01101100011()1()1,xxxxxxlxxxxxlxxx 01()(1)P xyy 18/602022-11-7二次插值问题二次插值问题 x x0 x1 x2 y y0 y1 y2已知函数表已知函数表求函数求函数 P(x)=a0+a1x+a2 x2 满足满足:P(x0)=y0,P(x1)=y1,P(x2)=y2P(x)=l0(x)y0+l1(x)y1+l2(x)y2l0(x)100l1(x)010l2(x)0 01 xx0 x1 x2001210121012012()00111()00()00P xyyyP xyyyP xyyyyyy 19/602022-11-7)()()(2010210 xxxxxxxxxl 二次插值函数二次插值函数:P(x)=l0(x)y0+l1(x)y1+l2(x)y2l0(x)100l1(x)010l2(x)0 01 xx0 x1 x2)()()(2101201xxxxxxxxxl )()()(1202102xxxxxxxxxl 20/602022-11-7拉格朗日方法拉格朗日方法插值条件插值条件:P(xk)=yk (k=0,1,n)0011()()()()nnP xlx ylx ylx y )()()()()()()(110110nkkkkkknkkkxxxxxxxxxxxxxxxxxl 其中第其中第k 个个插值基函数插值基函数 0()()()njkjkjjkxxlxxx 或或0()()()njkkjkjj kxxP xyxx 21/602022-11-7例例1求插值于点求插值于点(-2,-56),(-1,-16),(0,-2),(1,-2),(3,4)的次数小于等于的次数小于等于4的拉格朗日插值多项式。的拉格朗日插值多项式。(1)(0)(1)(3)(2)(0)(1)(3)()5616(21)(20)(21)(23)(12)(10)(11)(13)(2)(1)(1)(3)(2)(1)(0)(3)22(02)(01)(01)(03)(12)(11)(10)(13)xxxxxxxxP xxxxxxxxx 23(2)(1)(0)(1)4(32)(31)(30)(31)=2572xxxxxxx 22/602022-11-7程序片段程序片段1:Matlab Code:拉格朗日插值多项式拉格朗日插值多项式function v=polyinterp(x,y,u)%POLYINTERP Polynomial interpolation.%v=POLYINTERP(x,y,u)computes v(j)=P(u(j)where P is the%polynomial of degree d=length(x)-1 with P(x(i)=y(i).%Use Lagrangian representation.%Evaluate at all elements of u simultaneously.n=length(x);v=zeros(size(u);for k=1:n w=ones(size(u);for j=1:k-1 k+1:n w=(u-x(j)./(x(k)-x(j).*w;end v=v+w*y(k);end0()()()njkkkjjj kxxP xyxx 23/602022-11-7Demo2x=0:3;y=-5-6-1 16;u=-.25:.01:3.25;v=polyinterp(x,y,u);plot(x,y,o,u,v,-)syms u;L=polyinterp(x,y,u);L=simplify(L);24/602022-11-7则满足插值条件则满足插值条件 P(xk)=yk(k=0,1,n)的的 次数小次数小于等于于等于n次次的插值多项式的插值多项式 P(x)=a0+a1x+anxn存在而且唯一存在而且唯一。nnnnnnnnnyxaxaayxaxaayxaxaa101111000010证明证明:由插值条件由插值条件P(x0)=y0P(x1)=y1P(xn)=yn定理定理 若插值结点若插值结点x0,x1,xn 是是(n+1)个互异点个互异点,25/602022-11-700111 1det()1nnnnnxxxxAxx 回顾回顾1:非齐次方程组有唯一解的充分必要条件是系数非齐次方程组有唯一解的充分必要条件是系数矩阵行列式不等于零。矩阵行列式不等于零。系数矩阵行列式不等于零系数矩阵行列式不等于零,则方程组有唯一解。因此则方程组有唯一解。因此插值多项式插值多项式 P(x)存在且唯一。存在且唯一。0()0ijij nxx 回顾回顾2:范德蒙范德蒙(Vandermonde)矩阵矩阵0011102002131101211 1det()()()()()()()()nnnnnijnnj i nnnxxxxAxxAxxxxxxxxxxxxxxxx 则则()=()=21()()nnnnxxxx26/602022-11-7The voyage of discovery is not in seeking new landscapes but in having new eyes.-Marcel Proust27/602022-11-7给定给定x0,x1和和 x2,求二次函数求二次函数 P(x)=a0+a1(x x0)+a2(x x0)(x x1)满足条件满足条件 P(x0)=y0,P(x1)=y1,P(x2)=y2 00011010120220212()()()()ayaaxxyaaxxaxxxxy 满足插值条件的关于满足插值条件的关于a0,a1和和a2方程方程牛顿差商方法牛顿差商方法28/602022-11-7110110yyyxx 120220yyyxx 221221yyyxx 解下三角方程组过程中引入符号解下三角方程组过程中引入符号牛顿插值多项式牛顿插值多项式:000yy 012001122,ay ay ay 001201210()()()()yyP xxxxxyxx 29/602022-11-7定义定义 若已知函数若已知函数 f(x)在点在点 x0,x1,xn 处的值处的值 y0,y1,yn。如果如果 i j,则各阶差商则各阶差商(divided difference)定义如下定义如下(j=1,n)一阶差商一阶差商n阶差商阶差商010jjjyyyxx 二阶差商二阶差商121jjjyyyxx 三阶差商三阶差商(j=2,n)22232jjjyyyxx (j=3,n)1111nnnnnnnnyyyxx 30/602022-11-7更加一般地考虑更加一般地考虑000110101001 ()()()()nnnnnnayaa xxyaa xxaxxxxy 牛顿插值多项式牛顿插值多项式求解该方程组可得待定系数如下求解该方程组可得待定系数如下:001012010112()()()()+()()()nnnyyyyP xxxxxxxxxxxxx 001010121()()()()+()()()nnP xxxxxxxxxxxaaaxax 012001122,nnnay ay ayay 31/602022-11-7例例2求插值于点求插值于点(-2,-56),(-1,-16),(0,-2),(1,-2),(3,4)的的次数小于等于次数小于等于4的牛顿插值多项式。的牛顿插值多项式。2 561 16 0 2 1 23 4 4027 18 12 13 11 -7 2 20()56(2)(2)(21)(2)(1)()4013P xxxxxxx =56+(2)(+(1)(401)2+3xxx 32/602022-11-7例例3求插值于点求插值于点(-2,-56),(-1,-16),(0,-2),(1,-2),(3,4)的次数小于等于的次数小于等于4的拉格朗日插值多项式。的拉格朗日插值多项式。(1)(0)(1)(3)(2)(0)(1)(3)()5616(21)(20)(21)(23)(12)(10)(11)(13)(2)(1)(1)(3)(2)(1)(0)(3)22(02)(01)(01)(03)(12)(11)(10)(13)xxxxxxxxP xxxxxxxxx (2)(1)(0)(1)4(32)(31)(30)(31)xxxx 33/602022-11-7例例3求插值于点求插值于点(-2,-56),(-1,-16),(0,-2),(1,-2),(3,4),(2,-4)的的次数小于等于次数小于等于5的牛顿插值多项式。的牛顿插值多项式。2 561 16 0 2 1 23 4 2 4 0 0 加入一个新的点到拉格朗日插值多项式所需要额加入一个新的点到拉格朗日插值多项式所需要额外工作与牛顿插值多项式进行比较是很有趣的。牛顿外工作与牛顿插值多项式进行比较是很有趣的。牛顿差商方法具有拉格朗日方法所缺少的差商方法具有拉格朗日方法所缺少的”实时更新实时更新”性质。性质。4027 18 121 3 13 11 7 9 2 2 2 034/602022-11-7例例4 1 12 53 144 305 55 413/2 29/327/2 5/2 17/619/62222123(1)(21)/6nn nn证明证明:1/3 1/35/2()1(1)(1)(42)(1)(2)(1/33)S nnnnnnn 035/602022-11-7例例5求插值于点求插值于点(0,2),(1,1),(3,-1)的的次数小于等次数小于等于于2的插值多项式的插值多项式(拉格朗日方法和牛顿方法拉格朗日方法和牛顿方法)。(1)(3)(0)(3)(0)(1)()211(01)(03)(10)(13)(30)(31)=2 xxxxxxP xx 问题问题2：是否有多个经过这三个数据点的：是否有多个经过这三个数据点的次数小于等于次数小于等于2次的多项式？次的多项式？问题问题1：拉格朗日方法和牛顿方法的多项：拉格朗日方法和牛顿方法的多项式是否相同？式是否相同？36/602022-11-7问题问题3：是否有多个经过这三个数据点的次：是否有多个经过这三个数据点的次数大于数大于2次的多项式？次的多项式？1()=2P xx 加入第四个点加入第四个点(2,0),得到的多项式是得到的多项式是加入第四个点加入第四个点(2,3),得到的多项式是得到的多项式是322()=()(0)(1)(3)P xP xxxx 2(1)(3)(2)(0)(3)(2)()21(01)(03)(02)(10)(13)(12)(0)(1)(2)(0)(1)(3)13(30)(31)(32)(20)(21)(23)xxxxxxP xxxxxxx 37/602022-11-7压缩的概念压缩的概念:观测的离散数据可以想象成现实中无穷多观测的离散数据可以想象成现实中无穷多信息的代表。通过给定数据求出插值函数意味信息的代表。通过给定数据求出插值函数意味着用简单的规则代替无穷多信息。尽管期待这着用简单的规则代替无穷多信息。尽管期待这种简单规则精确地反映实际情况是不现实的种简单规则精确地反映实际情况是不现实的,但是它但是它可以充分接近实际可以充分接近实际。这一类压缩是这一类压缩是有损的压缩有损的压缩,即它会产生误差。即它会产生误差。用简单规则代替无穷多信息时会产生多大的误用简单规则代替无穷多信息时会产生多大的误差差,这是我们下面研究的内容。这是我们下面研究的内容。38/602022-11-701101010()xxxxP xyyxxxx 两点线性插值两点线性插值插值误差余项插值误差余项:R(x)=f(x)P(x)R(x)=?39/602022-11-7Rolle引理引理(),(,),()(),(,)()=0f xa ba bf af ba bf 设设函函数数在在闭闭区区间间连连续续 在在开开区区间间可可微微且且则则至至少少存存在在一一点点 属属于于使使得得。40/602022-11-7P(x)是满足是满足P(xk)=f(xk)的的n次插值多项式次插值多项式,则对则对任何任何xa,b,在在(a,b)内存在一点内存在一点 使得使得(1)1()()()()1)()!nnnRxf xPxxfxn 余余项项101()()()()nnxxxxxxx 。其中其中定理定理 设设 f(x)在在a,b连续且在连续且在(a,b)具有具有n+1阶阶导数导数,x0,x1,xn是是a,b内互不相同的节点内互不相同的节点。41/602022-11-7辅助函数构造辅助函数构造:(1)将欲证结论中将欲证结论中的改成的改成t;(2)将式子写成容易去掉导数符号的形式将式子写成容易去掉导数符号的形式(即容即容易积分的形式易积分的形式);(3)去掉导数符号去掉导数符号(即作积分即作积分),移项使等式一端移项使等式一端为零为零,另一端即为新作辅助函数另一端即为新作辅助函数F(t)。例如例如,拉格朗日中值定理和柯西中值定理。拉格朗日中值定理和柯西中值定理。42/602022-11-7证明证明:记记 n+1(x)=(x x0)(x xn)(0,),()(),(,)kkkxxknP xf xa b 首首先先注注意意到到如如果果则则在在内内任任意意选选择择 均均成成立立。(0,),kxxkna bt如如果果则则考考虑虑在在内内定定义义关关于于 的的函函数数如如下下。00()()()()()()()()()nntxtxg tf tP tf xP xxxxx00()()()()()()()()()()()=0()()00kkkkknkkkknt xxxxxxxg xf xP xf xP xxxxxxxf xP x对对=有有43/602022-11-700()()()()()()()()()=()()()()0nnt xxxxxg xf xP xf xP xxxxxf xP xf xP x对对=有有(1)()2Rolle,()(,)1,()(,)1ng tng ta bngta b 所所以以至至少少有有+个个互互异异的的零零点点。根根据据定定理理 在在内内至至少少存存在在+个个相相异异的的零零点点。以以此此类类推推在在内内至至少少存存在在 个个零零点点。(1)(1)0(1)0()(+1)!0()()()()()()()()()()()(+1)!nnnnnxngff xP xxxxxff xP xxxxxn 故故存存在在使使得得整整理理即即得得44/602022-11-7(1)0(1)10()()Taylor()=Taylor()(+1)!()()(+1)!)=nnnnR xfxxxxnfxnR xx 插插值值余余项项与与展展式式余余项项的的异异同同插插值值余余项项展展式式余余项项注释注释:45/602022-11-7Runge反例反例(rungeinterp)f(x)=1/(x2+1),(-5=x=5)x=-5:5;y=1./(x.2+1);u=-5:.01:5;v=polyinterp(x,y,u);plot(x,y,o,u,v,-)一般认为一般认为Pn(x)的次数的次数n越高逼近越高逼近f(x)的精度越好？的精度越好？46/602022-11-7例例 设设 y=f(x)在区间在区间 a,b上连续上连续,且且 f(x)在在(a,b)内具有内具有2阶导数。如果当阶导数。如果当x(a,b)时时|f(x)|M,则则21|()|()/8R xM ba证明证明 由由Lagrange插值误差公式插值误差公式11()()()()()()2fR xf xP xxaxb 令令h(x)=|(x a)(x b)|4)()2()(max2abbahxhbxa 21|()|()/8R xM ba47/602022-11-7分段线性插值分段线性插值插值节点满足插值节点满足:x0 x1xn 已知已知f(xk)=yk (k=0,1,2,n)11()()kkkkkkkkyyP xyxxysxx 先设定区间序号先设定区间序号k,使得使得 xk x xk+1再定义局部变量再定义局部变量s为为x-xk记一阶差商为记一阶差商为11kkkkkyyxx 插值函数可写为插值函数可写为:该式是通过该式是通过(xk,yk)和和(xk+1,yk+1)两点的线性函数。两点的线性函数。48/602022-11-7程序片段程序片段1:Matlab Code:分段线性插值分段线性插值function v=piecelin(x,y,u)%PIECELIN Piecewise linear interpolation.%v=piecelin(x,y,u)finds the piecewise linear L(x)%with L(x(j)=y(j)and returns v(k)=L(u(k).%First divided differencedelta=diff(y)./diff(x);%Find subinterval indices k so that x(k)=u x(k+1)n=length(x);k=ones(size(u);for j=2:n-1 k(x(j)=u)=j;end%Evaluate interpolants=u-x(k);v=y(k)+s.*delta(k);49/602022-11-7Demo x=-5:5;y=1./(x.2+1);u=-5:.01:5;v1=polyinterp(x,y,u);plot(x,y,o,u,v1,-)hold on,v2=piecelin(x,y,u);plot(u,v2,r-)分段线性插值保持了局部单调性分段线性插值保持了局部单调性(给定数据所描述的给定数据所描述的形状形状)。50/602022-11-7二维插值二维插值“If you dont know how to solve a problem,there must be a related but easier problem you know how to solve.See if youcan reduce the problem to the easier one.”G.Polya工程师的智慧工程师的智慧:分而治之分而治之divide and conquer 二维插值可以分解为两个一维插值。二维插值可以分解为两个一维插值。51/602022-11-752/602022-11-753/602022-11-7rowcolabcda ab ba ab bc cd dc cd dzero-orderfirst-order a (a+b)/2 b(a+c)/2 (a+b+c+d)/4 (b+d)/2 c (c+d)/2 d双线性插值双线性插值插值条件插值条件:P(x1,y1)=Q11,P(x2,y1)=Q21,P(x1,y2)=Q12,P(x2,y2)=Q22.21111211221(,)xxxxP x yQQxxxx 21212221221(,)xxxxP x yQQxxxx 21121221(,)(,)(,)yyyyP x yP x yP x yyyyy 56/602022-11-7 一维数据的插值一维数据的插值 interp1x=0:2:24;y=22,21,19,18,20,24,27,32,31,28,26,23,22;xi=0:0.1:24;yi=interp1(x,y,xi);plot(xi,yi,-,x,y,o);二维网格数据的插值二维网格数据的插值 interp2x,y=meshgrid(1:5,1:3);z=82 81 80 82 84;79 63 61 65 81;84 84 82 85 86;mesh(x,y,z)xi,yi=meshgrid(1:0.2:5,1:0.2:3);zi=interp2(x,y,z,xi,yi);mesh(xi,yi,zi);57/602022-11-7 散乱数据的插值散乱数据的插值(griddata和和ndgrid)rng(0,twister)x=rand(100,1)*4-2;y=rand(100,1)*4-2;z=x.*exp(-x.2-y.2);plot3(x,y,z,o),ti=-2:.25:2;xq,yq=meshgrid(ti,ti);zq=griddata(x,y,z,xq,yq);hold on,mesh(xq,yq,zq)58/602022-11-7x=double(imread(fl_orig.pgm);figure,subplot(131),imshow(x,);%help interp2(or imresize)%zero-order interpolation(replication)y0=interp2(x,nearest);subplot(132),imshow(y0,);%first-order interpolation(bilinear interpolation)y1=interp2(x,linear);subplot(133),imshow(y1,);59/602022-11-7%color imagex=double(imread(fl_orig.ppm);figure,subplot(131),imshow(x/255,);for i=1:3%zero-order interpolation(replication)y0(:,:,i)=interp2(x(:,:,i),nearest);%first-order interpolation(bilinear interpolation)y1(:,:,i)=interp2(x(:,:,i),linear);subplot(132),imshow(y0/255,);subplot(133),imshow(y1/255,);60/602022-11-7%How to implement image interpolation by yourself?figure,subplot(131),imshow(x,);%zero-order interpolation(replication)M,N=size(x);z0=zeros(2*M,2*N);z0(1:2:2*M,1:2:2*N)=x;h=1 1;1 1;z0=filter2(h,z0);subplot(132),imshow(z0,);%first-order interpolation(bilinear interpolation)z1=zeros(2*M,2*N);z1(1:2:2*M,1:2:2*N)=x;h=1 2 1;2 4 2;1 2 1/4;z1=filter2(h,z1);subplot(133),imshow(z1,);61/602022-11-7压缩的概念压缩的概念:观测的离散数据可以想象成现实中无穷多观测的离散数据可以想象成现实中无穷多信息的代表。通过给定数据求出插值函数意味信息的代表。通过给定数据求出插值函数意味着用简单的规则代替无穷多信息。尽管期待这着用简单的规则代替无穷多信息。尽管期待这种简单规则精确地反映实际情况是不现实的种简单规则精确地反映实际情况是不现实的,但是它但是它可以充分接近实际可以充分接近实际。这一类压缩是这一类压缩是有损的压缩有损的压缩,即它会产生误差。即它会产生误差。用简单规则代替无穷多信息时会产生多大的误用简单规则代替无穷多信息时会产生多大的误差差,这是我们下面研究的内容。这是我们下面研究的内容。62/602022-11-7 插值方法可以用来预测未来吗？插值方法可以用来预测未来吗？