数学实验讲义4

1数学实验——

基于MATLAB电子教案(四)——综合实验制作：郑学东zxd@2拟合与插值在生产与科学实验中,有时函数y=f(x)的解析式未知,只知该函数在若干观测点的函数值或导数值;然而却希望知道观测点以外的函数值,需要估计函数在该点的值。解决思路：根据观测点的值,构造一较简单函数

(x),使其在观测点处的值等于或总体最接近f(x)的已知值。

(x)构造方法一般有两类：①测量值与真实值有误差,不要求

(x)必须过已知点,一般用拟合法。②测量值准确无误,要求

(x)必须过已知点,一般用插值法。3插值在已知的离散点的数据，构造出一个能连接所有已知点的函数,该过程叫插值。一维插值指令：yi=interp1(X,Y,xi,method)yi=interp1(X,Y,xi,method,'extrap')yi=interp1(X,Y,xi,method,extrapval)其中：'extrap'表示当xi超过X的范围时,执行外插算法;extrapval一般取0或nan,表示当xi超过范围时,外插值取extrapval;method可取'nearest','linear','cubic','spline'等,指示函数在相邻已知点间如何取值;X,Y为已知数据点的坐标;xi为需计算插值的点。4二维插值指令：ZI=interp2(X,Y,Z,XI,YI,'method');例:(网格插值peaks函数.)[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])

例:(由时间、服务年限,插值查询工资.)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二维不规则数据的插值zi=griddata(x,y,z,xi,yi,'v4'); griddata可进行三角插值，两坐标向量x,y不必单调;一般xi,yi常为用meshgrid生成的平面网格点坐标。用griddata方法处理此类不规则数据最方便.6船在该海域会搁浅吗？在某海域测得一些点(x,y)处的水深z(单位:英尺)由下表给出，水深数据是在低潮时测得的。船的吃水深度为5英尺，问在矩形(75,200)×(-50,150)里的哪些地方船要避免进入。问题分析:

假设：该海域海底是平滑的。由于测量点是散乱分布的，先在平面上作出测量点的分布图，在利用二维插值方法补充一些点的水深，然后作出海底曲面图和等高线图，并求出水深小于5的海域范围。7水道水深测量数据(单位:英尺)问题求解步骤

1、作出测量点的分布图

2、作出海底地貌图

3、危险区域海底地貌图

4、危险区域平面图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曲线拟合已知离散点上的数据集[(x1,y1),(x2,y2),…(xn,yn)],求得一解析函数y=f(x),使f(x)在xi上的尽可能接近给定的yi值,这一过程叫曲线拟合。曲线拟合过程中f(x)的类型是由实际背景确定或用户指定的,即类型已知的,但其中的系数参量需要确定。确定f(x)系数的原则是使f(x)尽可能接近已知数据点。完成曲线拟合最常用的方法是最小二乘法拟合,即找出使最小的f(x)。由微积分知识可解出f(x)的系数参量a。9涉及拟合的常用MATLAB指令[a,SS,Re]=lsqcurvefit(fun,a0,x,y);

----最小二乘曲线拟合(拟合一元函数f(x),a是其系数向量,Re是残差,SS残差平方和)[a,Re,Jacobian]=nlinfit(X,y,fun,a0);

----非线性最小二乘拟合(拟合多元函数f(x),X的每一行是对应一个y的各自变量取值.)p=polyfit(x,y,n)

