matlab實(shí)現(xiàn)插值法和曲線擬合

1、精選優(yōu)質(zhì)文檔-傾情為你奉上插值法和曲線擬合電子科技大學(xué) 摘要：理解拉格朗日多項(xiàng)式插值、分段線性插值、牛頓前插，曲線擬合，用matlab編程求解函數(shù)，用插值法和分段線性插值求解同一函數(shù)，比較插值余項(xiàng)；用牛頓前插公式計(jì)算函數(shù)，計(jì)算函數(shù)值；對(duì)于曲線擬合，用不同曲線擬合數(shù)據(jù)。關(guān)鍵字：拉格朗日插值多項(xiàng)式；分段線性插值；牛頓前插；曲線擬合引言：在數(shù)學(xué)物理方程中，當(dāng)給定數(shù)據(jù)是不同散點(diǎn)時(shí)，無法確定函數(shù)表達(dá)式，求解函數(shù)就需要很大的計(jì)算量，我們有多種方法對(duì)給定的表格函數(shù)進(jìn)行求解，我們這里，利用插值法和曲線擬合對(duì)函數(shù)進(jìn)行求解，進(jìn)一步了解函數(shù)性質(zhì)，兩種方法各有利弊，適合我們進(jìn)行不同的散點(diǎn)函數(shù)求解。正文：一、插值法和分

2、段線性插值1拉格朗日多項(xiàng)式原理對(duì)某個(gè)多項(xiàng)式，已知有給定的k + 1個(gè)取值點(diǎn)：其中對(duì)應(yīng)著的位置，而對(duì)應(yīng)著函數(shù)在這個(gè)位置的取值。假設(shè)任意兩個(gè)不同的xj都互不相同，那么應(yīng)用拉格朗日插值公式所得到的拉格朗日插值多項(xiàng)式為：其中每個(gè)為拉格朗日基本多項(xiàng)式（或稱插值基函數(shù)），其表達(dá)式為：拉格朗日基本多項(xiàng)式的特點(diǎn)是在 上取值為1，在其它的點(diǎn) 上取值為0。2分段線性插值原理給定區(qū)間a,b, 將其分割成a=x0 <x1 <<xn =b, 已知函數(shù)y= f(x) 在這些插值結(jié)點(diǎn)的函數(shù)值為yk =f(xk)(k=0,1,n)求一個(gè)分段函數(shù)Ih(x), 使其滿足：(1) Ih(xk )=yk ，(k=0

3、,1,n) ;(2) 在每個(gè)區(qū)間xk ,xk+1 上,Ih (x)是個(gè)一次函數(shù)。易知,Ih(x)是個(gè)折線函數(shù), 在每個(gè)區(qū)間xk ,xk+1 上,(k=0,1,n) ， 于是, Ih (x)在a,b上是連續(xù)的，但其一階導(dǎo)數(shù)是不連續(xù)的。3拉格朗日插值多項(xiàng)式算法輸入，令。對(duì)，計(jì)算4分段線性插值算法輸入（x,y）,k=0,1,n;計(jì)算5插值法和分段線性插值程序按下列數(shù)據(jù)分別作五次插值和分段線性插值，畫出兩條插值曲線以及給定數(shù)據(jù)點(diǎn)。求x1=0.32, x2=0.55, x3=0.68 時(shí)的函數(shù)近似值，并比較兩種方法的插值余項(xiàng)。0.300.420.500.580.660.721.044031.084621

