Matlab MathCleve Morler著 陈文斌（wbchenfudan.edu.cn) 复旦大学2002多项式和样条在平面上给定n个点(xk,yk),可以唯一确定一个最多n- 1次的多项式通过这些点，这个多项式叫插值多项式插值多项式P(xk ) = yk , k = 1,2,nLagrange插值形式插值多项式例子x =0:3;y = -5 6 1 16;disp(x; y)Monomial基Vandermonde矩阵x =0:3; y = -5 6 1 16;V = vander(x)c = VyVandermonde矩阵是非 奇异的，但条件数是非 常坏的Polyinterp(Lagrange插值形式)function v = polyinterp(x,y,u)n = length(x);for k = 1:nw = ones(size(u);for j = 1:k-1 k+1:nw = (u x(j)./(x(k) x(j).*w;endv = v + w*y(k);end -0.500.511.522.533.5-10-50510152025u = -.25:.01:3.25; v = polyinterp(x,y,u);plot(x,y,o,u,v,-)Polyinterp(符号运算)symx =sym(x)P = polyinterp(x,y,symx)pretty(P)P = simplify(P)P = x3-2*x-5012345676810121416182022Polyinterp(另外的例子)x = 1:6; y = 16 18 21 17 15 12;u = .75:.05:6.25; v = polyinterp(x,y,u);plot(x,y,o,u,v,-);0123456710121416182022分片线性插值x = 1:6; y = 16 18 21 17 15 12;plot(x,y,o,u,v,-);function v,sigma = piecelin(x,y,u)d = diff(y)./diff(x); % First divided difference% Find subinterval indices, x(k) = x(j) = j;end% Evaluate interpolants = u - x(k);v = y(k) + s.*d(k);分片三次插值设s = x-xk , h = hkHermite插值（osculatory插值）pchiptx.msplinetx.mpchiptxpchip: piecewise cubic Hermit interpolating polynomial, Matlab中pchip算法基于Fritsch和Carlson1、如果左右导数是相反符号或有0，则dk=02、如果同号且区间长度相等，则是调和平均3、如果同号且区间长度不相等，则是加权调和平均0123456710121416182022三次样条三次样条也是分片三次插值函数。物理上的样条在满足插值限制的前提下，最小化势能。数学上的样条必须满足二次导数连续，且满足插值限制。参考文献：A Practical Guide to Spline, Carl de Boor。他也是Matlab的spline函数和spline工具箱的作者。三次样条二阶导 数连续等距“Not-a-knot”在边界上，把两段合成一段，在等距的情况下：这样，可以导致一个线性代数方程组：三次样条然后我们可以计算出各个节点的一阶导数0246810121416182022分片线性插值02468510152025full degree多项式插值0246810121416182022保形状Hermit插值0246810121416182022样条插值分 析 上面的图描述了光滑性和局部单调性（保形状）的一种 折衷。分片线性插值：保持单调性的，但光滑性比较差。Full degree多项式插值：无限可微，但不保持形状，特 别是在端点的地方。Pchip和spline插值在这两个极端之间。样条比pchip光滑 ，样条的两阶导数连续，而pchip一阶导数连续。不连续 的两阶导数隐含着不连续的曲率。人的眼睛可以检测出 图形上曲率的不连续。另一方面，pchip是保形状的，而 样条不一定保形状。pchiptx 和splinetx pchiptx和splinetx都是基于分片三次Hermite插值，在每个区间用monomial基 ：pchiptx.m,splinetx.m请比较一下与matlab中的pchip 与spline的区别123456-0.5-0.4-0.3-0.2-0.100.10.20.30.40.5Interpolatationinterpgui