數(shù)學(xué)實(shí)驗(yàn)講義4

1數(shù)學(xué)實(shí)驗(yàn)——

基于MATLAB電子教案(四)——綜合實(shí)驗(yàn)制作：鄭學(xué)東zxd@2擬合與插值在生產(chǎn)與科學(xué)實(shí)驗(yàn)中,有時(shí)函數(shù)y=f(x)的解析式未知,只知該函數(shù)在若干觀測(cè)點(diǎn)的函數(shù)值或?qū)?shù)值;然而卻希望知道觀測(cè)點(diǎn)以外的函數(shù)值,需要估計(jì)函數(shù)在該點(diǎn)的值。解決思路：根據(jù)觀測(cè)點(diǎn)的值,構(gòu)造一較簡(jiǎn)單函數(shù)

(x),使其在觀測(cè)點(diǎn)處的值等于或總體最接近f(x)的已知值。

(x)構(gòu)造方法一般有兩類：①測(cè)量值與真實(shí)值有誤差,不要求

(x)必須過(guò)已知點(diǎn),一般用擬合法。②測(cè)量值準(zhǔn)確無(wú)誤,要求

(x)必須過(guò)已知點(diǎn),一般用插值法。3插值在已知的離散點(diǎn)的數(shù)據(jù)，構(gòu)造出一個(gè)能連接所有已知點(diǎn)的函數(shù),該過(guò)程叫插值。一維插值指令：yi=interp1(X,Y,xi,method)yi=interp1(X,Y,xi,method,'extrap')yi=interp1(X,Y,xi,method,extrapval)其中：'extrap'表示當(dāng)xi超過(guò)X的范圍時(shí),執(zhí)行外插算法;extrapval一般取0或nan,表示當(dāng)xi超過(guò)范圍時(shí),外插值取extrapval;method可取'nearest','linear','cubic','spline'等,指示函數(shù)在相鄰已知點(diǎn)間如何取值;X,Y為已知數(shù)據(jù)點(diǎn)的坐標(biāo);xi為需計(jì)算插值的點(diǎn)。4二維插值指令：ZI=interp2(X,Y,Z,XI,YI,'method');例:(網(wǎng)格插值peaks函數(shù).)[X,Y]=meshgrid(-3:.25:3);Z=peaks(X,Y);[XI,YI]=meshgrid(-3:.125:3);ZI=interp2(X,Y,Z,XI,YI);mesh(X,Y,Z),hold,mesh(XI,YI,ZI+15),holdoffaxis([-33-33-520])

例:(由時(shí)間、服務(wù)年限,插值查詢工資.)years=1950:10:1990;service=10:10:30;wage=[150.697199.592187.625179.323195.072250.287203.212179.092322.767226.505153.706426.730249.633120.281598.243];w=interp2(service,years,wage,15,1975)5二維不規(guī)則數(shù)據(jù)的插值z(mì)i=griddata(x,y,z,xi,yi,'v4'); griddata可進(jìn)行三角插值，兩坐標(biāo)向量x,y不必單調(diào);一般xi,yi常為用meshgrid生成的平面網(wǎng)格點(diǎn)坐標(biāo)。用griddata方法處理此類不規(guī)則數(shù)據(jù)最方便.6船在該海域會(huì)擱淺嗎？在某海域測(cè)得一些點(diǎn)(x,y)處的水深z(單位:英尺)由下表給出，水深數(shù)據(jù)是在低潮時(shí)測(cè)得的。船的吃水深度為5英尺，問(wèn)在矩形(75,200)×(-50,150)里的哪些地方船要避免進(jìn)入。問(wèn)題分析:

假設(shè)：該海域海底是平滑的。由于測(cè)量點(diǎn)是散亂分布的，先在平面上作出測(cè)量點(diǎn)的分布圖，在利用二維插值方法補(bǔ)充一些點(diǎn)的水深，然后作出海底曲面圖和等高線圖，并求出水深小于5的海域范圍。7水道水深測(cè)量數(shù)據(jù)(單位:英尺)問(wèn)題求解步驟

1、作出測(cè)量點(diǎn)的分布圖

2、作出海底地貌圖

3、危險(xiǎn)區(qū)域海底地貌圖

4、危險(xiǎn)區(qū)域平面圖x129.0140.0103.588.0185.5195.0105.5Y7.5141.523.0147.022.5137.585.5Z4868688X157.5107.577.081.0162.0162.0117.5Y-6.5-81.03.056.5-66.584.0-33.5Z99889498曲線擬合已知離散點(diǎn)上的數(shù)據(jù)集[(x1,y1),(x2,y2),…(xn,yn)],求得一解析函數(shù)y=f(x),使f(x)在xi上的盡可能接近給定的yi值,這一過(guò)程叫曲線擬合。曲線擬合過(guò)程中f(x)的類型是由實(shí)際背景確定或用戶指定的,即類型已知的,但其中的系數(shù)參量需要確定。確定f(x)系數(shù)的原則是使f(x)盡可能接近已知數(shù)據(jù)點(diǎn)。完成曲線擬合最常用的方法是最小二乘法擬合,即找出使最小的f(x)。由微積分知識(shí)可解出f(x)的系數(shù)參量a。9涉及擬合的常用MATLAB指令[a,SS,Re]=lsqcurvefit(fun,a0,x,y);

----最小二乘曲線擬合(擬合一元函數(shù)f(x),a是其系數(shù)向量,Re是殘差,SS殘差平方和)[a,Re,Jacobian]=nlinfit(X,y,fun,a0);

----非線性最小二乘擬合(擬合多元函數(shù)f(x),X的每一行是對(duì)應(yīng)一個(gè)y的各自變量取值.)p=polyfit(x,y,n)

