復(fù)習(xí)建模培訓(xùn)matlab應(yīng)用

MATLAB

軟件及其應(yīng)用王林君江蘇大學(xué)理學(xué)院

ApplicationofMatlabLanguage1

MATLAB應(yīng)用--解線性規(guī)劃

2用MATLAB優(yōu)化工具箱解線性規(guī)劃minz=cX

1、模型：命令：x=linprog〔c，A，b〕2、模型：minz=cX

命令：x=linprog〔c，A，b，Aeq,beq〕注意：若沒(méi)有不等式：存在，則令A(yù)=[]，b=[].33、模型：minz=cX

VLB≤X≤VUB命令：[1]x=linprog〔c，A，b，Aeq,beq,VLB，VUB〕[2]x=linprog〔c，A，b，Aeq,beq,VLB，VUB,X0〕注意：[1]若沒(méi)有等式約束:,則令A(yù)eq=[],beq=[].[2]其中X0表示初始點(diǎn)4、命令：[x,fval]=linprog(…)返回最優(yōu)解ｘ及ｘ處的目標(biāo)函數(shù)值fval.4解編寫(xiě)M文件xxgh1.m如下：c=[-0.4-0.28-0.32-0.72-0.64-0.6];A=[0.010.010.010.030.030.03;0.02000.0500;00.02000.050;000.03000.08];b=[850;700;100;900];Aeq=[];beq=[];vlb=[0;0;0;0;0;0];vub=[];[x,fval]=linprog(c,A,b,Aeq,beq,vlb,vub)

5解:編寫(xiě)M文件xxgh2.m如下：

c=[634];A=[010];b=[50];Aeq=[111];beq=[120];vlb=[30,0,20];vub=[];[x,fval]=linprog(c,A,b,Aeq,beq,vlb,vub)6

MATLAB應(yīng)用—

解非線性規(guī)劃

7一、二次規(guī)劃(QuadraticProgram)概念8二、Matlab中求解二次規(guī)劃910轉(zhuǎn)化為matlab求解格式：1112定義如果目標(biāo)函數(shù)或約束條件中至少有一個(gè)是非線性函數(shù),那么最優(yōu)化問(wèn)題就叫做非線性規(guī)劃問(wèn)題．四、非線性規(guī)劃的根本概念一般形式:

（1）其中，是定義在Rn

上的實(shí)值函數(shù)()n

TnRxxxX?=,,,21L()()???íì===3.,...,2,1

0

m;1,2,...,

0..

ljXhiXgtsji13五、非線性規(guī)劃的根本解法SUTM外點(diǎn)法SUTM內(nèi)點(diǎn)法〔障礙罰函數(shù)法〕1．罰函數(shù)法2．近似線性規(guī)劃法14

1、罰函數(shù)法罰函數(shù)法根本思想是通過(guò)構(gòu)造罰函數(shù)把約束問(wèn)題轉(zhuǎn)化為一系列無(wú)約束最優(yōu)化問(wèn)題，進(jìn)而用無(wú)約束最優(yōu)化方法去求解．這類(lèi)方法稱為序列無(wú)約束最小化方法．簡(jiǎn)稱為SUMT法．其一為SUMT外點(diǎn)法，其二為SUMT內(nèi)點(diǎn)法．15Matlab求解非線性規(guī)劃問(wèn)題

其中X為n維變?cè)蛄?，G(X)與Ceq(X)均為非線性函數(shù)組成的向量。161．首先建立M文件fun.m,用來(lái)定義目標(biāo)函數(shù)F〔X〕:functionf=fun(X);f=F(X);MATLAB求解上述問(wèn)題，根本步驟分三步173．建立主程序.求解非線性規(guī)劃的函數(shù)是fmincon,命令的根本格式如下：(1)x=fmincon(‘fun’,X0,A,b)(2)x=fmincon(‘fun’,X0,A,b,Aeq,beq)(3)x=fmincon(‘fun’,X0,A,b,Aeq,beq,VLB,VUB)(4)x=fmincon(‘fun’,X0,A,b,Aeq,beq,VLB,VUB,’nonlcon’)(5)x=fmincon(‘fun’,X0,A,b,Aeq,beq,VLB,VUB,’nonlcon’,options)

(6)[x,fval]=fmincon(…)(7)[x,fval,exitflag]=fmincon(…)(8)[x,fval,exitflag,output]=fmincon(…)輸出極值點(diǎn)M文件迭代的初值參數(shù)說(shuō)明變量上下限fmincon函數(shù)可能會(huì)給出局部最優(yōu)解，這與初值X0的選取有關(guān)．181．寫(xiě)成標(biāo)準(zhǔn)形式：

