1、MatLab &數(shù)學(xué)建模授課：唐靜波（九江學(xué)院理學(xué)院）第五講數(shù)值計(jì)算（二）r=一、線性優(yōu)化min CTx, x g Rns.t.Ax bvlb x vub用命令 x=lp(C,A,b,vlb,vub)例最小值線性優(yōu)化f(x)=-5x -4x -6x 123x-x +x W20 1233x +2x +4x W 423x +2x W30(0 W x0 W x2,0 W x3)First, enter the coefficients:f = -5; -4; -6A = 1 -1 13 2 43 2 0;b = 20; 42; 30;lb = 0,0,0;% x 的最小值0,0,0ub = inf,

2、inf,inf;Next, call a linear programming routine: x= lp(f,A,b,lb,ub);Entering x x =0.000015.00003.0000實(shí)際此命令改為：x = linprog(f,A,b,Aeq,beq)x = linprog(f,A,b,Aeq,beq,lb,ub) 對(duì)以上的問題可做如下的操作： First, enter the coefficients: f = -5; -4; -6;A = 1 -1 13 2 43 2 0;b = 20; 42; 30;lb = zeros(3,1);Next, call a linear

3、 programming routine:x,fval,exitflag,output,lambda = linprog(f,A,b,lb)； x =0.000015.00003.0000fval =-78.0000exitflag =1output =iterations: 6cgiterations: 0 algorithm: lipsollambda =ineqlin: 3x1 double eqlin: 0 x1 double upper: 3x1 double lower: 3x1 double例線性優(yōu)化Min -400 x1-1000 x2-300 x3+200 x4-2x2 +

4、x3 + x02x1 +3x2=163x1 +4x2=0; x3=5c=-400,-1000,-300,200;% 目標(biāo)函數(shù)系數(shù)A=0 -2 1 1; 2 3 0 0; 3 4 0 0; % 約束條件系數(shù)b=0; 16; 24;xLB=0,0,0,0;% x取值范圍的最小值xUB=inf,inf,5,inf;% x取值范圍的最大值x0=0,0,0,0;% x取迭代初始值nEq=1;%約束條件中只有一個(gè)二號(hào)，其余為二x=lp(c,A,b,xLB,xUB,x0,nEq)disp(最優(yōu)值為：，num2str(c*x)結(jié)果：x =3.44833.03455.00001.0690最優(yōu)值為-5700二、非

5、線性優(yōu)化min f (x), x g Rns.t.g(x) 0用命令 x=constr(f ,x0)。例最小值非線性優(yōu)化Min f(x)=-x1x2x3,-x1-2x2-2x3 W0,x1+2x2+2x3 W72,初值：x = 10; 10; 10第一步：編寫M文件myfun.m function f,g=myfun(x) f=-x(1)*x(2)*x(3);g(1)=-x(1)-2*x(2)-2*x(3);g(2)=x(1)+2*x(2)+2*x(3)-72;第二步：求解在MATLAB工作窗中鍵入 x0=10,10,10;x=constr(myfun,x0)艮口可.x =24.000012.

6、000012.0000例非線性優(yōu)化Min f(x)=-x1x2(x1+ x2)x3=0; x3=2;第一步：編寫M文件fxxgh.mfunction F,G=fxxgh(x)F=-x(1)*x(2);G(1) 二 (x(1)+x(2)*x(3)-120;第二步：求解在MATLAB工作窗中鍵入x=1,1,1；%x取迭代初始值options(13)=0;%約束條件中有0個(gè)二XL=0,0,2;%x取值范圍的最小值XU=inf;inf;inf;%x取值范圍的最大值x,options=constr(fxxgh,x,options,XL,XU)options(8)%輸出最小值x號(hào)，其余為二ans = -9

7、00.0000 x =30.000030.00002.0000三、曲線擬合與插值曲線擬合和插值函數(shù)polyfit(x, y, n)interp1(x, y, xo)interp1(x, y, xo, spline )interp1(x, y, xo, cubic )interp2(x, y, Z, xi, yi)interp2(x, y, Z, xi, yi, cubic ) interp2(x, y, Z, xi, yi, nearest )對(duì)描述n階多項(xiàng)式y(tǒng)=f(x)的數(shù)據(jù)進(jìn)行最小二乘曲線擬合1維線性插值1維3次樣條插值1維3次插值2維線性插值2維3次插值2維最近鄰插值在大量的應(yīng)用領(lǐng)域中，

8、人們經(jīng)常面臨用一個(gè)解析函數(shù)描述數(shù)撼通常是測量值)的任務(wù)。對(duì)這 個(gè)問題有兩種方法。在插值法里，數(shù)據(jù)假定是正確的，要求以某種方法描述數(shù)據(jù)點(diǎn)之間所發(fā) 生的情況。曲線擬合或回歸是人們?cè)O(shè)法找出某條光滑曲線，它最佳地?cái)M合數(shù)據(jù)，但不必要經(jīng) 過任何數(shù)據(jù)點(diǎn)。圖1說明了這兩種方法。標(biāo)有0的是數(shù)據(jù)點(diǎn)；連接數(shù)據(jù)點(diǎn)的實(shí)線描繪了線 性內(nèi)插，虛線是數(shù)據(jù)的最佳擬合。曲線擬合曲線擬合涉及回答兩個(gè)基本問題:最佳擬合意味著什么？應(yīng)該用什么樣的曲線？可用許 多不同的方法定義最佳擬合，并存在無窮數(shù)目的曲線。所以，從這里開始，我們走向何方？ 正如它證實(shí)的那樣，當(dāng)最佳擬合被解釋為在數(shù)據(jù)點(diǎn)的最小誤差平方和，且所用的曲線限定為 多項(xiàng)式時(shí)，那么