4、.118031.156031.198171,23223拉格朗日插值程序：function lagrintxi=0.32,0.55,0.68;%xi=0.2:0.001:0.8;x=0.3,0.42,0.50,0.58,0.66,0.72;y=1.04403,1.08462,1.11803,1.15603,1.19817,1.23223;L=zeros(size(y);m=length(xi);for i=1:m dxi=xi(i)-x; L(1)=prod(dxi(2:6)/prod(x(1)-x(2:6); L(6)=prod(dxi(1:6-1)/prod(x(6)-x(1:6-1); f

5、or j=2:6-1 num=prod(dxi(1:j-1)*prod(dxi(j+1:6); den=prod(x(j)-x(1:j-1)*prod(x(j)-x(j+1:6); L(j)=num/den; end yi(i)=sum(y.*L); fprintf('x=%f,y=%fn',xi(i),yi(i); endplot(xi,yi,'r');axis(0.2 0.8 1.03 1.24);hold on plot(x,y,'b.','markersize',20)grid on分段線性插值算法程序：function

6、y=div %xi=0.3:0.001:0.72;x0=0.3,0.42,0.50,0.58,0.66,0.72;y0=1.04403,1.08462,1.11803,1.15603,1.19817,1.23223;k=1;xi=0.32,0.55,0.68;for j=1:3 for i=1:5 if xi(j)>=x0(i) && xi(j)<=x0(i+1) && k<=3 lx(1)=(xi(j)-x0(i+1)/(x0(i)-x0(i+1); lx(2)=(xi(j)-x0(i)/(x0(i+1)-x0(i); y(k)=lx(1)*

7、y0(i)+lx(2)*y0(i+1); k=k+1; end endendplot(xi,y,'r');axis(0.2 0.8 1.03 1.24);hold on plot(x0,y0,'b.','markersize',20)grid on6運(yùn)算結(jié)果拉格朗日插值結(jié)果x=0.,y=1.x=0.,y=1.x=0.,y=1.拉格朗日插值余項(xiàng)：分段插值結(jié)果ans =1.0508 1.1418 1.2095分段線性插值余項(xiàng)：由于拉格朗日插值的余項(xiàng)比分段線性插值的余項(xiàng)要求更為嚴(yán)格，點(diǎn)少、區(qū)間小的時(shí)候，拉格朗日插值要更好。但在區(qū)間較大、節(jié)點(diǎn)較多的時(shí)候，

8、分段線性插值要更好。二、牛頓前插1牛頓前插原理次牛頓前插公式：插值余項(xiàng)：，階差分記作。 階差商是差分和差商之間的關(guān)系是2牛頓前插算法輸入。對(duì)，計(jì)算各階差分計(jì)算函數(shù)值3牛頓前插程序：編寫一個(gè)用牛頓前插公式計(jì)算函數(shù)值的程序，要求先輸出差分表，再計(jì)算x點(diǎn)的函數(shù)值0.1250.2500.3750.5000.6250.7500.7960.7730.7440.7040.6560.602分別求x=0.158和x=0.636的三次插值的值，并比較二者的插值余項(xiàng)。這里以x=0.636為例function P=newtonchax0=0.636;X=0.125 0.250 0.375 0.500 0.625 0.

9、750;Y=0.796 0.773 0.744 0.704 0.656 0.602;h=abs(X(2)-X(1);n=find(abs(x0-X)<3*h);X=X(n(1):n(end);Y=Y(n(1):n(end);w=length(X);R=zeros(w,w);R(:,1)=Y(:);for k=2:w for j=k:w R(j,k)=R(j,k-1)-R(j-1,k-1); endendt=(x0-X(1)/h;T=1; for m=1:w-1 T=T*(t-m+1); N(m)=R(m+1,m+1)*T/factorial(m); endP=R(1,1)+sum(N);

10、4運(yùn)行結(jié)果：差分表0.0000000000.0000-0.000000000.0000-0.0000-0.00000000.0000-0.0000-0.0000-0.0000000.0000-0.0000-0.00000.00000.000000.0000-0.0000-0.00000.0000-0.0000-0.0000X=0.636時(shí)ans =0.4000 x=0,158時(shí)ans =0.0000三、曲線擬合1曲線擬合原理：給定數(shù)據(jù)。記擬合函數(shù)的形式為（1.1），其中為已知的線性無關(guān)函數(shù)。求系數(shù)使得（1.2）取最小值。稱（1.3）為擬合函數(shù)或經(jīng)驗(yàn)公式。如果，則（1.3）為次最小二乘擬合多項(xiàng)式

11、2曲線擬合算法：已知數(shù)據(jù)對(duì)，求多項(xiàng)式，使得為最小。注意到此時(shí)，多項(xiàng)式系數(shù)滿足下面的線性方程組：其中，然后只要調(diào)用線性方程組的函數(shù)程序即可3曲線擬合程序：試分別用拋物線y=a+bx2和指數(shù)曲線y=aebx擬合下列數(shù)據(jù)12.53.543.81.5026.033.0畫出數(shù)據(jù)點(diǎn)和兩條擬合曲線，并通過計(jì)算2個(gè)擬合函數(shù)殘差向量的2范數(shù)來比較擬合優(yōu)劣。用拋物線y=a+bx擬合程序：function ZXEx=1 2.52 3.52 42;y=3.8 1.50 26.0 33.0;m=1;S=zeros(1,2*m+1);T=zeros(m+1,1);for k=1:2*m+1 S(k)=sum(x.(k-1

12、);endfor k=1:m+1 T(k)=sum(x.(k-1).*y);endA=zeros(m+1,m+1);a=zeros(m+1,1);for i=1:m+1 for j=1:m+1 A(i,j)=S(i+j-1); endenda=AT;for k=1:m+1 fprintf('a%d=%fn',k,a(k);endp=polyfit(x,y,1);u=polyval(p,x);plot(sqrt(x),u,'b')hold onplot(sqrt(x),y,'b.')grid on指數(shù)曲線y=aebx擬合程序：function ZX

13、E2x=1 2.5 3.5 4;y=3.8 1.50 26.0 33.0;y=log(y);m=1;S=zeros(1,2*m+1);T=zeros(m+1,1);for k=1:2*m+1 S(k)=sum(x.(k-1);endfor k=1:m+1 T(k)=sum(x.(k-1).*y);endA=zeros(m+1,m+1);a=zeros(m+1,1);for i=1:m+1 for j=1:m+1 A(i,j)=S(i+j-1); endenda=AT;for k=1:m+1 fprintf('a%d=%fn',k,a(k);endp=polyfit(x,y,1);u=polyval(