s.t.

2x1+3x26

s.t.

x1+4x25

x1,x20例192．先建立M-文件fun3．m:

functionf=fun3(x);f=-x(1)-2*x(2)+(1/2)*x(1)^2+(1/2)*x(2)^23．再建立主程序youh2．m：

x0=[1;1];A=[23;14];b=[6;5];Aeq=[];beq=[];VLB=[0;0];VUB=[];[x,fval]=fmincon('fun3',x0,A,b,Aeq,beq,VLB,VUB)4．運(yùn)算結(jié)果為：

x=0．76471．0588fval=-2．0294201．先建立M文件fun4．m定義目標(biāo)函數(shù):

functionf=fun4(x);f=exp(x(1))*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+2*x(2)+1);

x1+x2=0s.t.1.5+x1x2-x1-x20

-x1x2–10

0例

2．再建立M文件mycon．m定義非線性約束：

function[g,ceq]=mycon(x)g=[x(1)+x(2);1．5+x(1)*x(2)-x(1)-x(2);-x(1)*x(2)-10];213．主程序youh3．m為:x0=[-1;1];A=[];b=[];Aeq=[11];beq=[0];vlb=[];vub=[];[x,fval]=fmincon('fun4',x0,A,b,Aeq,beq,vlb,vub,'mycon')4．運(yùn)算結(jié)果為：

x=-1．22501．2250fval=1．895122例1．先建立M文件fun．m定義目標(biāo)函數(shù):functionf=fun(x);f=-2*x(1)-x(2);2．再建立M文件mycon2．m定義非線性約束：function[g,ceq]=mycon2(x)g=[x(1)^2+x(2)^2-25;x(1)^2-x(2)^2-7];

233．主程序fxx．m為:

x0=[3;2．5];VLB=[00];VUB=[510];[x,fval,exitflag,output]=fmincon('fun',x0,[],[],[],[],VLB,VUB,'mycon2')24

MATLAB應(yīng)用—

微分方程解法

25要求目的主要內(nèi)容2、學(xué)會(huì)用Matlab求微分方程的數(shù)值解.1、學(xué)會(huì)用Matlab求簡(jiǎn)單微分方程的解析解.1、求簡(jiǎn)單微分方程的解析解.2、求微分方程的數(shù)值解.26求微分方程的數(shù)值解〔一〕常微分方程數(shù)值解的定義〔二〕建立數(shù)值解法的一些途徑〔三〕用Matlab軟件求常微分方程的數(shù)值解返回27微分方程的解析解求微分方程（組）的解析解命令:dsolve(‘方程1’,‘方程2’,…‘方程n’,‘初始條件’,‘自變量’)結(jié)果：u=tg(t-c)28解輸入命令:y=dsolve('D2y+4*Dy+29*y=0','y(0)=0,Dy(0)=15','x')結(jié)果為:y=3e-2xsin〔5x〕29解輸入命令：

[x,y,z]=dsolve('Dx=2*x-3*y+3*z','Dy=4*x-5*y+3*z','Dz=4*x-4*y+2*z','t')；

x=simple(x)%將x化簡(jiǎn)

y=simple(y)z=simple(z)結(jié)果為：x=(c1-c2+c3+c2e-3t-c3e-3t)e2t

y=-c1e-4t+c2e-4t+c2e-3t-c3e-3t+c1-c2+c3)e2tz=(-c1e-4t+c2e-4t+c1-c2+c3)e2t

