陳雄達-數(shù)學軟件介紹-2016-數(shù)模講座100

數(shù)學軟件介紹-Matlab

同濟大學陳雄達2016年9月4日提綱MatlaB數(shù)據(jù)分析工具MATLAB微分方程求解Matlab圖論工具Matlab最優(yōu)化問題求解MATLAB平臺上的算法實現(xiàn)2Matlab基礎(chǔ)知識矩陣功能(高維矩陣)計算可視化(圖形工具箱)超多工具箱(方程、統(tǒng)計、優(yōu)化等)和其他編程或軟件工具協(xié)同(C語言,Excel等)3Matlab數(shù)據(jù)分析工具數(shù)據(jù)文件的讀寫純文本數(shù)據(jù)文件Fprintf(fid,format,Data)DATA=FSCANF(FID,FORMAT,SIZE)相關(guān)命令：fopen,fread,fcloseMatlab內(nèi)置數(shù)據(jù)格式LOADfnamevarSAVEfnamevar–mat/-ascii4Matlab數(shù)據(jù)分析工具數(shù)據(jù)文件的讀寫EXcel數(shù)據(jù)文件DATA=XLSREAD(FNAME,SHEET,range)Xlswrite(fname,data,sheet,range)Range’A1:C5’例題2014D儲藥柜的設(shè)計2012A葡萄酒的評價5Matlab數(shù)據(jù)分析工具數(shù)據(jù)文件的讀寫圖片數(shù)據(jù)文件A=IMREAD(FNAME,ForMaT)IMWRITE(A,FNAME,ForMaT)例題2013B碎紙片的拼接復原2011B交巡警服務平臺的設(shè)置與調(diào)度2008A數(shù)碼相機定位6Matlab數(shù)據(jù)分析工具數(shù)據(jù)文件的讀寫視頻數(shù)據(jù)文件MOV=AVIREAD(FnAME)movie例題2013A車道被占用對城市道路通行能力的影響（2014B創(chuàng)意平板折疊桌）7Matlab數(shù)據(jù)分析工具數(shù)據(jù)處理插值interp1/2/3,

YI=INTERP1(X,Y,XI,method)griddata,

[XI,YI,ZI]=GRIDDATA(X,Y,Z,XI,YI)Spline,yy=spline(x,y,xx);例題2011A城市表層土壤重金屬污染分析8Matlab數(shù)據(jù)分析工具數(shù)據(jù)處理擬合P=Polyfit(x,Y,n)Beta=NLINFIT(X,Y,FUN,BETA0)functionyhat=fun(beta,x)

b1=beta(1);

b2=beta(2);

yhat=1./(1+exp(b1+b2*x));例題2008A數(shù)碼相機定位

9Matlab數(shù)據(jù)分析工具數(shù)據(jù)統(tǒng)計分析描述集中趨勢的統(tǒng)計量均值geomean,mean,median,trimean描述離散趨勢的統(tǒng)計量內(nèi)四分極值iqr,極差range,方差var例題2012A葡萄酒的評價10Matlab數(shù)據(jù)分析工具數(shù)據(jù)統(tǒng)計分析分組數(shù)據(jù)描述Means=grpstats(x,group)包含缺失數(shù)據(jù)的樣本描述Nan+max,min,mean,median,std,sum百分位Y=prctile(x,p)11Matlab數(shù)據(jù)分析工具其他數(shù)據(jù)統(tǒng)計工具中心矩M=moment(x,order)相關(guān)系數(shù)R=corrcoef(X)R=corrcoef(X,Y)峰度kurtosis,偏度skewness,

頻數(shù)tabulate,

列聯(lián)表crosstab12MATLAB微分方程求解13精確求解dsolveS=dsolve('Dx=y','Dy=-x','x(0)=0','y(0)=1')y=dsolve('D2y=sin(y)');pretty(y)MATLAB微分方程求解14近似求解Ode45,ode23s[t,y]=ode45(@fun,[020],2);plot(t,y);functiondydx=fun(t,y)

M=1e11;dydx=10*y*(1-y/M)

MATLAB微分方程求解15

MATLAB微分方程求解16兩點邊值問題近似求解例如：dydx

=

@(x,y)[y(2);(2*y(1)-4)/(1+x)];

%邊值微分方程

res

=

@(ya,yb)[yb(2)-yb(1)/2;ya(1)];

%邊界條件

solinit

=

bvpinit(linspace(0,1,10),[1

0]);%x及猜測值

sol

=

bvp4c(dydx,res,solinit);%求解

xint

=

linspace(0,1,50);

Sxint

=

deval_r(sol,xint);%畫圖

