复旦大学 Matlab课件 .pptVIP

Matlab Math Cleve Morler著 陈文斌（wbchen@fudan.edu.cn) 复旦大学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 = V\y’ Vandermonde矩阵是非奇异的，但条件数是非常坏的 Polyinterp(Lagrange插值形式) 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 -0.5 0 0.5 1 1.5 2 2.5 3 3.5 -10 -5 0 5 10 15 20 25 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 = x^3-2*x-5 0 1 2 3 4 5 6 7 6 8 10 12 14 16 18 20 22 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,’-’); 0 1 2 3 4 5 6 7 10 12 14 16 18 20 22 分片线性插值 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) = u u(k+1) n = length(x); k = ones(size(u)); for j = 2:n-1 k(u = x(j)) = 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、如果同号且区间长度不相等，则是加权调和平均 0 1 2 3 4 5 6 7 10 12 14 16 18 20 22 三次样条 三次样条也是分片三次插值函数。 物理上的样条在满足插值限制的前提下，最小化势能。 数学上的样条必须满足二次导数连续，且满足插值限制。 参考文献： A Practical Guide to Spline, Carl de Boor。 他也是Matlab的spline函数和spline工具箱的作者。 三次样条 二阶导数连续 等距 *