9、曲線擬合是相當(dāng)簡捷的。數(shù)學(xué)上，稱為多項(xiàng)式的最小二乘曲線擬合。如果這 種描述使你混淆，再研究圖1。虛線和標(biāo)志的數(shù)據(jù)點(diǎn)之間的垂直距離是在該點(diǎn)的誤差。對(duì)各 數(shù)據(jù)點(diǎn)距離求平方，并把平方距離全加起來，就是誤差平方和。這條虛線是使誤差平方和盡 可能小的曲線，即是最佳擬合。最小二乘這個(gè)術(shù)語僅僅是使誤差平方和最小的省略說法。在MATLAB中，函數(shù)polyfit求解最小二乘曲線擬合問題。為了闡述這個(gè)函數(shù)的用法， 讓我們以上面圖1中的數(shù)據(jù)開始。x=0 .1 .2 .3 .4 .5 .6 .7 .8 .9 1;y=-.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30

10、 11.2;為了用polyfit，我們必須給函數(shù)賦予上面的數(shù)據(jù)和我們希望最佳擬合數(shù)據(jù)的多項(xiàng)式的階 次或度。如果我們選擇n=1作為階次，得到最簡單的線性近似。通常稱為線性回歸。相反， 如果我們選擇n=2作為階次，得到一個(gè)2階多項(xiàng)式。現(xiàn)在，我們選擇一個(gè)2階多項(xiàng)式。 n=2;% polynomial order p=polyfit(x, y, n)P =-9.810820.1293-0.0317polyfit的輸出是一個(gè)多項(xiàng)式系數(shù)的行向量。其解是y = 9.8108x2 + 20.1293x 一 0.0317。為了將曲線擬合解與數(shù)據(jù)點(diǎn)比較，讓我們把二者都繪成圖。ezplot(-9.8108*x*x+

11、20.1293*x-0.0317) xi=linspace(0, 1, 100);% x-axis data for plotting z=polyval(p, xi);為了計(jì)算在xi數(shù)據(jù)點(diǎn)的多項(xiàng)式值，調(diào)用MATLAB的函數(shù)polyval。 plot(x, y, o , x, y, xi, z,:)畫出了原始數(shù)據(jù)x和y，用o標(biāo)出該數(shù)據(jù)點(diǎn)，在數(shù)據(jù)點(diǎn)之間，再用直線重畫原始數(shù)據(jù)， 并用點(diǎn)：線，畫出多項(xiàng)式數(shù)據(jù)xi和z。 xlabel( x ), ylabel( y=f(x) ), title( Second Order Curve Fitting )將圖作標(biāo)志。這些步驟的結(jié)果表示于前面的圖1中。多項(xiàng)式