MATLAB圖論工具17最短路徑求解、最大流、最小生成樹DIST=GRAPHSHORTESTPATH(G,S,D)graphallshortestpaths(G)[M,F,C]=graphmaxflow(G2,s,D)[ST,pred]=graphminspantree(G)MATLAB圖論工具18最短路徑求解及最大流R

=

[112412335734567];C

=

[233345566677888];W

=

[286711513429863];G

=

sparse(R,C,W);view(biograph(G,[],'ShowW','ON'));[distpath]=graphshortestpath(G,1,8);graphallshortestpaths(G);h=view(biograph(G,[],'showW','on'));edges=getedgesbynodeid(h,get(h.Nodes(path),'ID'));set(h.Nodes(path),'color',[100]);

set(edges,'LineColor',[100]);

set(edges,'LineWidth',1.5);[M,F,C]=graphmaxflow(G2,1,8);m=view(biograph(F,[],'showW','ON'));Matlab最優(yōu)化問題求解19無約束優(yōu)化問題x=fminbnd(fun,x1,x2)%單變量x=fminunc(fun,x0)%多變量x=fminsearch(fun,x0)

%多變量算法參數(shù)設(shè)置optimsetoptimgetMatlab最優(yōu)化問題求解20約束優(yōu)化問題線性規(guī)劃x=linprog(f,A,b,Aeq,beq,lb,ub,x0)二次規(guī)劃x=

quadprog(H,f,A,b,Aeq,beq,lb,ub,x0)例題：2000B鋼管訂購和運輸非線性規(guī)劃x=fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon)Matlab最優(yōu)化問題求解21特殊形式的約束優(yōu)化問題最大最小問題x=fminimax(fun,x0,A,B,Aeq,beq)最小二乘問題x=lsqnonneg(C,d,x0)x=lsqlin(C,D,A,B,Aeq,beq)非線性最小二乘x=lsqnonlin(fun,x0,lb,ub)非線性曲線擬合x=lsqcurvefit(fun,x0,xdata,ydata)Matlab最優(yōu)化問題求解22

Matlab最優(yōu)化問題求解23例：預測田鼠和貓頭鷹的種群數(shù)量變化。每兩個月作一次田間調(diào)查，已知它們種群的大致數(shù)目如下田鼠29.733.132.569.1134.2236.0269.6162.369.6貓頭鷹12810488968810414417619239.834.020.721.737.657.6124.6215.8272.7195.795.016815212012096728810412818419241.925.710.922.633.648.192.5183.3268.5230.6111.116813614411296807288104152184Matlab最優(yōu)化問題求解

24Matlab最優(yōu)化問題求解

25Matlab最優(yōu)化問題求解functionvolterrax=[29.733.132.569.1134.2236.0269.6162.369.639.834.0...20.721.737.657.6124.6215.8272.7195.795.041.925.7...10.922.633.648.192.5183.3268.5230.6111.1]’;y=[128104889688104144176192168152...

120120967288104128184192168136...14411296807288104152184]’;n=length(x);te=2;%timeelapseplot(te*(0:n-1),x,'ro',te*(0:n-1),y,'b*');legend('mice','owl');holdon;dt=te*ones(n-1,1);Sy=te*(y(1:end-1)+y(2:end))/2;Sx=te*(x(1:end-1)+x(2:end))/2;p1=[dt-Sy]\[diff(log(x))];p2=[-dtSx]\[diff(log(y))];26Matlab最優(yōu)化問題求解27delta=p1(1);mu=p2(1);alpha=p1(2);beta=p2(2);[t,y]=ode45(@vollot,[0:0.5:60],[x(1)y(1)]',odeset,...delta,mu,alpha,beta);pauseplot(t,y(:,1),'r-',t,y(:,2),'b-','linewidth',2);functiondydt=vollot(t,y,delta,mu,alpha,beta)dydt=[(delta-alpha*y(2))*y(1)(-mu+beta*y(1))*y(2)];Matlab最優(yōu)化問題求解28如果分析不到位，直接處理整個流程？（效率很低）

functiondydt=vollot(t,y,delta,mu,alpha,beta)dydt=[(delta-alpha*y(2))*y(1)(-mu+beta*y(1))*y(2)];Matlab最優(yōu)化問題求解29如果分析不到位，直接處理整個流程？（效率很低）functionerr=pray(para)delta=para(1);alpha=para(2);mu=para(3);beta=para(4);x0=para(5);y0=para(6);x=[29.733.132.569.1134.2236.0269.6162.369.639.834.0...20.721.737.657.6124.6215.8272.7195.795.041.925.7...10.922.633.648.192.5183.3268.5230.6111.1]';y=[128104889688104144176192168152...120120967288104128184192168136...