----多項(xiàng)式擬合(擬合函數(shù)f(x)是n次多項(xiàng)式,x要單調(diào),p為降冪多項(xiàng)式的系數(shù)向量.)曲線擬合示例x=0:.1:5;y=sin(x);p=polyfit(x,y,5)xf=0:0.05:8;yf=polyval(p,xf);plot(xf,sin(xf),'ro',xf,yf,'b-')---------------------[a,ss]=lsqcurvefit(@myfit,[1,1,1,1],x,y)yfit=myfit(a,xf);plot(xf,sin(xf),'ro',xf,yfit,'b-')-----myfit.m-----functiony=myfit(a,x)y=a(1)*cos(a(2)*x+a(3))+a(4);11種群的分布格局實(shí)驗(yàn)植物種群的分布規(guī)律在很大程度上由該種群的生物學(xué)與生態(tài)學(xué)特性決定的。由于不同環(huán)境適合于具有不同適應(yīng)性和需求特性的生物生長(zhǎng),于是在長(zhǎng)期自然選擇和自然協(xié)同過(guò)程中,便形成了各種群特有的空間分布格局。通過(guò)調(diào)查可以得到種群在空間分布各樣方中的種群數(shù),以此可進(jìn)一步推測(cè)、估計(jì)整個(gè)種群的分布格局。為此我們首先要由各樣方中的種群數(shù)統(tǒng)計(jì)出頻數(shù)、頻率密度,然后再估計(jì)出頻率密度函數(shù),進(jìn)而計(jì)算達(dá)到一定種群數(shù)的概率。12分布格局實(shí)驗(yàn)要用的MATLAB函數(shù)findSize、lengthbarpolyfitsumdoubledispformatshortg13各樣方中的種群數(shù)據(jù)a=[541100020000000101...00041000001000102...010000101000020000...01000001002000000...010011210000101001...02000000001000021...000011000200000220...20000011043101120...002000001001201012...10000001332000000...1112040000011000000...00010011302010004...000000000020000001...00000000000000000];n=length(a);fori=min(a):max(a)fn(i+1)=sum(a==i);endf=fn/n;p=[(0:5)',fn',f']bar(p(:,1),p(:,2))bar(0:5,fn)bar(0:5,f)[pc,s]=polyfit(0:5,f,5)symsxxx=[x^5,x^4,x^3,x^2,x,1]fx=sum(pc.*xx)vpa(int(fx,0,4.8),4)14Richards彈性生長(zhǎng)模型1959年Richards提出了一個(gè)生長(zhǎng)模型：上式中y(t)為t時(shí)刻的總生長(zhǎng)量,參數(shù)a為總生長(zhǎng)量的極限值,b為初始值參數(shù),k為生長(zhǎng)速率參數(shù),m是曲線形狀參數(shù)。它的圖形是以y=a為漸近線的S形曲線。當(dāng)m=0時(shí)為Mitscherlich模型當(dāng)m→1時(shí)近似為Compertz模型當(dāng)m=2時(shí)為L(zhǎng)ogistic模型所以Richards模型具有較廣泛的適應(yīng)性。15Richards曲線形狀的計(jì)算機(jī)模擬改變Richards模型的曲線形狀參數(shù)m和初值b(假設(shè)k=0.07,a=100),模擬出其曲線形狀。clearall,clcm=input('注意：m不能為1,取m=');b=input('注意：-1<=b<=1,取b=');a=100;k=0.07;t=0:100;y=a*(1-b*exp(-k*t)).^(1/(1-m));plot(t,y),xlabel('時(shí)間'),ylabel('生長(zhǎng)量')16Richards生長(zhǎng)曲線的拐點(diǎn)Richards生長(zhǎng)曲線的拐點(diǎn)處就是最大生長(zhǎng)速率時(shí)刻。clearall,clcsymsYABkxm;Y=A*(1-B*exp(-k*x))^(1/(1-m))equl=diff(Y,x,2);x=solve(equl,'x')%解符號(hào)方程equl=0,未知量為xY=subs(Y)%符號(hào)表達(dá)式Y(jié)中的x用方程根x代入17用變形蟲(chóng),水稻葉數(shù)據(jù)擬合Richards模型變形蟲(chóng)細(xì)胞重量生長(zhǎng)數(shù)據(jù)時(shí)間01.252.503.755.006.257.508.7510.0011.2512.50重量10.8511.3112.3013.4413.6314.1915.1815.6115.9016.9817.38時(shí)間13.7515.0016.2517.5018.7520.0021.2522.5023.7525.00重量17.7818.6619.1918.7819.2119.1419.7419.9620.0619.91水稻葉伸長(zhǎng)生長(zhǎng)數(shù)據(jù)時(shí)間11.82.63.44.14.85.46.16.87.48.1長(zhǎng)度0.30.50.91.42.53.24.37.610.114.418.5時(shí)間8.89.410.110.811.712.413.114.415.115.7長(zhǎng)度23.025.230.433.738.841.743.744.845.545.3采用曲線擬合方法,估計(jì)出Richards模型中各參數(shù),并進(jìn)行預(yù)測(cè)。18變形蟲(chóng)數(shù)據(jù)擬合的MATLAB程序t=[0,1.25,2.50,3.75,5.00,6.25,7.50,8.75,10.00,11.25,12.50,... 13.75,15.00,16.25,17.50,18.75,20.00,21.25,22.50,23.75,25.00];y=[10.85,11.31,12.30,13.44,13.63,14.19,15.18,15.61,15.90,...16.98,17.38,17.78,18.66,19.19,18.78,...19.21,19.14,19.74,19.96,20.06,19.91];scatter(t,y,5,'r','filled'),holdon;Rechards=inline('b(1)*(1-b(2)*exp(-b(3)*t)).^(1/(1-b(4)))','b','t');b0=[20,-15,0.1,4];[beta,Res,Re]=lsqcurvefit(Rechards,b0,t,y);parameters=betass=sum((y-mean(y)).^2);R=(ss-Res)/sssymsbt;b=beta;y=subs(b(1)*(1-b(2)*exp(-b(3)*t)).^(1/(1-b(4))));ezplot(y,[0,26]),xlabel('時(shí)間h'),ylabel('重量ug')title('變形蟲(chóng)細(xì)胞重量生長(zhǎng)的擬合曲線')19季節(jié)性生物種群消長(zhǎng)模型實(shí)驗(yàn)生物種群實(shí)際增長(zhǎng)過(guò)程可用下面公式描述：

其中N(t)為t時(shí)刻種群密度,K為環(huán)境最大容量,r為種群內(nèi)稟增長(zhǎng)率過(guò)程,a是由K和種群初始密度決定的常數(shù)。f(t)0時(shí)即為L(zhǎng)ogistic模型。對(duì)季節(jié)性生物,其生長(zhǎng)過(guò)程受季節(jié)性生態(tài)因子周期變化影響。因此影響因素f(t)可設(shè)為其中h為影響幅度,

為影響的初相位,為生態(tài)因子變化的角頻率。對(duì)較復(fù)雜情況,h和可以是t的函數(shù)。所以,季節(jié)性生物種群總的增長(zhǎng)規(guī)律可以表示為：20綿蚜蟲(chóng)種群季節(jié)性變化的模擬實(shí)驗(yàn)綿蚜蟲(chóng)種群變化的實(shí)測(cè)數(shù)據(jù)t1234567891011y3940518662085420307433841706110071031620根據(jù)季節(jié)性生物種群總的增長(zhǎng)規(guī)律模型,采用曲線擬合方法,估計(jì)出模型中各參數(shù);并進(jìn)行預(yù)測(cè)。21綿蚜蟲(chóng)數(shù)據(jù)擬合的MATLAB程序t=1:11;y=[3940,5186,6208,5420,3074,3384,1706,1100,710,316,20];scatter(t,y,5,'r','filled'),holdon;Jijie=inline('b(1)*(1+b(2)*cos(b(3)+b(4)*t))./(1+exp(b(5)-b(6)*t))','b','t');b0=[4200,1,1,0.1,0.1,0.20];[beta,Res,Re]=lsqcurvefit(Jijie,b0,t,y);parameters=beta,predictions=y+Ress=sum((y-mean(y)).^2);R=(ss-Res)/sssymsbt;b=beta;y=subs(b(1)*(1+b(2)*cos(b(3)+b(4)*t))./(1+exp(b(5)-b(6)*t)));ezplot(y,[0,11]),xlabel('時(shí)間'),ylabel('數(shù)量')title('棉蚜蟲(chóng)種群消長(zhǎng)實(shí)驗(yàn)數(shù)據(jù)擬合曲線')22線性規(guī)劃線性規(guī)劃模型：23求解線性規(guī)劃的指令[x,fval,exitflag,output,lambda]

=linprog(c,A,b,Aeq,beq,lb,ub,x0,options)注：凸體字的部分是可選項(xiàng)。exitflag結(jié)束標(biāo)志：">0"表示解是收斂的;"=0"表示因到達(dá)最大迭代次數(shù)而結(jié)束;"<0"表示解不收斂！output結(jié)構(gòu)變量：返回迭代次數(shù)、優(yōu)化算法等信息。lambda結(jié)構(gòu)變量：返回Lagrange乘子信息。fval：返回最優(yōu)值。不需要某項(xiàng)輸入?yún)?shù)時(shí)，可用空方括號(hào)[]占據(jù)其原有位置。24應(yīng)用舉例解：

