多項式與插值

多項式與插值第一頁，共六十頁，2022年，8月28日

補(bǔ)充：將多項式的數(shù)值形式轉(zhuǎn)化為字符串或符號表達(dá)式的形式poly2str(p,’s’)將多項式p轉(zhuǎn)化為以s為變量的字符串形式poly2sym(p)將多項式p轉(zhuǎn)化為符號表達(dá)式poly2sym(p,’s’)將多項式p轉(zhuǎn)化為以s為變量的符號表達(dá)式>>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第二頁，共六十頁，2022年，8月28日2、用多項式的全部根X來生成

poly(X)該函數(shù)返回以X為全部根的一個多項式P（首項系數(shù)為1），當(dāng)X是一個長度為m的向量時，P是一個長度為m+1的向量。3、用polyfit生成多項式的系數(shù)給定n+1個點(diǎn)可以唯一確定一個n階多項式調(diào)用格式為：

p=polyfit(x,y,n)

其中x,y是同維向量，代表n+1個數(shù)據(jù)點(diǎn)的橫、縱坐標(biāo)，n

是多項式的階數(shù)。第三頁，共六十頁，2022年，8月28日二、多項式的運(yùn)算1.、求根運(yùn)算

roots(P)P是多項式p(x)的系數(shù)向量，該函數(shù)返回p(x)=0的全部根（含重根，復(fù)根）2、求值運(yùn)算

polyval(P,x)

求多項式P在某點(diǎn)或某些點(diǎn)的函數(shù)值；若x為一數(shù)值，則求多項式P在該點(diǎn)處的值；若x為向量或矩陣，則求多項式P在向量或矩陣中的每個元素處的值第四頁，共六十頁，2022年，8月28日

polyvalm(P,x)此處x必須為方陣，它以方陣為自變量求多項式P的值第五頁，共六十頁，2022年，8月28日例已知一個多項式(1)計算f(x)=0的全部根(2)由方程f(x)=0的根構(gòu)造一個多項式g(x)，并與f(x)

進(jìn)行對比(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)的區(qū)別是首項系數(shù)為1X0=[5,7.8,9.6,12.3];f=polyval(P,X0)%求多項式f(x)在給定點(diǎn)的值第六頁，共六十頁，2022年，8月28日3.多項式的四則運(yùn)算(1)多項式的加減法當(dāng)兩多項式向量的大小相等，直接兩向量相加（減）當(dāng)兩多項式次數(shù)不同時，低階多項式必須用首0填補(bǔ)，使其與高階多項式次數(shù)相同>>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第七頁，共六十頁，2022年，8月28日(2)多項式的乘法

conv(P1,P2)求多項式P1和P2的乘積(3)多項式的除法

[Q,r]=deconv(P1,P2)求P1/P2；其中Q為多項式P1除以P2的商式，r為P1除以P2的余式。這里，Q和r仍是多項式系數(shù)向量。deconv是conv的逆函數(shù)，即有P1=conv(P2,Q)+r第八頁，共六十頁，2022年，8月28日例：已知三個一次多項式的零點(diǎn)分別為0.4、0.8、1.2，求三個一次多項式的乘積。法1：>>X=[0.4,0.8,1.2];c=poly(X),P=poly2sym(c)運(yùn)行后的結(jié)果：c=1.0000-2.40001.7600-0.3840P=x^3-12/5*x^2+44/25*x-48/125第九頁，共六十頁，2022年，8月28日法2：>>P1=poly(0.4);P2=poly(0.8);P3=poly(1.2);c=conv(conv(P1,P2),P3),P=poly2sym(c)運(yùn)行后的結(jié)果：c=1.0000-2.40001.7600-0.3840P=x^3-12/5*x^2+44/25*x-48/125第十頁，共六十頁，2022年，8月28日例設(shè)有兩個多項式，計算：(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。第十一頁，共六十頁，2022年，8月28日4、多項式的微分與積分（1）對多項式求導(dǎo)數(shù)：

p=polyder(P)

求多項式P的導(dǎo)函數(shù)

p=polyder(P,Q)

求P*Q的導(dǎo)函數(shù)

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

求P/Q的導(dǎo)函數(shù)，導(dǎo)函數(shù)的分子存入p，分母存入q。（2）對多項式求積分：

d=polyint(P)d是多項式P積分后的系數(shù)，但不包括積分常數(shù)第十二頁，共六十頁，2022年，8月28日例求有理分式的導(dǎo)數(shù)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)第十三頁，共六十頁，2022年，8月28日§4.2MATLAB插值常遇到，對于y=f(x),它是存在的，可并不知道它的具體函數(shù)表達(dá)式，只知道離散數(shù)據(jù)希望找到一個簡單函數(shù)且滿足被插函數(shù)插值函數(shù)（常取多項式）插值節(jié)點(diǎn)插值條件：插值區(qū)間：的余項或誤差：第十四頁，共六十頁，2022年，8月28日若找到P(x)=a0+a1x+···+anxn，則稱相應(yīng)的插值為代數(shù)插值（或多項式插值）定理：若插值結(jié)點(diǎn)x1,x2,…,xn+1

