多项式与插值

多项式与插值第1页，共60页，2023年，2月20日，星期一

补充：将多项式的数值形式转化为字符串或符号表达式的形式poly2str(p,’s’)将多项式p转化为以s为变量的字符串形式poly2sym(p)将多项式p转化为符号表达式poly2sym(p,’s’)将多项式p转化为以s为变量的符号表达式>>p=[2,5,0,4,1,4];

>>A=poly2sym(p)>>B=poly2str(p,'x')

A=4+2x5+5x4+4x2+xB=2x^5+5x^4+4x^2+x+4第2页，共60页，2023年，2月20日，星期一2、用多项式的全部根X来生成

poly(X)该函数返回以X为全部根的一个多项式P（首项系数为1），当X是一个长度为m的向量时，P是一个长度为m+1的向量。3、用polyfit生成多项式的系数给定n+1个点可以唯一确定一个n阶多项式调用格式为：

p=polyfit(x,y,n)

其中x,y是同维向量，代表n+1个数据点的横、纵坐标，n

是多项式的阶数。第3页，共60页，2023年，2月20日，星期一二、多项式的运算1.、求根运算

roots(P)P是多项式p(x)的系数向量，该函数返回p(x)=0的全部根（含重根，复根）2、求值运算

polyval(P,x)

求多项式P在某点或某些点的函数值；若x为一数值，则求多项式P在该点处的值；若x为向量或矩阵，则求多项式P在向量或矩阵中的每个元素处的值第4页，共60页，2023年，2月20日，星期一