c=[-5,-4,-6] A=[1-11 324 320]; b=[20;42;30]; lb=zeros(3,1); [x,fval,exitflag,output,lambda]=linprog(c,A,b,[],[],lb)Min25飼料配方的線性規(guī)劃實(shí)驗(yàn)背景：為了保證畜禽生產(chǎn)的質(zhì)量和高產(chǎn),應(yīng)配制完善而平衡的飼料。一般每種飼養(yǎng)動(dòng)物都有飼養(yǎng)標(biāo)準(zhǔn)，規(guī)定動(dòng)物飼料中必須達(dá)到的營(yíng)養(yǎng)標(biāo)準(zhǔn);通過(guò)調(diào)查也可知道各種飼料原料的價(jià)格和所含營(yíng)養(yǎng)成分量。問(wèn)題：如何在滿足飼料的營(yíng)養(yǎng)標(biāo)準(zhǔn)前提下,配制出成本最低的飼料。方法：利用線性規(guī)劃解決該類配方問(wèn)題。26雞飼料營(yíng)養(yǎng)標(biāo)準(zhǔn)和原料限制情況雞飼料營(yíng)養(yǎng)標(biāo)準(zhǔn)代謝能兆卡/千克粗蛋白克/千克量比蛋白質(zhì)能克/兆卡粗纖維克/千克賴氨酸克/千克蛋氨酸克/千克鈣克/千克磷克/千克食鹽克/千克2.7~2.8135~14551<505.62.523~404.6~6.53.7每1000千克飼料中,各原料用量限制(適口性)玉米小麥麥米糠豆餅菜籽餅魚(yú)粉槐葉粉碳酸鈣400~500100~150100~20015010030~505030~502527雞飼料原料的營(yíng)養(yǎng)成分及價(jià)格編號(hào)(千克)品名單價(jià)元/千克代謝能千卡/千克粗蛋白克/千克粗纖維克/千克賴氨酸克/千克蛋氨酸克/千克鈣克/千克磷克/千克食鹽克/千克X1玉米0.2423.3078162.31.20.70.3X2小麥0.2723.08114223.41.70.60.34X3麥0.101.78142956.02.30.310.0X4米糠0.102.10117726.52.71.013.0X5豆餅0.172.504024924.15.13.25.0X6菜籽餅0.1341.623601138.17.15.38.4X7魚(yú)粉0.822.00450029.111.86327X8槐葉粉0.201.6117010810.62.24.04.0X9DL蛋氨酸12.00980X10骨粉0.192300140X11碳酸鈣0.52400X12食鹽0.18100028配制1千千克飼料的線性規(guī)劃模型(目標(biāo)函數(shù))(總重量約束)(代謝能約束)(粗蛋白約束)(粗纖維約束)(賴氨酸約束)(蛋氨酸約束)(鈣約束)(磷約束)(原料用量約束)DEA評(píng)價(jià)模型簡(jiǎn)介DEA是對(duì)一些同類型的部門(mén)或單位(稱為決策單元DMU)進(jìn)行相對(duì)有效性評(píng)價(jià)的一種方法.一個(gè)決策單元的有效性是用該單元的多指標(biāo)輸出的加權(quán)和與多指標(biāo)輸入的加權(quán)和之比來(lái)定義的.對(duì)n個(gè)決策單元

，向量xi=[x1i,x2i,…,xmi]T和向量yi=[y1i,y2i,…,ysi]T,分別表示第i個(gè)DMU的m項(xiàng)輸入和s項(xiàng)輸出,X=[x1,x2,…,xn]和Y=[y1,y2,…,yn]分別稱為多指標(biāo)輸入矩陣和輸出矩陣.2930DEA模型C2R12…

…n

…x1x2

…

…xn y1y2

…

…yn

…

評(píng)價(jià)DMU－j0的DEA模型（C2R）為分式規(guī)劃（h0為DMU－j0的效率指數(shù)）：分式規(guī)劃(其中x0=xj0,y0=yj0,1≤j0≤n．31

用1962年Charnes和Cooper對(duì)于分式規(guī)劃的Charnes－Cooper變換（稱為C2－變換）:32Matlab線性規(guī)劃標(biāo)準(zhǔn)形式DEA評(píng)價(jià)示例Dea程序%五個(gè)決策單元，3個(gè)投入性指標(biāo)X=[3060554070;2540703090;13015012070180]%五個(gè)決策單元，2個(gè)產(chǎn)出性指標(biāo)Y=[3543765263;6080534271]k=5%計(jì)算第k個(gè)單元的效率s=size(X,1);t=size(Y,1);n=size(X,2);%s=3t=2n=5c=[zeros(1,s),-Y(:,k)'];A=[-X',Y'];b=zeros(n,1);Aeq=[X(:,k)',zeros(1,t)];beq=1;lb=zeros(s+t,1);[x,f]=linprog(c,A,b,Aeq,beq,lb);w=x(1:s),v=x(s+1:end),h=-f%上述程序可設(shè)計(jì)為一個(gè)求解Dea的指令函數(shù)：%[w,v,h]=dea(X,Y,k)34DEA交叉評(píng)價(jià)由于前面提到的DEA模型的決策單元效率是利用對(duì)其最有利的權(quán)重計(jì)算出來(lái)的(自我評(píng)價(jià))。在實(shí)際問(wèn)題中，往往有較多的決策單元都能取到最大的效率值1，從而不能區(qū)分這些決策單元的優(yōu)劣。為了解決這個(gè)問(wèn)題，我們引入交叉評(píng)價(jià)機(jī)制.用每一個(gè)DMUi的最佳權(quán)重去計(jì)算其它DMUk的效率值35由于線性規(guī)劃的最優(yōu)解ui*和vi*不唯一，得出的交叉評(píng)價(jià)值Eik具有不確定性。對(duì)抗型交叉評(píng)價(jià)(1)先用自我評(píng)價(jià)模型得出每一個(gè)DMUi的自我評(píng)價(jià)值Eii;(2)在保證DMUi得到最大值Eii的前提下，使其它的DMUk得到盡可能小的交叉評(píng)價(jià)值Eik。對(duì)抗型交叉評(píng)價(jià)的實(shí)質(zhì)是:每一個(gè)DMUi,在盡可能抬高自己的前提下，盡可能地貶低其它DMUk。36對(duì)抗型交叉評(píng)價(jià)模型算法37將E的第i列的平均值ei作為衡量DMUi的優(yōu)劣的一項(xiàng)指標(biāo);ei可視為諸決策單元對(duì)DMUi的總的評(píng)價(jià),越大越好。對(duì)抗型交叉評(píng)價(jià)模型程序functionE=deaEij(X,Y)%輸入指標(biāo)矩陣X，輸出指標(biāo)矩陣Y，ei=mean(E)s=size(X,1);t=size(Y,1);n=size(X,2);E=zeros(n);A=[-X',Y'];b=zeros(n,1);lb=zeros(n,1);options=optimset('Largescale','off');fori=1:nc=[zeros(1,s),-Y(:,i)'];Aeq=[X(:,i)',zeros(1,t)];beq=1;[x,f]=linprog(c,A,b,Aeq,beq,lb);w=x(1:s);v=x(s+1:end);E(i,i)=-f;x0=x;fork=1:nifk~=ic=[zeros(1,s),Y(:,k)'];Aeq=[E(i,i)*X(:,i)',-Y(:,i)';X(:,k)',zeros(1,t)];beq=[0;1];[x,E(i,k),ex]=linprog(c,A,b,Aeq,beq,lb,[],x0,options);endendend3839實(shí)驗(yàn)指標(biāo)銀行員工人數(shù)((投入)機(jī)構(gòu)總數(shù)(投入)營(yíng)業(yè)費(fèi)用率(投入)總資產(chǎn)(億)(投入)人均利潤(rùn)率(萬(wàn)元)(產(chǎn)出)資產(chǎn)收益率(產(chǎn)出)工商銀行i=1381713164760.55466586842.8830.1572120.010126農(nóng)業(yè)銀行i=2447519244520.78563160501.277.275890.00208中國(guó)銀行i=3237379101450.46610159955.5337.8950960.01095建設(shè)銀行i=4298868134480.54189665981.7733.7326180.011479浦發(fā)銀行i=5141284080.5840729149.803574.148790.006855招商銀行i=6289715790.4086613105.5272.596730.013582華夏銀行i=793902870.7296145923.382740.690030.010898興業(yè)銀行i=8118513900.5081238513.352792.061180.016909浙商銀行i=91372210.217699575.2279940.1203350.007526上海銀行i=1047681880.5139033089.860179.6233010.010262上海農(nóng)村商業(yè)銀行i=1141963270.2685721574.727517.189590.002014xi1,xi2,xi3,xi4yi1,yi2綜合作業(yè)1作業(yè)1、第34頁(yè)的DEA程序改造為用戶函數(shù)： [w,v,h]=dea(X,Y,k)輸入?yún)?shù)為：多指標(biāo)輸入矩陣Xmxn和輸出矩陣Ysxn以及單元標(biāo)號(hào)k。輸出參數(shù)為：輸入指標(biāo)的權(quán)重w和輸出指標(biāo)的權(quán)重v以及第k單元的效率h.作業(yè)2、用對(duì)抗性交叉評(píng)價(jià)程序(第39頁(yè))來(lái)測(cè)算第十五頁(yè)上各單位的效率高低。41非線性規(guī)劃數(shù)學(xué)模型：42求解非線性規(guī)劃的指令[x,fval,exitflag,output,lambda,grad,hessian]=fmincon('fun',x0,A,b,Aeq,beq,lb,ub,'nonlcon',options,P1,P2,...)其中nonlcon為非線性約束及導(dǎo)數(shù)：[C,Ceq,gC,gCeq]=nonlcon(x)；