是(n+1)個互異點(diǎn)，則滿足插值條件P(xk)=yk(k=0,1,…,n)的n次插值多項式

P(x)=a0+a1x+……+anxn存在而且是唯一的此定理表明：在n+1個互異節(jié)點(diǎn)構(gòu)造出來的n次多項式是唯一的，也就是說，后面將要學(xué)到的冪級數(shù)插值多項式、lagrange插值多項式、牛頓插值多項式是恒等的第十五頁，共六十頁，2022年，8月28日證明:由插值條件P(x1)=y1P(x2)=y2··············P(xn+1)=yn+1系數(shù)行列式故方程組有唯一解.

從而插值多項式P(x)存在而且是唯一的.第十六頁，共六十頁，2022年，8月28日例已知誤差函數(shù)在四個點(diǎn)處函數(shù)值

x0 0.6000 1.2000 1.8000Erf(x)

0 0.6039 0.9103 0.9891構(gòu)造3次多項式P(x)逼近Erf(x)設(shè)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第十七頁，共六十頁，2022年，8月28日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)第十八頁，共六十頁，2022年，8月28日一、線性插值(一次多項式插值）線性插值是兩個數(shù)據(jù)點(diǎn)的直線擬合或誤差估計：第十九頁，共六十頁，2022年，8月28日在MATLAB中，命令interp1可做線性插值，

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

x表示已知數(shù)據(jù)點(diǎn)橫坐標(biāo)的向量；y表示已知數(shù)據(jù)點(diǎn)縱坐標(biāo)的向量；xi為插值點(diǎn)的橫坐標(biāo)，輸出的yi為插值點(diǎn)的函數(shù)值第二十頁，共六十頁，2022年，8月28日例已知數(shù)據(jù)表如下，分別求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第二十一頁，共六十頁，2022年，8月28日另外，interp1命令有5種可選參數(shù)yi=interp1(x,y,xi,’linear’)

線性插值（缺?。﹜i=interp1(x,y,xi,’nearst’)

最近鄰點(diǎn)插值yi=interp1(x,y,xi,’pchip’)分段三次Hermite插值yi=interp1(x,y,xi,’cubic’)同上yi=interp1(x,y,xi,’spline’)分段三次樣條插值（建議使用三次樣條插值）若已知的數(shù)據(jù)點(diǎn)x是等距節(jié)點(diǎn)，插值速度可提高，此時，填寫最后一個插值方法時，以星號引導(dǎo)，即為‘*linear’,’*spline’,’*cubic’,’*pchip’

,’*nearst’

第二十二頁，共六十頁，2022年，8月28日例：假設(shè)已知的數(shù)據(jù)點(diǎn)來自函數(shù)試根據(jù)生成的數(shù)據(jù)進(jìn)行插值,得到較光滑的曲線解：>>x=0:0.12:1;>>y=(x.^2-3*x+5).*exp(-5*x).*sin(x);>>plot(x,y,'-o')第二十三頁，共六十頁，2022年，8月28日可看出，由這樣的數(shù)據(jù)直接連線繪制出來的曲線十分粗糙，可再選擇一組插值點(diǎn)，然后用interp1插值第二十四頁，共六十頁，2022年，8月28日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:')第二十五頁，共六十頁，2022年，8月28日

interp1默認(rèn)的linear插值得到的曲線和plot畫出來的曲線一樣粗糙,而’nearest’選項得到的插值效果就更差了。而采用’cubic’和’spline’選項得到的插值更接近與理論值.第二十六頁，共六十頁，2022年，8月28日可看出：

interp1默認(rèn)的線性插值得到的曲線和plot畫出來的曲線一樣粗糙，因為該方法是對各個點(diǎn)直線連接。而’nearest’選項得到的插值效果就更差了。而采用’cubic’和’spline’選項得到的插值更接近與理論值。事實上，采用樣條插值算法得到的插值十分逼近理論值，甚至用肉眼難以分辨。第二十七頁，共六十頁，2022年，8月28日二、用冪級數(shù)做多項式插值給定n+1個數(shù)據(jù)點(diǎn)：過n+1個點(diǎn)的n

階多項式可寫為冪級數(shù)形式：注：過n＋1個數(shù)據(jù)點(diǎn)的n階插值多項式是唯一的。第二十八頁，共六十頁，2022年，8月28日對n＋1個數(shù)據(jù)點(diǎn)，設(shè)，則得到n+1個線性方程，可以表示為矩陣形式求解該方程組可確定系數(shù)（或用polyfit(x,y,n)確定）第二十九頁，共六十頁，2022年，8月28日例求用下列數(shù)據(jù)點(diǎn)擬合出的多項式的系數(shù)，并求當(dāng)x取2.101、4.234時y的值，并畫出已知數(shù)據(jù)點(diǎn)和擬合曲線。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第三十頁，共六十頁，2022年，8月28日第三十一頁，共六十頁，2022年，8月28日三、Lagrange插值多項式1、給定n+1個數(shù)據(jù)點(diǎn)：插值基函數(shù)：滿足插值條件：的次數(shù)不超過n的lagrange多項式g(x)為：第三十二頁，共六十頁，2022年，8月28日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）:調(diào)用格式：

y=lagrange1(x0,y0,x)第三十三頁，共六十頁，2022年，8月28日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)例：寫出擬合下面四個數(shù)據(jù)點(diǎn)的Lagrange插值公式，并計算x＝2.101、4.234時y