----多项式拟合(拟合函数f(x)是n次多项式,x要单调,p为降幂多项式的系数向量.)曲线拟合示例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种群的分布格局实验植物种群的分布规律在很大程度上由该种群的生物学与生态学特性决定的。由于不同环境适合于具有不同适应性和需求特性的生物生长,于是在长期自然选择和自然协同过程中,便形成了各种群特有的空间分布格局。通过调查可以得到种群在空间分布各样方中的种群数,以此可进一步推测、估计整个种群的分布格局。为此我们首先要由各样方中的种群数统计出频数、频率密度,然后再估计出频率密度函数,进而计算达到一定种群数的概率。12分布格局实验要用的MATLAB函数findSize、lengthbarpolyfitsumdoubledispformatshortg13各样方中的种群数据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弹性生长模型1959年Richards提出了一个生长模型：上式中y(t)为t时刻的总生长量,参数a为总生长量的极限值,b为初始值参数,k为生长速率参数,m是曲线形状参数。它的图形是以y=a为渐近线的S形曲线。当m=0时为Mitscherlich模型当m→1时近似为Compertz模型当m=2时为Logistic模型所以Richards模型具有较广泛的适应性。15Richards曲线形状的计算机模拟改变Richards模型的曲线形状参数m和初值b(假设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('时间'),ylabel('生长量')16Richards生长曲线的拐点Richards生长曲线的拐点处就是最大生长速率时刻。clearall,clcsymsYABkxm;Y=A*(1-B*exp(-k*x))^(1/(1-m))equl=diff(Y,x,2);x=solve(equl,'x')%解符号方程equl=0,未知量为xY=subs(Y)%符号表达式Y中的x用方程根x代入17用变形虫,水稻叶数据拟合Richards模型变形虫细胞重量生长数据时间01.252.503.755.006.257.508.7510.0011.2512.50重量10.8511.3112.3013.4413.6314.1915.1815.6115.9016.9817.38时间13.7515.0016.2517.5018.7520.0021.2522.5023.7525.00重量17.7818.6619.1918.7819.2119.1419.7419.9620.0619.91水稻叶伸长生长数据时间11.82.63.44.14.85.46.16.87.48.1长度0.30.50.91.42.53.24.37.610.114.418.5时间8.89.410.110.811.712.413.114.415.115.7长度23.025.230.433.738.841.743.744.845.545.3采用曲线拟合方法,估计出Richards模型中各参数,并进行预测。18变形虫数据拟合的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('时间h'),ylabel('重量ug')title('变形虫细胞重量生长的拟合曲线')19季节性生物种群消长模型实验生物种群实际增长过程可用下面公式描述：

其中N(t)为t时刻种群密度,K为环境最大容量,r为种群内禀增长率过程,a是由K和种群初始密度决定的常数。f(t)0时即为Logistic模型。对季节性生物,其生长过程受季节性生态因子周期变化影响。因此影响因素f(t)可设为其中h为影响幅度,

为影响的初相位,为生态因子变化的角频率。对较复杂情况,h和可以是t的函数。所以,季节性生物种群总的增长规律可以表示为：20绵蚜虫种群季节性变化的模拟实验绵蚜虫种群变化的实测数据t1234567891011y3940518662085420307433841706110071031620根据季节性生物种群总的增长规律模型,采用曲线拟合方法,估计出模型中各参数;并进行预测。21绵蚜虫数据拟合的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('时间'),ylabel('数量')title('棉蚜虫种群消长实验数据拟合曲线')22线性规划线性规划模型：23求解线性规划的指令[x,fval,exitflag,output,lambda]

=linprog(c,A,b,Aeq,beq,lb,ub,x0,options)注：凸体字的部分是可选项。exitflag结束标志：">0"表示解是收敛的;"=0"表示因到达最大迭代次数而结束;"<0"表示解不收敛！output结构变量：返回迭代次数、优化算法等信息。lambda结构变量：返回Lagrange乘子信息。fval：返回最优值。不需要某项输入参数时，可用空方括号[]占据其原有位置。24应用举例解：

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饲料配方的线性规划实验背景：为了保证畜禽生产的质量和高产,应配制完善而平衡的饲料。一般每种饲养动物都有饲养标准，规定动物饲料中必须达到的营养标准;通过调查也可知道各种饲料原料的价格和所含营养成分量。问题：如何在满足饲料的营养标准前提下,配制出成本最低的饲料。方法：利用线性规划解决该类配方问题。26鸡饲料营养标准和原料限制情况鸡饲料营养标准代谢能兆卡/千克粗蛋白克/千克量比蛋白质能克/兆卡粗纤维克/千克赖氨酸克/千克蛋氨酸克/千克钙克/千克磷克/千克食盐克/千克2.7~2.8135~14551<505.62.523~404.6~6.53.7每1000千克饲料中,各原料用量限制(适口性)玉米小麦麦米糠豆饼菜籽饼鱼粉槐叶粉碳酸钙400~500100~150100~20015010030~505030~502527鸡饲料原料的营养成分及价格编号(千克)品名单价元/千克代谢能千卡/千克粗蛋白克/千克粗纤维克/千克赖氨酸克/千克蛋氨酸克/千克钙克/千克磷克/千克食盐克/千克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鱼粉0.822.00450029.111.86327X8槐叶粉0.201.6117010810.62.24.04.0X9DL蛋氨酸12.00980X10骨粉0.192300140X11碳酸钙0.52400X12食盐0.18100028配制1千千克饲料的线性规划模型(目标函数)(总重量约束)(代谢能约束)(粗蛋白约束)(粗纤维约束)(赖氨酸约束)(蛋氨酸约束)(钙约束)(磷约束)(原料用量约束)DEA评价模型简介DEA是对一些同类型的部门或单位(称为决策单元DMU)进行相对有效性评价的一种方法.一个决策单元的有效性是用该单元的多指标输出的加权和与多指标输入的加权和之比来定义的.对n个决策单元

，向量xi=[x1i,x2i,…,xmi]T和向量yi=[y1i,y2i,…,ysi]T,分别表示第i个DMU的m项输入和s项输出,X=[x1,x2,…,xn]和Y=[y1,y2,…,yn]分别称为多指标输入矩阵和输出矩阵.2930DEA模型C2R12…

…n

…x1x2

…

…xn y1y2

…

…yn

…

评价DMU－j0的DEA模型（C2R）为分式规划（h0为DMU－j0的效率指数）：分式规划(其中x0=xj0,y0=yj0,1≤j0≤n．31

用1962年Charnes和Cooper对于分式规划的Charnes－Cooper变换（称为C2－变换）:32Matlab线性规划标准形式DEA评价示例Dea程序%五个决策单元，3个投入性指标X=[3060554070;2540703090;13015012070180]%五个决策单元，2个产出性指标Y=[3543765263;6080534271]k=5%计算第k个单元的效率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%上述程序可设计为一个求解Dea的指令函数：%[w,v,h]=dea(X,Y,k)34DEA交叉评价由于前面提到的DEA模型的决策单元效率是利用对其最有利的权重计算出来的(自我评价)。在实际问题中，往往有较多的决策单元都能取到最大的效率值1，从而不能区分这些决策单元的优劣。为了解决这个问题，我们引入交叉评价机制.用每一个DMUi的最佳权重去计算其它DMUk的效率值35由于线性规划的最优解ui*和vi*不唯一，得出的交叉评价值Eik具有不确定性。对抗型交叉评价(1)先用自我评价模型得出每一个DMUi的自我评价值Eii;(2)在保证DMUi得到最大值Eii的前提下，使其它的DMUk得到尽可能小的交叉评价值Eik。对抗型交叉评价的实质是:每一个DMUi,在尽可能抬高自己的前提下，尽可能地贬低其它DMUk。36对抗型交叉评价模型算法37将E的第i列的平均值ei作为衡量DMUi的优劣的一项指标;ei可视为诸决策单元对DMUi的总的评价,越大越好。对抗型交叉评价模型程序functionE=deaEij(X,Y)%输入指标矩阵X，输出指标矩阵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实验指标银行员工人数((投入)机构总数(投入)营业费用率(投入)总资产(亿)(投入)人均利润率(万元)(产出)资产收益率(产出)工商银行i=1381713164760.55466586842.8830.1572120.010126农业银行i=2447519244520.78563160501.277.275890.00208中国银行i=3237379101450.46610159955.5337.8950960.01095建设银行i=4298868134480.54189665981.7733.7326180.011479浦发银行i=5141284080.5840729149.803574.148790.006855招商银行i=6289715790.4086613105.5272.596730.013582华夏银行i=793902870.7296145923.382740.690030.010898兴业银行i=8118513900.5081238513.352792.061180.016909浙商银行i=91372210.217699575.2279940.1203350.007526上海银行i=1047681880.5139033089.860179.6233010.010262上海农村商业银行i=1141963270.2685721574.727517.189590.002014xi1,xi2,xi3,xi4yi1,yi2综合作业1作业1、第34页的DEA程序改造为用户函数： [w,v,h]=dea(X,Y,k)输入参数为：多指标输入矩阵Xmxn和输出矩阵Ysxn以及单元标号k。输出参数为：输入指标的权重w和输出指标的权重v以及第k单元的效率h.作业2、用对抗性交叉评价程序(第39页)来测算第十五页上各单位的效率高低。41非线性规划数学模型：42求解非线性规划的指令[x,fval,exitflag,output,lambda,grad,hessian]=fmincon('fun',x0,A,b,Aeq,beq,lb,ub,'nonlcon',options,P1,P2,...)其中nonlcon为非线性约束及导数：[C,Ceq,gC,gCeq]=nonlcon(x)；

而fun则为目标函数及导数和Hess矩阵：

[f,g,H]=myfun(x)；fval返回最优目标值。注：凸体字的部分是可选项；输入空矩阵[]则表示该项不需要输入或用缺省值。43应用举例

function[f,g,H]=myoptfun(x)%目标函数(包括函数值,梯度,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优化选项optionsoptions.DerivativeCheck[on|{off}]检查梯度输入options.Diagnostics[on|{off}]

options.Display[off|iter|{final}]结果等信息何时显示★options.DiffMaxChange[正数|{1e-1}]导数的最大改变options.DiffMinChange[正数|{1e-8}]导数的最小改变options.GradConstr[on|{off}]约束求梯度options.GradObj[on|{off}]目标求梯度options.LargeScale[{on}|off]大规模算法★options.MaxFunEvals[正整数]

options.MaxIter[正整数]最大迭代次数options.TolFun[正数]函数值终止误差标准★options.Tolx[正数]x的终止误差标准★46用optimset函数设置options参数options=optimset('param1',value1,'param2',value2,…)

%建立options结构中的各参数。options=optimset(oldopts,'param1',value1',…)

%修改oldopts中的参数param1等。例：options=optimset('Display','iter','TolFun',1e-8)无条件极值fminunc

/fminsearch问题：求minf(x)?求解指令：[x,fval]=fminunc(fun,x0,options)[x,fval]=fminsearch(fun,x0,options)47二次规划模型指令

[x,fval,exitflag,output,lambda]

=quadprog(H,f,A,b,Aeq,beq,lb,ub,x0,options)应用举例解: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组合投资计算设一组投资由m个投资项目组成，它们的单项期望回报率为(r1,r2,…,rm)，该m个项目的投资比例为(x1,x2,…,xm)，则该组投资的总回报率的期望值R：投资组合的总回报率的期望值描述了多项投资的总体平均回报水平。同样地。仅汉只用总回报率期望值来描述投资组合的效果是不够的，还需描述总回报率的离散趋势，也就是整个投资组合风险的大小。一组投资的总回报率的风险（或离散趋势）的常用测度是总回报率的方差和标准方差。总回报率R的方差总回报率R的方差的估算公式其中：为该m个项目的投资比例，是该m个项目的各单项回报率的标准差，为第i个投资项目以与第j个投资项目的相关系数，V为各单项回报率间的协方差矩阵。投资者的目标是获得大的投资回报和承担小的投资风险。投资组合优化模型就是确定一组投资项目的最优投资比例,在该投资组合的总回报率的方差不超过某个可接受的值的约束下(即在可接受的风险水平下)，使得总回报率的期望值最大(即投资回报水平最高)；或者在投资组合的回报率的期望值不低于某个所要求的值的约束下(即在所要求的投资回报水平)，使得总回报率的方差最小(即投资风险最小)。投资组合优化模型多目标优化模型：总回报率最大模型：总风险最小模型：有效前沿Markowitz资产组合理论就是寻找一个能使在同样收益水平下,总风险水平最低的有效组合。不同的收益与相应的最小风险构成有效前沿(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))有效前沿计算示例例考虑一个三资产组合,分别为资产1、资产2与资产3，其预期收益率分别为0.2、0.1、0.15。协方差如下表资产1资产2资产3资产10.01-0.00610．0042资产2-0.00610.04-0.0252资产30.0042-0.02520.0225求该资产组合有效前沿。示例程序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)练习投资三种资产的收益率相关系数矩阵、预期回报、标准差如表所示：资产A资产B资产C资产A10.80.4资产B0.810.3资产C0.40.31预期回报0.10.150.20各资产标准差0.20.250.18试给出有效前沿。相关阵到协方差阵的转换指令corr2cov函数用法：

Covariances=corr2cov(STDs,Correlations);57非线性方程组求数值解MATLAB求‘Fun(x)=0’解的指令一般格式：

[x,fval,exitflag,output]=fsolve(Fun,x0,options,p1,p2,...)例：求解方程组sinx+y+z2ex-4=0,x+yx=0,xyz=0的具体步骤如下 创建函数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指令补充说明返回参数exitflag为求解结束状态标志：

>0----解收敛;=0----达到迭代次数;<0----解发散。输入参数Fun可为向量函数(表示方程组)，其定义形式为：

functionF=Fun(x)

F=...%Computefunctionvaluesatx

或

function[F,J]=Fun(x)

F=...%objectivefunctionvaluesatx

ifnargout>1%twooutputarguments

J=...%Jacobianofthefunctionevaluatedatx

endKMV信用监控模型KMV模型是美国旧金山市KMV公司用来估计借款企业违约概率的方法。该模型认为，在债务到期日，如果公司资产的市场价值高于公司债务值（违约点），则公司股权价值为公司资产市场价值与债务值之间的差额；如果此时公司资产价值低于公司债务值，则公司变卖所有资产用以偿还债务，股权价值变为零。59KMV模型的运用步骤首先,利用Black-Scholes期权定价公式,根据企业资产的市场价值、资产价值的波动性、到期时间、无风险借贷利率及负债的账面价值估计出企业股权的市场价值及其波动性。其次根据公司的负债计算出公司的违约实施点(defaultexercisepoint,为企业1年以下短期债务的价值加上未清偿长期债务账面价值的一半),计算借款人的违约距离。最后,根据企业的违约距离与预期违约率(EDF)之间的对应关系，求出企业的预期违约率。60违约距离DD计算61式中：DP(DefaultPoint)为违约点值DP=短期负债(STD)+0.5长期负债(LTD)式中:E为企业股权市场价值,V为企业资产市场价值,D为企业债务面值,r为无风险收益率,τ为债务偿还期限,N(.)为标准累积正态分布函数,

V为企业资产价值波动率,

E为企业股权市场价值波动率。其中V和

V是未知量,解非线性方程组可得。据违约距离DD的定义可知,理论上发生违约的概率为1-N(DD).而基于违约数据库,依据违约距离可以映射出公司实际的期望违约频率EDF。由于我国当前还没有公开的违约数据库可以使用，以违约距离DD作为上市公司信用评价的依据。上市公司股权市场价值

中国证券市场发展历史较为特殊，上市公司股票被人为分割为上市流通股票和暂不上市流通股票两种。在计算上市公司股权市场价值时需要考虑以什么样的价格来计算非流通股市场价值/由于非流通股没有市场交易价格，因此如何给非流通股定价是一件困难的事情。参考上市公司股票全流通研究中非流通股定价，以每股净资产计算非流通股的价格。流通股市场价值=12月份周平均收盘价格×流通股股数非流通股市场价值=每股净资产×非流通股股数上市公司股权市场价值=流通股市场价值+非流通股市场价值公司债务面值D为公司财务年报中总负债面值。考虑到数据和工作量的限制，我们设定违约距离的计算时间为一年，即τ=1。无风险利率使用中国人民银行公布的一年期定期整存整取的存款利率。KMV测算示例某公司流动负债为1亿元，长期负债为5000万元，根据上市公司的股价行情表，可以统计出公司股权价值E和公司股权价值波动率

E。试计算公司的违约概率。月份总市值(元)收益率(%)月份总市值(元)收益率(%)1129523558-13.6581409724643.31214988546213.5891484050955.013142316387-5.3210144898861-2.4241494409124.77111449046090.005147924524-1.0312130292794-11.216130439432-13.40均值14127642771363130244.31标准差8.35(1)基本参数计算公司股价波动率为8.35%，公司股权价值波动率公司股权价值E=141276427(2)求解非线性方程组，解出Va和

a

。KMV测算示例(续-程序)KMV计算示例的Matlab程序表示联列方程组的向量函数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测算示例(续-程序执行)求解的指令：%KMV计算示例1%r:无风险利率r=0.0225;%T:债务期限T=1;%输入年数%DP:违约点t%SD:短期债务,LD:长期债务SD=1e8;%输入LD=50000000;%输入%计算违约点DP=SD+0.5*LD;%D:债务的市场价值D=DP;%债务的市场价值,可改为公司年报公布的总负债面值%sigma:波动率sigma=0.0835;%(输入）%Esigma:总波动率Esigma=sigma*sqrt(12);%E:公司股权价值E=141276427;%开始计算VaandAsigma;若exitflag<0,重新设定x0再算一遍！x0=[100000000;1];%x(1)的初值应接近E，x(2)的初值应接近Esigma[x,f,exitflag]=fsolve(@KMVe1fun,x0,[],E,D,Esigma,r,T)Va=x(1),Asigma=x(2)%计算违约距离DD=(Va-DP)/(Va*Asigma)%计算违约率EDF=normcdf(-DD)混沌分形简介与实验函数迭代与不动点给定一函数f(x)以及初始点x0，定义数列

xk+1=f(xk)

,k=0,1,2,…

{xk}称为函数f(x)的迭代序列。函数f的不动点：满足f(x)=x的点x称为f的不动点。吸引点与排斥点：附近的点经函数f迭代后都趋向f的某一不动点x，则称x为吸引点；附近的点经函数f迭代后都远离不动点x，则称它是排斥点。周期点与预周期点K-周期点：若f(u1)=u2,f(u2)=u3,…,f(uk)=u1;则点集u1,u2,u3,…,uk形成一个k循环;u1称为k-周期点;k称为周期。类似地，周期点也可以分吸引点与排斥点。预周期点:如果点x最终归宿于某个循环中，则称它为预周期点。如1是x2-1的预周期点。注：迭代序列{xk}的收敛与发散性质不仅与函数f有关，而且与初值的选择有关。例如，对于迭代xk+1=2xk+1,当初值x0=-1时，迭代序列收敛，否则发散。吸引与排斥点判断//迭代函数y=ax(1-x)定理:设x是f(x)的不动点,如果在x附近有|f'(x)|<1,则x是f(x)的吸引不动点;否则,x是f(x)的排斥不动点。对迭代函数f(x)=ax(1-x)(0<a4)而言：不动点：x=ax(1-x)x=0或x=(a-1)/a吸引/排斥点：f’(0)=a,f’((a-1)/a)=2-a，故当0<a<1时，0为吸引点，(a-1)/a为排斥点。当1<a<3,0为排斥点，(a-1)/a为吸引点。2-周期点：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模拟迭代引起的混沌现象对迭代函数xk+1=axk(1-xk)进行模拟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('混沌现象观察');xlabel('迭代次数'),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=['混沌现象观察:','a=',num2str(a),',x0=',num2str(x0)];t=1:200;scatter(t,x,3,'r','filled');title(str);xlabel('迭代次数'),ylabel('序列{x(n)}')---------------x=itfun(0.5,1.4);思考与实验(y=ax(1-x))1〉对几组不同的参数值a(如a=0.5,a=1.4以及不同的初值x0，观察迭代是否收敛。2〉取参数a=0.8，用不同的初值x0做迭代。你能找到一个吸引的不动点吗？一个排斥的不动点吗？哪些初值收敛到吸引的不动点？哪些初值使序列发散？取不动的参数a=1,1.6,2,2.5回答同样的问题。3〉找出一个参数a,使它对应的迭代具有2周期点。这种性质依赖于初值吗?4〉对任意的整数k，你能找到一个a值使得它对应的迭代具有k周期点吗？对哪些k值能给出k周期点？在每种情况下，结果是否依赖于初值？(对3.4a

3.6和3.6a

4的值进行验证)72混沌实验初探迭代公式在方程求解等数值计算中,有着重要的地位。下面对一个迭代公式(也可看作差分方程),讨论其解的稳定性。对迭代公式：yk+1=ayk(1-yk)1)当1<a<3时,方程有稳定的平衡点,y收敛到1-1/a。2)当3<a<3.499时,有2个平衡点,出现2倍周期现象。从k的连续两期看,y的变化是稳定的。3)当3.499<a<3.544时,有4个平衡点,出现4倍周期。4)当a>3.57时,便不出现任何2n倍周期现象,{yn}的变化趋势是一片混乱,这就是常说的混沌现象。混沌的特性:①对初值敏感.初值虽近,迭代后散开.②混沌不是随机.不同初值,x值落在某区间频率不定.73混沌现象模拟的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('混沌现象观察');xlabel('迭代次数'),ylabel('序列{x(n)}')74分枝与混沌(采用系列不同参数a观测)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)产生的分枝与混沌图')xlabel('参数a');ylabel('迭代式的根x');end;holdoff;75迭代结果的分布情况x0=input('请输入初值(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('迭代次数n');ylabel('x值');disp(['x0=',num2str(x0),'p=n1/n=',num2str(n1/n)]);disp(['x值落在0~1/2中的次数n1=',num2str(n1)]);disp(['x值落在0~1/2外的次数n2=',num2str(n2)]);76分形几何与分形艺术什么是分形几何？通俗的说法就是研究无限复杂但具有一定意义下自相似图形和结构的几何学。"分形"译于英文Fractal,词本身具有"破碎"、"不规则"等含义。分形几何创始人Mandelbrot最精彩的部分研究是1980年发现了Mandelbrot集合，他发现整个宇宙以一种出人意料的方式构成自相似的结构。Mandelbrot集合图形的边界处，具有无限复杂和精细的结构。当你放大某个区域，它的结构就在变化，展现出新的结构元素。用数学方法对放大区域进行着色处理，这些区域就变成一幅幅精美的艺术图案，这些艺术图案人们称之为"分形艺术"。

77复平面中的神奇迭代-Julia集合在复平面上,x轴代表实数,y轴代表虚数。现在复平面上任意取一个点z,将其代入下面方程中进行反复迭代运算：

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

①z值没有界限增加（趋向无穷）

②z值衰减（趋向于零）

③z值是变化的,即非1或非2趋向无穷和趋向于零的点叫定常吸引子，很多点在定常吸引子处结束，被定常吸引子所吸引。非趋向无穷和趋向于零的点是"Julia集合"部分，也叫混沌吸引子。Mandelbrot集合是所有的Julia集合的合并。78得到"Julia集合"的近似算法一般按下述算法近似计算"Julia集合"：n=0;while((n++<Nmax)&&((Real(Z)^2+Imag(Z)^2)<Rmax)){Z=Z*Z+C;}其中：Nmax为最大迭代次数,Rmax为逃离界限。由(Real(Z)^2+Imag(Z)^2)>=Rmax退出循环时,相当于"①z值没有界限增加(趋向无穷)",为定常吸引子,把这些区域着成白色。由(n>=Nmax)退出循环时,相当于"②z值衰减(趋向于零)"或"③z值是变化的",把这些区域着成黑色。黑色区域图形的边界处即为"Julia集合"。"Julia集合"有着极其复杂的形态和精细的结构。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;%设置对参数c的缺省值x=0;y=0;L=2;Nmax=100;%设置最大迭代次数M=100;%设置一个大数，表示无穷大[x1,y1]=meshgrid(linspace(x-L,x+L,512),linspace(y-L,y+L,512));%计算采样点坐标值N=ones(size(x1))*Nmax;%设置N的初值z=x1+y1*i;%生成复数坐标点fork=1:Nmax;z=z.^2+c;%迭代计算新的z值

N(abs(z)>M)=k;%找出模值超过M的点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;%最大迭代次数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');迭代函数系统迭代函数系统（iteratedfunctionsystem，IFS）是根据仿射变换而来的，仿射变换的数学表达式为或其中x和y是变换前的坐标，x'和y'是变换后的坐标，a、b、c、d、e和f是变换的系数。在变换表达式里面体现了缩放、旋转和平移。绘制复杂的图形，需要多个仿射变换，即仿射变换集{

n}。对应于每个仿射变换在绘图中都有一个固定的被选中几率，被选中则被调用。如：绘制Sierpinski垫片需要建立3个仿射变换：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个部分作递归而得到的，所以需要建立3个仿射变换。同时3个部分形状是完全一样的，只是边长是前一层边长的一半，所以T为[0.5,0;0,0.5];而3个部分的左下角位置可以这样表示：(0,0)、(0,1)和(0.5,sqrt(3)/2);3个部分出现的几率相等。IFS绘制的实现functionIFS_draw(M,p);%迭代函数法生成分形图形的通用函数%M是矩阵系数；p是对应的几率。比如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;%迭代次数fork=1:length(p);%生成映射矩阵a1,a2,...eval(['a',num2str(k),'=reshape(M(',num2str(k),',:),2,3);']);endxy=zeros(2,N);%设置迭代矩阵初值pp=meshgrid(p);%生成网格矩阵pp=tril(pp);%取上三角矩阵pp=sum(pp,2);%对行向量求和,从而得到几率的累加结果fork=1:N-1;a=rand-pp;%计算随机数和几率向量的差值

