數(shù)學(xué)建模

插值拉格朗日插值分段線性插值三次樣條插值一維插值一、插值的定義二、插值的方法三、用Matlab解插值問題返回返回二維插值一、二維插值定義二、網(wǎng)格節(jié)點插值法三、用Matlab解插值問題最鄰近插值分片線性插值雙線性插值網(wǎng)格節(jié)點數(shù)據(jù)的插值散點數(shù)據(jù)的插值一維插值的定義已知n+1個節(jié)點其中互不相同，不妨設(shè)求任一插值點處的插值節(jié)點可視為由產(chǎn)生,，表達式復(fù)雜,，或無封閉形式,，或未知.。

構(gòu)造一個(相對簡單的)函數(shù)通過全部節(jié)點,即再用計算插值，即返回

稱為拉格朗日插值基函數(shù)。

已知函數(shù)f(x)在n+1個點x0,x1,…,xn處的函數(shù)值為y0,y1,…,yn

。求一n次多項式函數(shù)Pn(x)，使其滿足：

Pn(xi)=yi,i=0,1,…,n.

解決此問題的拉格朗日插值多項式公式如下其中Li(x)為n次多項式：拉格朗日(Lagrange)插值拉格朗日(Lagrange)插值特別地:兩點一次(線性)插值多項式:三點二次(拋物)插值多項式:用Matlab作Lagrange插值

Matlab中沒有現(xiàn)成lagrange插值函數(shù)，必須編寫一個M文件實現(xiàn)lagrange插值。設(shè)n個節(jié)點數(shù)據(jù)以數(shù)組x0,y0輸入，m個插值點以數(shù)組x輸入，輸出數(shù)組y為m個插值。編寫一個名為lagrange.m

的M件：程序如下：functiony=lagrange(x0,y0,x)n=length(x0);m=length(x);fori=1:mz=x(i);s=0.0;fork=1:np=1.0;forj=1:nifj~=kp=p*(z-x0(j))/(x0(k)-x0(j));endends=p*y0(k)+s;endy(i)=s;end

拉格朗日多項式插值的這種振蕩現(xiàn)象叫Runge現(xiàn)象

采用拉格朗日多項式插值：選取不同插值節(jié)點個數(shù)n+1，其中n為插值多項式的次數(shù)，當n分別取2,4,6,8,10時，繪出插值結(jié)果圖形.例返回分段線性插值

分段線性插值通常有較好的收斂性和穩(wěn)定性，算法簡單，能夠克服拉格朗日多項式插值中出現(xiàn)的龍格現(xiàn)象，但其插值函數(shù)不如拉格朗日多項式光滑。分段線性插值n越大，誤差越小.xjxj-1xj+1x0xnxoy返回例用分段線性插值法求插值,并觀察插值誤差.1.在[-6,6]中平均選取5個點作插值(xch11)4.在[-6,6]中平均選取41個點作插值(xch14)2.在[-6,6]中平均選取11個點作插值(xch12)3.在[-6,6]中平均選取21個點作插值(xch13)

樣條插值

分段線性插值函數(shù)在結(jié)點的一階導(dǎo)數(shù)一般不存在，光滑性不高，這就導(dǎo)致了樣條插值的出現(xiàn)。主要應(yīng)用于工程技術(shù)上問題上，如飛機的機翼外形，內(nèi)燃機的進、排氣門的凸輪線。比分段線性插值更光滑。xyxi-1xiab

在數(shù)學(xué)上，光滑程度的定量描述是：函數(shù)(曲線)的k階導(dǎo)數(shù)存在且連續(xù)，則稱該曲線具有k階光滑性。光滑性的階次越高，則越光滑。是否存在較低次的分段多項式達到較高階光滑性的方法？三次樣條插值就是一個很好的例子。三次樣條插值

三次樣條插值g(x)為被插值函數(shù)。例用三次樣條插值選取11個基點計算插值(ych)返回用MATLAB作插值計算一維插值函數(shù)：yi=interp1(x，y，xi，'method')插值方法被插值點插值節(jié)點xi處的插值結(jié)果‘nearest’

：最鄰近插值‘linear’

：線性插值；‘spline’

：三次樣條插值；‘cubic’

：立方插值。缺省時：分段線性插值。注意：所有的插值方法都要求x是單調(diào)的，并且xi不能夠超過x的范圍。

例：在1-12的11小時內(nèi)，每隔1小時測量一次溫度，測得的溫度依次為：5，8，9，15，25，29，31，30，22，25，27，24。試估計每隔1/10小時的溫度值。hours=1:12;temps=[589152529313022252724];h=1:0.1:12;t=interp1(hours,temps,h,'spline');(直接輸出數(shù)據(jù)將是很多的)plot(hours,temps,'+‘)holdonplot(h,t,hours,temps,'r:')%作圖xlabel('Hour'),ylabel('DegreesCelsius’)xy機翼下輪廓線例已知飛機下輪廓線上數(shù)據(jù)如下，求x每改變0.1時的y值。返回二維插值的定義xyO第一種（網(wǎng)格節(jié)點）：

