捕食者与被捕食者问题微分方程与matlab
数学建模狐狸野兔问题
狐狸野兔问题摘要：封闭自然环境中的狐狸和野兔存在捕食与被捕食关系，本题旨在通过对自然状态下两物种数量变化规律的分析，推测加入人类活动（即人工捕获）时两物种数量的变化，进而得出人类活动对自然物种的影响，为人类活动提供参考，使其在自然允许的范围内，促进人与自然和谐相处。
对于问题一，首先建立微分方程，描述两物种数量随时间变化的Volterra 模型()0,0,0,021212211>>>>⎪⎪⎩⎪⎪⎨⎧+-=-=r r k k xyr y k dtdy xy r x k dtdx并用解析法求得狐狸与野兔数量的关系 ()()2211k r xk r yxeyec --=为直观反映两物种数量随时间的变化规律，选取三组有代表性的初值，利用Matlab 软件绘图。
在狐狸和野兔随时间的变化图像中，大致得出其数量呈周期变化，为进一步检验周期性，再用Matlab 绘图做出狐狸与野兔数量的关系图，得到封闭曲线，因此分析结果为：狐狸和野兔的数量都呈现周期性的变化，但不在同一时刻达到峰值。
对于问题二，利用数值解法，令模型中两式皆为0，即求得狐狸和野兔数量的平衡状态。
且由问题一中狐狸与野兔数量的关系图知野兔和狐狸的平衡量恰为他们在一个周期内的平均值。
对于问题三，在Volterra 模型基础上引入人工捕获系数。
只捕获野兔时，野兔的自然增长率降低，狐狸自然死亡率增加，改进后模型同问题二处理方式一样，求得平衡状态，得出结论：捕获野兔时，狐狸数量减少，野兔数量反而增加，即Volterra 原理：为了减少强者，只需捕获弱者。
只捕获狐狸时，分析方法与只捕获野兔时相同，并得出野兔狐狸数量皆增加的结论。
问题三为自然界人类捕获生物提供了新的思路，即可以在正常允许范围内，为了达到减少某一种群数量的目的，相应的捕获其食饵，或适度地捕获捕食者使捕食者与被捕食者的数量都有所增加。
关键词：Volterra 模型Matlab 软件解析法周期性一、问题重述在一个封闭的大草原里生长着狐狸和野兔。
数学建模实验二：微分方程模型Matlab求解与分析
实验二： 微分方程模型Matlab 求解与分析一、实验目的[1] 掌握解析、数值解法，并学会用图形观察解的形态和进行解的定性分析； [2] 熟悉MATLAB 软件关于微分方程求解的各种命令；[3] 通过范例学习建立微分方程方面的数学模型以及求解全过程； [4] 熟悉离散 Logistic 模型的求解与混沌的产生过程。
二、实验原理1. 微分方程模型与MATLAB 求解解析解用MATLAB 命令dsolve(‘eqn1’,’eqn2’, ...) 求常微分方程（组）的解析解。
其中‘eqni'表示第i 个微分方程，Dny 表示y 的n 阶导数,默认的自变量为t 。
（1） 微分方程 例1 求解一阶微分方程 21y dxdy+= (1) 求通解 输入：dsolve('Dy=1+y^2')输出：ans =tan(t+C1)（2）求特解 输入：dsolve('Dy=1+y^2','y(0)=1','x')指定初值为1，自变量为x 输出：ans =tan(x+1/4*pi)例2 求解二阶微分方程 221()04(/2)2(/2)2/x y xy x y y y πππ'''++-=='=-原方程两边都除以2x ，得211(1)04y y y x x'''++-= 输入:dsolve('D2y+(1/x)*Dy+(1-1/4/x^2)*y=0','y(pi/2)=2,Dy(pi/2)=-2/pi','x')ans =- (exp(x*i)*(pi/2)^(1/2)*i)/x^(1/2) +(exp(x*i)*exp(-x*2*i)*(pi/2)^(3/2)*2*i)/(pi*x^(1/2))试试能不用用simplify 函数化简 输入: simplify(ans)ans =2^(1/2)*pi^(1/2)/x^(1/2)*sin(x) （2）微分方程组例3 求解 d f /d x =3f +4g ; d g /d x =-4f +3g 。
捕食者-被捕食者模型中参数估计的新方法
) )
5 6
其中，k (1 k 6)为模型的待定参数。
通过对这个生态系统的观测可以得到若干组捕
食者和被捕食者在 t 时刻各自数目的观测数据（具
体可见2006年全国研究生数学建模竞赛B题中数据
文件DATA1.TXT、DATA2.TXT、DATA3.TXT和
DATA4.TXT）。利用有关数据，解决以下参数估计
ln复相关系数决定系数调整的决定系数回归方程的估计标准误差1000100010000000095718251误差来源离差平方和自由度方差f统计量概率p值回归平方和残差平方和总平方和1142910000114291380970000415810方差分析表变量回归系数取值回归系数的标准误差t统计量概率95的置信区间置信下限置信上限常数项13820000051201000001381913821200000001078100000200020001200000005852100000120001200110000000561110000010001000各参数的显著性检验表1表明回归分析的拟合优度为1000
观测值为参数 5,6 的估计值。对于问题1，取：
5 10, 6 60
综合以上分析结果得到所有参数的高精度估计值：
(1, 2 , 3, 4 , 5, 6 )T (2, 0.2,12, 1,10, 60)T
问题2.观测数据无误差，但2 未知时确定需要观
测数据的最少组数
图1 残差图
从残差图来看，残差围绕着0随机波动，说明使 用该模型比较合理。
由观测数据的散点图看出，捕食者-被捕食者模 型的解具有明显的周期性。即存在常数 T 0 ，使
得t，x(t T ) x(t)且 y(t T ) y(t) 。只要观测序列
认证考试数学建模试验谜底
实验一：食饵与捕食者的相互依存与制约分析 练习1：解微分方程组：,⎥⎦⎤⎢⎣⎡==⎥⎦⎤⎢⎣⎡-+⎥⎦⎤⎢⎣⎡⎥⎦⎤⎢⎣⎡--=⎥⎦⎤⎢⎣⎡3)0(2)0(,)sin (cos 999sin 299999812''v u x x x v u v u Input :function,f=f(x,y)f=[-2,1;998,-999]*y+[2*sin(x);999*(cos(x)-sin(x))];[x,y]=ode23('f',[0,10],[2,3]);xpauseyplot(x,y)Result :,练习2：解二阶微分方程4|,0)0(,0202022===++=t dt dg g g dt dg a dt g d ω Input :●function,f=zuni(t,y);global,a,wo;A=[0,1;,,-wo^2,-2*a];f=A*y;●global,a,wo;a=5;wo=1;[t,y]=ode45('zuni',[0::2*pi],[0,4])yplot(t,y(:,1),'b')hold,onplot(t,y(:,2),'g')output:,2．利用matlab实现食饵与捕食者系统的仿真：V olterra食饵与捕食者模型：⎪⎩⎪⎨⎧+-=+-=-=-=bxydy bx d y t y axy rx ay r x t x )()()()(.. 取2,25,02.0,1.0,5.0,100======y x b a d rInput:●function,dx=shier(t,x)global,r,d,a,b;dx=zeros(2,1);dx(1)=(r -a*x(2))*x(1)dx(2)=(-d+b*x(1))*x(2)●global,r,d,a,b;r=1;d=;a=;b=;ts=0::15;x0=[25,2];[t,x]=ode45('shier',ts,x0);xplot(x(:,1),x(:,2))3. 利用matlab 实现两种群共生系统的仿真：仿照上例编写程序。
Matlab上机实验题及参考解答
Matlab上机实验题及参考解答目录实验一Matlab初步实验 (2)一matlab基本功能介绍 (2)二Matlab扩展功能 (2)三练习 (2)四练习题参考解答 (3)实验二概率模型实验 (5)一复习 (5)二事件的响应 (5)三Matlab中随机数字的生成与处理 (5)四练习 (5)五练习题参考解答 (5)实验三插值与拟合 (7)实验四线性规划与非线性规划 (8)4.1 实验目的 (8)4.2 实验内容 (9)4.3 综合练习 (10)4.4 课外作业 (11)实验五数值计算 (12)5.1 实验目的 (12)5.2 实验内容 (12)4.3 综合练习 (15)4.4 课外作业 (15)实验六计算机图像处理 (16)6.1 实验目的 (16)6.2 实验内容 (16)6.3 综合练习 (17)6.4 课外作业 (19)实验七综合练习 (19)7.1 实验目的 (19)7.2 实验内容 (19)7.3 综合练习 (20)7.4 课外作业 (21)实验一 Matlab 初步实验 一 matlab 基本功能介绍1 编程环境2语法规范：for … end; if …else if …end; 3 矩阵运算 4 图形绘制二 Matlab 扩展功能1 编程练习：(1) 绘出序列kk x x r r 0(1),0.2083=+=；(2) 绘出曲线rtx t x e t 0(),0=>2 扩展功能(1) 矩阵中全部数据、部分数据的截取、更改； (2) 矩阵的初始化与赋值如：A=zeros(5,5); A(2:2:)=[1,2 3 4 5] 3 微积分基础(见实验4) 符号计算三 练习(课上编程完成下列练习，课后上机验证) 1 求和S=1+2+3+…+100; 2 求和e 1111!2!10!1...=++++3求和S 1112310!1...=++++4设A 234576138⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦, 求A 的逆、特征值和特征向量；验证Ax=λx 5 画函数图()011mrtm x x t x e x -=⎛⎫+- ⎪⎝⎭6 展开 (x-1)(x-2)…(x-100)7 因式分解 x 8—y 8; 因数分解200520068 求极限312lim +∞→⎪⎭⎫⎝⎛++n n n n9 )](sin[cos 22x x y += 求dxdy10 求积分x xdx 10ln ⎰11 求积分3⎰并且画出所求的平面区域12 设x+2y=1, 2x+3y=6, y=2x 2, 画出各个方程图形，求出曲线交点.四 练习题参考解答%MatlabTrain1.m clear all % 2nd e=1; temp=1; for I=1:1:10temp=temp*I; e=e+1/temp; end e%%%%%%%%%%% clear all % 3nd S=0; temp=1;for I=1:1:100temp=temp*I; endfor J=1:1:temp S=S+1/J; end S%%%%%%%%%%%%%% clear all % 11ndx=linspace(0,4); y=1./sqrt(x.^5+1); plot(x,y) for t=1:0.1:3yt=1./sqrt(t.^5+1);hold online([t,t],[0,yt]);end%fill(t,yt,'b') %%%%%%%%%%%%% clear all% 12ndx=linspace(-2,2);y=[0.5-0.5*x; 2-2/3.*x; 2*x.^2]; plot(x,y)grid实验二概率模型实验一复习1 小结上次编程练习中存在的问题，讲述部分习题答案2 画图命令介绍：line二事件的响应(1) 获取鼠标的位置%MatlabTrain2.mclear all% 鼠标响应p=ginput(3)plot(p(:,1),p(:,2),'r*')(2) 键盘输入相应t=input('How many apples? t=');m=t+3三Matlab中随机数字的生成与处理1 随机数的生成2 产生随机数字3 产生某区间的整数4 生日模拟问题的Montecaro法设计技术、思路学生尝试编程四练习(1) 编程验证人数在不同年龄段的生日的概率计算(2) 编程实现游戏”聪明伶俐100分”(3) 编程实现两家电影院的座位数问题(4) 编程实现某图形面积的计算五练习题参考解答(1) 生日问题程序示例：%birthPro.mn=0;nStudents=30;for I=1:1000 %how many times testy=0;x=1+floor(365*rand(1,nStudents));%get nStudents random numbersfor J=1:nStudents-1for K=J+1:nStudentsif x(J)==x(K)y=1;break;endendendn=n+y;%count, n times of that there are two people's dirthday in the same dayendfreq=n/I % caculating the frequently(2) 编程实现游戏”聪明伶俐100分”参考答案%MatlabTrain2.mclear all% 鼠标响应x=floor(10*rand(1,4))t=input('填入四个数字[n1 n2 n3 n4]=');flag=0;A=0;B=0;for I=1:1:8flag=flag+1;A=0;B=0;if t==xswitch flagcase 1disp('聪明绝顶！');case 2disp('聪明！');case 3disp('有点聪明！');case 4disp('还可以！');case 5disp('聪明伶俐100分！');case 6disp('聪明伶俐90分！');case 7disp('聪明伶俐85分！');case 8disp('聪明伶俐80分！');otherwisedisp('赫赫！');endbreak;endfor J=1:1:4for K=1:1:4if x(J)==t(K) & J==KA=A+1;else if x(J)==t(K) & J~=KB=B+1;endendendends='AABB';s(1)=INT2STR(A);s(3)=INT2STR(B);disp(s);t=input('不重复填入四个数字[n1 n2 n3 n4]=');endif flag>0disp('太烂了! 正确答案是:');xend实验三插值与拟合一复习讲述聪明伶俐100分的编程中的问题二插值三拟合课堂练习2 某之股票价格from 2003 09 01 to 2004 01 02，试进行插值、拟合%TimerS.m%from 2003 09 01 to 2003 01 02clear all;dataST=[15.09 14.7514.95 14.722.88 21.8619.82 19.09];plot(dataST)四课外练习112）进行多项式拟合，求出拟合多项式，并求出多项式在t=4, 5处的值.实验四线性规划与非线性规划4.1 实验目的1 用Matlab求解线性规划2 用Matlab求解非线性规划4.2 实验内容4.2.1 线性规划求解实用格式：x=lp(c, A, b, xLB,xUB,x0,nEq)可以求解下列线性规划模型：min f=c’xs.t. Ax=<=b(其中前nEq个约束为等式约束,即等式约束的个数，其余是不等式约束<=) xLB<=x<=xUB函数中x0参数是算法迭代的初始点，任意取值例1 求解下列线性规划1）123123123123min2..360210200,1,2,3jz x x xs t x x xx x xx x xx j=--+⎧⎪++≤⎪⎪-+≤⎨⎪+-≤⎪≥=⎪⎩，2）1235635623416367min..3621060,1,,7jz x x x x xs t x x xx x xx xx x xx j=-++-⎧⎪++=⎪⎪+-=⎪⎨-+=⎪⎪++=⎪≥=⎪⎩例1求解示例c=[-2 -1 1]';%book page 72 Number 16-1A=[3 1 1;1 -1 2;1 1 -1];b=[60 10 20]';xlb=[0 0 0]';xub=[inf inf inf]';x0=[0 0 0]'; x=lp(c,A,b,xlb,xub,x0,0)% x=(15 5 0)'例2 求解示例c2=[1 -1 1 0 1 -1 0]';%book page 72 Number 16-3A2=[0 0 3 0 1 1 0;...0 1 2 -1 0 0 0;...-1 0 0 0 0 1 0;...0 0 1 0 0 1 1];b2=[6 10 0 6]';xlb2=[0 0 0 0 0 0 0]';xub2=[inf inf inf inf inf inf inf]';x02=[0 0 0 0 0 0 0]';x2=lp(c2,A2,b2,xlb2,xub2,x02,4)% unbounded4.2.2 非线性规划1）命令格式1：[X, OPTIONS]=constr(‘FUN’, X, OPTIONS,VLB,VUB)2）命令格式2：X=FMINCON(FUN,X0,A,B,Aeq,Beq)% minimizes FUN subject to the linear equalities% Aeq*X = Beq as well as A*X <= B. (Set A=[] and B=[] if no inequalities exist.)例2 求解非线性规划y x x x x s t x3211221min22 ..1=++-≤-求解示例%unconop.mfunction y=unconop(x)y=x(1).^3+2*x(1).*x(2)+2*x(2).^2;%book page 148 ex.7-1 后建立调用函数xx=fmincon('unconop',[0 0]',[-1 0],-1,[],[])%book page 148 ex.7-1 4.3 综合练习学生独立编写程序，求解一个含有2个变量的线性规划问题，要求：1)编写程序，把可行域画上阴影；2)求出最优解，在可行域上标出最优解；3)求出基本解，并在上图中表示出来；4)求出基本可行解，观察单纯形方法迭代时，顶点的变化.可行域画图与表出阴影示例：syms x y[u(1),v(1)]=solve('y=x+2','y=2*x');%求出交点坐标[u(2),v(2)]=solve('y=-x+2','y=2*x');[u(3),v(3)]=solve('y=x+2','y=-x+2');x=linspace(0,3,5); %直线作图y=[2*x;-x+2;x+2];line(x,y); gridpatch(double(u),double(v),'b'); 运行结果：4.4 课外作业1 求解线性规划131223min ..250.530,1,2,3i x x s t x x x x x i +⎧⎪+≤⎪⎨+=⎪⎪≥=⎩ （1） 求解线性规划；x *=()（2） 目标函数中c 1由1变为（-1.25）时求最优解；（3） 目标函数中c 1由1变为（-1.25），c 3由1变为2时求最优解；（4） 约束条件中53b ⎛⎫= ⎪⎝⎭变为21b -⎛⎫'= ⎪⎝⎭时，求解；（5） 约束条件中53b ⎛⎫= ⎪⎝⎭变为23b ⎛⎫'= ⎪⎝⎭时，求解[刁在筠，运筹学（第二版），高等教育出版社，2004，01 p74第20题]2 求解非线性规划y x x x x x x x 3221122233min 2223=++++ 注：无约束非线性规划问题, 命令：fminunc子函数% unconop.mfunction y=unconop(x)y=x(1).^2+2*x(1).*x(2)+2*x(2).^2+2*x(2).*x(3)+3*x(3).^2;%book page 148 ex.7-1 主函数：xx=fminunc('unconop',[0.1 0.1 1]')思考：绘出两个变量的线性规划问题的可行域、标出可行的整数解和求出可行解；演示单纯形方法的迭代过程，如j z x x s t x x x x x j 121212min 2..360200,1,2=--⎧⎪+≤⎪⎪+≤⎨⎪⎪≥=⎪⎩实验五 数值计算5.1 实验目的1 掌握代数数值计算2 掌握常微分方程数值计算5.2 实验内容5.2.1 关于多项式设多项式1110()n n n n p x a x a x a x a --=++++表示为110[,,,,]n n p a a a a -=1）求多项式的根 roots(p) %求出p(x)=0的解。
食饵捕食者研究
具有干扰因素的食饵-捕食者模型分析目录目录摘要…………………………………………………………………第一部分前言………………………………………………………1.1 生态数学的的研究背景及发展…………………………………1.2 基础知识…………………………………………………………第二部分 Lotka-Volterra模型的改进及其稳定性的研究…………2.1Lotka-Volterra模型………………………………………………2.2模型的研究对象及改进…………………………………………2.3 模型的稳定性的研究……………………………………………第三部分数值模拟3.1利用matlab对模型进行了数值模拟……………………………3.2模型缺陷…………………………………………………………第四部分总结………………………………………………………致谢…………………………………………………………………参考文献………………………………………………………………第一部分 前言1.1 生态数学的的研究背景及发展生态系统具有稳定性、可测性和可控性三大属性，是多层次的、多因子的、多变量的系统，只用常规的定性描述和一般的数理统计，搞不清楚它的内在规律，运用数学模型对生态系统实行管理、预测和调控，使其持续稳定发展是现代生态学研究的重要领域。
种群动力学是生态学的一个重要分支.它广泛地利用数学思想加积分方程、差分方程、泛函微分方程、偏微分方程、算子理论等数学学科中的理论和方法，通过数学建模研究生物种群的生存条件、生物种群与环境之间相互作用的过程、生物种群的演变和发展趋势.揭示生物种群的变化规律，合理利用资源，促进生态平衡这是迄今为止数学在生态学中应用深入，发展最为系统和成熟的分支，种群 动 力 学的研究有着悠久的历史.早在1798年，Malthus 在研究人类的增长时，他引入数学方法，建立了最早的连续确定模型一一Malth 。
模型)(/)(t rN dt t dN =这是一个单种群模型.它反应了人类数量的变化，在t 不很长时是比较符合实际的，但当+∞→t 时种群规模将无限增长是不合实际的，究其原因在于它没有考虑到有限的资源对种群增长的制约作用.针对这个模型，后人不断分析各种因素的影响，完善和改进这一模型，使之能较好地反应人口(单种群)的变化规律，如P. F.Verhulst(1938年)建立的Logistic 模型)/)(1)((/)(k t N t rN dt t dN -=E.M .W right(1945年)建立的有确定时滞的Logistic 模型)/)(1)((/)(k r t N t rN dt t dN --=P. M. Nisbet 和W. S. C. Gurney(1984年)建立的具有生理阶段结构(stagestructure) 模型以及H. 1. Freedman 研究的具有斑块迁移的单种群模型等，无一不是对Malthus 模型的完善和扩展，极大地推动了种群动力学的发展现 实 世 界中种群不可能单独存在，它必与相关种群相互作用，相互依存.Lotka-Volterra 模型是种群动力学中最为经典和重要的两种群相互作用的动力学模型，该模型分别由意大利数学家Volterra(1923年)解释鱼群变化规律和美国种群学家Lotka(1921年)在研究化学反应时提出。
微分方程数值解
微分方程数值解4.1当常微分方程能解析求解时，可利用Matlab符号工具箱中的功能找到精确解. 见下例求解方程,,,. 键入： yyy，，,20syms x y %定义符号变量diff_equ= ‘D2y+2*Dy-y=0’; %D2y表示,,,Dy= y,yy=dsolve (diff_equ, ‘x’) %定义x为自变量 y=cl*exp ((2^(1/2)-1)*x+c2*exp (-(2^(1/2)+1)*x)%表达式中含c1与c2，表示通解.%初始条件为y (0)=0,,y(0)=1时，按如下方式调用 y=dsolve (diff_equ,‘y (0)=0’, ‘Dy (0)=1’, ‘x’) y=1/4*2^(1/2)*exp ((2^(1/2)-1)*x)—1/4*2^(1/2)*exp (-(2^(1/2)+1)*x)%画出函数y=y (x)的图形ezplot (y,[-2,2])图形具体形式请上机试之.在方程无法获得解析解的情况下，可方便地获得数值解. 下面的例子说明用Matlab求数值解的方法及应注意的问题.[例1] 求解范德堡（vander pol）方程2dxdx2,,，,,(1)0xx 2dtdt求解高阶方程，必须等价地变换为一阶微分方程组，对本例，通过定义两个新的变量，实现这一变换yxydxdt1,2/,, 则令 dydty1/2,2dydtyyy2/(11)*21,,,,编写求解程序分为两部分，第一部分为待求解的方程，存盘的文件名为,待求解方程的函数名.m,，第二部分为求解主程序，本例中取名为main1.m.首先编写待求解方程的文件. 文件存盘名为“vdpol.m”. M,function yprime=vdpol(,)tyyprime (1)=y (2);; (1(1)^2)*(2)(1),,yyymu=2 yprime=[yprime (1);yprime (2)]; yprime (2)=mu*说明函数yprime=vdpol中. 定义为自变量，的形式取决于求解方程的阶数，本(,)tyyt 例中，,为解向量，为导数向量. yprime， y(2)yyyy,[(1),(2)],(1)(1)(1),y yprime,,，函数返回vander pol方程的导数列向量. 因为所求结果为方程数值解，(2)(1),y所以各向量维数只有在主程序求解时定下精度后才能确定.主程序定名为main1.m，你可用你所喜欢的其它名子，但vdpol.m除外. clear functions%调试程序时，放置这一语句是必要的. 它清除前边已编译的存在于内存中的废弃程序[]=ode23 (‘vdpol’,[0,30],[1,0]); ty,y1=y (:,1); %解曲线.y2=y (:,2); %解曲线的导数.polt ( ‘_ _’) tyty,1,,2,说明龙格_库塔的2阶与4阶改进型求解公式的实现，其指令分别为：[]=ode23 (‘f’,tsx,0,options) tx,[]=ode45 (‘f’,tsx,0,options) tx,其中可由系统依据精度要求自动设定，亦可由使用者依据实际需要自己确定，分别说明之. ts（1）若令tstttf,[0,1,,]，则输出在指定时刻tttf0,1,,给出，当tstktf,0::时，输出在区间[0,]ttf的等分点上给出，为步长. k（2）若tsttft,[0,],0为自变量初值，tf为终值，此时，options决定自变量的维数，t中的时间点不是等间隔的，这是为了保证所需的相对精度，积分算法改变了步长. 用于t,3,6设定误差限的参数options可缺省，此时系统设定相对误差为，绝对误差为，若1010自行设定误差限，可用如下语句：options=odeset (‘reltol’,, ‘abstol’,) rtat这里的与分别为设定的相对与绝对误差. rtat须注意的是无论用哪种方法确定ttf0,的取值方式，必须由使用者确定且应与相匹配. x0t,y0[1,0],ts,[0,30]y(0)1,y(0)0,为初始条件，本例中，因为，这意味着解曲线，，x0一般说，当解nnn个未知函数的方程组时，为维向量，共含有个初始条件. x0两个输出参数是列向量xx与矩阵，它们具有相同的行数，而矩阵的列数等于方程t组的个数，本例中y(:,1)y(:,2)的列数为2，其中，为自变量上各点函数值，为上各ytt点导数值.最后，提请读者注意的是：ode45也不总是比ode23好，在很多时候，低阶算法更有效，有关微分方程数值解法的更进一步信息，请参考数值分析方面的书籍. 有些参考书提供了一些关于算法选择和如何处理那些时间常数变化范围大的病态方程的非常实用的算法. 4.2 -设有一阶方程与初始条件,yfxy,(,), （4.1） ,yxy(),00,其中适当光滑，关于满足Lipschitz条件，即存在使 fxy(,)LyfxyfxyLyy(,)(,),,,1212则（4.1）式的解存在且惟一.关于yyx,()的解析解一般难以求到或根本无解析解，因而，实际问题中，通常，采用差分的方法. 在一系列离散点xxx,,,,yyy,,,,上寻求其数值近似解. 12n12n相邻两个节点间的间距xxnhn,，,,1,2,hxx,,称为步长，一般地取等步长，则hn0nnn，11、欧拉方法在区间[,]xx上用差商 nn，1yxyx()(),nn，1 h代替（4.1）式中,[,]xxxxxy，对fxy(,)中在上取值还是，而形成向前欧拉公式nn，1nn，1与向后欧拉公式.（1）向前欧拉公式xfxy(,)取左端点，得如下公式 nyxyxhfxyx()()(,()),，（4.2） nnnn，1从yxy(),x点出发，由初值代入（4.2）求得 000yyhfxy,，(,) （4.3） 1000反复利用（4.2），有yyhfxyn,，,(,) 0,1,2, （4.4） nnnn，1其几何意义如图4.1所示.y 图中yyx,()为方程（4.1）的精确 P P43P 2解曲线，其上任意点(,)xy处切线斜率为误差 P 1yyx,() 32Pxy(,)fxy(,). 从初值点出发，用该 P000 0y 0点斜率fxy(,)xx,作一直线段，在 001yyx,() yx() 3处得到Pxy(,)y，由（4.2）式确定， 1111y 3再从Pxy(,)fxy(,)出发，以为斜率 11111作直线段，在xx,Pxy(,)处得到， 2222xxxxx x O03124PPP， 012作为积分曲线yx()的近似，用图4.1 yyx,()n，1这一过程继续下去，形成折线表示在xy处的精确值，为解的近似值，不难得到 n，1n，12h32,,()()()()yxyyxOhOh,,，, nnn，，112P，1这一误差称为局部截断误差. 若一种算法局部截断误差为Oh()，则称该算法具有阶P精度，所以向前欧拉公式具有1阶精度.（2）向后欧拉公式若[,]xxxx中取中的，则有如下公式： fxy(,)nn，1n，1yyhfxyn,，,(,) 0,1,2, （4.5） nnnn，，，111称式（4.5）为向后欧拉公式，因为此式中y未知，故称其为隐式公式，无法用其直n，1接计算y，一般用向前欧拉公式产生初值. n，1(0)yyhfxyn,，,(,) 0,1,2, 11nnnn，，再按下式迭代(1)()kk，yyhfxykn,，,,(,),0,1,,0,1, nnnn，，，111其误差估计如下2h32,,()()()()yxyyxohoh,,，, nnn，，112精度亦为1阶，将向前欧拉公式（4.4）与向后欧拉公式（4.5）及它们的误差的几何说明作一对比，是十分有益的，见图4.2.y 为讨论局部截断误差，在图4.2中设点APxy(,)落在积分曲线yyx,()上，按式 nnnyyx,() （4.4）及式（4.5）分别得 ,P点为与， ABn，1 B且P AB,yyx,()点一定在积分曲线上相应 n点的上、下两边,所以将式（4.4）与（4.5） , 平均之，一定能得到更好的结果.xxx （3）梯形公式 nn，1 将向前与向后欧拉公式加以平均得到所图4.2 谓梯形公式hyyfxyfxyn,，，,[(,)(,)] 0,1,2, （4.6） nnnnnn，，，11123其局部截断误差为Oh()，具有2阶精度.（4）改进的欧拉公式为使计算简单，又免去迭代的繁复，将公式（4.6）简化为两步yyhfxy,，(,)nnnn，1h （4.7） yyfxyfxyn,，，,[(,)(,)], 0,1,2,nnnnn，，11n，12或写为h,yykk,，，()nn，112,2,1nn （4.8） n,0,1,2,kfxy,(,),,211nn，kfxyhk,，(,),,最后指出，上述欧拉方法可推广至微分方程组，如,yfxyz,(,,),,,zgxyz,(,,), ,yxy(),00,,zxz(),,00向前欧拉公式为yyhfxyz,，(,,),nnnnn，1 n,0,1,2, ,zzhgxyz,，(,,),nnnnn，12、龙格_库塔方法由微分中值定理,[()()]/(),01yxyxhyxh,,，,,,, nnn，1又因为,,yxhfxhyxh()(,())，,，，,,,yfxy,(,)，所以 nnn从而有yxyxhfxhyxh()()(),(),，，，,, （4.9） nnnn，1令[,]xx，称其为区间上的平均斜率，由（4.9）可知，给kfxhyxh,，，(,()),,nn，1nn出一种平均斜率，可相应导出一种算法. 向前欧拉公式中，精度低. 改进欧kfxy,(,)nn1拉公式中取[,]xxkfxyfxy,，((,)(,))，精度提高，下面，我们在区间内nn，1nnn，1n，12多取几个点，将其斜率加权平均，就能构造出精度更高的计算公式，公式的推导不再具体给出，只开列具体结果.（1）2阶龙格_库塔公式yyhkk,，，(),,,nn，11122,kfxy,(,) （4.10） 1nn,,21nnkfxahyhka,，，,,(,),0,1,,,1,其中,,,，,,,1,,1a，由于4个未知数只有3个方程，所以解不惟一，若令1222a1，即得改进的欧拉公式，具有2阶精度. ,,,，,,,,1a122（3）4阶龙格_库塔公式只给出精典格式中最常用的一种.h,yykkkk,，，，，(22)nn，11234,6,kfxy,(,),1nn,hh, （4.11） kfxyk,，，(,),21nn22,hh,kfxyk,，，(,)32nn,22,kfxhyhk,，，(,),43nn,其计算精度为4阶4.31、模型与问题[例2] 单摆运动图4.3中一根长的细线，一端固定，另一端悬挂质量为 lm的小球，在重力作用下，小球处于竖直的平衡位置. 现使小球偏离平衡位置一个小的角度，然后使其自由运动，在不 ,考虑空气阻力情形下，小球将沿弧线作周期一定的简谐运动.为平衡位置，在小球摆动过程中，当与平衡位置夹 ,,0角为,mgsin,时，小球所受重力在其动运轨迹的分量为 , , l（负号表示力的方向使减少），由牛顿第二定律可得微分方 ,程,,mltmg,,()sin,, （4.12）设小球初始偏离角度为,，且初速为0，式（4.12）的初 0,始条件为,,,,(0),(0)0,, （4.13） 0 mg 当,不大时，，式（4.12）化为线性常系数微sin,,,0分方程图4.3g,,,,，,0 （4.14） l解得g,,()costt, （4.15） 0ll简谐运动的周期为. T,2,g现在的问题是：当,较大时，仍用近似，误差太大，式（4.12）又无解析解，,sin,0试用数值方法在,,30,10,,两种情况下求解，画出的图形，与近似解（4.15）,()t00比较，这里设. l,25cm[例3] 捕食与被捕食当鲨鱼捕食小鱼，简记为乙捕食甲，在时刻，小鱼的数量为，鲨鱼的数量为，xt()yt()t当甲独立生存时它的（相对）增长率与种群数量成正比，即有,r，为增长率，xtrxt()(),而乙的存在使甲的增长率r减少，设减少率与乙的数量成正比，而得微分方程,xtxtraytrxaxy()()(()),,,, （4.16）比例系数a反映捕食者掠取食饵的能力.乙离开甲无法生存，设乙独自存在时死亡率为,，ytdyt()(),,，甲为乙提供食物，d使乙的死亡率降低，而促其数量增长，这一作用与甲的数量成正比，于是yt()满足 d,ytydbxtdybxy()(()),,，,,，（4.17）比例系数反映甲对乙的供养能力，设若甲，乙的初始数量分别为 bxxyy(0);(0),, （4.18） 00则微分方程（4.16），（4.17）及初始条件（4.18）确定了甲，乙数量xt()、yt()随时间变化而演变的过程，但该方程无解析解，试用数值解讨论以下问题：（1）设rdabxy,,,,,,1,0.5,0.1,0.02,25,2，求方程（4.16），（4.17）在00条件（4.18）下的数值解，画出xtyt(),()的图形及相图(,)xy，观察解xtyt(),()的周期变化，近似确定解的周期和的最大、小值，近似计算在一个周期内的平均值. xy,xy,（2）从式（4.16）和（4.17）消去得到 dtdxxray(),, （4.19） dyydbx(),，解方程（4.19），得到的解即为相轨线，说明这是封闭曲线，即解确为周期函数.（3）将方程（4.17）改写为,1yxtd()[],，（4.20） by在一个周期内积分，得到xt()yt()一周期内的平均值，类似可得一周期内的平均值，将近似计算的结果与理论值比较.2、实验（1）方程（单摆问题）,,mlmg,,,,sin, ,,,,,(0),(0)0,,0,无解析解，为求其数值解，先将其化为等价的一阶方程组. 令,xx,,,,,，方程化为 12,,xx,12,g,,21,,xxsin ,l,102,,xx(0),(0)0,,,其中，gl,,9.8,25,,为（弧度）与（弧度）两种情况，具100.1745,300.5236,0体编程如下：先以danbai.m为文件名存放待解方程. 键入： function xdot =danbai (t,x) %x=[x (1),x (2)],g=9.8;1=25; %x (1)为解向量, x (2)是导数 xdot (1)=x (2); %xdot(1)=,,(1)= x,xdot (2)=-g/1*sin (x (1)); %xdot (2)=,, ,xdot=[xdot (1);xdot (2)]; %必须将导数向量以列向量形式给出.再以主程序(求数值解)调用待求方程，主程序用main2.m为文件名存盘，其代码如下：clear functions[t,x]=ode23 (‘danbai’,[0,10],[0.1745,0])%只计算,,100.1745,（弧度）的情形.01g,,()cos,twtw,,%对近似解，周期 T2,,01gw=sqrt (9.8/25);y=0.1745*cos (w*t);[t,x (:,1),y] %显示数据，无分号. hold on %欲在同一幅图中画近似解. plot (t,x (:,1), ‘b’) %画数值解，绿色plot (t,y, ‘g*’) %用*号，红色画近似解. Hold off（2）食饵_捕食者方程,xtrxaxy(),,,ytdybxy(),,，可化为如下形式，，xrayx,0，，，，,,, ,,,,,,0,，dbxy，，，，y，，初始条件xxyy(0),(0),,表示为 00xx(0)，，，，0 ,,,,,yy(0)，，，，0以shier.m存盘如下代码function xdot=shier (t,x)r=1;d=0.5;a=0.1;b=0.02;xdot=diag ([r-a*x (2),-d+b*x (1)])*x;，，xxx(1)，，，，%x=,,,xdot= ,,,,,,,xy(2)，，，，y，，以main3.m存盘如下代码.clear functionsts=0:0.1:15;x0=[25,2];[t,x]=ode45 (‘shier’,ts,x0);[t,x] %显示数据t,x,y plot (t,x)gird %加网格线gtext (‘x(t)’),gtext (‘y(t)’), %用点鼠标方式pause,figure (2) %将x(t),x(t)放至指定点 12plot (x (:,1),x (:,2)) %以x (1)与x (2)为坐标点画相图xlabel (‘x’),ylabel (‘y’)可以猜测，xt()xt()(,)xx与是周期函数，相图是封闭曲线，从数值解可近似定出2121周期为10.7,x2.0,x的最大和最小值分别为与的最大和最小值分别为和，99.328.42.012为求xx与在一个周期的平均值，只需键入： 12y1=x (1:108,1); %x周期为10.7. 1x1p=trapz (y1)*0.1/10.7, %trapz (y1)返回按 y2=x (1:108,2); %梯形法对y1的积107yiyi1()1(1)，，x2p=trapz (y2)*0.1/10.7, %分值,x ,i12i,可得xx,,25,10 12对方程dxxray(),,，,dbxray,dxdy,化为 dyydbx(),，xy积分得解dbxray,,()()xeyec, （4.21）即为原方程组的相轨线，其中c由初始条件确定. 为说明上述相轨线是封闭的，令dbxray,,fxxegyye();(),,设其最大值分别为fg,xy,，若满足 mm00fxffyg(),(),, 00mmdr则有,,xy,fgc,,（令，可解出），又当时，相轨线（4.21）fg,,0,0xy,,,00mm00ba有意义. 当fgc,0,,cfgPxy(,)时，相轨线退化为一个点，对时，相轨线如mmmm00图4.4（c），而图（a），（b）分别为fx()与gy()的图像，下面讨论如何由（a），（b）画（c）.fxf(),nm y gyg(),gv()nm fx() y2 gm fm Pxy(,)qq 00 1 yP q02 q3 y1x y x x xxxy xyxy012001221 （a）（b）（c）图4.4设cpgpf,,,(0)yy,，xx,fxp(),，若令，则有，由（a）知，，使mm012qxy(,)xxx,,qxy(,)fxfxp()(),,，且于是相轨线应过与，对11010222012 fxgxpggyg()()(),,,xxx,(,)fxp(),gyq(),，，由，令，又由（b）知，mm12 存在yyy,,qxy(,)yy,gygyq()(),,qxy(,)使，且，于是相轨线又过与10231121242yyyyy,(),,xqq,两点，所以对间每个，轨线总要通过纵坐标为的两点，于是1210212相轨线是一条封闭曲线.相轨线封闭等价于xtyt(),()是周期函数，记周期为，为求其在一个周期内的平均T值(,)xy，由1yxtd()(),， by两边在一个周期内积分有((0)())yyT,：T11ln()ln(0)yTydTd,，，xxtdt,,，,() ,,,0TTbbb，，同理ry, a从而xxyy,,,00即xy的平均值恰为相轨线中心点的坐标，提请读者注意的是：这里的与与xtyt(),()p00初始条件中的xy,不是一件事. 00注意到在生态学上的意义，上述结果表明，捕食者的数量与食饵的增长率成rdab,,,正比，与它掠取食饵的能力成反比，食饵的数量与捕食者死亡率成正比，与它供养捕食者的能力成反比.3、练习内容（1）编写改进欧拉公式求微分方程数值解的程序，并用其与ode23求下列微分方程数值解，对二者作出比较.a）22,yxyy,,,,(0)0或y(0)1,.222 b),,,xyxyxny，，,,()02,,,21 yxsin,n,,（Bessel方程，这里令，其精确解为）. yy()2,(),,,x222,c),,,yyxyy，,,,cos0,(0)1,(0)0.（2）倒圆锥形容器，上底面直径为1.2m，容器的高亦为1.2m，在锥尖的地方开有一直径为3cm的小孔，容器装满水后，下方小孔开启，由水利学知识可知当水面高度为时，h水从小孔中流出的速度为为重力加速度，若孔口收缩系数为0.6（即若一个面vghg,2,积单位的小孔向外出水时，水柱截面积为0.6），问水从小孔中流完需多少时间？2分钟时，水面高度是多少？（3）一只小船渡过宽为的河流，目标是起点正对着的另一岸上点，已知河水ABd流速vv与船在静水中的速度之比为. k12（a）建立小船航线的方程，求其解析解.（b）设dvv,,,100m,1m/s,2m/s，用数值解法求渡河所需时间，任意时刻小12 船的位置及航行曲线，作图并与解析解比较.（c）若流速v为0,0.5,2 (m/s)，结果将如何？ 1（4）研究种群竞争模型. 当甲、乙两个种群各自生存时，数量演变服从下面规律xyxtrxytry()(1),()(1),,,, 12nn12其中，rr,nn,xtyt(),()分别为时刻甲，乙两个种群的数量，为其固有增长率，为它t1212们的最大容量，而当这两个种群在同一环境中生存时，由于乙消耗有限资源而对甲的增长产生影响，将甲的方程修改为xyxrxs,,,(1) （4.22） 11nn12这里s的含意是：对于供养甲的资源而言，单位数量乙（相对n）的消耗率为单位数12量甲（相对n）消耗的s倍，类似地，甲的存在亦影响乙的增长，乙的方程应改为 11xyytrys()(1),,, （4.23） 22nn12给定种群的初始值为xxyy(0),(0),, （4.24） 00及参数rrssnn,,,,,后，方程（4.22）与（4.23）确定了两种群的变化规律，因其解析121212解不存在，试用数值解法研究以下问题：（a）设rrnnssxy,,,,,,,,1,100,0.5,2,10，计算，画出xtyt(),()12121200 它们的图形及相图(,)xy，说明时间充分大以后xtyt(),()的变化趋势（人们今天看到的t已经是自然界长期演变的结果）.（b）改变rrxy,,,，但s与s不变（保持ss,,1,1），分析所得结果，若12001212ss,,,,1.5(1),0.7(1)，再分析结果，由此你得到什么结论，请用各参数生态学上的12含义作出解释.（c）试验当ss,,,,0.8(1),0.7(1)ss,,,,1.5(1),1.7(1)时会有什么结果；当1212时又会出现什么结果，能解释这些结果吗？。
数学建模实验三 Lorenz模型与食饵模型
数学建模实验三 Lorenz模型与食饵模型一、实验目的1、学习用Mathematica求常微分方程的解析解和数值解，并进行定性分析；2、学习用MATLAB求常微分方程的解析解和数值解，并进行定性分析。
二、实验材料问题图是著名的洛仑兹混沌吸引子，洛仑兹吸引子已成为混沌理论的徽标，好比行星轨道图代表着哥白尼、开普勒理论一样。
洛仑兹是学数学出身的，1948年起在美国麻省理工学院（MIT）作动力气象学博士后工作，1963年他在《大气科学杂志》上发表的论文《确定性非周期流》是混沌研究史上光辉的著作。
以前科学家们不自觉地认为微分方程的解只有那么几类：1)发散轨道；2)不动点；3)极限环；4)极限环面。
除此以外，大概没有新的运动类型了，这是人们的一种主观猜测，谁也没有给出证明。
事实上这种想法是非常错误的。
1963年美国麻省理工学院气象科学家洛仑兹给出一个具体模型，就是著名的Lorenz模型，清楚地展示了一种新型运动体制：混沌运动，轨道既不收敛到极限环上也不跑掉。
而今Lorenz 模型在科学与工程计算中经常运用的问题。
例如，数据加密中。
我们能否绘制出洛仑兹吸引子呢图洛仑兹混沌吸引子假设狐狸和兔子共同生活在同一个有限区域内，有足够多的食物供兔子享用，而狐狸仅以兔子为食物.x为兔子数量,y表狐狸数量。
假定在没有狐狸的情况下，兔子增长率为400％。
如果没有兔子，狐狸将被饿死，死亡率为90％。
狐狸与兔子相互作用的关系是，狐狸的存在使兔子受到威胁，且狐狸越多兔子增长受到阻碍越大，设增长的减小与狐狸总数成正比，比例系数为。
而兔子的存在又为狐狸提供食物，设狐狸在单位时间的死亡率的减少与兔子的数量成正比，设比例系数为。
建立数学模型，并说明这个简单的生态系统是如何变化的。
预备知识1、求解常微分方程的Euler 折线法求初值问题⎩⎨⎧=='00)(),,(y x y y x f y （） 在区间],[0n x x 上的数值解，并在区间插入了结点)()(110n n x x x x <<<<- 。
食饵捕食者模型进一步研究(matlab)
一、食饵-捕食者模型的进一步研究
1）在食饵—捕食者模型（231页1，2式）中研究参数及初始值的变化对食饵和捕食者数量的周期、最大（小）值的影响。
【注：给出不同参数画图即可】
解：取三组不同初值分别为①x0=25,y0=2②x0=25,y0=4③x0=100,y0=2，在matlab中绘图如下
.
分析：
第一组作为对照组，第二组与第一组相比，捕食者初始数量增加，由图象可以看出，捕食者和食饵数量周期均缩短，最大值均变小，最小值均变大；第三组与第一组相比，捕食者和食饵数量周期均增长，最大值均变大，最小值均变小。
2）在上述模型中引入Logistic项（235页16，17式），分析相轨线及参数的影响。
【注：给出不同参数画图即可】
解：取五组不同的参数
①r1=1;r2=0.3;σ1=2; σ2=8;N1=3000;N2=400;（对照）
②r1=2;r2=0.3;σ1=2; σ2=8;N1=3000;N2=400;（食饵增长率变高）
③r1=1;r2=0.6;σ1=2; σ2=8;N1=3000;N2=400;（捕食者增长率高）
④r1=1;r2=0.3;σ1=4; σ2=8;N1=3000;N2=400;（σ1即捕食者掠取食饵能力变大）
⑤r1=1;r2=0.3;σ1=2; σ2=12;N1=3000;N2=400;（σ2即食饵对捕食者的供养能力变大）
在matlab中绘图如下：
从图象可以看出：
⑴改变食饵增长率r1和捕食者增长率r2不会改变最后的稳定点，即最终的稳定点与食饵和捕食者的增长率无关。
⑵第一二三组比较，。
Matlab微分方程的解法
-0.5
-0.55
-0.6
-0.65
-0.7
-0.75
-0.8
-0.85
-0.9
-0.95
-1
0
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
1
time t0=0,tt=1
图3 给定新的初始数据，由函数xprim2定义的ODE解的图形
（d） 求解下面方程组并不难：
x x x x ì ' = - 0.1
在下面的初值问题中，有两个未知函数：x1(t)和x2(t),并用以下式子表达其微... 页码，1/11
Matlab关于微分方程的解法
MATLAB使用龙格-库塔-芬尔格（Runge-Kutta-Fehlberg）方法来解ODE问题。在有限点内计算求解。而 这些点的间距有解的本身来决定。当解比较平滑时，区间内使用的点数少一些，在解变化很快时，区间内应使 用较多的点。 为了得到更多的有关何时使用哪种解法和算法的信息，推荐使用helpdesk。所有求解方程通用的语法或句法在 命令集中头两行给出。时间间隔将以向量t=[t0,tt]给出。 命令ode23可以求解(2,3)阶的常微分方程组，函数ode45使用(4,5)阶的龙格-库塔-芬尔格方法。注意，在这种情 况下x’是x的微分不是x的转置。 在命令集中solver将被诸如ode45函数所取代。
0.6
0.55
0.5
0
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
1
time t0=0,tt=1
图1 由函数xprim1定义的ODE解的图形
（b） 解下面的ODE过程是等价的：
ïíìx' = x2
ïîx(0) = 1
基于MATLAB的捕食模型研究
ｍａ ａ ．Ａｎ ｅ ｓ ｂ ｌ ｙ ａ ｄ ａ ｐ ｉａ ｉｎ ｆｔ ｅ ｓ ｓｅ ｒ ｉｃ ｓｅ ．Ｔ ｅ ｒ ｓ ｌ ｈ ｗ ｔ ａ ， ｆ ｂ ｓ ｆ ａ ｅ ｄｂ ｄ ｔ ｔ ｉ ｔ ｎ ｐ ｌ ｔ ｓ ｏ ｙ ｔｍｓ ａ ｅ ｄ ｓ ｕ ｓ ｄ ｈ ｅ ｕ ｔ ｓ ｏ ｈ ｔ ｍａ ａ ｏ ｔ ｒ ｈ ａ ｉ ｃ ｏ ｈ ｓ ｌ ｗ ｃ ｎ ｗ ｌ ａ ｐ ｉｄ ｔ ｈ ｏ ｅ ｆｐ ｅ ａ ｉｎ ａ ｅｌ ｐ ｌ ｏ ｔ ｅ ｍ ｄ ｌｏ ｒ ｄ ｔ ．Ａｎ ｈ ｏ ｙ ｏ ｓ ｅ ｉｓ ａ ｄ ａ ｒｃ ｌ ｒ ｌｐ ｏ ｕ ｔ ｎ ｐ ｓ ｃ ｎｒ ｌ ｎ ｅ ｏ ｄ ｔ ｅ ｒ ｆｆ ｈ ｒ ｎ ｇ ｉ ｕｔ ａ ｒ ｄ ｃｉ ｅ ｔ ｏ ｔ ｌ ｇ ｉ ｅ ｕ ｏ ｏ ｉ
ｃ ｎ ｂ ｂｔｉ ｅ ａ ｅ ｏ ａ ｎ ｄ，ａ ｄ ｕ ｅ ｏ ｇ ｉ ｅ ｔ ｒ ｄ ｃｉ ｎ ｎ ｓ ｄ ｔ ｕ ｄ ｈｅ ｐ ｏ ｕ ｔ ． ｏ Ｋｅ ｙ ｗｏｒ ｓ： ｒ ｄ ｔ ｎ ｍｏ ｌ ｄ Ｐ ｅ ａ ｉ ｄｅ ；Ｎｏ ｌｎ ａ ｙ ｔｍｓ ｈ ｓ ｉｇ ａ ；Ｄｉ ｅ ｅ ｔａ ｑ ａ ｏ ｓ ｏ ｎ ｉ ｅ ｒｓ ｓｅ ；Ｐ ａ ｅ ｄ ａ ｍ ｒ ｆ ｒ ｎ ｉ ｌｅ ｕ ｔ ｎ ；Ｍａ ｌｂ ｉ ｔ ａ
众 所周 知 ， 捕食 现象 是生态 学研究 的基本 现象 之一 ， 捕食 关 系是 种群 间相互 作 用 的基 本关 系 之一 ， 在
生态学 和数学 上被广 泛地 进行 了研究 ， 中 Ｖ ｈｒ 其 ｏ ｅ ａ模 型是刻 画捕食 者 一 饵 系统 的最 简单 的数 学模 型 。 ｒ 食 在 该模 型 中考 虑系统 自身 阻滞作 用和外 界扰动 等 因素 ， 由此 演 变 出多个 复杂 捕食 系 统 的动 态过 程 和稳 定
MATLAB计算方法与实现
（1）：恢复窗口：在Desktop 中下拉式菜单中的Desktop Layout,选择Default 来恢复。
（2）：在同一坐标系中，画出函数y=x^3-x-1和y=abs(x)*sin5x 的图像。
x=-1:0.1:2;y1=x.^3-x-1; y2=abs(x).*sin(5*x); plot(x,y1,'k',x,y2,':ro')legend('y1=x.^3-x-1','y2=abs(x).*sin(5*x)'),xlabel('x'),ylabel('y'),title('y1,y2画在同一坐标系中')-1-0.500.51 1.52xyy1,y2画在同一坐标系中（3）：根据数据建立一个人口增长模型。
（百万）的函数并绘制出这一函数图形。
根据数学相关理论，用3,4阶多项式拟合这一函数，拟合时不计2000年的数据对，而是将这对数据用来检验并确定模型。
最后用确定的模型预测2010年美国人口。
在Command window 中输入： t=1850:10:1990;p=[23.2,31.4,38.6,50.2,62.9,75.995,91.972,105.711,123.203,131.699,150.697,179.323,203.212,226.505,249.633]; %读取数据plot(t,p,’o ’);axis([1850 2020 0 400]); title(‘Population of the U.s.1850-1990’);ylabel(‘Millions ’);%绘制出数据的函数图形并加以修饰f1=polyfit(t,p,3);f2=polyfit(t,p,4);%对数据做3,4阶多项式拟合，结果分别为f1和f2 v=[polyval(f1,2000),polyval(f2,2000)];%计算当t=2000时多项式f1,f2的值 abs(v-251.422) %计算两个模型与2000年人口数的绝对误差。
Matlab关于微分方程的解法
Matlab关于微分方程的解法MATLAB使用龙格-库塔-芬尔格（Runge-Kutta-Fehlberg）方法来解ODE问题。
在有限点内计算求解。
而这些点的间距有解的本身来决定。
当解比较平滑时，区间内使用的点数少一些，在解变化很快时，区间内应使用较多的点。
为了得到更多的有关何时使用哪种解法和算法的信息，推荐使用helpdesk。
所有求解方程通用的语法或句法在命令集中头两行给出。
时间间隔将以向量t=[t0,tt]给出。
命令ode23可以求解(2,3)阶的常微分方程组，函数ode45使用(4,5)阶的龙格-库塔-芬尔格方法。
注意，在这种情况下x’是x的微分不是x的转置。
在命令集中solver将被诸如ode45函数所取代。
命令集龙格-库塔-芬尔格方法[time,x]=solver(str,t,x0)计算ODE或由字符串str给定的ODE的值，部分解已在向量time中给出。
在向量time中给出部分解，包含的是时间值。
还有部分解在矩阵x中给出，x的列向量是每个方程在这些值下的解。
对于标量问题，方程的解将在向量x中给出。
这些解在时间区间t(1)到t(2)上计算得到。
其初始值是x0即x(t(1)).此方程组有str指定的M文件中函数表示出。
这个函数需要两个参数：标量t和向量x,应该返回向量x’(即x的导数)。
因为对标量ODE来说，x和x’都是标量。
在M文件中输入odefile可得到更多信息。
同时可以用命令numjac来计算Jacobi函数。
[t,x]=solver(str,t,x0,val)此方程的求解过程同上，结构val包含用户给solver的命令。
参见odeset和表1，可得到更多信息。
Ode45此方法被推荐为首选方法。
Ode23这是一个比ode45低阶的方法。
Ode113用于更高阶或大的标量计算。
Ode23t用于解决难度适中的问题。
Ode23s用于解决难度较大的微分方程组。
对于系统中存在常量矩阵的情况也有用。
两种群相互作用规律及其简单应用
两种群相互作用规律及其简单应用摘要：本文主要研究两种生物种群的变化规律及其简单的应用。
我们先建立两种群相互作用的伏特拉模型，分别讨论了相互竞争型、互利互惠型、捕食与被捕食型三种基本情况的相互作用规律，用方程组的形式给出相应的解答，并利用图像对两种群相互作用规律进行直观的分析，接着又在此基础上进行了推广，对三种群的系统进行分析，确定了3种不同状态下的平衡点。
最后，我们利用已经创建的伏特拉模型解释了实际中的问题。
关键词：种群伏特拉模型捕食关系三种群平衡点目录1.问题重述2.基本假设3.符号说明4.问题分析5.模型的建立与求解5.1两种群相互作用规律5.1.1两种群线性作用规律5.1.2应用问题的解释5.2线性作用规律的推广，三种群平衡点分析6.模型评价与体会7.参考文献1.问题重述两种群相互作用存在三种基本类型：相互竞争型、互利互惠型、捕食与被捕食型，试建立合适的数学模型来描述两种生物种群相互作用的规律，并解释为什么在战争期间捕鱼的总数减少而掠肉鱼的比例却会上升这一现象。
2．基本假设1.两种群的自身发展规律和相互制约规律函数都是线性的;2.只考虑到种内自身的发展规律和种间相互作用的影响两个方面;3.种内自身发展规律和种间相互作用均为连续过程，不是在某一时刻突然发生;4.两种生物之间的影响系数，由它们之间的相互关系确定；5.两个种群都受其内部因素制约，只用制约因数来统一描述，制约系数为负，不再使用出生率、死亡率等。
3符号说明由假设可知，两个种群要分别受到种群自身内部和种间相互作用的影响，且自身发展规律和相互制约规律函数都是线性的，故我们可以采用伏特拉模型，通过建立生物种群增长率与自身以及另一个相关种群的数量的关系式来分析问题。
再分别讨论相互竞争型、互利互惠型、捕食与被捕食型三种基本情况下，两个生物种群之间的相互作用关系，求解出两种群最终稳定时的种群数量，从而解决问题。
5.模型的建立与求解5.1两种群相互作用规律5.1.1两种群线性作用规律应用伏特拉模型知两种群相互作用可表示为dN1 / dt = r1 N1（1 - N1 / K1 +αN2 / K1）————（1）其中，N/K可以理解为该种群已经利用的空间（称为“已利用空间项”），（1-N/K）可以理解为尚未利用的空间（称为“未利用空间项”）。
捕食者-被捕食者模型稳定性分析
被捕食者—捕食者模型稳定性分析【摘要】自然界中不同种群之间还存在着一种非常有趣的既有相互依存、又有相互制约的生活方式：种群甲靠丰富的天然资源生存，种群乙靠捕食甲为生，形成食饵-捕食者系统，如食用鱼和鲨鱼，美洲兔和山猫，害虫和益虫等。
本文是基于食饵—捕食者之间的有关规律，建立具有自身阻滞作用的两种群食饵—捕食者模型，分析平衡点的稳定性，进行相轨线分析，并用数值模拟方法验证理论分析的正确性。
【关键词】食饵—捕食者模型相轨线平衡点稳定性一、问题重述在自然界中，存在这种食饵—捕食者关系模型的物种很多。
下面讨论具有自身阻滞作用的两种群食饵-捕食者模型，首先根据该两种群的相互关系建立模型，解释参数的意义，然后进行稳定性分析，解释平衡点稳定的实际意义，对模型进行相轨线分析来验证理论分析的正确性。
二、问题分析本文选择渔场中的食饵(食用鱼)和捕食者(鲨鱼)为研究对象，建立微分方程，并利用数学软件MATLAB 求出微分方程的数值解，通过对数值结果和图形的观察，猜测出它的解析解构造。
然后，从理论上研究其平衡点及相轨线的形状，验证前面的猜测。
三、模型假设1.假设捕食者（鲨鱼）离开食饵无法生存；2.假设大海中资源丰富，食饵独立生存时以指数规律增长；四、符号说明)(t x /)(1t x ——食饵(食用鱼)在时刻t 的数量；)(t y /)(2t x ——捕食者(鲨鱼)在时刻t 的数量；1r ——食饵(食用鱼)的相对增长率；2r ——捕食者(鲨鱼)的相对增长率；1N ——大海中能容纳的食饵(食用鱼)的最大容量；2N ——大海中能容纳的捕食者(鲨鱼)的罪的容量；1σ——单位数量捕食者（相对于2N ）提供的供养食饵的实物量为单位数量捕食者(相对于1N )消耗的供养甲实物量的1σ倍；2σ——单位数量食饵（相对于1N ）提供的供养捕食者的实物量为单位数量捕食者(相对于2N )消耗的供养食饵实物量的2σ倍；d ——捕食者离开食饵独立生存时的死亡率。
捕食者—被捕食者、竞争、共生三种模型的参数估计问题
捕食者—被捕食者、竞争、共生三种模型的参数估计问题 章栋恩为了说明标题中出现的数学模型的重要性，我在这里首先引用MCM 评阅人的一段话。
2009 MCM Judges’ Commentary—Problem BBy Marie Vanisko, Carroll College, Helena, MontanaGeneral Remarks………..The Problem and Selected Modeling Approaches…………Interesting models were constructed for the transitional phase of the cell phone “takeover.” Some teams considered the spread of cell phones as the spread of a disease and used the Verhulst model for logistic growth , using the population of the United States as the carrying capacity and estimating the rate of growth of cell phones from published reports on the growth of cell phone use in the United States. Other teams generalized this to an SIR model or used the Lotka Volterra p redator‐prey model , with cell phones as the predators and landline phones as the prey. A few used the competing species model . The judges looked very favorably upon models for which sufficient rationale was given as to why that model might be appropriate in this circumstance. Interpretation of the parameters and solutions as they applied to the problem at hand was essential.因为捕食者—被捕食者、共生、竞争三种模型都属微分方程组建模问题，在美国微积分教科书和数学建模教材中都有研究。
基于MatLab的三种群Volterra模型数值求解
ｉⅢ８弹．ＫＤＤ ２∞４：５４５－５５０
【卸ｏｈｍＤｎ，Ｎｅｉｌ Ｅ ａｎｄ Ｊａｊｏｄｋ Ｓｕｓｈｌｌ．Ｓｔｅｆ，ａｎｏｇａｐｂｙ：Ｓｅｅｉｎｇ
Ｕｎｍ，ＩＥＥＥ ｔｈｅ
Ｃｏｍｐｕｔｅｒ，ＰＰ．２６—３４，Ｆｃｂｎ“ｙ １９９８
［３］Ｊ＇ｏｈｎｓｅｎ，Ｎｅｉｌ Ｆ．ａｎｄ Ｊａｊｏｄｈ Ｓｕｓｈｉ］．Ｓｔｅｇａｎａｌｙｍ ０ｆ Ｉｍａ弹
Ａｂｅｔｒａｃｔ：ｈｑ”螂－“ｑ掣“ｙｔｋ矗ｌ耐ｔＩｌⅢｈＩ目ｅｄ硼ａ妇ｃｈａｒｅｃｔｅ血ｆｉｃ ０ｆ ｉｍａｇｅ，．ｈ扭ｅ１９０ｒｉｔｈｍ ｉｓ ｂａｓｅｄ佣
ｍｔ瑚ｕｌｔ． ｅｅｌｅｔ－ｄｕｔｒａｃｔｅｍｔｉｅ ０ｆ ２４ ｂｉｔ ＢＭＰ ｉｎｌ噼．ＩｎＩ砧ｅ≮ｆ曲ｔｕｒｅ靶咖ｗＭ“岫砌０１１ ｃｏｌｏｒ ｃｈｌｌ铲ｄ
力。
方程（１）、（３）、（４）构成植物、哺乳动物、爬行动
物三者依存、制约现象的数学模型．即
机 物独立生存时以指数规律增长。相对增长率为ｎ，即
＾
墓
ｘｍ－－－＂ｒｌｘ，，而哺乳动物的存在使植物的增长率减小．设
（５）
ｊＩ旧ｌ（ｎ—＾一ｄ 尊 减小的程度与捕食者数量成正比．于是植物的模型为：
二
＾
＾
七
收罄日期：２００７－０６－１８謦罄日期：２００７－０８—１１
（８）
Ｎ：是哺乳动物的最大容量，仃Ｉ的意义是：单位数
量的哺乳动物（相对Ｎ２而言）掠取毋倍的单位植物
量（相对Ｎ。）。 哺乳动物没有植物的存在会灭亡．设其死亡率
为．则哺乳动物独立存在时有：
五（Ｉ）—ｍ
（９）
植物为哺乳动物提供食物．于是（９）式右端应加 上植物对哺乳动物的促进作用．哺乳动物的增长又受
到自身的阻滞作用．于是有：
系统动力学_“捕食者-被捕食者”Vensim建模
“捕食者-被捕食者”——基于Vensim的模型模拟14307130034光电信息科学与工程毛臻岑-目录-一、模型背景 (3)二、建模过程 (4)2.1新建model (4)2.2创建变量、设置方程 (4)2.3绘制因果图流图 (6)三、运行和调试模型： (6)3.1以原始值运行模拟 (7)3.2调整参数，出现平衡 (14)引用及参考： (21)一、模型背景本文参考捕食者-被捕食者模型（Predator-Prey Model），参考《社会系统动力学》第58页所示流图。
二、建模过程2.1新建model开启Vensim 6.4b PLE新建model，参数设置如下：2.2创建变量、设置方程2.3绘制因果图流图三、运行和调试模型：3.1以原始值运行模拟prey：predator：Time(Month)predator Runs:predator 0Current10113217322428535644755869987 10109 111373.2调整参数，出现平衡设置predator死亡速率为1/16，则两个种群的数量基本上达到平衡：prey：Time(Month)prey Runs:prey 0Current1000 11300 21620 31878 41878 51878 61878 71878 81878 91878 101878 111878predator：引用及参考：1.李旭著：社会系统动力学：政策研究的原理、方法和应用[M].上海：复旦大学出版社，ISBN：978-7-309-06360-82.朱老师的ppt。
地中海鲨鱼问题
地中海鲨鱼问题摘要自然界中存在捕食者与被捕食者系统。
这表示两个物种或多个物种之间的既相互制约又相互依存的自然法则。
本文首先建立了鲨鱼数量在战争前期和战争期间的增长比例模型。
考虑到外界因素不同，模型分为鲨鱼（捕食者）与食用鱼（食饵）在自然条件下战争前、战争期间的比较模型，还有鲨鱼（捕食者）、食用鱼（食饵）在人工捕获条件下战争前、战争期间比较模型。
采用解微分方程组的方法分别求出战争前、战争中两种条件下鲨鱼（捕食者）与食用鱼（鱼饵）的增长比例之间的关系。
由意大利著名生物学家（volterra）建立的系列数学模型。
食用鱼一多，鲨鱼容易得到食物，鲨鱼的数量就会增加。
而由于鲨鱼的数量增加，会吃掉大量的食用鱼，食用鱼的数量减少，鲨鱼缺少食物，鲨鱼的数量就会减少，而食用鱼天敌数量减少，食用鱼的数量就会增加，就这样，食用鱼鲨鱼的数量交替增减，无休止的循环，遂形成生态的动态平衡。
由于战争中人对食用鱼的捕获量减少，鲨鱼食物增多，战争中鲨鱼的比例比战争前要高。
关键词：捕获量百分比，Lotka-Volterra模型，非线性模型。
一、 问题重述意大利生物学家Ancona 曾致力于鱼类种群相互制约关系的研究，他从第一次世界大战期间,地中海各港口捕获的几种鱼类捕获量百分比的资料中，发现鲨鱼等的比例有明显增加（见下表）。
而供其捕食的食用鱼的百分比却明显下降.显然战争使捕鱼量下降，食用鱼增加，鲨鱼等也随之增加，但为何鲨鱼的比例大幅增加呢？他无法解释这个现象，于是求助于著名的意大利数学家V.Volterra ，希望建立一个食饵—捕食系统的数学模型定量地回答这个问题.二、 模型假设食饵由于捕食者的存在使增长率降低，假设降低的程度与捕食者数量成正比； 捕食者由于食饵为它提供食物的作用使其死亡率降低或使之增长，假定增长的程度与食饵数量成正比。
三、 符号说明1()x t —食饵在t 时刻的数量；2()x t —捕食者在t 时刻的数量；1γ—食饵独立生长的增长率；2γ—捕食者独立存在的死亡率；1λ—捕食者掠取食饵的能力；2λ—食饵对捕食者的供养能力；e —捕获能力系数；四、 问题分析战争给人们的生活带来了很多不可磨灭的记忆，同时也会给我们的生活带来不便，战争的到来，征召了许多身强体壮的劳动力去打仗，没有业余的工夫去打鱼，这就造成了鱼类的大量繁殖和生长，同时，捕食者也因为有充足的食饵而得到繁殖和生长。