而fun則為目標(biāo)函數(shù)及導(dǎo)數(shù)和Hess矩陣：

[f,g,H]=myfun(x)；fval返回最優(yōu)目標(biāo)值。注：凸體字的部分是可選項(xiàng)；輸入空矩陣[]則表示該項(xiàng)不需要輸入或用缺省值。43應(yīng)用舉例

function[f,g,H]=myoptfun(x)%目標(biāo)函數(shù)(包括函數(shù)值,梯度,HESS矩陣)f=exp(x(1))*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+2*x(2)+1);ifnargout>1%funcalledwithtwooutputargumentsg1=exp(x(1))*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+8*x(1)+6*x(2)+1);g2=exp(x(1))*(4*x(2)+4*x(1)+2);g=[g1,g2];%Gradientofthefunctionevaluatedatxifnargout>2h11=exp(x(1))*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+16*x(1)+10*x(2)+9);h12=exp(x(1))*(4*x(2)+4*x(1)+6);h21=h12;h22=4*exp(x(1));H=[h11,h12;h21,h22];%Hessianevaluatedatxendend44約束條件、MATLAB的解法function[c,ceq,GC,GCeq]=mycon(x)c1=1.5+x(1)*x(2)-x(1)-x(2);c2=-x(1)*x(2)-10;c=[c1;c2];%Nonlinearinequalitiesatxceq=[];%Nonlinearequalitiesatxifnargout>2%nonlconcalledwith4outputsGC1=[x(2)-1,x(1)-1];GC2=[-x(2),-x(1)];GC=[GC1;GC2];%GradientsoftheinequalitiesGCeq=[];%Gradientsoftheequalitiesend---------------------------------------------------------------------------------------------------------------------x0=[10;10];lb=[0;0];ub=[inf;inf];options=optimset('Algorithm','active-set','Display','off');[x,fval]=fmincon(@myoptfun,x0,[],[],[],[],lb,ub,@mycon,options)45優(yōu)化選項(xiàng)optionsoptions.DerivativeCheck[on|{off}]檢查梯度輸入options.Diagnostics[on|{off}]

options.Display[off|iter|{final}]結(jié)果等信息何時(shí)顯示★options.DiffMaxChange[正數(shù)|{1e-1}]導(dǎo)數(shù)的最大改變options.DiffMinChange[正數(shù)|{1e-8}]導(dǎo)數(shù)的最小改變options.GradConstr[on|{off}]約束求梯度options.GradObj[on|{off}]目標(biāo)求梯度options.LargeScale[{on}|off]大規(guī)模算法★options.MaxFunEvals[正整數(shù)]

options.MaxIter[正整數(shù)]最大迭代次數(shù)options.TolFun[正數(shù)]函數(shù)值終止誤差標(biāo)準(zhǔn)★options.Tolx[正數(shù)]x的終止誤差標(biāo)準(zhǔn)★46用optimset函數(shù)設(shè)置options參數(shù)options=optimset('param1',value1,'param2',value2,…)

%建立options結(jié)構(gòu)中的各參數(shù)。options=optimset(oldopts,'param1',value1',…)

%修改oldopts中的參數(shù)param1等。例：options=optimset('Display','iter','TolFun',1e-8)無(wú)條件極值fminunc

/fminsearch問(wèn)題：求minf(x)?求解指令：[x,fval]=fminunc(fun,x0,options)[x,fval]=fminsearch(fun,x0,options)47二次規(guī)劃模型指令

[x,fval,exitflag,output,lambda]

=quadprog(H,f,A,b,Aeq,beq,lb,ub,x0,options)應(yīng)用舉例解:H=[2,0;0,2];c=[-4,0];A=[-1,1;1,-1];b=[2;-1];

[x,fval]=quadprog(H,c,A,b,[],[],zeros(2,1)),mf=fval+4組合投資計(jì)算設(shè)一組投資由m個(gè)投資項(xiàng)目組成，它們的單項(xiàng)期望回報(bào)率為(r1,r2,…,rm)，該m個(gè)項(xiàng)目的投資比例為(x1,x2,…,xm)，則該組投資的總回報(bào)率的期望值R：投資組合的總回報(bào)率的期望值描述了多項(xiàng)投資的總體平均回報(bào)水平。同樣地。僅漢只用總回報(bào)率期望值來(lái)描述投資組合的效果是不夠的，還需描述總回報(bào)率的離散趨勢(shì)，也就是整個(gè)投資組合風(fēng)險(xiǎn)的大小。一組投資的總回報(bào)率的風(fēng)險(xiǎn)（或離散趨勢(shì)）的常用測(cè)度是總回報(bào)率的方差和標(biāo)準(zhǔn)方差?？偦貓?bào)率R的方差總回報(bào)率R的方差的估算公式其中：為該m個(gè)項(xiàng)目的投資比例，是該m個(gè)項(xiàng)目的各單項(xiàng)回報(bào)率的標(biāo)準(zhǔn)差，為第i個(gè)投資項(xiàng)目以與第j個(gè)投資項(xiàng)目的相關(guān)系數(shù)，V為各單項(xiàng)回報(bào)率間的協(xié)方差矩陣。投資者的目標(biāo)是獲得大的投資回報(bào)和承擔(dān)小的投資風(fēng)險(xiǎn)。投資組合優(yōu)化模型就是確定一組投資項(xiàng)目的最優(yōu)投資比例,在該投資組合的總回報(bào)率的方差不超過(guò)某個(gè)可接受的值的約束下(即在可接受的風(fēng)險(xiǎn)水平下)，使得總回報(bào)率的期望值最大(即投資回報(bào)水平最高)；或者在投資組合的回報(bào)率的期望值不低于某個(gè)所要求的值的約束下(即在所要求的投資回報(bào)水平)，使得總回報(bào)率的方差最小(即投資風(fēng)險(xiǎn)最小)。投資組合優(yōu)化模型多目標(biāo)優(yōu)化模型：總回報(bào)率最大模型：總風(fēng)險(xiǎn)最小模型：有效前沿Markowitz資產(chǎn)組合理論就是尋找一個(gè)能使在同樣收益水平下,總風(fēng)險(xiǎn)水平最低的有效組合。不同的收益與相應(yīng)的最小風(fēng)險(xiǎn)構(gòu)成有效前沿(R,

2).H=2*V;n=length(r);c=zeros(1,n);Aeq=[r;ones(1,n)];beq=[R;1];[x,fval]=quadprog(H,c,[],[],Aeq,beq,zeros(n,1))有效前沿計(jì)算示例例考慮一個(gè)三資產(chǎn)組合,分別為資產(chǎn)1、資產(chǎn)2與資產(chǎn)3，其預(yù)期收益率分別為0.2、0.1、0.15。協(xié)方差如下表資產(chǎn)1資產(chǎn)2資產(chǎn)3資產(chǎn)10.01-0.00610．0042資產(chǎn)2-0.00610.04-0.0252資產(chǎn)30.0042-0.02520.0225求該資產(chǎn)組合有效前沿。示例程序clearall;clcR=0.1:0.01:0.2;V=[0.01,-0.0061,0.0042;-0.0061,0.04,-0.0252;0.0042,-0.0252,0.0225];r=[0.2,0.1,0.15];H=2*V;n=length(r);c=zeros(1,n);Aeq=[r;ones(1,n)];lb=zeros(n,1);options=optimset('Display','off','Algorithm','active-set');fork=1:length(R)beq=[R(k);1];[X(:,k),s2(k),ex(k)]=quadprog(H,c,[],[],Aeq,beq,lb,[],[],options);endplot(R,s2)練習(xí)投資三種資產(chǎn)的收益率相關(guān)系數(shù)矩陣、預(yù)期回報(bào)、標(biāo)準(zhǔn)差如表所示：資產(chǎn)A資產(chǎn)B資產(chǎn)C資產(chǎn)A10.80.4資產(chǎn)B0.810.3資產(chǎn)C0.40.31預(yù)期回報(bào)0.10.150.20各資產(chǎn)標(biāo)準(zhǔn)差0.20.250.18試給出有效前沿。相關(guān)陣到協(xié)方差陣的轉(zhuǎn)換指令corr2cov函數(shù)用法：