polyvalm(P,x)此处x必须为方阵，它以方阵为自变量求多项式P的值第5页，共60页，2023年，2月20日，星期一例已知一个多项式(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),与f(x)的区别是首项系数为1X0=[5,7.8,9.6,12.3];f=polyval(P,X0)%求多项式f(x)在给定点的值第6页，共60页，2023年，2月20日，星期一3.多项式的四则运算(1)多项式的加减法当两多项式向量的大小相等，直接两向量相加（减）当两多项式次数不同时，低阶多项式必须用首0填补，使其与高阶多项式次数相同>>P1=[2,5,0,4,1,4];>>P2=[0,0,5,1,3,2];>>P3=P1+P2;>>P=poly2str(P3,'x')P=2x^5+5x^4+5x^3+5x^2+4x+6第7页，共60页，2023年，2月20日，星期一(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第8页，共60页，2023年，2月20日，星期一例：已知三个一次多项式的零点分别为0.4、0.8、1.2，求三个一次多项式的乘积。法1：>>X=[0.4,0.8,1.2];c=poly(X),P=poly2sym(c)运行后的结果：c=1.0000-2.40001.7600-0.3840P=x^3-12/5*x^2+44/25*x-48/125第9页，共60页，2023年，2月20日，星期一法2：>>P1=poly(0.4);P2=poly(0.8);P3=poly(1.2);c=conv(conv(P1,P2),P3),P=poly2sym(c)运行后的结果：c=1.0000-2.40001.7600-0.3840P=x^3-12/5*x^2+44/25*x-48/125第10页，共60页，2023年，2月20日，星期一例设有两个多项式，计算：(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];g1=[0,0,0,3,5,-3];f+g1%求f(x)+g(x)f-g1

%求f(x)-g(x)conv(f,g)%求f(x)*g(x)[Q,r]=deconv(f,g)%求f(x)/g(x)，商式送Q，余式送r。第11页，共60页，2023年，2月20日，星期一4、多项式的微分与积分（1）对多项式求导数：

p=polyder(P)

求多项式P的导函数

p=polyder(P,Q)

求P*Q的导函数

[p,q]=polyder(P,Q)

求P/Q的导函数，导函数的分子存入p，分母存入q。（2）对多项式求积分：

d=polyint(P)d是多项式P积分后的系数，但不包括积分常数第12页，共60页，2023年，2月20日，星期一例求有理分式的导数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=polyint(P)第13页，共60页，2023年，2月20日，星期一§4.2MATLAB插值常遇到，对于y=f(x),它是存在的，可并不知道它的具体函数表达式，只知道离散数据希望找到一个简单函数且满足被插函数插值函数（常取多项式）插值节点插值条件：插值区间：的余项或误差：第14页，共60页，2023年，2月20日，星期一若找到P(x)=a0+a1x+···+anxn，则称相应的插值为代数插值（或多项式插值）定理：若插值结点x1,x2,…,xn+1

是(n+1)个互异点，则满足插值条件P(xk)=yk(k=0,1,…,n)的n次插值多项式

P(x)=a0+a1x+……+anxn存在而且是唯一的此定理表明：在n+1个互异节点构造出来的n次多项式是唯一的，也就是说，后面将要学到的幂级数插值多项式、lagrange插值多项式、牛顿插值多项式是恒等的第15页，共60页，2023年，2月20日，星期一证明:由插值条件P(x1)=y1P(x2)=y2··············P(xn+1)=yn+1系数行列式故方程组有唯一解.

从而插值多项式P(x)存在而且是唯一的.第16页，共60页，2023年，2月20日，星期一例已知误差函数在四个点处函数值

x0 0.6000 1.2000 1.8000Erf(x)

0 0.6039 0.9103 0.9891构造3次多项式P(x)逼近Erf(x)设P(x)=a0+a1x+a2x2+a3x3,令P(xi)=Erf(xi)得求解，得a0=0,a1=1.293,a2=-0.5099,a3=0.0538所以,P(x)=1.293x–0.5099x2+0.0538x3第17页，共60页，2023年，2月20日，星期一MATLAB计算程序:x=0:0.6:1.8;y=erf(x);x=x';y=y';A=[ones(4,1)xx.^2x.^3];p=A\y;a0=p(1);a1=p(2);a2=p(3);a3=p(4);t=0:0.2:2;u=a0+a1*t+a2*t.^2+a3*t.^3;plot(x,y,'o',t,u)第18页，共60页，2023年，2月20日，星期一一、线性插值(一次多项式插值）线性插值是两个数据点的直线拟合或误差估计：第19页，共60页，2023年，2月20日，星期一在MATLAB中，命令interp1可做线性插值，

yi＝interp1(x,y,xi)其中：

x表示已知数据点横坐标的向量；y表示已知数据点纵坐标的向量；xi为插值点的横坐标，输出的yi为插值点的函数值第20页，共60页，2023年，2月20日，星期一例已知数据表如下，分别求x＝0.9，0.7，0.6，0.5

处y的值。xy0.912600.81090.250.69310.500.55960.750.40551.00x=[0.9126,0.8109,0.6931,0.5596,0.4055];y=[0.0,0.25,0.5,0.75,1.0];

xi=[0.9,0.7,0.6,0.5]';yi=interp1(x,y,xi,'linear');[xi,yi]ans＝0.90000.03100.70000.48540.60000.67430.50000.8467第21页，共60页，2023年，2月20日，星期一另外，interp1命令有5种可选参数yi=interp1(x,y,xi,’linear’)

线性插值（缺省）yi=interp1(x,y,xi,’nearst’)

最近邻点插值yi=interp1(x,y,xi,’pchip’)分段三次Hermite插值yi=interp1(x,y,xi,’cubic’)同上yi=interp1(x,y,xi,’spline’)分段三次样条插值（建议使用三次样条插值）若已知的数据点x是等距节点，插值速度可提高，此时，填写最后一个插值方法时，以星号引导，即为‘*linear’,’*spline’,’*cubic’,’*pchip’

,’*nearst’

第22页，共60页，2023年，2月20日，星期一例：假设已知的数据点来自函数试根据生成的数据进行插值,得到较光滑的曲线解：>>x=0:0.12:1;>>y=(x.^2-3*x+5).*exp(-5*x).*sin(x);>>plot(x,y,'-o')第23页，共60页，2023年，2月20日，星期一可看出，由这样的数据直接连线绘制出来的曲线十分粗糙，可再选择一组插值点，然后用interp1插值第24页，共60页，2023年，2月20日，星期一x=0:0.12:1;y=(x.^2-3*x+5).*exp(-5*x).*sin(x);x0=0:0.02:1;y0=(x0.^2-3*x0+5).*exp(-5*x0).*sin(x0);y01=interp1(x,y,x0);y02=interp1(x,y,x0,'cubic');y03=interp1(x,y,x0,'spline');y04=interp1(x,y,x0,'nearest');plot(x,y,'mo',x0,y0,'k',x0,y01,'r:',x0,y02,'g:',...x0,y02,'g:',x0,y03,'m:',x0,y04,'b:')第25页，共60页，2023年，2月20日，星期一

interp1默认的linear插值得到的曲线和plot画出来的曲线一样粗糙,而’nearest’选项得到的插值效果就更差了。而采用’cubic’和’spline’选项得到的插值更接近与理论值.第26页，共60页，2023年，2月20日，星期一可看出：

interp1默认的线性插值得到的曲线和plot画出来的曲线一样粗糙，因为该方法是对各个点直线连接。而’nearest’选项得到的插值效果就更差了。而采用’cubic’和’spline’选项得到的插值更接近与理论值。事实上，采用样条插值算法得到的插值十分逼近理论值，甚至用肉眼难以分辨。第27页，共60页，2023年，2月20日，星期一二、用幂级数做多项式插值给定n+1个数据点：过n+1个点的n

阶多项式可写为幂级数形式：注：过n＋1个数据点的n阶插值多项式是唯一的。第28页，共60页，2023年，2月20日，星期一对n＋1个数据点，设，则得到n+1个线性方程，可以表示为矩阵形式求解该方程组可确定系数（或用polyfit(x,y,n)确定）第29页，共60页，2023年，2月20日，星期一例求用下列数据点拟合出的多项式的系数，并求当x取2.101、4.234时y的值，并画出已知数据点和拟合曲线。clear,clf,holdoffx=[1.1,2.3,3.9,5.1]';y=[3.887,4.276,4.651,2.117]';n=length(x)-1;c=polyfit(x,y,n);xi=[2.101,4.234];yi=polyval(c,xi)plot(x,y,'ro')holdonxp=1.1:0.05:5.1;yp=polyval(c,xp);pauseplot(xp,yp)yi=4.14574.3007第30页，共60页，2023年，2月20日，星期一第31页，共60页，2023年，2月20日，星期一三、Lagrange插值多项式1、给定n+1个数据点：插值基函数：满足插值条件：的次数不超过n的lagrange多项式g(x)为：第32页，共60页，2023年，2月20日，星期一functiony=lagrange1(x0,y0,x)np=length(x0);y=zeros(size(x));fori=1:npz=ones(size(x));forj=1:npifj~=iz=z.*(x-x0(j))/(x0(i)-x0(j));endendy=y+y0(i)*zend

2、MATLAB程序（法1）:调用格式：

y=lagrange1(x0,y0,x)第33页，共60页，2023年，2月20日，星期一clearx0=[1.1,2.3,3.9,5.1];y0=[3.887,4.276,4.651,2.117];x=[2.101,4.234];y=lagrange1(x0,y0,x)例：写出拟合下面四个数据点的Lagrange插值公式，并计算x＝2.101、4.234时y

的值。x0y01.13.8872.34.2763.94.6515.12.117结果为y=4.14574.3007第34页，共60页，2023年，2月20日，星期一由于lagrange插值基函数li(x)在节点xk(k=1,2,…,n)处满足又由于过n+1个数据点的n次插值多项式是唯一的,故li(x)也可看做是拟合下列数据点的n次多项式可编写程序求出

2、MATLAB程序实现（法2）:进而求出lagrange插值公式第35页，共60页，2023年，2月20日，星期一functionp=lag_base(x)np=length(x);forj=1:npy=zeros(1,np);y(j)=1;p(j,:)=polyfit(x,y,np-1);end其调用格式为：p=lag_base(x)其中:

x是数据点的横坐标数组，p是一个矩阵，它的第i行为li(x)

的幂系数。先编写程序求出第36页，共60页，2023年，2月20日，星期一如前例也可求解如下：x=[1.1,2.3,3.9,5.1];y=[3.887,4.276,4.651,2.117];xi=[2.101,4.234];p=lag_base(x);np=length(x);s=0;fori=1:nps=s+p(i,:).*y(i);endsyi=polyval(s,xi)结果为：yi＝

4.14574.3007第37页，共60页，2023年，2月20日，星期一3、Lagrange插值公式的微分插值公式微分为计算Lagrange插值多项式的一阶导数，可用polyder函数将p的每一行转换为一阶导数的系数数组。第38页，共60页，2023年，2月20日，星期一例：求出拟合下面四个数据点的Lagrange插值多项式在x＝2.101、4.234处的一阶导数值。x0y01.13.8872.34.2763.94.6515.12.117x=[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=lag_base(x);s=0;fori=1:nps=s+polyder(p(i,:)).*y(i);endyi=polyval(s,xi)结果为：yi＝

0.6292-1.4004第39页，共60页，2023年，2月20日，星期一例考虑一个著名的例子假设已知其中一些点的坐标，则可以采用下面的命令进行Lagrange插值，得出相应的插值曲线x0=-1:0.2:1;y0=1./(1+25*x0.^2);x=-1:0.01:1;y=lagrange1(x0,y0,x);ya=1./(1+25*x.^2);plot(x,y,'b:p',x,ya,'r')4、Runge现象第40页，共60页，2023年，2月20日，星期一可看到：由lagrange插值得出的效果和精确值相差甚远，这种多项式阶次越高越发散的现象称为Runge现象所以对此例来说，传统的lagrange算法失效第41页，共60页，2023年，2月20日，星期一下面考虑用interp1来解决此问题x0=-1:0.2:1;y0=1./(1+25*x0.^2);x=-1:0.01:1;ya=1./(1+25*x.^2);y1=interp1(x0,y0,x,'cubic');y2=interp1(x0,y0,x,'spline');plot(x,ya,x,y1,':',x,y2,'--')第42页，共60页，2023年，2月20日，星期一interp1函数效果可见：Matlab中提供的算法不存在Runge现象，可大胆放心使用第43页，共60页，2023年，2月20日，星期一总结：（1）尽可能在小区间上使用多项式插值；（2）只能在一定范围内依靠增加插值点个数提高插值精度，如果插值点个数过多往往会适得其反。第44页，共60页，2023年，2月20日，星期一定义:

若已知函数f(x)在点x1,

x2,···,xn+1

处的值f(x1),f(x2),···,f(xn+1),如果

i≠j,则n阶均差(j=1,2,…,n)一阶均差二阶均差(j=1,2,…,n-1

)四、牛顿插值（均差插值）第45页，共60页，2023年，2月20日，星期一x -2 -1 0 1 3y -56 -16 -2 -2 4例由函数表求各阶均差解:按公式计算各阶差商如下x f(x)一阶二阶三阶四阶差商差商差商差商-2 -56

-1 -16

40

0 -214 -13

1 -20 -7 23 4 3 1 20第46页，共60页，2023年，2月20日，星期一程序：x=[-2-1013];y=[-56-16-2-24];n=length(x);A=zeros(n,n);A(:,1)=y';forj=2:nfori=j:nA(i,j)=(A(i,j-1)-A(i-1,j-1))/(x(i)-x(i-j+1));endendA-560000-1640000-214-1300-20-72043120第47页，共60页，2023年，2月20日，星期一有了均差，就能得到牛顿插值公式（均差插值公式）：第48页，共60页，2023年，2月20日，星期一以3阶均差插值公式为例，改写形式第49页，共60页，2023年，2月20日，星期一function[A,C,N]=newton_poly(x,y)%x为已知点的横坐标，y为已知点的纵坐标%A为均差矩阵，C为牛顿插值多项式的系数向量%N为牛顿插值多项式%先产生均差矩阵An=length(x);A=zeros(n,n);A(:,1)=y';forj=2:nfori=j:nA(i,j)=(A(i,j-1)-A(i-1,j-1))/(x(i)-x(i-j+1));endend牛顿插值的Matlab算法：第50页，共60页，2023年，2月20日，星期一%下面求牛顿插值多项式系数向量CC=A(n,n);fork=(n-1):-1:1C=conv(C,poly(x(k)));d=length(C);C(d)=C(d)+A(k,k);end%最后求牛顿插值多项式NN=poly2sym(C);第51页，共60页，2023年，2月20日，星期一例：给出节点x=[-2.15-1.000.011.022.033.25]，y=[17.037.241.052.0317.0623.05]，作五阶牛顿插值多项式和差商。>>x=[-2.15-1.000.011.022.033.25];y=[17.037.241.052.0317.0623.05];[A,C,N]=newton_poly(x,y)第52页，共60页，2023年，2月20日，星期一结果：A=17.0300000007.2400-8.513000001.0500-6.12871.10390002.03000.97033.51440.76040017.060014.88126.88661.11290.0843023.