MATLAB数学建模14个范例

1、精选优质文档-倾情为你奉上1.专心-专注-专业1.整数规划的蒙特卡洛解法2015-06-10已知非线性整数规划为： 如果用显枚举试探，共需要计算1005=1010个点，其计算量非常大。然而应用蒙特卡洛去随机模拟计算106个点，便可以找到满意解，那么这种方法的可信度究竟怎么样呢？ 下面就分析随机采样106个点计算时，应用概率理论估计下可信度。 不是一般性，假设一个整数规划的最优点不是孤立的奇点。 假设目标函数落在高值区的概率分别为0.01，0.00001，则当计算106个点后，有任一个点落在高值区的概率分别为：1-0.99=0.99.99(100多位)1-0.99999=0. 解 （1）首先编写

2、M文件 mengte.m定义目标函数f和约束向量g,程序如下：function f,g=mengte(x);f=x(1)2+x(2)2+3*x(3)2+4*x(4)2+2*x(5)-8*x(1)-2*x(2)-3*x(3)-.x(4)-2*x(5);g=sum(x)-400x(1)+2*x(2)+2*x(3)+x(4)+6*x(5)-8002*x(1)+x(2)+6*x(3)-200x(3)*x(3)+x(4)+5*x(5)-200; （2）编写M文件mainint.m如下求问题的解：rand('state',sum(clock);p0=0;tic for i=1:105 x=

3、99*rand(5,1); x1=floor(x);%向下取整 x2=ceil(x);%向上取整 f,g=mengte(x1); if sum(g<=0)=4 if p0<=f x0=x1; p0=f; end end f,g=mengte(x2); if sum(g<=0)=4 if p0<=f x0=x2; p0=f; end endendx0,p0Matlab求解整数规划祥见第二章（优秀教材）2. 罚函数法 2015-06-11 利用罚函数法，可将非线性规划问题的求解，转化为求解一系列无约束极值问题，因而也称这种方法为系列无约束最小化技术，简记为SUMT。 罚函数

4、法求解非线性规划问题的思想是，利用问题中的约束函数作出适当的罚函数，由此构造出带参数的增广目标函数，把问题转化为无约束非线性规划问题。主要有两种形式，一种叫外函数法，另一种叫内函数法，下面介绍外函数法。 考虑问题：取一个充分大的数M>0,构造函数： 增广目标函数P(x,M)为目标函数的无约束极值问题 minP(x,M)的最优解x也是原问题的最优解。Eg:求解下列非线性规划问题解 编写M文件test1.mfunction g=test1(x);M=50000;f=x(1)2+x(2)2+8;g=f-M*min(x(1),0)-M*min(x(2),0)-M*min(x(1)2-x(2),0

5、)+M*abs(-x(1)-x(2)2+2);然后命令窗口输入x,y=fminunc('test1',rand(2,1)扩展： 将有约束的极值问题转化为无约束的极值问题，这样便可以很好的扩展matlab神经网络43个案例分析里面求无约束的函数极值问题。Matlab求解非线性规划祥见第三章（优秀教材）3. 层次分析 2015-06-12 层次分析祥见第八章 准则层b=1 1 1 4 1 1/2 ;方案层cb1=1 1/4 1/2; cb2=1 1/4 1/5;cb3=1 3 1/3;cb4=1 1/3 5;cb5=1 1 7; cb6=1 7 9; 准则B1B2B3B4B5B6总

6、权值准则层权值0.16000.16000.16000.04000.16000.3200工作10.14290.10000.23080.23810.46670.79750.4152工作20.57140.40000.07690.71430.46670.11390.3074工作30.28570.50000.69230.04760.06670.08860.2774由总权值可以知道工作1最合适。 权重的计算weight1.m 传入各因素重要性 例如b=1 3 9就可以计算相应的权重了weigth1(b)。Weight.m是近似算法和weight1.m 用法相同。主要还是用 weight1.m weight

7、1.m function w,cr=weight1(b)%b=1 1/3 1/9;%模糊比较 n1=length(b);ri=0 0 0.58 0.90 1.12 1.24 1.32 1.41 .145;%一致性指标a=zeros(n1,n1);for i=1:n1 for j=1:n1 a(i,j)=b(j)/b(i); endendx,y=eig(a);lamda=max(diag(y);num=find(diag(y)=lamda);w=x(:,num)/sum(x(:,num);cr=(lamda-n1)/(n1-1)/ri(n1);%一致性检验 若cr<0.1,认为矩阵的一致性

8、可以接受主程序 score3.mb=1 1 1 4 1 1/2 ;cb1=1 1/4 1/2; cb2=1 1/4 1/5;cb3=1 3 1/3;cb4=1 1/3 5;cb5=1 1 7; cb6=1 7 9;w0=weight1(b);w1=weight1(cb1);w2=weight1(cb2);w3=weight1(cb3);w4=weight1(cb4);w5=weight1(cb5);w6=weight1(cb6);ww=w1 w2 w3 w4 w5 w6;score3=ww*w0 %计算总的权值4. 粒子群优化算法的寻优算法-非线性函数极值寻优 2015-06-13 PSO算法

9、首先在可解空间中初始化一群粒子，每个粒子都代表极值最优问题的一个潜在最优解。用位置、速度和适应度三项指标表示该粒子的特征，适应度值由适应度函数计算得到，其值的好坏表示例子的优劣。粒子在解空间中运动，通过跟踪个体极值Pbest和群体极值Gbest跟新个体位置；个体极值Pbest是指个体所经历位置中计算得到的适应度最优解的位置，群体的极值Gbest是指种群中的所有粒子搜索到的适应度最优位置。粒子每更新一次位置，就计算一次适应度值，并通过比较粒子的适应度值和个体极值、群体极值的适应度值来跟新个体极值Pbest和群体极值Gbest位置。 假设在一个D维的搜索空间中有n个粒子组成的种群，其中第i个粒子表

10、示为一个D维的向量,代表第i个粒子在D为搜索空间中的位置，亦代表问题的一个潜在解。根据目标函数即可以计算出每个粒子位置对应的适应度值。第i个粒子的速度，其中个体极值为,种群的全局极值为。 在每一次迭代过程中，粒子通过个体极值和全局极值更新自身的速度和位置，更新公式如下：速度：位置：式中，为惯性权重；为非负的常数，称为加速因子；为分布于0,1之间的随机数。为了防止粒子的盲目搜索，一般建议将奇位置和速度限制在一定的区间。例子：求非线性函数的最小值 Min f(x)如果是带约束的可以通过罚函数法转化为无约束的。程序如下所求函数的最小值，函数如下：function y = fun(x)%函数用于计算粒

11、子适应度值%x input 输入粒子 %y output 粒子适应度值 y=-20*exp(-0.2*sqrt(x(1)2+x(2)2)/2)-exp(cos(2*pi*x(1)+cos(2*pi*x(2)/2)+20+exp(1)*(x(1)<10)*(x(2)<10)+x(3);%y=x(1)2+x(2)2+7+sin(x(1)+x(2);%y=x(1)2+x(2)2+7+sin(x(1)+x(2)*(x(1)+x(2)>10);%y=x(1)2-10*cos(2*pi*x(1)+10+x(2)2-10*cos(2*pi*x(2)+10;%y=-x(1)2-x(2)+x(

12、3)*x(2)*x(1);主程序： PSO.mclcclear % 参数初始化%粒子群算法中的两个参数c1 = 1.49445;c2 = 1.49445;%c1=20;%c2=20;maxgen=500; % 进化次数 sizepop=50; %种群规模 Vmax=1;%速度的范围Vmin=-1;popmax=5;%变量的取值范围popmin=-5; % 产生初始粒子和速度for i=1:sizepop %随机产生一个种群 pop(i,:)=5*rands(1,3); %初始种群 变量的维数3 V(i,:)=rands(1,3); %初始化速度 速度的维数和变量的维数相同 %计算适应度 fit

13、ness(i)=fun(pop(i,:); %染色体的适应度end % 个体极值和群体极值bestfitness bestindex=min(fitness);zbest=pop(bestindex,:); %全局最佳gbest=pop; %个体最佳fitnessgbest=fitness; %个体最佳适应度值fitnesszbest=bestfitness; %全局最佳适应度值 % 迭代寻优for i=1:maxgen for j=1:sizepop %速度更新 V(j,:) = V(j,:) + c1*rand*(gbest(j,:) - pop(j,:) + c2*rand*(zbest

14、 - pop(j,:); V(j,find(V(j,:)>Vmax)=Vmax; V(j,find(V(j,:)<Vmin)=Vmin; %种群更新 pop(j,:)=pop(j,:)+0.5*V(j,:); pop(j,find(pop(j,:)>popmax)=popmax; pop(j,find(pop(j,:)<popmin)=popmin; %适应度值 fitness(j)=fun(pop(j,:); end for j=1:sizepop %个体最优更新 if fitness(j) < fitnessgbest(j) gbest(j,:) = pop(

15、j,:); fitnessgbest(j) = fitness(j); end %群体最优更新 if fitness(j) < fitnesszbest zbest = pop(j,:);%最优值的位置 fitnesszbest = fitness(j);%最小的函数值 end end yy(i)=fitnesszbest; end% 结果分析plot(yy)title('最优个体适应度','fontsize',12);xlabel('进化代数','fontsize',12);ylabel('适应度',

16、9;fontsize',12);axis(0 500 0 0.4)增加自适应变异变异操作扩展了在迭代中不断缩小的种群空间，使粒子能够跳出先前搜索到的最优值位置，在更大的空间中开展搜索，同时保持了种群多样性，提高了算法寻找到更优值的可能性。代码：PSOMutation.m % 该代码为基于变异粒子群算法的函数极值寻优算法% 清空环境clcclear % 参数初始化%粒子群算法中的两个参数%c1 = 1.49445;%c2 = 1.49445;c1=200;c2=100;maxgen=500; % 进化次数 sizepop=50; %种群规模 Vmax=1;Vmin=-1;popmax=5

17、;popmin=-5; % 产生初始粒子和速度for i=1:sizepop %随机产生一个种群 pop(i,:)=5*rands(1,3); %初始种群 V(i,:)=rands(1,3); %初始化速度 %计算适应度 fitness(i)=fun(pop(i,:); %染色体的适应度 %目标函数end % 个体极值和群体极值bestfitness bestindex=min(fitness);zbest=pop(bestindex,:); %全局最佳gbest=pop; %个体最佳fitnessgbest=fitness; %个体最佳适应度值fitnesszbest=bestfitness

18、; %全局最佳适应度值 % 迭代寻优for i=1:maxgen %maxgen最大的迭代次数 for j=1:sizepop omg=1;%惯性权重 %速度更新 V(j,:) =omg* V(j,:) + c1*rand*(gbest(j,:) - pop(j,:) + c2*rand*(zbest - pop(j,:); V(j,find(V(j,:)>Vmax)=Vmax; V(j,find(V(j,:)<Vmin)=Vmin; %种群更新 pop(j,:)=pop(j,:)+0.5*V(j,:); pop(j,find(pop(j,:)>popmax)=popmax;

19、 pop(j,find(pop(j,:)<popmin)=popmin; if rand>0.98 pop(j,:)=rands(1,2); %变异操作 end %适应度值 fitness(j)=fun(pop(j,:); end for j=1:sizepop %个体最优更新 if fitness(j) < fitnessgbest(j) gbest(j,:) = pop(j,:); fitnessgbest(j) = fitness(j); end %群体最优更新 if fitness(j) < fitnesszbest zbest = pop(j,:); fitn

20、esszbest = fitness(j); end end yy(i)=fitnesszbest; end% 结果分析plot(yy)title('最优个体适应度','fontsize',12);xlabel('进化代数','fontsize',12);ylabel('适应度','fontsize',12);axis(0 500 0 0.4)惯性权重的选择：PSOMutationOmg.m % 该代码为基于变异粒子群算法的函数极值寻优算法% 清空环境clcclear % 参数初始化%粒子群算法中的

21、两个参数%c1 = 1.49445;%c2 = 1.49445;c1=200;c2=100;maxgen=500; % 进化次数 sizepop=50; %种群规模omgstar=0.9;%初始惯性权重omgend=0.4;%迭代最大次数时的惯性权重 Vmax=1;Vmin=-1;popmax=5;popmin=-5; % 产生初始粒子和速度for i=1:sizepop %随机产生一个种群 pop(i,:)=5*rands(1,3); %初始种群 V(i,:)=rands(1,3); %初始化速度 %计算适应度 fitness(i)=fun(pop(i,:); %染色体的适应度 %目标函数e

22、nd % 个体极值和群体极值bestfitness bestindex=min(fitness);zbest=pop(bestindex,:); %全局最佳gbest=pop; %个体最佳fitnessgbest=fitness; %个体最佳适应度值fitnesszbest=bestfitness; %全局最佳适应度值 % 迭代寻优for i=1:maxgen %maxgen最大的迭代次数 for j=1:sizepop omg=omgstar-(omgstar-omgend)*i/maxgen;%惯性权重 %速度更新 V(j,:) =omg* V(j,:) + c1*rand*(gbest(

23、j,:) - pop(j,:) + c2*rand*(zbest - pop(j,:); V(j,find(V(j,:)>Vmax)=Vmax; V(j,find(V(j,:)<Vmin)=Vmin; %种群更新 pop(j,:)=pop(j,:)+0.5*V(j,:); pop(j,find(pop(j,:)>popmax)=popmax; pop(j,find(pop(j,:)<popmin)=popmin; if rand>0.98 pop(j,:)=rands(1,2); %变异操作 end %适应度值 fitness(j)=fun(pop(j,:); e

24、nd for j=1:sizepop %个体最优更新 if fitness(j) < fitnessgbest(j) gbest(j,:) = pop(j,:); fitnessgbest(j) = fitness(j); end %群体最优更新 if fitness(j) < fitnesszbest zbest = pop(j,:); fitnesszbest = fitness(j); end end yy(i)=fitnesszbest; end% 结果分析plot(yy)title('最优个体适应度','fontsize',12);xlab

25、el('进化代数','fontsize',12);ylabel('适应度','fontsize',12);axis(0 500 0 0.4)主要的文字说明见MATLAB神经网络43个案例分析35章5有约束函数极值APSO寻优 2015-06-14目标函数：约束条件：Bestsolution = 0.2124 3.3917 8.8909 0.2125fmin = 1.7501使用者只需要改前三部分% 改进的快速粒子群优化算法 (APSO):function apso Lb=0.1 0.1 0.1 0.1; %下边界Ub=2.0 10

26、.0 10.0 2.0; %上边界% 默认参数para=25 150 0.95; %粒子数，迭代次数，gama参数 % APSO 优化求解函数gbest,fmin=pso_mincon(cost,constraint,Lb,Ub,para); % 输出结果Bestsolution=gbest % 全局最优个体fmin % 目标函数function f=cost(x)f=1.10471*x(1)2*x(2)+0.04811*x(3)*x(4)*(14.0+x(2); % 非线性约束function g,geq=constraint(x)% 不等式限制条件Q=6000*(14+x(2)/2);D=

27、sqrt(x(2)2/4+(x(1)+x(3)2/4);J=2*(x(1)*x(2)*sqrt(2)*(x(2)2/12+(x(1)+x(3)2/4);alpha=6000/(sqrt(2)*x(1)*x(2);beta=Q*D/J;tau=sqrt(alpha2+2*alpha*beta*x(2)/(2*D)+beta2);sigma=/(x(4)*x(3)2);delta=/(30*106*x(4)*x(3)3);F=4.013*(30*106)/196*sqrt(x(3)2*x(4)6/36)*(1-x(3)*sqrt(30/48)/28); g(1)=tau-13600;g(2)=si

28、gma-30000;g(3)=x(1)-x(4);g(4)=0.10471*x(1)2+0.04811*x(3)*x(4)*(14+x(2)-5.0;g(5)=0.125-x(1);g(6)=delta-0.25;g(7)=6000-F;% 如果没有等式约束，则置geq=；geq=; % APSO Solverfunction gbest,fbest=pso_mincon(fhandle,fnonlin,Lb,Ub,para)if nargin<=4, para=20 150 0.95;endn=para(1);% 粒子种群大小time=para(2); %时间步长，迭代次数gamma=

29、para(3); %gama参数 scale=abs(Ub-Lb); %取值区间% 验证约束条件是否合乎条件if abs(length(Lb)-length(Ub)>0, disp('Constraints must have equal size'); returnend alpha=0.2; % alpha=0,1粒子随机衰减因子 beta=0.5; % 收敛速度(0->1)=(slow->fast); % 初始化粒子群best=init_pso(n,Lb,Ub); fbest=1.0e+100;% 迭代开始for t=1:time, %寻找全局最优个体

30、for i=1:n, fval=Fun(fhandle,fnonlin,best(i,:); % 更新最有个体 if fval<=fbest, gbest=best(i,:); fbest=fval; end end % 随机性衰减因子 alpha=newPara(alpha,gamma); % 更新粒子位置 best=pso_move(best,gbest,alpha,beta,Lb,Ub); % 结果显示 str=strcat('Best estimates: gbest=',num2str(gbest); str=strcat(str,' iteration

31、='); str=strcat(str,num2str(t); disp(str); fitness1(t)=fbest; plot(fitness1,'r','Linewidth',2) grid on hold on title('适应度')end % 初始化粒子函数function guess=init_pso(n,Lb,Ub)ndim=length(Lb);for i=1:n, guess(i,1:ndim)=Lb+rand(1,ndim).*(Ub-Lb); end %更新所有的粒子 toward (xo,yo)function

32、 ns=pso_move(best,gbest,alpha,beta,Lb,Ub)% 增加粒子在上下边界区间内的随机性n=size(best,1); ndim=size(best,2);scale=(Ub-Lb);for i=1:n, ns(i,:)=best(i,:)+beta*(gbest-best(i,:)+alpha.*randn(1,ndim).*scale;endns=findrange(ns,Lb,Ub); % 边界函数function ns=findrange(ns,Lb,Ub)n=length(ns);for i=1:n, % 下边界约束 ns_tmp=ns(i,:); I=

33、ns_tmp<Lb; ns_tmp(I)=Lb(I); % 上边界约束 J=ns_tmp>Ub; ns_tmp(J)=Ub(J); %更新粒子 ns(i,:)=ns_tmp; end % 随机性衰减因子function alpha=newPara(alpha,gamma);alpha=alpha*gamma; % 带约束的d维目标函数的求解function z=Fun(fhandle,fnonlin,u)% 目标z=fhandle(u); z=z+getconstraints(fnonlin,u); % 非线性约束 function Z=getconstraints(fnonlin

34、,u)% 罚常数 >> 1PEN=1015;lam=PEN; lameq=PEN; Z=0;% 非线性约束g,geq=fnonlin(u); %通过不等式约束建立罚函数for k=1:length(g), Z=Z+ lam*g(k)2*getH(g(k);end% 等式条件约束for k=1:length(geq), Z=Z+lameq*geq(k)2*geteqH(geq(k);end % Test if inequalities function H=getH(g)if g<=0, H=0; else H=1; end % Test if equalities holdf

35、unction H=geteqH(g)if g=0, H=0;else H=1; end6.模拟退火算法 TSP问题2015-06-15 已知敌方100个目标的经度、纬度如下表所示： 经度和纬度数据表53.7121 15.304651.1758 0.032246.3253 28.275330.3313 6.934856.5432 21.418810.8198 16.252922.7891 23.104510.1584 12.481920.1050 15.45621.9451 0.205726.4951 22.122131.4847 8.964026.2418 18.176044.0356 13

36、.540128.9836 25.987938.4722 20.173128.2694 29.001132.1910 5.869936.4863 29.72840.9718 28.14778.9586 24.663516.5618 23.614310.5597 15.117850.2111 10.29448.1519 9.532522.1075 18.55690.1215 18.872648.2077 16.888931.9499 17.63090.7732 0.465647.4134 23.778341.8671 3.566743.5474 3.906153.3524 26.725630.81

37、65 13.459527.7133 5.070623.9222 7.630651.9612 22.851112.7938 15.73074.9568 8.366921.5051 24.090915.2548 27.21116.2070 5.144249.2430 16.704417.1168 20.035434.1688 22.75719.4402 3.920011.5812 14.567752.1181 0.40889.5559 11.421924.4509 6.563426.7213 28.566737.5848 16.847435.6619 9.933324.4654 3.16440.7

38、775 6.957614.4703 13.636819.8660 15.12243.1616 4.242818.5245 14.359858.6849 27.148539.5168 16.937156.5089 13.709052.5211 15.795738.4300 8.464851.8181 23.01598.9983 23.644050.1156 23.781613.7909 1.951034.0574 23.396023.0624 8.431919.9857 5.790240.8801 14.297858.8289 14.522918.6635 6.743652.8423 27.28

39、8039.9494 29.511447.5099 24.066410.1121 27.266228.7812 27.66598.0831 27.67059.1556 14.130453.7989 0.219933.6490 0.39801.3496 16.835949.9816 6.082819.3635 17.662236.9545 23.026515.7320 19.569711.5118 17.388444.0398 16.263539.7139 28.42036.9909 23.180438.3392 19.995024.6543 19.605736.9980 24.39924.159

40、1 3.185340.1400 20.303023.9876 9.403041.1084 27.7149 我方有一个基地，经度和纬度为（70，40）。假设我方飞机的速度为1000公里/小时。我方派出一架飞机从基地出发，侦察完敌方所有目标，再返回原来的基地。在敌方每一目标点的侦察时间不计，求该架飞机所花费的时间。 这是一个旅行商问题。我们依次给出基地编号为1，敌方目标依次编号为2，3，.,101,最后我方基地再重复编号102.距离矩阵，其中表示i,j两点的距离，i,j=1,2.,102，这里D为实对称矩阵。则问题是求从一个点1出发，走遍所有中间点，达到点102的一个最短路径。 上面问题中给出的是

41、地理坐标（经度和纬度），我们必须求两点间的实际距离。设A,B两点的地理坐标分别为（x1,y1）,(x2,y2)，过A,B两点的大圆的劣弧即为两点的实际距离。以地心为坐标原点O,以赤道平面为XOY平面，以0度经线圈所在的平面为XOZ平面建立三维直角坐标系。则A,B两点的直角坐标分别为：A（Rcosx1cosy1,Rsinx1siny1）B（Rcosx2cosy2,Rsinx2siny2）其中R=6370为地球半径。A,B两点的实际距离为化简得：d=Rarccoscos(x1-x2)cosy1cosy2+siny1siny2求解模拟退火的步骤见第二十三章 现代优化算法求解程序 monituihuo

42、.mclc, clearsj0=load('sj.txt'); %加载100个目标的数据，数据按照表格中的位置保存在纯文本文件sj.txt中x=sj0(:,1:2:8);x=x(:);y=sj0(:,2:2:8);y=y(:);sj=x y; d1=70,40; sj=d1;sj;d1; sj=sj*pi/180; %角度化成弧度d=zeros(102); %距离矩阵d初始化%计算距离矩阵for i=1:101 for j=i+1:102d(i,j)=6370*acos(cos(sj(i,1)-sj(j,1)*cos(sj(i,2)*cos(sj(j,2)+sin(sj(i,2

43、)*sin(sj(j,2); endendd=d+d'path=;long=inf; %巡航路径及长度初始化rand('state',sum(clock); %初始化随机数发生器%求较好的初始解for j=1:1000 path0=1 1+randperm(100),102; %随机产生一个路径 %起始点出发走一圈回到起始点 故路径为102个点 temp=0; for i=1:101 temp=temp+d(path0(i),path0(i+1); end if temp<long path=path0; long=temp; endende=0.130;%终止温

44、度L=20000;at=0.999;%降温系数T=1;%初始温度for k=1:L %退火过程c=2+floor(100*rand(1,2); %产生新解c=sort(c); c1=c(1);c2=c(2); %计算代价函数值的增量df=d(path(c1-1),path(c2)+d(path(c1),path(c2+1)-d(path(c1-1),path(c1)-d(path(c2),path(c2+1); if df<0 %接受准则 path=path(1:c1-1),path(c2:-1:c1),path(c2+1:102); long=long+df; elseif exp(-

45、df/T)>rand %以概率接受新解 path=path(1:c1-1),path(c2:-1:c1),path(c2+1:102); long=long+df; end T=T*at; if T<e break; 退火温度 endendpath, long % 输出巡航路径及路径长度xx=sj(path,1);yy=sj(path,2);plot(xx,yy,'-*') %画出巡航路径7. 右端步连续微分方程求解 2015-06-16 考虑下面的具体例子：程序如下：dvp2.mfunction dy=dvp2(t,y)dy=zeros(4,1);p1=.8;c1

46、=400;p2=.9;c2=64;r1=1;r2=0.05;b=1/1000;a1=0.1;a2=0.25;m=1;gam=0.8;dy(1)=y(1)*(r1-b*y(1)-a1*y(2)/(m*y(2)gam+y(1)-y(3)*y(1)*(y(1)>492);dy(2)=y(2)*(-r2+a2*y(1)/(m*y(2)gam+y(1)-y(4)*y(2)*(y(2)>70);dy(3)=0.8*y(3)*(p1*y(1)-c1)*(y(1)>492);dy(4)=0.8*y(4)*(p2*y(2)-c2)*(y(2)>70);主函数Solve2.mt,y=ode

47、45('dvp2',0,20,100 50 0.1 0.1);ca=input('请输入画:x y E1 E2用1 2 3 4 表示:');if ca=1plot(t,y(:,1),'r')xlabel('t');ylabel('x(t)');title('种群x随时间的变化图')elseif ca=2plot(t,y(:,2),'r')xlabel('t');ylabel('y(t)');title('种群y随时间的变化图')else

48、if ca=3plot(t,y(:,3)-0.1,'r')xlabel('t');ylabel('E1(t)');title('E1随时间的变化图')else plot(t,y(:,4),'r') xlabel('t'); ylabel('E2(t)'); title('E2随时间的变化图')End8. 多元方差分析 2015-06-17 某养鸡场的鸡蛋育成期的配合饲料主要由5种成分（玉米、麦皮、豆饼、鱼粉和食盐）组成，分别记为A,B,C,D,E，为了研究饲料配方对鸡产蛋效果的影响，对哥成分均选取3个水平，进行了5因素3水平的正交实验，通过试验找出饲料的最佳配方，试验要求考虑交互作用AB、AC和AE。5个因素的水平表如下：因素ABCDE161.56.56302668950.1370.6141590.25对这样的5个因素3水平试验的正交实验表如下：试验号ABCED产蛋量y11111156921122255431133363741212356651223156561231264871313258181321356891332153510211115931121222615122133