的值。x0y01.13.8872.34.2763.94.6515.12.117結(jié)果為y=4.14574.3007第三十四頁，共六十頁，2022年，8月28日由于lagrange插值基函數(shù)li(x)在節(jié)點(diǎn)xk(k=1,2,…,n)處滿足又由于過n+1個數(shù)據(jù)點(diǎn)的n次插值多項式是唯一的,故li(x)也可看做是擬合下列數(shù)據(jù)點(diǎn)的n次多項式可編寫程序求出

2、MATLAB程序?qū)崿F(xiàn)（法2）:進(jìn)而求出lagrange插值公式第三十五頁，共六十頁，2022年，8月28日functionp=lag_base(x)np=length(x);forj=1:npy=zeros(1,np);y(j)=1;p(j,:)=polyfit(x,y,np-1);end其調(diào)用格式為：p=lag_base(x)其中:

x是數(shù)據(jù)點(diǎn)的橫坐標(biāo)數(shù)組，p是一個矩陣，它的第i行為li(x)

的冪系數(shù)。先編寫程序求出第三十六頁，共六十頁，2022年，8月28日如前例也可求解如下：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)結(jié)果為：yi＝

4.14574.3007第三十七頁，共六十頁，2022年，8月28日3、Lagrange插值公式的微分插值公式微分為計算Lagrange插值多項式的一階導(dǎo)數(shù)，可用polyder函數(shù)將p的每一行轉(zhuǎn)換為一階導(dǎo)數(shù)的系數(shù)數(shù)組。第三十八頁，共六十頁，2022年，8月28日例：求出擬合下面四個數(shù)據(jù)點(diǎn)的Lagrange插值多項式在x＝2.101、4.234處的一階導(dǎo)數(shù)值。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)結(jié)果為：yi＝

0.6292-1.4004第三十九頁，共六十頁，2022年，8月28日例考慮一個著名的例子假設(shè)已知其中一些點(diǎn)的坐標(biāo)，則可以采用下面的命令進(jìn)行Lagrange插值，得出相應(yīng)的插值曲線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現(xiàn)象第四十頁，共六十頁，2022年，8月28日可看到：由lagrange插值得出的效果和精確值相差甚遠(yuǎn)，這種多項式階次越高越發(fā)散的現(xiàn)象稱為Runge現(xiàn)象所以對此例來說，傳統(tǒng)的lagrange算法失效第四十一頁，共六十頁，2022年，8月28日下面考慮用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,'--')第四十二頁，共六十頁，2022年，8月28日interp1函數(shù)效果可見：Matlab中提供的算法不存在Runge現(xiàn)象，可大膽放心使用第四十三頁，共六十頁，2022年，8月28日總結(jié)：（1）盡可能在小區(qū)間上使用多項式插值；（2）只能在一定范圍內(nèi)依靠增加插值點(diǎn)個數(shù)提高插值精度，如果插值點(diǎn)個數(shù)過多往往會適得其反。第四十四頁，共六十頁，2022年，8月28日定義:

若已知函數(shù)f(x)在點(diǎn)x1,

x2,···,xn+1

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

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

)四、牛頓插值（均差插值）第四十五頁，共六十頁，2022年，8月28日x -2 -1 0 1 3y -56 -16 -2 -2 4例由函數(shù)表求各階均差解:按公式計算各階差商如下x f(x)一階二階三階四階差商差商差商差商-2 -56

-1 -16

40

0 -214 -13

1 -20 -7 23 4 3 1 20第四十六頁，共六十頁，2022年，8月28日程序：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第四十七頁，共六十頁，2022年，8月28日有了均差，就能得到牛頓插值公式（均差插值公式）：第四十八頁，共六十頁，2022年，8月28日以3階均差插值公式為例，改寫形式第四十九頁，共六十頁，2022年，8月28日function[A,C,N]=newton_poly(x,y)%x為已知點(diǎn)的橫坐標(biāo)，y為已知點(diǎn)的縱坐標(biāo)%A為均差矩陣，C為牛頓插值多項式的系數(shù)向量%N為牛頓插值多項式%先產(chǎ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算法：第五十頁，共六十頁，2022年，8月28日%下面求牛頓插值多項式系數(shù)向量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);第五十一頁，共六十頁，2022年，8月28日例：給出節(jié)點(diǎn)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)第五十二頁，共六十頁，2022年，8月28日結(jié)果：A=17.0300000007.2400-8.513000001.0500-6.12871.10390002.03000.97033.51440.76040017.060014.88126.88661.11290.08430