12、階次的選擇是有點(diǎn)任意的。兩點(diǎn)決定一直線或一階多項(xiàng)式。三點(diǎn)決定一個(gè)平方或 2階多項(xiàng)式。按此進(jìn)行，n+1數(shù)據(jù)點(diǎn)唯一地確定n階多項(xiàng)式。于是，在上面的情況下，有11 個(gè)數(shù)據(jù)點(diǎn)，我們可選一個(gè)高達(dá)10階的多項(xiàng)式。然而，高階多項(xiàng)式給出很差的數(shù)值特性，人 們不應(yīng)選擇比所需的階次高的多項(xiàng)式。此外，隨著多項(xiàng)式階次的提高，近似變得不夠光滑， 因?yàn)檩^高階次多項(xiàng)式在變零前，可多次求導(dǎo)。例如，選一個(gè)10階多項(xiàng)式 pp=polyfit(x, y, 10); format short e % change display format pp.% display polynomial coefficients as a col

13、umnans =-4.6436e+0052.2965e+006-4.8773e+0065.8233e+006-4.2948e+0062.0211e+006-6.0322e+0051.0896e+005-1.0626e+0044.3599e+002-4.4700e-001要注意在現(xiàn)在情況下，多項(xiàng)式系數(shù)的規(guī)模與前面的2階擬合的比較。還要注意在最小 (-4.4700e-001)和最大(5.8233e+006)系數(shù)之間有7個(gè)數(shù)量級(jí)的幅度差。將這個(gè)解作圖，并把此 圖與原始數(shù)據(jù)及2階曲線擬合相比較，結(jié)果如何呢？ zz=polyval(pp, xi);% evaluate 10th order polyno

14、mial plot(x, y, o , xi, z, : , xi, zz) % plot data xlabel( x ), ylabel( y=f(x) ), title( 2nd and 10th Order curve Fitting )在下面的圖11.2中，原始數(shù)據(jù)標(biāo)以o，2階曲線擬合是虛線，10階擬合是實(shí)線。注意， 在10階擬合中，在左邊和右邊的極值處，數(shù)據(jù)點(diǎn)之間出現(xiàn)大的紋波。當(dāng)企圖進(jìn)行高階曲線 擬合時(shí)，這種紋波現(xiàn)象經(jīng)常發(fā)生。根據(jù)圖2,顯然，越多就越好的觀念在這里不適用。圖2 2階和10階曲線擬合維插值正如在前對(duì)曲線擬合所描述的那樣，插值定義為對(duì)數(shù)據(jù)點(diǎn)之間函數(shù)的估值方法，這些數(shù) 據(jù)

15、點(diǎn)是由某些集合給定。當(dāng)人們不能很快地求出所需中間點(diǎn)的函數(shù)值時(shí)，插值是一個(gè)有價(jià)值 的工具。例如，當(dāng)數(shù)據(jù)點(diǎn)是某些實(shí)驗(yàn)測量的結(jié)果或是過長的計(jì)算過程時(shí)，就有這種情況。或許最簡單插值的例子是MATLAB的作圖。按缺省，MATLAB用直線連接所用的數(shù)據(jù) 點(diǎn)以作圖。這個(gè)線性插值猜測中間值落在數(shù)據(jù)點(diǎn)之間的直線上。當(dāng)然，當(dāng)數(shù)據(jù)點(diǎn)個(gè)數(shù)的增 加和它們之間距離的減小時(shí)，線性插值就更精確。例如， x1=linspace(0, 2*pi, 60); x2=linspace(0, 2*pi, 6); plot(x1, sin(x1), x2, sin(x2),-) xlabel( x ), ylabel( sin(x) )

