第四章多项式与插值课件.ppt

,第四章 多项式与插值,4.1 MATLAB多项式,一、 多项式的建立1. MATLAB中多项式用行向量表示，其元素为 多项式的系数，且从左至右按降幂排列；2. 已知一个多项式的全部根 X求多项式系数的函 数是poly(X)，该函数返回以 X为全部根的一个 多项式 P（首项系数为1），当 X是一个长度为 m的向量时，P是一个长度为 m+1的向量。,3. 给定 n+1个点可以唯一确定一个 n 阶多项式，利 用polyfit可以确定多项式的系数。 调用格式为： p=polyfit(x, y, n) 其中 x, y是同维向量，代表数据点的横、纵坐标， n 是多项式的阶数。,二、 多项式计算1. 多项式求根 求多项式 p(x)的根的函数是roots(P)，这里，P是 p(x)的系数向量，该函数返回方程 p(x)=0 的全部根 (含重根，复根)。2. 多项式求值 求多项式 p(x)在某点或某些点的函数值的函数是polyval(P, x)。若x为一数值，则求多项式在该点的值；若x为向量或矩阵，则对向量或矩阵中的每个元素求其多项式的值。,例1 已知一个多项式 (1)计算f (x)=0 的全部根。(2)由方程f (x)=0的根构造一个多项式 g(x)，并与 f (x) 进行对比。(3)计算f (5)、f (7.8)、f (9.6)、f (12.3)的值。,P=3,0,4,-5,-7.2,5; X=roots(P) %求方程f(x)=0的根 G=poly(X) %求多项式g(x) X0=5,7.8,9.6,12.3; f=polyval(P,X0) %求多项式f(x)在给定点的值,3. 多项式的四则运算(1)多项式的加减法,注：多项式求值还有一个函数是polyvalm，其调用格式与polyval相同，但含义不同。polyvalm函数要求x为方阵，它以方阵为自变量求多项式的值。,function p3 = poly_add(p1,p2)n1=length(p1); n2 = length(p2);if n1=n2 p3 = p1 + p2; endif n1n2， p3 = p1 + zeros(1,n1-n2) ,p2;endif n1n2 ，p3 = zeros(1,n2-n1) ,p1 + p2; end,加法：c=poly_add(a,b) 减法: c=poly_add(a,-b),(2)多项式的乘法 函数conv(P1,P2)用于求多项式P1和P2的乘积。(3)多项式的除法 函数Q,r=deconv(P1,P2)用于对多项式P1和P2作除法运算。其中Q返回多项式P1除以P2的商式，r返回P1除以P2的余式。这里，Q和r仍是多项式系数向量。deconv是conv的逆函数，即有P1=conv(P2,Q)+r。,例2 设有两个多项式，计算：(1)求f(x)+g(x)、f(x)-g(x)。(2)求f(x)g(x)、f(x)/g(x)。,f=3,-5,2,-7,5,6;g=3,5,-3;poly_add(f, g) %求f(x)+g(x)poly_add(f, -g) %求f(x)-g(x)conv(f, g) %求f(x)*g(x)Q,r=deconv(f,g) %求f(x)/g(x)，商式送Q，余式送r。,4. 多项式的微分与积分（1）对多项式求导数的函数是： p=polyder(P) 求多项式P的导函数 p=polyder(P,Q) 求P*Q的导函数 p,q=polyder(P,Q) 求P/Q的导函数，导函数的 分子存入p，分母存入q。（2）对多项式的积分函数： d=poly_itg(c) d是多项式c积分后的系数，但 不包括积分常数,function py = poly_itg(p)n=length(p); py = p.*n:-1:1.(-1),0;,例3 求有理分式的导数。,P=3,5,0,-8,1,-5;Q=10,5,0,0,6,0,0,7,-1,0,-100;p,q=polyder(P,Q),若求多项式P的积分：,c=poly_itg(P),4.2 MATLAB插值,通常取 为多项式函数代数插值（多项式插值）,注:一次多项式插值 - 过两点直线; 二次多项式插值 - 过三点抛物线; 不用待定系数法 - (1)计算量大; (2)不易讨论误差;,几何意义：两条曲线有交点（公共点）,一、线性插值,线性插值是两个数据点的直线拟合,或,误差估计：,在MATLAB中，命令 interp1可做线性插值，调用格式为： yiinterp1( x, y, xi),其中 x 表示数据点横坐标的列数组，y 表示数据纵坐标的列数组（可以有多列）。,另外，interp1命令有三种可选参数,yi = interp1(x,y,xi, linear ) 线性插值（缺省）yi = interp1(x,y,xi, spline ) 三次样条yi = interp1(x,y,xi, cubic ) 三次插值,例3 已知数据表如下，分别求 y0.9，0.7，0.6，0.5 处 x 的值。,x = 0.0,0.25,0.5,0.75,1.0; y = 0.9126,0.8109,0.6931,0.5596,0.4055;yi = 0.9,0.7,0.6,0.5xi = interp1(y,x,yi,linear);yi , xi,ans 0.9000 0.0385 0.7000 0.4854 0.6000 0.6743 0.5000 0.8467,二、用幂级数做多项式插值,给定 n+1 个数据点：,过 n+1个点的 n 阶多项式可写为幂级数形式：,注：过 n1个数据点的 n 阶插值多项式是惟一的。,对 n1个数据点，设 ，则得到 n+1个线性方程，可以表示为矩阵形式,求解该方程组可确定系数（或用 polyfit(x,y,n)确定）,例3 求下列数据点拟合多项式的系数，并求当x2.101和4.234处 y 的值，并画出数据点和曲线。,x=1.1, 2.3, 3.9, 5.1;y=3.887, 4.276, 4.651, 2.117 ;n=length(x)-1 ;a(:,n+1)=ones(size(x);a(:,n)=x;for j=n-1:-1:1 a(:,j)=a(:,j+1).*x;endcoeff=ay;,xi=2.101, 4.234;yi=zeros(size(xi);for k=1:n+1 yi = yi + coeff(k)*xi.(n+1-k);endxp=1.1:0.05:5.1;yp=zeros(size(xp);for k=1:n+1 yp=yp + coeff(k)*xp.(n+1-k);endplot(xp,yp, x,y, ro),三、Lagrange插值多项式,1.插值基函数,形状函数的图形如下（n8）,2.Lagrange插值多项式,function fi = Lagran_(x, f, xi)fi=zeros(size(xi);np1=length(f);for i=1:np1 z=ones(size(xi); for j=1:np1 if i=j z = z.*(xi - x(j)/(x(i)-x(j); end end fi=fi+z*f(i);endreturn,3. MATLAB程序实现:,调用格式： yi =Lagran_( x, y, xi),clearx = 1.1, 2.3, 3.9, 5.1;y =3.887, 4.276, 4.651, 2.117;xi = 2.101, 4.234;yi = Lagran_(x, y, xi);,例4：写出拟合下面三个数据点的Lagrange插值公式，并计算 x2.101、4.234时 y 的值。,4. 截断误差：,例5 用 的5个等距点对函数进行插值估计。,插值结果及误差分布如图,可见，误差峰值出现在端点附近的区间里，这是由于 的局部峰值在端点附近。,减小误差的方案：（1）减小插值区域，即 ba；（2）增加数据点个数；（3）使用可变间距的数据点（Chebyshev点）。,总结：（1）尽可能在小区间上使用多项式插值；（2）只能在一定范围内依靠增加插值点个数提高插值精度，如果插值点个数过多往往会适得其反。,5. Lagrange 插值公式的微分与积分,插值公式,微分,实际上， 是一个拟合如下数据点的 n 阶多项式,一般地， 是拟合数据点的多项式,所以，多项式 可通过拟合待定数据点的 n 阶多项式表示为幂级数形式。,对所有 i， 的幂级数形式可用函数 shape_pw 计算,function p = shape_pw(x)np = length(x);for j=1:np y = zeros(1,np); y(j) = 1; p(j,:)=polyfit(x,y,np-1);end,如例4也可求解如下：,x = 1.1, 2.3, 3.9, 5.1;y = 3.887, 4.276, 4.651, 2.117;xi = 2.101 ,4.234;np = length(x);p=shape_pw(x);s=0;for i=1:np s = s+p(i,:).*y(i); endsyi=polyval(s,xi),结果为：yi 4.1457 4.3007,为计算Lagrange插值多项式的一阶导数，可用polyder函数将 p 的每一行转换为一阶导数的系数数组。,x = 1.1, 2.3, 3.9, 5.1;y = 3.887, 4.276, 4.651, 2.117;xi = 2.101 ,4.234;np = length(x);p=shape_pw(x);s=0;for i=1:np s = s+polyder(p(i,:).*y(i); endyi=polyval(s,xi),结果为：yi 0.6292 -1.4004,四、Chebyshev多项式,等距Lagrange插值的误差在中间部分最小，越靠近端点处误差越大，降低误差的一种方法是改变插值点的分布，即增大定义域中间部分的插值步长，同时减少定义域两端的插值步长，但是最优的插值点分布取决于插值多项式的目的。,若目的是近似一个函数，则最优插值点是Chebyshev多项式的零点，因为（1）此时L(x)分布是最平坦的分布；（2）误差值不会像等距插值那样随不同插值区间 而变化。,Chebyshev多项式：,当 时，Chebyshev多项式的图形如下,幂级数形式Chebyshev多项式的系数可由函数Cheby_ pw计算,function pn = Cheby_pw(n)pbb=1; if n=0, pn=pbb; break; endpb=1 0; if n=1, pn=pb; break; endfor i=2:n; pn= 2*pb,0 - 0, 0, pbb ; pbb=pb; pb=pn;end,调用格式：p=Cheby_pw(n)n为多项式阶数，p是系数的行数组,Chebyshev多项式的根系数可由 roots命令计算 roots(Cheby_ pw(n),若令 ，则由,得,有 k 个根，都在-1,1，也可用下面公式计算：,若插值区间为a,b，可建立映射,下图比较了区间 上9个Chebyshev点和等距点条件下L(x)的分布,五、Lobatto多项式,Lobatto点是Chebyshev多项式的零点加上 x -1和 x1对区间a,b，k +1个Lobatto点可由 k 阶Chebyshev多项式得到,六、Legendre多项式,Legendre多项式：,n 阶Legendre多项式的系数可由函数Legen_pw计算：,function pn = Legen_pw(n)pbb=1; if n=0, pn=pbb; break; endpb=1 0; if n=1, pn=pb; break; endfor i=2:n; pn= ( (2*i-1)*pb,0 - (i-1)*0, 0, pbb )/i; pbb=pb; pb=pn;end,调用：p=Legen_pw(n),n 阶Legendre多项式的根可由roots计算。,七、三次Hermite多项式,可以同时拟合函数值与导数值的多项式称为Hermite插值多项式或密切多项式。,三次Hermite多项式：,可以拟合两点函数值和导数值,已知两点 和 及它们的函数值和导数值，可建立方程组：,若用有限差分近似表示导数边界条件,可通过Lagrange插值确定该方程组的系数。,当y(x)、x(y)含有奇点时，常需引入一个参数 s，并将 x 和 y 表示为 x(s) 和 y(s)，假设点A处 s =0，B处 s =1，则 x，y 表示为关于 s 的三次多项式：,x(s) 和 y(s)的边界条件可写为：,其中a、b是任意参数，在一定程度上影响曲线形状。,例6 确定一条经过点A、B的曲线，使满足,解：引入 s，按上述方法可得到两方程组,分别求解：,a=3;b=3;c=0 0 0 1;0 0 1 0;1 1 1 1;3 2 1 01;a;4;0d=0 0 0 1;0 0 1 0;1 1 1 1;3 2 1 01;0;2;bs=0:0.01:1;x=polyval(c,s)y=polyval(d,s);plot(x,y),结果：c=-3 3 3 1 d=1 0 0 1,若用差分近似，设置四个数据点：,取 ，求解,z=0.01; a=3; b=3;s(1) = 0; x(1) = 1; y(1) = 1;s(2) = z; x(2) = 1+z*a; y(2) = 1;s(3) = 1 - z; x(3) = 4; y(3) = 2 - z*b;s(4) = 1; x(4) = 4; y(4) = 2;,c=polyfit(s,x,length(s)-1);d=polyfit(s,y,length(s)-1);ss=0:0.1:1;xp = polyval(c,ss);yp = polyval(d,ss);plot(xp,yp)ylabel(y)xlabel(x),结果是：c=-3.9021 3.1231 2.9691 1d=1.0307 -0.0309 0.0002 1,八、样条插值,问题：结点增多，多项式次数增高，逼近精度越 好？未必！多结点高次插值往往在局部误 差更大Runge（龙格）现象。实用：采用分段低次插值 有分段线性，分段二次插值等，几何意义,缺点：分段插值函数只能保证连续性， 不能保证 光滑性。,分段插值可以得到整体连续函数，但在连接点处一般不光滑，而 Hermite插值虽然在连接点处一阶光滑，但整体插值由于结点多次数高而有可能发生龙格现象。,既想分段插值，又想在结点处保持光滑，甚至二阶光滑三次样条。,希望：,定义：设有 n +1个点的点列,若函数 满足：,(3) 在整个区间 内具有二阶连续导数。,(2) 在每个小区间,上是三次多项式；,(1) 此时 叫插值函数；,则称 为点列的三次样条插值函数或三次样条 多项式，简称三次样条。,则样条函数可写为,考虑区间,令,列方程组,求解得：,所以,则在区间 上,同理，在区间 上,由在节点 处连续，式、相等，消去 ，得：,式除了对第一个点和最后一个点不成立外，其它各点均成立，所以共有 n-2 个方程，而待定的 个数为 n，所以还需两个额外条件，可由下述边界条件提供。,常见到的边界条件有三种：,（a）在端点处指定 ，若已知 ，则方程组为：,该方程组为含 n-2 个未知量的 n-2 个方程其系数矩阵为三对角阵，可用追赶法求解。,大多数情形中，可令 ，等价于在几何上趋于端点的曲线变为直线。,（b）从内部外推 ， 利用 对 进行外推,将其代入方程组的第一个方程中，整理得：,类似，利用 对 进行外推,将其代入方程组的最后一个方程，整理得：,此时，同样得到一三对角方程组，求解即可。,（c）循环边界条件 适用于封闭曲线，即第一个点和最后一个点相 同，它们的导数值也相同，此时相当于有 n-1个点， 相应有 n-1个方程，可以求解。,在MATLAB中，三次样条插值运算实现如下： yi=interp1( x, y, xi, spline)或 yi=spline( x, y, xi )其中 x，y 都是向量形式的点，xi是进行插值的点的横坐标向量，yi 为插值函数值。,例7 求三次样条插值并作图比较。,%以s为参数，分别做x，y的样条函数xx = -1 -0.866 -0.5 0 0.5 0.866 1 1 1.0402 1.1500 1.3. 1.54 1.8280 2.1736 2.5883 3.0860;yy = 0 -0.25 -0.433 -0.5 -0.433 -0.25 0 0 0.15 0.2598 0.3. 0.3 0.3 0.3 0.3 0.3;s=1:length(xx);sp=1:length(xx)/100:length(xx);xp=spline(s,xx,sp);yp=spline(s,yy,sp);plot(xp,yp); hold onplot(xx,yy, ro);xlabel(x); ylabel(y);axis(-1 3.5 -1 1),若不取s为参数，直接以x，y做样条函数，曲线不光滑。,xi=-1:1/100:3.5;yi=spline(xx,yy,xi);plot(xi,yi,xx,yy, ro),九、二维插值,双线性插值 用两次线性插值对二维函数的数据表格做插值，二维表格是矩形网格点 上的函数值阵列,设需估计如图所示矩形区域 上一些点的函数值,在 y 方向上做线性插值，则E、F处的值分别为：,然后做 的线性插值,将这两步结合为一个方程，即为双线性插值：,在MATLAB中，双线性插值命令为table2 g = table2(tab, x , y)其中tab为一二维数据表格，第一列是 值的数组，第一行是 值的数组，按升序排列，剩余行列为 ，x，y是待插值点的横、纵坐标值，可以是标量、向量或矩阵。,2. 双Lagrange插值 在两个维度上利用两次Lagrange插值的方法。,其中形状函数,当MN2时，方程化简为双线性插值表达式。,作业： 4.5，4.6，4.9，4.12，4.16，4.20，4.22，4.26 9.1，9.2要求：编写程序，上机运行后保存结果。,