已知mn個節(jié)點其中互不相同，不妨設(shè)

構(gòu)造一個二元函數(shù)通過全部已知節(jié)點,即再用計算插值，即第二種（散亂節(jié)點）：yx0已知n個節(jié)點其中互不相同，

構(gòu)造一個二元函數(shù)通過全部已知節(jié)點,即再用計算插值，即返回3種二維插值方法：（1）、最鄰近插值‘nearest’

（2）、雙線性插值‘linear’

（3）、雙三次插值‘cubic’

要求x0,y0單調(diào)；x，y可取為矩陣，或x取行向量，y取為列向量，x,y的值分別不能超出x0,y0的范圍。z=interp2(x0,y0,z0,x,y,’method’)被插值點插值方法用MATLAB作網(wǎng)格節(jié)點數(shù)據(jù)的插值插值節(jié)點被插值點的函數(shù)值‘nearest’

最鄰近插值‘linear’

雙線性插值‘cubic’

雙三次插值缺省時,雙線性插值例：測得平板表面3*5網(wǎng)格點處的溫度分別為：828180828479636165818484828586試作出平板表面的溫度分布曲面z=f(x,y)的圖形。輸入以下命令：x=1:5;y=1:3;temps=[8281808284;7963616581;8484828586];mesh(x,y,temps)1.先在三維坐標畫出原始數(shù)據(jù)，畫出粗糙的溫度分布曲圖.2．以平滑數(shù)據(jù),在x、y方向上每隔0.2個單位的地方進行插值.再輸入以下命令:xi=1:0.2:5;yi=1:0.2:3;zi=interp2(x,y,temps,xi',yi,'cubic');mesh(xi,yi,zi)畫出插值后的溫度分布曲面圖.

通過此例對最近鄰點插值、雙線性插值方法和雙三次插值方法的插值效果進行比較。返回

插值函數(shù)griddata格式為:

cz

=griddata（x，y，z，cx，cy，‘method’）用MATLAB作散點數(shù)據(jù)的插值計算

要求cx取行向量，cy取為列向量。被插值點插值方法插值節(jié)點被插值點的函數(shù)值‘nearest’

最鄰近插值‘linear’

雙線性插值‘cubic’

雙三次插值'v4'-Matlab提供的插值方法缺省時,雙線性插值

例在某海域測得一些點(x,y)處的水深z由下表給出，船的吃水深度為5英尺，在矩形區(qū)域（75，200）*（-50，150）里的哪些地方船要避免進入。返回4.作出水深小于5的海域范圍,即z=5的等高線.數(shù)學(xué)建模

擬合擬合2.擬合的基本原理1.擬合問題引例擬合問題引例1溫度t(0C)20.532.751.073.095.7電阻R()7658268739421032已知熱敏電阻數(shù)據(jù)：求600C時的電阻R。

設(shè)

R=at+ba,b為待定系數(shù)擬合與插值的關(guān)系

函數(shù)插值與曲線擬合都是要根據(jù)一組數(shù)據(jù)構(gòu)造一個函數(shù)作為近似，由于近似的要求不同，二者的數(shù)學(xué)方法上是完全不同的。

實例：下面數(shù)據(jù)是某次實驗所得，希望得到X和f之間的關(guān)系？問題：給定一批數(shù)據(jù)點，需確定滿足特定要求的曲線或曲面解決方案：若不要求曲線（面）通過所有數(shù)據(jù)點，而是要求它反映對象整體的變化趨勢，這就是數(shù)據(jù)擬合，又稱曲線擬合或曲面擬合。若要求所求曲線（面）通過所給所有數(shù)據(jù)點，就是插值問題；最臨近插值、線性插值、樣條插值與曲線擬合結(jié)果：曲線擬合問題最常用的解法——線性最小二乘法的基本思路第一步:先選定一組函數(shù)

r1(x),r2(x),…rm(x),m<n,

令

f(x)=a1r1(x)+a2r2(x)+…+amrm(x)（1）其中

a1,a2,…am

為待定系數(shù)。

第二步:確定a1,a2,…am

的準則（最小二乘準則）：使n個點（xi,yi)與曲線y=f(x)的距離i的平方和最小

。記

問題歸結(jié)為，求

a1,a2,…am

使

J(a1,a2,…am)

最小。線性最小二乘法的求解：預(yù)備知識超定方程組：方程個數(shù)大于未知量個數(shù)的方程組即Ra=y其中超定方程一般是不存在解的矛盾方程組。

如果有向量a使得達到最小，則稱a為上述超定方程的最小二乘解。常用多項式擬合的一般方法Step1：根據(jù)所給數(shù)據(jù)畫出散點圖，確定擬合多項式的次數(shù)n;Step2：利用matlab