16、, title( Linear Interpolation )圖3線性插值圖3是sine函數(shù)的兩個(gè)圖，一個(gè)在數(shù)據(jù)點(diǎn)之間用60個(gè)點(diǎn)，它比另一個(gè)只用6個(gè)點(diǎn)更光 滑和更精確。如曲線擬合一樣，插值要作決策。根據(jù)所作的假設(shè)，有多種插值。而且，可以在一維以 上空間中進(jìn)行插值。即如果有反映兩個(gè)變量函數(shù)的插值，z=f(x, y)，那么就可在x之間和在 y之間，找出z的中間值進(jìn)行插值。MATLAB在一維函數(shù)interpl和在二維函數(shù)interp2中， 提供了許多的插值選擇。其中的每個(gè)函數(shù)將在下面闡述。為了說明一維插值，考慮下列問題，12小時(shí)內(nèi)，一小時(shí)測量一次室外溫度。數(shù)據(jù)存儲(chǔ)在兩個(gè)MATLAB變量中。 hour

17、s=1:12;% index for hour data was recorded temps=5 8 9 15 25 29 31 30 22 25 27 24; % recorded temperatures plot(hours, temps, hours, temps, + )% view temperatures title( Temperature ) xlabel( Hour ), ylabel( Degrees Celsius )圖4在線性插值下室外溫度曲線正如圖.4看到的，MATLAB畫出了數(shù)據(jù)點(diǎn)線性插值的直線。為了計(jì)算在任意給定時(shí)間 的溫度，人們可試著對(duì)可視的圖作解釋。另外一

18、種方法，可用函數(shù)interp。 t=interp1(hours, temps, 9.3) t =22.9000% estimate temperature at hour=9.3 t=interp1(hours, temps, 4.7) t =22% estimate temperature at hour=4.7 t=interp1(hours, temps, 3.2 6.5 7.1 11.7) t =% find temp at many points!10.200030.000030.900024.9000interpl的缺省用法是由interp1(x, y, xo)來描述，這里x是獨(dú)立

19、變量(橫坐標(biāo))，y是應(yīng)變 量(縱坐標(biāo))，xo是進(jìn)行插值的一個(gè)數(shù)值數(shù)組。另外，該缺省的使用假定為線性插值。若不采用直線連接數(shù)據(jù)點(diǎn)，我們可采用某些更光滑的曲線來擬合數(shù)據(jù)點(diǎn)。最常用的方法 是用一個(gè)3階多項(xiàng)式，即3次多項(xiàng)式，來對(duì)相繼數(shù)據(jù)點(diǎn)之間的各段建模，每個(gè)3次多項(xiàng)式的 頭兩個(gè)導(dǎo)數(shù)與該數(shù)據(jù)點(diǎn)相一致。這種類型的插值被稱為3次樣條或簡稱為樣條。函數(shù)interpl 也能執(zhí)行3次樣條插值。 t=interp1(hours, temps, 9.3, spline )% estimate temperature at hour=9.3t =21.8577 t=interp1(hours, temps, 4.7,

20、spline )% estimate temperature at hour=4.7t =22.3143 t=interp1(hours, temps, 3.2 6.5 7.1 11.7, spline )t =9.673430.042731.175525.3820注意，樣條插值得到的結(jié)果，與上面所示的線性插值的結(jié)果不同。因?yàn)椴逯凳且粋€(gè)估計(jì) 或猜測的過程，其意義在于，應(yīng)用不同的估計(jì)規(guī)則導(dǎo)致不同的結(jié)果。一個(gè)最常用的樣條插值是對(duì)數(shù)據(jù)平滑。也就是，給定一組數(shù)據(jù)，使用樣條插值在更細(xì)的 間隔求值。例如， h=1:0.1:12;% estimate temperature every 1/10 hour

