數(shù)學(xué)建模插值擬合

1、插值與擬合數(shù)學(xué)建模插值求一個(gè)函數(shù)，經(jīng)過(guò)已知的數(shù)據(jù)點(diǎn)。常用的插值方法有：多項(xiàng)式插值，Hermite插值，分段線性插值，三次樣條插值，多項(xiàng)式插值（拉格朗日多項(xiàng)式插值，牛頓插值）：可找到一個(gè)唯一的n次多項(xiàng)式經(jīng)過(guò)已知的n+1個(gè)數(shù)據(jù)點(diǎn)；但高次多項(xiàng)式插值往往會(huì)產(chǎn)生嚴(yán)重的振蕩現(xiàn)象（龍格現(xiàn)象）常用分段低次插值代替。如分段線性插值，三次樣條插值等。Hermite插值要求在節(jié)點(diǎn)處與函數(shù)有相同的一階、二階甚至更高階的導(dǎo)數(shù)值，用Matlab插值y=interp1(x0,y0,x) :分段線性插值，x0，y0為節(jié)點(diǎn)數(shù)組，x為插值點(diǎn)數(shù)組，y為插值結(jié)果。y=interp1(x0,y0,x,spline): 三次樣條插值。y

2、=spline(x0,y0,x): 三次樣條插值。或： y=ppval(spline(x0,y0),x)y=pchip(x0,y0,x): 分段三次Hermite插值。 或：y=ppval(pchip(x0,y0,x)三次樣條插值和分段三次Hermite插值的區(qū)別：三次樣條插值更光滑，而分段三次Hermite插值振蕩少，更平坦。見 help pchip二維插值插值節(jié)點(diǎn)為網(wǎng)格節(jié)點(diǎn)z=interp2(x0,y0,z0,x,y,method) method有nearest，linear，cubic，splinepp=csape(x0,y0,z0,conds)，z=fnval(pp,x,y) cond

3、s是邊界條件例題在一丘陵地帶測(cè)量高程，x 和 y 方向每隔100米測(cè)一個(gè)點(diǎn)，得高程如下表，試插值一曲面，確定合適的模型，并由此找出最高點(diǎn)和該點(diǎn)的高程。Matlab求解二維網(wǎng)格節(jié)點(diǎn)插值編寫程序如下：x=100:100:500;y=100:100:400;z=636 697 624 478 450 698 712 630 478 420 680 674 598 412 400 662 626 552 334 310;pp=csape(x,y,z)xi=100:10:500;yi=100:10:400;cz=fnval(pp,xi,yi);i,j=find(cz=max(max(cz) %找最高點(diǎn)的

4、地址xh=xi(i),yh=yi(j),zmax=cz(i,j) %求最高點(diǎn)的坐標(biāo)xx,yy=meshgrid(xi,yi); x0,y0=meshgrid(x,y);mesh(xx, yy, cz), hold on, plot3(x0,y0,z,ko)二維插值插值節(jié)點(diǎn)為散亂節(jié)點(diǎn)Matlab 中提供了插值函數(shù) griddata，其格式為： z=griddata(x0,y0,z0,x,y,method) method有nearest，linear，cubic例題在某海域測(cè)得一些點(diǎn)（x,y）處的水深 z 由下表給出，在矩形區(qū)域（75,200）(-50,150) 內(nèi)畫出海底曲面的圖形。Matlab

5、求解二維網(wǎng)格節(jié)點(diǎn)插值編寫程序如下：clc, clearx=129,140,103.5,88,185.5,195,105,157.5,107.5,77,81,162,162,117.5;y=7.5,141.5,23,147,22.5,137.5,85.5,-6.5,-81,3,56.5,-66.5,84,-33.5;z=-4,8,6,8,6,8,8,9,9,8,8,9,4,9;xmm=minmax(x) %求x的最小值和最大值ymm=minmax(y) %求y的最小值和最大值xi=xmm(1):xmm(2);yi=ymm(1):ymm(2);zi1=griddata(x,y,z,xi,yi,cu

6、bic); %立方插值z(mì)i2=griddata(x,y,z,xi,yi,nearest); %最近點(diǎn)插值z(mì)i=zi1; %立方插值和最近點(diǎn)插值的混合插值的初始值z(mì)i(isnan(zi1)=zi2(isnan(zi1); %把立方插值中的不確定值換成最近點(diǎn)插值的結(jié)果subplot(1,2,1), plot(x,y,*)subplot(1,2,2), mesh(xi,yi,zi)已知一組（二維）數(shù)據(jù)，即平面上 n個(gè)點(diǎn)（xi, yi) i=1,n, 尋求一個(gè)函數(shù)（曲線）y=f(x), 使 f(x) 在某種準(zhǔn)則下與所有數(shù)據(jù)點(diǎn)最為接近，即曲線擬合得最好。 +xyy=f(x)(xi,yi)ii 為點(diǎn)（xi

7、, yi) 與曲線 y=f(x) 的距離曲線擬合第一步:先選定一組函數(shù) r1(x), r2(x), rm(x), mn, 令 f(x)=a1r1(x)+a2r2(x)+ +amrm(x) （1）其中 a1,a2, ,am 為待定系數(shù)。 第二步: 確定a1,a2, am 的準(zhǔn)則（最小二乘準(zhǔn)則）：使n個(gè)已知點(diǎn)（xi,yi) 與曲線 y=f(x) 的距離i 的平方和最小 。記 問題歸結(jié)為，求 a1,a2, am 使 J(a1,a2, am) 最小。線性最小二乘解超定方程組：方程個(gè)數(shù)大于未知量個(gè)數(shù)的方程組即 Ra=y其中超定方程一般是不存在解的矛盾方程組。 如果有向量a使得 達(dá)到最小，則稱a為上述超定

8、方程的最小二乘解。 線性最小二乘法的求解：預(yù)備知識(shí) 定理：當(dāng)RTR可逆時(shí)，超定方程組（3）存在最小二乘解，且即為方程組 RTRa=RTy的解：a=(RTR)-1RTy 所以，曲線擬合的最小二乘法要解決的問題，實(shí)際上就是求以下超定方程組的最小二乘解的問題。其中Ra=y （3）線性最小二乘法的求解線性最小二乘擬合中函數(shù)r1(x), rm(x)的選取 1. 通過(guò)機(jī)理分析建立數(shù)學(xué)模型來(lái)確定 f(x)；+f=a1+a2xf=a1+a2x+a3x2f=a1+a2x+a3x2f=a1+a2/xf=aebxf=ae-bx 2. 將數(shù)據(jù) (xi, yi) i=1, n 作圖，通過(guò)直觀判斷確定 f(x)：用MAT

9、LAB解擬合問題1、線性最小二乘擬合2、非線性最小二乘擬合1. 作多項(xiàng)式 f (x)=a1xm+ +amx+am+1擬合,可利用已有程序:a=polyfit(x, y, m)2. 對(duì)超定方程組可得最小二乘意義下的解。，用2. 多項(xiàng)式在 x 處的值 y 可用以下命令計(jì)算： y=polyval(a, x)輸出擬合多項(xiàng)式系數(shù)a=a1, ,am , am+1 (數(shù)組)輸入同長(zhǎng)度的數(shù)組X，Y擬合多項(xiàng)式次數(shù)用MATLAB作線性最小二乘擬合即要求 出二次多項(xiàng)式:中 的使得:例 對(duì)下面一組數(shù)據(jù)作二次多項(xiàng)式擬合1）輸入以下命令：x=0:0.1:1; y=-0.447 1.978 3.28 6.16 7.08 7

10、.34 7.66 9.56 9.48 9.30 11.2; R=(x.2) x ones(11,1); A=Ry解法1用解超定方程組的方法2）計(jì)算結(jié)果： = -9.8108 20.1293 -0.0317擬合效果1）輸入以下命令： x0=0:0.1:1; y0=-0.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2; A=polyfit(x0,y0,2) x=0:0.01:1; y=polyval(A,x); plot(x0,y0,k+,x,y,r) %作出數(shù)據(jù)點(diǎn)和擬合曲線的圖形2）計(jì)算結(jié)果： = -9.8108 20.1293 -0

11、.0317解法2用多項(xiàng)式擬合的命令1. lsqcurvefit已知數(shù)據(jù)點(diǎn)： xdata=（xdata1，xdata2，xdatan）， ydata=（ydata1，ydata2，ydatan） 用MATLAB作非線性最小二乘擬合 Matlab提供了兩個(gè)求非線性最小二乘擬合的函數(shù)：lsqcurvefit 和 lsqnonlin。兩個(gè)命令都要先建立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)

12、x = lsqcurvefit (fun,x0,xdata,ydata); (2) x =lsqcurvefit (fun,x0,xdata,ydata,Lb,Ub); (3) x = lsqcurvefit (fun,x0,xdata,ydata,Lb,Ub,options); fun是一個(gè)事先建立的定義函數(shù)F(x,xdata) 的M-文件, 自變量為x和xdata說(shuō)明：x = lsqcurvefit (fun,x0,xdata,ydata,Lb,Ub,options);輸入輸出格式迭代初值已知數(shù)據(jù)上下界控制選項(xiàng)2. lsqnonlin已知數(shù)據(jù)點(diǎn)： xdata=(xdata1，xdata2，

13、xdatan) ydata=( ydata1，ydata2，ydatan)其中 fi(x)=f(x，xdatai，ydatai） =F(x,xdatai)-ydatai 最小。lsqnonlin用以求含參量x（向量）的向量值函數(shù) f (x)=(f1(x),f2(x),fn(x)T 中的參量x，使得 1） x=lsqnonlin(fun, x0)； 2） x= lsqnonlin (fun, x0, Lb, Ub)； 3） x= lsqnonlin (fun, x0, Lb, Ub, options)； 說(shuō)明：x= lsqnonlin (fun, x0, Lb, Ub, options)；fun

14、是一個(gè)事先建立的定義函數(shù)f(x)的M-文件，自變量為x輸入輸出格式迭代初值上下界控制選項(xiàng)例題1用下面的語(yǔ)句生成一組數(shù)據(jù)，其中 ai 為待定系數(shù)x=0:.1:10; y=0.12*exp(-0.213*x)+0.54*exp(-0.17*x).*sin(1.23*x);顯然該組數(shù)據(jù)滿足函數(shù)采用最小二乘曲線擬合獲得這些待定系數(shù)解法1：用lsqcurvefit函數(shù)使用內(nèi)聯(lián)函數(shù)f=inline(a(1)*exp(-a(2)*x)+a(3)* exp(-a(4)*x).*sin(a(5)*x),a,x);xx,res=lsqcurvefit(f,1,1,1,1,1,x,y); xx, resOptimi

15、zation terminated successfully: Relative function value changing by less than OPTIONS.TolFunans = 0.1197 0.2125 0.5404 0.1702 1.2300res = 7.1637e-007修改最優(yōu)化選項(xiàng)ff=optimset; ff.TolFun=1e-20; ff.TolX=1e-15; % 修改精度限制xx,res=lsqcurvefit(f,1,1,1,1,1,x,y,ff); xx,resOptimization terminated successfully: Relativ

16、e function value changing by less than OPTIONS.TolFunans = 0.1200 0.2130 0.5400 0.1700 1.2300res = 9.5035e-0211）編寫M-文件 curvefun1.mfunction f=curvefun1(a, x)f=a(1)*exp(-a(2)*x)+a(3)* exp(-a(4)*x).*sin(a(5)*x); %其中 a(1)=a1; a(2)=a2; a(3)=a3; a(4)=a4; a(5)=a5;2）輸入命令x=0:.1:10;y=0.12*exp(-0.213*x)+0.54*e

17、xp(-0.17*x).*sin(1.23*x); a0=1, 1, 1, 1, 1; a, res=lsqcurvefit (curvefun1, a0, x, y)例題2用最小二乘法擬合 中的未知參數(shù) ，其中已知數(shù)據(jù)值xi, yi分別放在數(shù)據(jù)文件data2.mat中的x0和y0中Matlab程序如下：clc, clearload data2 %分別加載xi的觀測(cè)值x0，yi的觀測(cè)值y0mf=(cs,xdata)1/sqrt(2*pi)/cs(2)*exp(-(xdata-cs(1).2/cs(2)2/2);% yc=mf(2,1,1) %測(cè)試匿名函數(shù)cs=lsqcurvefit(mf,rand(2,1),x0,y0) %擬合參數(shù)的初始值是任意取的data2.m由如下程序產(chǎn)生：x0=-10:0.01:10;y0=normpdf(x0,2,5); %計(jì)算正態(tài)分布密度函數(shù)在