多項式擬合函數(shù)polyfit進行擬合，求出多項式系數(shù)矩陣A，通過函數(shù)polyval可求出擬合函數(shù)f(x)在x=x0處的函數(shù)值y。線性最小二乘擬合f(x)=a1r1(x)+…+amrm(x)中函數(shù){r1(x),…rm(x)}的選取

1.通過機理分析建立數(shù)學(xué)模型來確定f(x)；++++++++++++++++++++++++++++++f=a1+a2xf=a1+a2x+a3x2f=a1+a2x+a3x2f=a1+a2/xf=aebxf=ae-bx2.將數(shù)據(jù)(xi,yi)i=1,…n作圖，通過直觀判斷確定f(x)：用MATLAB解擬合問題1、線性最小二乘擬合2、非線性最小二乘擬合用MATLAB作線性最小二乘擬合1.作多項式f(x)=a1xm+…+amx+am+1擬合,可利用已有程序:a=polyfit(x,y,m)2.對超定方程組可得最小二乘意義下的解。，用3.多項式在x處的值y可用以下命令計算：

y=polyval（a，x）輸出擬合多項式系數(shù)a=[a1,…am,

am+1](數(shù)組)）輸入同長度的數(shù)組X，Y擬合多項式次數(shù)即要求出二次多項式:中的使得:例對下面一組數(shù)據(jù)作二次多項式擬合1）輸入以下命令：

x=0:0.1:1;y=[-0.4471.9783.286.167.087.347.669.569.489.3011.2];A=polyfit(x,y,2)z=polyval(A,x);

plot(x,y,'k+',x,z,'r')%作出數(shù)據(jù)點和擬合曲線的圖形2）計算結(jié)果：Ａ=-9.810820.1293-0.0317用matlab多項式擬合的命令1.lsqcurvefit已知數(shù)據(jù)點：xdata=（xdata1，xdata2，…，xdatan），

ydata=（ydata1，ydata2，…，ydatan）

用MATLAB作非線性最小二乘擬合

Matlab的提供了兩個求非線性最小二乘擬合的函數(shù)：lsqcurvefit和lsqnonlin。兩個命令都要先建立M-文件fun.m，在其中定義函數(shù)f(x)，但兩者定義f(x)的方式是不同的,可參考例題.

lsqcurvefit用以求含參量x（向量）的向量值函數(shù)F(x,xdata)=（F（x，xdata1），…，F(xiàn)（x，xdatan））T中的參變量x(向量),使得

輸入格式為:(1)x=lsqcurvefit(‘fun’,x0,xdata,ydata);(2)x=lsqcurvefit(‘fun’,x0,xdata,ydata,options);(3)x=lsqcurvefit(‘fun’,x0,xdata,ydata,options,’grad’);(4)[x,options]=lsqcurvefit(‘fun’,x0,xdata,ydata,…);(5)[x,options,funval]=lsqcurvefit(‘fun’,x0,xdata,ydata,…);(6)[x,options,funval,Jacob]=lsqcurvefit(‘fun’,x0,xdata,ydata,…);fun是一個事先建立的定義函數(shù)F(x,xdata)

的M-文件,自變量為x和xdata說明：x=lsqcurvefit(‘fun’,x0,xdata,ydata,options);迭代初值已知數(shù)據(jù)點選項見無約束優(yōu)化

lsqnonlin用以求含參量x（向量）的向量值函數(shù)

f(x)=(f1(x),f2(x),…,fn(x))T

中的參量x，使得

最小。其中fi（x）=f（x，xdatai，ydatai）

=F(x,xdatai)-ydatai

2.lsqnonlin已知數(shù)據(jù)點：xdata=（xdata1，xdata2，…，xdatan）

ydata=（ydata1，ydata2，…，ydatan）輸入格式為：

1）x=lsqnonlin（‘fun’，x0）；

2）x=lsqnonlin

（‘fun’，x0，options）；

3）x=lsqnonlin

（‘fun’，x0，options，‘grad’）；

4）[x，options]=lsqnonlin

（‘fun’，x0，…）；

5）[x，options，funval]=lsqnonlin

（‘fun’，x0，…）；說明：x=lsqnonlin

（‘fun’，x0，options）；fun是一個事先建立的定義函數(shù)f(x)的M-文件，自變量為x迭代初值選項見無約束優(yōu)化

例2用下面一組數(shù)據(jù)擬合中的參數(shù)a，b，k該問題即解最優(yōu)化問題：

1）編寫M-文件curvefun1.m

functionf=curvefun1(x,tdata)f=x(1)+x(2)*exp(-0.02*x(3)*tdata)%其中x(1)=a;x(2)=b；x(3)=k;2）輸入命令tdata=100:100:1000cdata=1e-03*[4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59];x0=[0.2,0.05,0.05];

x=lsqcurvefit('curvefun1',x0,tdata,cdata)f=curvefun1(x,tdata)

F(x，tdata)=，x=(a