Covariances=corr2cov(STDs,Correlations);57非線性方程組求數(shù)值解MATLAB求‘Fun(x)=0’解的指令一般格式：

[x,fval,exitflag,output]=fsolve(Fun,x0,options,p1,p2,...)例：求解方程組sinx+y+z2ex-4=0,x+yx=0,xyz=0的具體步驟如下 創(chuàng)建函數(shù)funs.m functionF=funs(x) F(1)=sin(x(1))+x(2)+x(3)^2*exp(x(1))-4; F(2)=x(1)+x(1)*x(2);F(3)=x(1)*x(2)*x(3);

解方程

[x,fv]=fsolve(@funs,[1,1,1])58fslove指令補(bǔ)充說(shuō)明返回參數(shù)exitflag為求解結(jié)束狀態(tài)標(biāo)志：

>0----解收斂;=0----達(dá)到迭代次數(shù);<0----解發(fā)散。輸入?yún)?shù)Fun可為向量函數(shù)(表示方程組)，其定義形式為：

functionF=Fun(x)

F=...%Computefunctionvaluesatx

或

function[F,J]=Fun(x)

F=...%objectivefunctionvaluesatx

ifnargout>1%twooutputarguments

J=...%Jacobianofthefunctionevaluatedatx

endKMV信用監(jiān)控模型KMV模型是美國(guó)舊金山市KMV公司用來(lái)估計(jì)借款企業(yè)違約概率的方法。該模型認(rèn)為，在債務(wù)到期日，如果公司資產(chǎn)的市場(chǎng)價(jià)值高于公司債務(wù)值（違約點(diǎn)），則公司股權(quán)價(jià)值為公司資產(chǎn)市場(chǎng)價(jià)值與債務(wù)值之間的差額；如果此時(shí)公司資產(chǎn)價(jià)值低于公司債務(wù)值，則公司變賣(mài)所有資產(chǎn)用以償還債務(wù)，股權(quán)價(jià)值變?yōu)榱恪?9KMV模型的運(yùn)用步驟首先,利用Black-Scholes期權(quán)定價(jià)公式,根據(jù)企業(yè)資產(chǎn)的市場(chǎng)價(jià)值、資產(chǎn)價(jià)值的波動(dòng)性、到期時(shí)間、無(wú)風(fēng)險(xiǎn)借貸利率及負(fù)債的賬面價(jià)值估計(jì)出企業(yè)股權(quán)的市場(chǎng)價(jià)值及其波動(dòng)性。其次根據(jù)公司的負(fù)債計(jì)算出公司的違約實(shí)施點(diǎn)(defaultexercisepoint,為企業(yè)1年以下短期債務(wù)的價(jià)值加上未清償長(zhǎng)期債務(wù)賬面價(jià)值的一半),計(jì)算借款人的違約距離。最后,根據(jù)企業(yè)的違約距離與預(yù)期違約率(EDF)之間的對(duì)應(yīng)關(guān)系，求出企業(yè)的預(yù)期違約率。60違約距離DD計(jì)算61式中：DP(DefaultPoint)為違約點(diǎn)值DP=短期負(fù)債(STD)+0.5長(zhǎng)期負(fù)債(LTD)式中:E為企業(yè)股權(quán)市場(chǎng)價(jià)值,V為企業(yè)資產(chǎn)市場(chǎng)價(jià)值,D為企業(yè)債務(wù)面值,r為無(wú)風(fēng)險(xiǎn)收益率,τ為債務(wù)償還期限,N(.)為標(biāo)準(zhǔn)累積正態(tài)分布函數(shù),

V為企業(yè)資產(chǎn)價(jià)值波動(dòng)率,

E為企業(yè)股權(quán)市場(chǎng)價(jià)值波動(dòng)率。其中V和

V是未知量,解非線性方程組可得。據(jù)違約距離DD的定義可知,理論上發(fā)生違約的概率為1-N(DD).而基于違約數(shù)據(jù)庫(kù),依據(jù)違約距離可以映射出公司實(shí)際的期望違約頻率EDF。由于我國(guó)當(dāng)前還沒(méi)有公開(kāi)的違約數(shù)據(jù)庫(kù)可以使用，以違約距離DD作為上市公司信用評(píng)價(jià)的依據(jù)。上市公司股權(quán)市場(chǎng)價(jià)值

中國(guó)證券市場(chǎng)發(fā)展歷史較為特殊，上市公司股票被人為分割為上市流通股票和暫不上市流通股票兩種。在計(jì)算上市公司股權(quán)市場(chǎng)價(jià)值時(shí)需要考慮以什么樣的價(jià)格來(lái)計(jì)算非流通股市場(chǎng)價(jià)值/由于非流通股沒(méi)有市場(chǎng)交易價(jià)格，因此如何給非流通股定價(jià)是一件困難的事情。參考上市公司股票全流通研究中非流通股定價(jià)，以每股凈資產(chǎn)計(jì)算非流通股的價(jià)格。流通股市場(chǎng)價(jià)值=12月份周平均收盤(pán)價(jià)格×流通股股數(shù)非流通股市場(chǎng)價(jià)值=每股凈資產(chǎn)×非流通股股數(shù)上市公司股權(quán)市場(chǎng)價(jià)值=流通股市場(chǎng)價(jià)值+非流通股市場(chǎng)價(jià)值公司債務(wù)面值D為公司財(cái)務(wù)年報(bào)中總負(fù)債面值?？紤]到數(shù)據(jù)和工作量的限制，我們?cè)O(shè)定違約距離的計(jì)算時(shí)間為一年，即τ=1。無(wú)風(fēng)險(xiǎn)利率使用中國(guó)人民銀行公布的一年期定期整存整取的存款利率。KMV測(cè)算示例某公司流動(dòng)負(fù)債為1億元，長(zhǎng)期負(fù)債為5000萬(wàn)元，根據(jù)上市公司的股價(jià)行情表，可以統(tǒng)計(jì)出公司股權(quán)價(jià)值E和公司股權(quán)價(jià)值波動(dòng)率

