复旦大学 Matlab课件 (2).ppt

1、Matlab Math,Cleve Morler著 陈文斌（) 复旦大学2002,多项式和样条,在平面上给定n个点(xk,yk),可以唯一确定一个最多n-1次的多项式通过这些点，这个多项式叫插值多项式,插值多项式,P(xk ) = yk , k = 1,2,n,Lagrange插值形式,插值多项式例子,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 = Vy,Vandermonde矩阵是非奇异的，但条件数是非常坏的,Polyinterp(Lagran

2、ge插值形式),function v = polyinterp(x,y,u) n = length(x); 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); end,u = -.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 =

3、 x3-2*x-5,Polyinterp(另外的例子),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,-);,分片线性插值,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

4、) = j; end % Evaluate interpolant s = u - x(k); v = y(k) + s.*d(k);,分片三次插值,设s = x-xk , h = hk,Hermite插值（osculatory插值）,pchiptx.m splinetx.m,pchiptx,pchip: piecewise cubic Hermit interpolating polynomial, Matlab中pchip算法基于Fritsch和Carlson,1、如果左右导数是相反符号或有0，则dk=0 2、如果同号且区间长度相等，则是调和平均,3、如果同号且区间长度不相等，则是加权调和

5、平均,三次样条,三次样条也是分片三次插值函数。 物理上的样条在满足插值限制的前提下，最小化势能。 数学上的样条必须满足二次导数连续，且满足插值限制。 参考文献： A Practical Guide to Spline, Carl de Boor。 他也是Matlab的spline函数和spline工具箱的作者。,三次样条,二阶导数连续,等距,“Not-a-knot”,在边界上，把两段合成一段，在等距的情况下：,这样，可以导致一个线性代数方程组：,三次样条,然后我们可以计算出各个节点的一阶导数,分 析,上面的图描述了光滑性和局部单调性（保形状）的一种折衷。 分片线性插值：保持单调性的，但光滑性比较差。 Full degree多项式插值：无限可微，但不保持形状，特别是在端点的地方。 Pchip和spline插值在这两个极端之间。样条比pchip光滑，样条的两阶导数连续，而pchip一阶导数连续。不连续的两阶导数隐含着不连续的曲率。人的眼睛可以检测出图形上曲率的不连续。另一方面，pchip是保形状的，而样条不一定保形状。,pchiptx 和splinetx,pchiptx和splinetx都是基于分片三次H