30微分方程的數(shù)值解〔一〕常微分方程數(shù)值解的定義在生產(chǎn)和科研中所處理的微分方程往往很復(fù)雜且大多得不出一般解。而在實(shí)際上對(duì)初值問(wèn)題，一般是要求得到解在假設(shè)干個(gè)點(diǎn)上滿足規(guī)定精確度的近似值，或者得到一個(gè)滿足精確度要求的便于計(jì)算的表達(dá)式。因此，研究常微分方程的數(shù)值解法是十分必要的。31〔二〕建立數(shù)值解法的一些途徑1、用差商代替導(dǎo)數(shù)若步長(zhǎng)h較小，則有故有公式：此即歐拉法。322、使用數(shù)值積分對(duì)方程y’=f(x,y),兩邊由xi到xi+1積分，并利用梯形公式，有：實(shí)際應(yīng)用時(shí)，與歐拉公式結(jié)合使用：此即改進(jìn)的歐拉法。故有公式：333、使用泰勒公式以此方法為根底，有龍格-庫(kù)塔法、線性多步法等方法。4、數(shù)值公式的精度當(dāng)一個(gè)數(shù)值公式的截?cái)嗾`差可表示為O〔hk+1〕時(shí)〔k為正整數(shù)，h為步長(zhǎng)〕，稱它是一個(gè)k階公式。k越大，那么數(shù)值公式的精度越高。歐拉法是一階公式，改進(jìn)的歐拉法是二階公式。龍格-庫(kù)塔法有二階公式和四階公式。線性多步法有四階阿達(dá)姆斯外插公式和內(nèi)插公式。34〔三〕用Matlab軟件求常微分方程的數(shù)值解[t，x]=solver〔’f’,ts,x0,options〕ode45ode23ode113ode15sode23s由待解方程寫(xiě)成的m-文件名ts=[t0，tf]，t0、tf為自變量的初值和終值函數(shù)的初值ode23：組合的2/3階龍格-庫(kù)塔-芬爾格算法ode45：運(yùn)用組合的4/5階龍格-庫(kù)塔-芬爾格算法自變量值函數(shù)值用于設(shè)定誤差限(缺省時(shí)設(shè)定相對(duì)誤差10-3,絕對(duì)誤差10-6),命令為：options=odeset（’reltol’,rt,’abstol’,at）,rt，at：分別為設(shè)定的相對(duì)誤差和絕對(duì)誤差.35

1、在解n個(gè)未知函數(shù)的方程組時(shí)，x0和x均為n維向量，m-文件中的待解方程組應(yīng)以x的分量形式寫(xiě)成.

2、使用Matlab軟件求數(shù)值解時(shí)，高階微分方程必須等價(jià)地變換成一階微分方程組.注意:36解:令y1=x，y2=y1’1、建立m-文件vdp1000.m如下：

functiondy=vdp1000(t,y)dy=zeros(2,1);dy(1)=y(2);dy(2)=1000*(1-y(1)^2)*y(2)-y(1);

2、取t0=0，tf=3000，輸入命令：

[T,Y]=ode15s('vdp1000',[03000],[20]);plot(T,Y(:,1),'-')3、結(jié)果如圖37解

1、建立m-文件rigid.m如下：

functiondy=rigid(t,y)dy=zeros(3,1);dy(1)=y(2)*y(3);dy(2)=-y(1)*y(3);dy(3)=-0.51*y(1)*y(2);2、取t0=0，tf=12，輸入命令：

[T,Y]=ode45('rigid',[012],[011]);plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'+')3、結(jié)果如圖圖中，y1的圖形為實(shí)線，y2的圖形為“*〞線，y3的圖形為“+〞線.38Matlab的應(yīng)用

-----------插值和擬合39拉格朗日插值分段線性插值三次樣條插值一維插值一、插值的定義二、插值的方法三、用Matlab解插值問(wèn)題40二維插值一、二維插值定義二、網(wǎng)格節(jié)點(diǎn)插值法三、用Matlab解插值問(wèn)題最鄰近插值分片線性插值雙線性插值網(wǎng)格節(jié)點(diǎn)數(shù)據(jù)的插值散點(diǎn)數(shù)據(jù)的插值41一維插值的定義已知n+1個(gè)節(jié)點(diǎn)其中互不相同，不妨設(shè)求任一插值點(diǎn)處的插值

節(jié)點(diǎn)可視為由產(chǎn)生,，表達(dá)式復(fù)雜,，或無(wú)封閉形式或未知

42構(gòu)造一個(gè)(相對(duì)簡(jiǎn)單的)函數(shù)通過(guò)全部節(jié)點(diǎn),即再用計(jì)算插值，即

43稱為拉格朗日插值基函數(shù)。函數(shù)f(x)在n+1個(gè)點(diǎn)x0,x1,…,xn處的函數(shù)值為y0,y1,…,yn。求一n次多項(xiàng)式函數(shù)Pn(x)，使其滿足：Pn(xi)=yi,i=0,1,…,n.解決此問(wèn)題的拉格朗日插值多項(xiàng)式公式如下其中Li(x)為n次多項(xiàng)式：拉格朗日(Lagrange)插值44拉格朗日(Lagrange)插值特別地:兩點(diǎn)一次(線性)插值多項(xiàng)式:三點(diǎn)二次(拋物)插值多項(xiàng)式:45分段線性插值計(jì)算量與n無(wú)關(guān);n越大，誤差越小.

xjxj-1xj+1x0xnxoy46比分段線性插值更光滑。

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

三次樣條插值48用MATLAB作插值計(jì)算一維插值函數(shù)：yi=interp1(x，y，xi，'method')插值方法被插值點(diǎn)插值節(jié)點(diǎn)xi處的插值結(jié)果‘nearest’

：最鄰近插值‘linear’

：線性插值；‘spline’

：三次樣條插值；‘cubic’

：立方插值。缺省時(shí)：分段線性插值。注意：所有的插值方法都要求x是單調(diào)的，并且xi不能夠超過(guò)x的范圍。49例：在1-12的11小時(shí)內(nèi)，每隔1小時(shí)測(cè)量一次溫度，測(cè)得的溫度依次為：5，8，9，15，25，29，31，30，22，25，27，24。試估計(jì)每隔1/10小時(shí)的溫度值。hours=1:12;temps=[589152529313022252724];h=1:0.1:12;t=interp1(hours,temps,h,'spline');(直接輸出數(shù)據(jù)將是很多的)plot(hours,temps,'+',h,t,hours,temps,'r:')%作圖xlabel('Hour'),ylabel('DegreesCelsius’)50

要求x0,y0單調(diào)；x，y可取為矩陣，或x取行向量，y取為列向量，x,y的值分別不能超出x0,y0的范圍。z=interp2(x0,y0,z0,x,y,’method’)被插值點(diǎn)插值方法用MATLAB作網(wǎng)格節(jié)點(diǎn)數(shù)據(jù)的插值插值節(jié)點(diǎn)被插值點(diǎn)的函數(shù)值‘nearest’最鄰近插值‘linear’雙線性插值‘cubic’雙三次插值缺省時(shí),雙線性插值51例：測(cè)得平板外表3*5網(wǎng)格點(diǎn)處的溫度分別為：828180828479636165818484828586試作出平板外表的溫度分布曲面z=f(x,y)的圖形。輸入以下命令：x=1:5;y=1:3;temps=[8281808284;7963616581;8484828586];mesh(x,y,temps)1.先在三維坐標(biāo)畫(huà)出原始數(shù)據(jù)，畫(huà)出粗糙的溫度分布曲圖.2．以平滑數(shù)據(jù),在x、y方向上每隔0.2個(gè)單位的地方進(jìn)行插值.52再輸入以下命令:xi=1:0.2:5;yi=1:0.2:3;zi=interp2(x,y,temps,xi',yi,'cubic');mesh(xi,yi,zi)畫(huà)出插值后的溫度分布曲面圖.53

插值函數(shù)griddata格式為:cz=griddata〔x，y，z，cx，cy，‘method’〕用MATLAB作散點(diǎn)數(shù)據(jù)的插值計(jì)算

要求cx取行向量，cy取為列向量。被插值點(diǎn)插值方法插值節(jié)點(diǎn)被插值點(diǎn)的函數(shù)值‘nearest’最鄰近插值‘linear’雙線性插值‘cubic’雙三次插值'v4'-Matlab提供的插值方法缺省時(shí),雙線性插值54擬合2.擬合的根本原理1.擬合問(wèn)題引例55擬合問(wèn)題引例1溫度t(0C)20.532.751.073.095.7電阻R()7658268739421032熱敏電阻數(shù)據(jù)：求600C時(shí)的電阻R。

設(shè)

R=at+ba,b為待定系數(shù)56擬合問(wèn)題引例2

t(h)0.250.511.523468c(g/ml)19.2118.1515.3614.1012.899.327.455.243.01已知一室模型快速靜脈注射下的血藥濃度數(shù)據(jù)(t=0注射300mg)求血藥濃度隨時(shí)間的變化規(guī)律c(t).作半對(duì)數(shù)坐標(biāo)系(semilogy)下的圖形57曲線擬合問(wèn)題的提法一組〔二維〕數(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,yi)與曲線y=f(x)的距離58擬合與插值的關(guān)系函數(shù)插值與曲線擬合都是要根據(jù)一組數(shù)據(jù)構(gòu)造一個(gè)函數(shù)作為近似，由于近似的要求不同，二者的數(shù)學(xué)方法上完全不同。實(shí)例：下面數(shù)據(jù)是某次實(shí)驗(yàn)所得，希望得到x和f之間的關(guān)系？問(wèn)題：給定一批數(shù)據(jù)點(diǎn)，需確定滿足特定要求的曲線或曲面解決方案：假設(shè)不要求曲線〔面〕通過(guò)所有數(shù)據(jù)點(diǎn)，而是要求它反映對(duì)象整體的變化趨勢(shì)，就是數(shù)據(jù)擬合，又稱曲線擬合或曲面擬合。假設(shè)要求所求曲線〔面〕通過(guò)所給所有數(shù)據(jù)點(diǎn)，就是插值問(wèn)題；59最臨近插值、線性插值、樣條插值與曲線擬合結(jié)果：60曲線擬合問(wèn)題最常用的解法——線性最小二乘法的根本思路第一步:先選定一組函數(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的準(zhǔn)那么〔最小二乘準(zhǔn)那么〕：使n個(gè)點(diǎn)〔xi,yi)與曲線y=f(x)的距離i的平方和最小。記

問(wèn)題歸結(jié)為，求

a1,a2,…am

使

J(a1,a2,…am)

最小。61線性最小二乘法的求解：預(yù)備知識(shí)超定方程組：方程個(gè)數(shù)大于未知量個(gè)數(shù)的方程組即Ra=y其中超定方程一般是不存在解的矛盾方程組。如果有向量a使得達(dá)到最小，則稱a為上述超定方程的最小二乘解。62線性最小二乘法的求解定理：當(dāng)RTR可逆時(shí)，超定方程組〔3〕存在最小二乘解，

且即為方程組RTRa=RTy------正那么〔正規(guī)〕方程組的解：a=(RTR)-1RTy所以，曲線擬合的最小二乘法要解決的問(wèn)題，實(shí)際上就是求以下超定方程組的最小二乘解的問(wèn)題。其中Ra=y（3）63用MATLAB解擬合問(wèn)題1、線性最小二乘擬合2、非線性最小二乘擬合64用MATLAB作線性最小二乘擬合1.作多項(xiàng)式f(x)=a1xm+…+amx+am+1擬合,可利用已有程序:a=polyfit(x,y,m)2.對(duì)超定方程組可得最小二乘意義下的解。，用3.多項(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ù)65即要求出二次多項(xiàng)式:中的使得:例對(duì)下面一組數(shù)據(jù)作二次多項(xiàng)式擬合661〕輸入命令：x=0:0.1:1;y=[-0.447,1.978,3.28,6.16,7.08,7.34,7.66,9.56,9.48,9.30,11.2];R=[(x.^2)',x',ones(11,1)]；A=R\y'解法1．解超定方程的方法2〕計(jì)算結(jié)果：Ａ=[-9.8108,20.1293,-0.0317]672〕計(jì)算結(jié)果：Ａ=[-9.8108,20.1293,-0.0317]解法2．用多項(xiàng)式擬合的命令MATLAB(zxec2)1〕輸入命令：x=0:0.1:1;y=[-0.447,1.978,3.28,6.16,7.08,7.34,7.66,9.56,9.48,9.30,11.2];A=polyfit(x,y,2)z=polyval(A,x);plot(x,y,'k+',x,z,'r')%作出數(shù)據(jù)點(diǎn)和擬合曲線的圖形681.lsqcurvefit數(shù)據(jù)點(diǎn)：xdata=〔xdata1，xdata2，…，xdatan〕ydata=〔ydata1，ydata2，…，ydatan〕用MATLAB作非線性最小二乘擬合兩個(gè)求非線性最小二乘擬合的函數(shù)：lsqcurvefit、lsqnonlin。相同點(diǎn)和不同點(diǎn)：兩個(gè)命令都要先建立M-文件fun.m，定義函數(shù)f(x)，但定義f(x)的方式不同，請(qǐng)參考例題。

lsqcurvefit用以求含參量x（向量）的向量值函數(shù)F(x,xdata)=（F（x，xdata1），…，F(xiàn)（x，xdatan））T中的參變量x(向量),使得69輸入格式:(1)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);(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是一個(gè)事先建立的定義函數(shù)F(x,xdata)

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

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

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

中的參量x，使得

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

=F(x,xdatai)-ydatai2.lsqnonlin數(shù)據(jù)點(diǎn)：xdata=〔xdata1，xdata2，…，xdatan〕ydata=〔ydata1，ydata2，…，ydatan〕71輸入格式：1〕x=lsqnonlin〔‘fun’，x0〕；2〕x=lsqnonlin〔‘fun’，x0，lb，ub〕；3〕x=lsqnonlin〔‘fun’，x0，，lb，ub，options〕；4〕[x，options]=lsqnonlin〔‘fun’，x0，…〕；5〕[x，options，funval]=lsqnonlin〔‘fun’，x0，…〕；說(shuō)明：x=lsqnonlin〔‘fun’，x0，options〕；fun是一個(gè)事先建立的定義函數(shù)f(x)的M-文件，自變量為x迭代初值選項(xiàng)見(jiàn)無(wú)約束優(yōu)化72例2用下面一組數(shù)據(jù)擬合中的參數(shù)a，b，k該問(wèn)題即解最優(yōu)化問(wèn)題：731〕編寫(xiě)M-文件curvefun1.m