E。試計(jì)算公司的違約概率。月份總市值(元)收益率(%)月份總市值(元)收益率(%)1129523558-13.6581409724643.31214988546213.5891484050955.013142316387-5.3210144898861-2.4241494409124.77111449046090.005147924524-1.0312130292794-11.216130439432-13.40均值14127642771363130244.31標(biāo)準(zhǔn)差8.35(1)基本參數(shù)計(jì)算公司股價(jià)波動(dòng)率為8.35%，公司股權(quán)價(jià)值波動(dòng)率公司股權(quán)價(jià)值E=141276427(2)求解非線性方程組，解出Va和

a

。KMV測(cè)算示例(續(xù)-程序)KMV計(jì)算示例的Matlab程序表示聯(lián)列方程組的向量函數(shù)KMVe1fun.m：functionF=KMVe1fun(x,E,D,Esigma,r,T)%Va=x(1),Asigma=x(2)d1=(log(x(1)/D)+(r+0.5*x(2)^2))/(x(2)*sqrt(T));d2=d1-x(2)*sqrt(T);f1=x(1)*normcdf(d1)-D*exp(-r*T)*normcdf(d2)-E;f2=normcdf(d1)*x(1)*x(2)-E*Esigma;F=[f1;f2];KMV測(cè)算示例(續(xù)-程序執(zhí)行)求解的指令：%KMV計(jì)算示例1%r:無(wú)風(fēng)險(xiǎn)利率r=0.0225;%T:債務(wù)期限T=1;%輸入年數(shù)%DP:違約點(diǎn)t%SD:短期債務(wù),LD:長(zhǎng)期債務(wù)SD=1e8;%輸入LD=50000000;%輸入%計(jì)算違約點(diǎn)DP=SD+0.5*LD;%D:債務(wù)的市場(chǎng)價(jià)值D=DP;%債務(wù)的市場(chǎng)價(jià)值,可改為公司年報(bào)公布的總負(fù)債面值%sigma:波動(dòng)率sigma=0.0835;%(輸入）%Esigma:總波動(dòng)率Esigma=sigma*sqrt(12);%E:公司股權(quán)價(jià)值E=141276427;%開(kāi)始計(jì)算VaandAsigma;若exitflag<0,重新設(shè)定x0再算一遍！x0=[100000000;1];%x(1)的初值應(yīng)接近E，x(2)的初值應(yīng)接近Esigma[x,f,exitflag]=fsolve(@KMVe1fun,x0,[],E,D,Esigma,r,T)Va=x(1),Asigma=x(2)%計(jì)算違約距離DD=(Va-DP)/(Va*Asigma)%計(jì)算違約率EDF=normcdf(-DD)混沌分形簡(jiǎn)介與實(shí)驗(yàn)函數(shù)迭代與不動(dòng)點(diǎn)給定一函數(shù)f(x)以及初始點(diǎn)x0，定義數(shù)列

xk+1=f(xk)

,k=0,1,2,…

{xk}稱為函數(shù)f(x)的迭代序列。函數(shù)f的不動(dòng)點(diǎn)：滿足f(x)=x的點(diǎn)x稱為f的不動(dòng)點(diǎn)。吸引點(diǎn)與排斥點(diǎn)：附近的點(diǎn)經(jīng)函數(shù)f迭代后都趨向f的某一不動(dòng)點(diǎn)x，則稱x為吸引點(diǎn)；附近的點(diǎn)經(jīng)函數(shù)f迭代后都遠(yuǎn)離不動(dòng)點(diǎn)x，則稱它是排斥點(diǎn)。周期點(diǎn)與預(yù)周期點(diǎn)K-周期點(diǎn)：若f(u1)=u2,f(u2)=u3,…,f(uk)=u1;則點(diǎn)集u1,u2,u3,…,uk形成一個(gè)k循環(huán);u1稱為k-周期點(diǎn);k稱為周期。類似地，周期點(diǎn)也可以分吸引點(diǎn)與排斥點(diǎn)。預(yù)周期點(diǎn):如果點(diǎn)x最終歸宿于某個(gè)循環(huán)中，則稱它為預(yù)周期點(diǎn)。如1是x2-1的預(yù)周期點(diǎn)。注：迭代序列{xk}的收斂與發(fā)散性質(zhì)不僅與函數(shù)f有關(guān)，而且與初值的選擇有關(guān)。例如，對(duì)于迭代xk+1=2xk+1,當(dāng)初值x0=-1時(shí)，迭代序列收斂，否則發(fā)散。吸引與排斥點(diǎn)判斷//迭代函數(shù)y=ax(1-x)定理:設(shè)x是f(x)的不動(dòng)點(diǎn),如果在x附近有|f'(x)|<1,則x是f(x)的吸引不動(dòng)點(diǎn);否則,x是f(x)的排斥不動(dòng)點(diǎn)。對(duì)迭代函數(shù)f(x)=ax(1-x)(0<a4)而言：不動(dòng)點(diǎn)：x=ax(1-x)x=0或x=(a-1)/a吸引/排斥點(diǎn)：f’(0)=a,f’((a-1)/a)=2-a，故當(dāng)0<a<1時(shí)，0為吸引點(diǎn)，(a-1)/a為排斥點(diǎn)。當(dāng)1<a<3,0為排斥點(diǎn)，(a-1)/a為吸引點(diǎn)。2-周期點(diǎn)：x2=f(x1),x1=f(x2);即x=f(f(x))的解。

得x(1)=0,x(2)=(a-1)/a

x(3,4)=(1+asqrt(a2-2a-3))/(2a),(a3)70Matlab模擬迭代引起的混沌現(xiàn)象對(duì)迭代函數(shù)xk+1=axk(1-xk)進(jìn)行模擬clearall,clca=input('a=');x(1)=input('x(1)=');fori=2:200x(i)=a*x(i-1)*(1-x(i-1));endt=1:200;scatter(t,x,3,'r','filled')%plot(x,'r.');title('混沌現(xiàn)象觀察');xlabel('迭代次數(shù)'),ylabel('序列{x(n)}')functionx=itfun(a,x0)x=zeros(1,200);x(1)=x0;fori=2:200x(i)=a*x(i-1)*(1-x(i-1));endstr=['混沌現(xiàn)象觀察:','a=',num2str(a),',x0=',num2str(x0)];t=1:200;scatter(t,x,3,'r','filled');title(str);xlabel('迭代次數(shù)'),ylabel('序列{x(n)}')---------------x=itfun(0.5,1.4);思考與實(shí)驗(yàn)(y=ax(1-x))1〉對(duì)幾組不同的參數(shù)值a(如a=0.5,a=1.4以及不同的初值x0，觀察迭代是否收斂。2〉取參數(shù)a=0.8，用不同的初值x0做迭代。你能找到一個(gè)吸引的不動(dòng)點(diǎn)嗎？一個(gè)排斥的不動(dòng)點(diǎn)嗎？哪些初值收斂到吸引的不動(dòng)點(diǎn)？哪些初值使序列發(fā)散？取不動(dòng)的參數(shù)a=1,1.6,2,2.5回答同樣的問(wèn)題。3〉找出一個(gè)參數(shù)a,使它對(duì)應(yīng)的迭代具有2周期點(diǎn)。這種性質(zhì)依賴于初值嗎?4〉對(duì)任意的整數(shù)k，你能找到一個(gè)a值使得它對(duì)應(yīng)的迭代具有k周期點(diǎn)嗎？對(duì)哪些k值能給出k周期點(diǎn)？在每種情況下，結(jié)果是否依賴于初值？(對(duì)3.4a

3.6和3.6a

4的值進(jìn)行驗(yàn)證)72混沌實(shí)驗(yàn)初探迭代公式在方程求解等數(shù)值計(jì)算中,有著重要的地位。下面對(duì)一個(gè)迭代公式(也可看作差分方程),討論其解的穩(wěn)定性。對(duì)迭代公式：yk+1=ayk(1-yk)1)當(dāng)1<a<3時(shí),方程有穩(wěn)定的平衡點(diǎn),y收斂到1-1/a。2)當(dāng)3<a<3.499時(shí),有2個(gè)平衡點(diǎn),出現(xiàn)2倍周期現(xiàn)象。從k的連續(xù)兩期看,y的變化是穩(wěn)定的。3)當(dāng)3.499<a<3.544時(shí),有4個(gè)平衡點(diǎn),出現(xiàn)4倍周期。4)當(dāng)a>3.57時(shí),便不出現(xiàn)任何2n倍周期現(xiàn)象,{yn}的變化趨勢(shì)是一片混亂,這就是常說(shuō)的混沌現(xiàn)象。混沌的特性:①對(duì)初值敏感.初值雖近,迭代后散開(kāi).②混沌不是隨機(jī).不同初值,x值落在某區(qū)間頻率不定.73混沌現(xiàn)象模擬的MATLAB程序clearall,clca=input('a=');x(1)=input('x(1)=');fori=2:200x(i)=a*x(i-1)*(1-x(i-1));endt=1:200;scatter(t,x,3,'r','filled')%plot(x,'r.');title('混沌現(xiàn)象觀察');xlabel('迭代次數(shù)'),ylabel('序列{x(n)}')74分枝與混沌(采用系列不同參數(shù)a觀測(cè))n0=50;n=100;x0=0.5;clf;holdon;fora=1.7:0.01:3.7x(1)=x0;fork=1:100x(k+1)=a*x(k)*(1-x(k));endplot(a*ones(1,n-n0+1),x(n0+1:n+1),'k.')%scatter(a*ones(1,n-n0+1),x(n0+1:n+1),2,'k','filled')title('用x=ax(1-x)產(chǎn)生的分枝與混沌圖')xlabel('參數(shù)a');ylabel('迭代式的根x');end;holdoff;75迭代結(jié)果的分布情況x0=input('請(qǐng)輸入初值(0<x0<1)：x0=');x=x0;n1=0;n2=0;n=200;fori=1:nx=4*x*(1-x);y1(i)=x;if(0<=x&x<1/2)n1=n1+1;elsen2=n2+1;endendplot(y1,'-k.');xlabel('迭代次數(shù)n');ylabel('x值');disp(['x0=',num2str(x0),'p=n1/n=',num2str(n1/n)]);disp(['x值落在0~1/2中的次數(shù)n1=',num2str(n1)]);disp(['x值落在0~1/2外的次數(shù)n2=',num2str(n2)]);76分形幾何與分形藝術(shù)什么是分形幾何？通俗的說(shuō)法就是研究無(wú)限復(fù)雜但具有一定意義下自相似圖形和結(jié)構(gòu)的幾何學(xué)。"分形"譯于英文Fractal,詞本身具有"破碎"、"不規(guī)則"等含義。分形幾何創(chuàng)始人Mandelbrot最精彩的部分研究是1980年發(fā)現(xiàn)了Mandelbrot集合，他發(fā)現(xiàn)整個(gè)宇宙以一種出人意料的方式構(gòu)成自相似的結(jié)構(gòu)。Mandelbrot集合圖形的邊界處，具有無(wú)限復(fù)雜和精細(xì)的結(jié)構(gòu)。當(dāng)你放大某個(gè)區(qū)域，它的結(jié)構(gòu)就在變化，展現(xiàn)出新的結(jié)構(gòu)元素。用數(shù)學(xué)方法對(duì)放大區(qū)域進(jìn)行著色處理，這些區(qū)域就變成一幅幅精美的藝術(shù)圖案，這些藝術(shù)圖案人們稱之為"分形藝術(shù)"。

