matlab實現插值法和曲線擬合

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

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

3、,1,n) ;(2) 在每個區(qū)間xk ,xk+1 上,Ih (x)是個一次函數。易知,Ih(x)是個折線函數, 在每個區(qū)間xk ,xk+1 上,(k=0,1,n) ， 于是, Ih (x)在a,b上是連續(xù)的，但其一階導數是不連續(xù)的。3拉格朗日插值多項式算法輸入，令。對，計算4分段線性插值算法輸入（x,y）,k=0,1,n;計算5插值法和分段線性插值程序按下列數據分別作五次插值和分段線性插值，畫出兩條插值曲線以及給定數據點。求x1=0.32, x2=0.55, x3=0.68 時的函數近似值，并比較兩種方法的插值余項。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運算結果拉格朗日插值結果x=0.,y=1.x=0.,y=1.x=0.,y=1.拉格朗日插值余項：分段插值結果ans =1.0508 1.1418 1.2095分段線性插值余項：由于拉格朗日插值的余項比分段線性插值的余項要求更為嚴格，點少、區(qū)間小的時候，拉格朗日插值要更好。但在區(qū)間較大、節(jié)點較多的時候，

8、分段線性插值要更好。二、牛頓前插1牛頓前插原理次牛頓前插公式：插值余項：，階差分記作。 階差商是差分和差商之間的關系是2牛頓前插算法輸入。對，計算各階差分計算函數值3牛頓前插程序：編寫一個用牛頓前插公式計算函數值的程序，要求先輸出差分表，再計算x點的函數值0.1250.2500.3750.5000.6250.7500.7960.7730.7440.7040.6560.602分別求x=0.158和x=0.636的三次插值的值，并比較二者的插值余項。這里以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運行結果：差分表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時ans =0.4000 x=0,158時ans =0.0000三、曲線擬合1曲線擬合原理：給定數據。記擬合函數的形式為（1.1），其中為已知的線性無關函數。求系數使得（1.2）取最小值。稱（1.3）為擬合函數或經驗公式。如果，則（1.3）為次最小二乘擬合多項式

11、2曲線擬合算法：已知數據對，求多項式，使得為最小。注意到此時，多項式系數滿足下面的線性方程組：其中，然后只要調用線性方程組的函數程序即可3曲線擬合程序：試分別用拋物線y=a+bx2和指數曲線y=aebx擬合下列數據12.53.543.81.5026.033.0畫出數據點和兩條擬合曲線，并通過計算2個擬合函數殘差向量的2范數來比較擬合優(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指數曲線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(