21、t=interp1(hours, temps, h, spline ); plot(hours, temps, - , hours, temps, + , h, t) % plot comparative results title( Springfield Temperature ) xlabel( Hour ), ylabel( Degrees Celsius )在圖5中，虛線是線性插值，實(shí)線是平滑的樣條插值，標(biāo)有+ 的是原始數(shù)據(jù)。如要求 在時(shí)間軸上有更細(xì)的分辨率，并使用樣條插值，我們有一個(gè)更平滑、但不一定更精確地對(duì)溫 度的估計(jì)。尤其應(yīng)注意，在數(shù)據(jù)點(diǎn)，樣條解的斜率不突然改變。作為這個(gè)平滑插

22、值的回報(bào)，3次樣條插值要求更大量的計(jì)算，因?yàn)楸仨氄业?次多項(xiàng)式以描述給定數(shù)據(jù)之間的特征。圖5在不同插值下室外溫度曲線在討論二維插值之前，了解interpl所強(qiáng)制的二個(gè)強(qiáng)約束是很重要的。首先，人們不能 要求有獨(dú)立變量范圍以外的結(jié)果，例如，interp1(hours, temps, 13.5)導(dǎo)致一個(gè)錯(cuò)誤，因?yàn)?hours在1到12之間變化。其次，獨(dú)立變量必須是單調(diào)的。即獨(dú)立變量在值上必須總是增 加的或總是減小的。在我們的例子里，hours是單調(diào)的。然而，如果我們已經(jīng)定義獨(dú)立變量 為一天的實(shí)際時(shí)間， time_of_day=7:12 1:6% start at 7AM,end at 6PMtime

23、_of_day =789101112123456則獨(dú)立變量將不是單調(diào)的，因?yàn)閠ime_of_day增加到12，然后跌到1，再然后增加。如果用 time_of_day代替interp1中的hours，將會(huì)返回一個(gè)錯(cuò)誤。同樣的理由，人們不能對(duì)temps 插值來找出產(chǎn)生某溫度的時(shí)間(小時(shí))，因?yàn)閠emps不是單調(diào)的。二維插值二維插值是基于與一維插值同樣的基本思想。然而，正如名字所隱含的，二維插值是對(duì) 兩變量的函數(shù)z=f(x, y)進(jìn)行插值。為了說明這個(gè)附加的維數(shù)，考慮一個(gè)問題。設(shè)人們對(duì)平 板上的溫度分布估計(jì)感興趣，給定的溫度值取自平板表面均勻分布的格柵。采集了下列的數(shù)據(jù)： depth=1:3;%in

24、dex for depth of plate (i,e,the y-dimension) temps=8281 808284; 79 63 61 65 81; 84 84 82 85 86%temperature datatemps =828180828479636165818484828586% index for width of plate (i.e.,the x-dimension)如同在標(biāo)引點(diǎn)上測量一樣，矩陣temps表示整個(gè)平板的溫度分布。temps的列與下標(biāo) depth或y-維相聯(lián)系，行與下標(biāo)width或x-維相聯(lián)系(見圖11.6)。為了估計(jì)在中間點(diǎn)的溫度， 我們必須對(duì)它們進(jìn)行辨

25、識(shí)。 width=1:5; wi=1:0.2:5;% estimate across width of plate d=2;% at a depth of 2 zlinear=interp2(width, depth, temps, wi, d) ;% linear interpolation zcubic=interp2(width, depth, temps, wi,d, cubic ) ;% cubic interpolation plot(wi, zlinear, - , wi, zcubic) % plot results xlabel( Width of Plate ), ylabel( Degrees Celsius ) title( Temperature at Depth = num2str(d)另一種方法，我們可以在兩個(gè)方向插值。先在三維坐標(biāo)畫出原始數(shù)據(jù)，看一下該數(shù)據(jù)的 粗糙程度(見圖7)。 mesh(width, depth, temps)% use mesh plot xlabel( Width of Plate ), ylab