77復(fù)平面中的神奇迭代-Julia集合在復(fù)平面上,x軸代表實(shí)數(shù),y軸代表虛數(shù)。現(xiàn)在復(fù)平面上任意取一個(gè)點(diǎn)z,將其代入下面方程中進(jìn)行反復(fù)迭代運(yùn)算：

zn+1=zn2+c最后得到的z值有3種可能性：

①z值沒(méi)有界限增加（趨向無(wú)窮）

②z值衰減（趨向于零）

③z值是變化的,即非1或非2趨向無(wú)窮和趨向于零的點(diǎn)叫定常吸引子，很多點(diǎn)在定常吸引子處結(jié)束，被定常吸引子所吸引。非趨向無(wú)窮和趨向于零的點(diǎn)是"Julia集合"部分，也叫混沌吸引子。Mandelbrot集合是所有的Julia集合的合并。78得到"Julia集合"的近似算法一般按下述算法近似計(jì)算"Julia集合"：n=0;while((n++<Nmax)&&((Real(Z)^2+Imag(Z)^2)<Rmax)){Z=Z*Z+C;}其中：Nmax為最大迭代次數(shù),Rmax為逃離界限。由(Real(Z)^2+Imag(Z)^2)>=Rmax退出循環(huán)時(shí),相當(dāng)于"①z值沒(méi)有界限增加(趨向無(wú)窮)",為定常吸引子,把這些區(qū)域著成白色。由(n>=Nmax)退出循環(huán)時(shí),相當(dāng)于"②z值衰減(趨向于零)"或"③z值是變化的",把這些區(qū)域著成黑色。黑色區(qū)域圖形的邊界處即為"Julia集合"。"Julia集合"有著極其復(fù)雜的形態(tài)和精細(xì)的結(jié)構(gòu)。79生成"Julia集合"的MATLAB程序M=zeros(256,256);r=linspace(-1,1,256);Nmax=100;Rmax=1000000000;c=-1.1;fora=1:length(r)forb=1:length(r)z=r(a)+r(b)*i;n=0;while(n<Nmax)&(abs(z)<Rmax)z=z*z+c;n=n+1;endifn>=NmaxM(a,b)=1;elseM(a,b)=2;endendendmap=[0,0,0;1,1,1];colormap(map);image([-1,1],[-1,1],M)生成"Julia集合"c=0.302-0.577i;%設(shè)置對(duì)參數(shù)c的缺省值x=0;y=0;L=2;Nmax=100;%設(shè)置最大迭代次數(shù)M=100;%設(shè)置一個(gè)大數(shù)，表示無(wú)窮大[x1,y1]=meshgrid(linspace(x-L,x+L,512),linspace(y-L,y+L,512));%計(jì)算采樣點(diǎn)坐標(biāo)值N=ones(size(x1))*Nmax;%設(shè)置N的初值z(mì)=x1+y1*i;%生成復(fù)數(shù)坐標(biāo)點(diǎn)fork=1:Nmax;z=z.^2+c;%迭代計(jì)算新的z值