d=find(a<=0);%找出满足几率条件的位置

xy(:,k+1)=eval(['a',num2str(d(1)),'(:,1:2)'])*xy(:,k)+...eval(['a',num2str(d(1)),'(:,3)']);%计算迭代点和映射矩阵的乘积endP=complex(xy(1,:),xy(2,:));%生成复数点plot(P,'k.','markersize',2);%画图axisequal;%设置坐标轴属性用IFS绘制两种树叶编号abcdef10.65000.650.20.420.55000.550.20.13530.38-0.280.280.380.30.440.380.28-0.280.380.30.1编号abcdefp10000.15000.0820.83-0.020.040.83020.830.20.22-0.30.3020.0964-0.140.240.310.2700.50.096p均为1/4L系统

L系统是美国生物学家Lindenmayer1968年为模拟生物形态而设计的描述植物形态与生长的方法。L系统实际上是字符串重写系统。即把字符串解释成图形，于是只要能生成字符串，也就等于生成了图形。从一个初始串（叫做公理）记为W开始，将生成规则P多次作用于其上，最后产生一个较长的命令串，用它来绘图。L系统的符号规定与解释：F：从当前位置向前移一步，步长为h，同时画线；G：从当前位置向前移一步，步长为h，但不画线；+：从当前方向逆时针转一个给定的角度δ；-：从当前方向顺时针转一个给定的角度δ；｜：原地转向180°；［：Push，将龟行图当前状态压进栈（stack）；］：Pop，将图形状态重置为栈顶的状态，并去掉该栈中的内容；A：记录状态的方向；z：记录当前的位置。

L系统模拟分形植物L系统用于植物结构绘制，比如一棵树，它是分支结构，即一根树干带大量的分枝，每个分枝都有一个终点，是一种一个起点多个终点的图形。这就意味着在某一运算中，当画到一个分枝的尽头时画笔必须退回来再画其它结构，即产生一种所谓进退操作。该操作符号是一对方括号[·]，方括号中是3个简单符号，即F,+,-。当执行完方括号中的指令后，画笔回到方括号“[”前的位置并保持原方向不变。设公理W：F；生成规则P：F→FF+[+F-F-F]-[-F+F+F]；角度增量a：22.5°。在公理中，从起点往上两步后，先后作出两个分枝，而每个分枝又分别右凸左凸，最后形成一棵风吹动着树的模样。运行：Ltree(5)L系统绘制分形树的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