N(abs(z)>M)=k;%找出模值超過(guò)M的點(diǎn)endI1=image([x-L,x+L],[y-L,y+L],N);xlabel('X');ylabel('Y');生成"Mandelbrot集合"x=0;y=0;L=2;%x=0.145;y=0.65;L=0.03;[x1,y1]=meshgrid(linspace(x-L,x+L,512),linspace(y-L,y+L,512));c=x1+i*y1;z=zeros(size(c));N=100;%最大迭代次數(shù)S=ones(size(c))*N;M=10;fork=1:N;z=z.^2+c;S(abs(z)>M)=k;c(abs(z)>M)=0;z(abs(z)>M)=0;endII=image([x-L,x+L],[y-L,y+L],S);xlabel('X');ylabel('Y');迭代函數(shù)系統(tǒng)迭代函數(shù)系統(tǒng)（iteratedfunctionsystem，IFS）是根據(jù)仿射變換而來(lái)的，仿射變換的數(shù)學(xué)表達(dá)式為或其中x和y是變換前的坐標(biāo)，x'和y'是變換后的坐標(biāo)，a、b、c、d、e和f是變換的系數(shù)。在變換表達(dá)式里面體現(xiàn)了縮放、旋轉(zhuǎn)和平移。繪制復(fù)雜的圖形，需要多個(gè)仿射變換，即仿射變換集{

n}。對(duì)應(yīng)于每個(gè)仿射變換在繪圖中都有一個(gè)固定的被選中幾率，被選中則被調(diào)用。如：繪制Sierpinski墊片需要建立3個(gè)仿射變換：ia(i)b(i)c(i)d(i)e(i)f(i)p(i)10.5000.5001/320.5000.5101/330.5000.50.5sin(pi/3)1/3Sierpinski墊片繪制M=[0.5,0,0,0.5,0,00.5,0,0,0.5,1,00.5,0,0,0.5,0.5,sin(pi/3)];p=[1/3,1/3,1/3];IFS_draw(M,p);注：在Sierpinski墊片中主要是由3個(gè)部分作遞歸而得到的，所以需要建立3個(gè)仿射變換。同時(shí)3個(gè)部分形狀是完全一樣的，只是邊長(zhǎng)是前一層邊長(zhǎng)的一半，所以T為[0.5,0;0,0.5];而3個(gè)部分的左下角位置可以這樣表示：(0,0)、(0,1)和(0.5,sqrt(3)/2);3個(gè)部分出現(xiàn)的幾率相等。IFS繪制的實(shí)現(xiàn)functionIFS_draw(M,p);%迭代函數(shù)法生成分形圖形的通用函數(shù)%M是矩陣系數(shù)；p是對(duì)應(yīng)的幾率。比如M的取值是:%i

a(i)b(i)c(i)d(i)e(i)f(i)%10.3330000.333000%20.1670-0.28900.28900.16700.33300%30.16700.2890-0.28900.16700.50000.2890%40.3330000.33300.66700%p的取值是：0.25000.25000.25000.2500N=300000;%迭代次數(shù)fork=1:length(p);%生成映射矩陣a1,a2,...eval(['a',num2str(k),'=reshape(M(',num2str(k),',:),2,3);']);endxy=zeros(2,N);%設(shè)置迭代矩陣初值pp=meshgrid(p);%生成網(wǎng)格矩陣pp=tril(pp);%取上三角矩陣pp=sum(pp,2);%對(duì)行向量求和,從而得到幾率的累加結(jié)果fork=1:N-1;a=rand-pp;%計(jì)算隨機(jī)數(shù)和幾率向量的差值

d=find(a<=0);%找出滿足幾率條件的位置

xy(:,k+1)=eval(['a',num2str(d(1)),'(:,1:2)'])*xy(:,k)+...eval(['a',num2str(d(1)),'(:,3)']);%計(jì)算迭代點(diǎn)和映射矩陣的乘積endP=complex(xy(1,:),xy(2,:));%生成復(fù)數(shù)點(diǎn)plot(P,'k.','markersize',2);%畫(huà)圖axisequal;%設(shè)置坐標(biāo)軸屬性用IFS繪制兩種樹(shù)葉編號(hào)abcdef10.65000.650.20.420.55000.550.20.13530.38-0.280.280.380.30.440.380.28-0.280.380.30.1編號(hào)abcdefp10000.15000.0820.83-0.020.040.83020.830.20.22-0.30.3020.0964-0.140.240.310.2700.50.096p均為1/4L系統(tǒng)

L系統(tǒng)是美國(guó)生物學(xué)家Lindenmayer1968年為模擬生物形態(tài)而設(shè)計(jì)的描述植物形態(tài)與生長(zhǎng)的方法。L系統(tǒng)實(shí)際上是字符串重寫(xiě)系統(tǒng)。即把字符串解釋成圖形，于是只要能生成字符串，也就等于生成了圖形。從一個(gè)初始串（叫做公理）記為W開(kāi)始，將生成規(guī)則P多次作用于其上，最后產(chǎn)生一個(gè)較長(zhǎng)的命令串，用它來(lái)繪圖。L系統(tǒng)的符號(hào)規(guī)定與解釋：F：從當(dāng)前位置向前移一步，步長(zhǎng)為h，同時(shí)畫(huà)線；G：從當(dāng)前位置向前移一步，步長(zhǎng)為h，但不畫(huà)線；+：從當(dāng)前方向逆時(shí)針轉(zhuǎn)一個(gè)給定的角度δ；-：從當(dāng)前方向順時(shí)針轉(zhuǎn)一個(gè)給定的角度δ；｜：原地轉(zhuǎn)向180°；［：Push，將龜行圖當(dāng)前狀態(tài)壓進(jìn)棧（stack）；］：Pop，將圖形狀態(tài)重置為棧頂?shù)臓顟B(tài)，并去掉該棧中的內(nèi)容；A：記錄狀態(tài)的方向；z：記錄當(dāng)前的位置。

L系統(tǒng)模擬分形植物L(fēng)系統(tǒng)用于植物結(jié)構(gòu)繪制，比如一棵樹(shù)，它是分支結(jié)構(gòu)，即一根樹(shù)干帶大量的分枝，每個(gè)分枝都有一個(gè)終點(diǎn)，是一種一個(gè)起點(diǎn)多個(gè)終點(diǎn)的圖形。這就意味著在某一運(yùn)算中，當(dāng)畫(huà)到一個(gè)分枝的盡頭時(shí)畫(huà)筆必須退回來(lái)再畫(huà)其它結(jié)構(gòu)，即產(chǎn)生一種所謂進(jìn)退操作。該操作符號(hào)是一對(duì)方括號(hào)[·]，方括號(hào)中是3個(gè)簡(jiǎn)單符號(hào)，即F,+,-。當(dāng)執(zhí)行完方括號(hào)中的指令后，畫(huà)筆回到方括號(hào)“[”前的位置并保持原方向不變。設(shè)公理W：F；生成規(guī)則P：F→FF+[+F-F-F]-[-F+F+F]；角度增量a：22.5°。在公理中，從起點(diǎn)往上兩步后，先后作出兩個(gè)分枝，而每個(gè)分枝又分別右凸左凸，最后形成一棵風(fēng)吹動(dòng)著樹(shù)的模樣。運(yùn)行：Ltree(5)L系統(tǒng)繪制分形樹(shù)的MATLAB程序functionLtree(n)S='F';a=pi/8;A=pi/2;z=0;zA=[0,pi/2];p='FF+[+F-F-F]-[-F+F+F]';fork=2:n;S=strrep(S,'F',p);endfigure;holdon;fork=1:length(S);switchS(k);case'F'plot([z,z+exp(i*A)],'linewidth',2);z=z+exp(i*A);case'+'A=A+a;case'-'A=A-a;case'['zA=[zA;[z,A