数学建模实验二：微分方程模型Matlab求解与分析
用matlab对微分方程求解实验报告.
o 《高等数学》上机作业(三一、上机目的1、学会用 M a t l a b 求简单微分方程的解析解。
2、学会用 M a t l a b 求微分方程的数值解。
二、上机内容1、求简单微分方程的解析解.2、求微分方程的数值解.3、数学建模实例.4、上机作业. 三、上机作业1. 求微分方程:在初值条件下的特解,并画出解函数的图形. 命令>> y =d s o l v e ('x *D y +y -e x p (x =0','y (1=2*e x p (1','x ' 运行结果:y = 1/x *e x p (x +1/x *e x p (1'xxy y e +-=12(y e =函数图象:2. 求微分方程的特解.22450(00,'(110d y dyy dx dx y y ?+-=???==?命令>> y=dsolve('D2y+4*Dy-5*y=0','y(0=0,Dy(1=10','x' 运行结果:y=10/(exp(1+5*exp(-5*exp(x-10/(exp(1+5*exp(-5*exp(-5*x3. 鱼雷追击问题一敌舰在某海域内沿着正北方向航行时,我方战舰恰好位于敌舰的正西方向 1 公里处.我舰向敌舰发射制导鱼雷,敌舰速度为0.42 公里/分,鱼雷速度为敌舰速度的2倍。
试问敌舰航行多远时将被击中?M文件x0=0; xf=0.9999999999999; [x,y]=ode15s('eq1',[x0 xf],[0 0]; plot(x,y(:,1,'b.'hold on;y=0:0.1:1;plot(1,y, '*'运行结果图像:结论:大概在y=0.67处击中敌方舰艇!(选做一个慢跑者在平面上沿椭圆以恒定的速率v=1跑步,设椭圆方程为:x=10+20cost, y=20+5sint. 突然有一只狗攻击他. 这只狗从原点出发,以恒定速率w跑向慢跑者,狗的运动方向始终指向慢跑者.分别求出w=20,w=5时狗的运动轨迹.W=20M文件代码function dy=eq3(t,ydy=zeros(2,1;dy(1=20*(10+20*cos(t-y(1/sqrt((10+20*cos(t-y(1^2+(2 0+15*sin(t-y(2^2;dy(2=20*(20+15*sin(t-y(2/sqrt((10+20*cos(t-y(1^2+(2 0+15*sin(t-y(2^2;运行命令t0=0;tf=10;[t,y]=ode45('eq3',[t0 tf],[0 0];T=0:0.1:2*pi;X=10+20*cos(T;Y=20+15*sin(T;plot(X,Y,'-'hold onplot(y(:,1,y(:,2,'r*'运行结果:利用二分法更改tf tf=5时tf=2.5时tf=3.15时:所以在t=3.15时刻恰好追上!W=5M文件代码function dy=eq4(t,ydy=zeros(2,1;dy(1=5*(10+20*cos(t-y(1/sqrt((10+20*cos(t-y(1^2+(20 +15*sin(t-y(2^2; dy(2=5*(20+15*sin(t-y(2/sqrt((10+20*cos(t-y(1^2+(20 +15*sin(t-y(2^2; 命令:t0=0;tf=10;[t,y]=ode45('eq4',[t0 tf],[0 0]; T=0:0.1:2*pi;X=10+20*cos(T; Y=20+15*sin(T; plot(X,Y,'-'hold onplot(y(:,1,y(:,2,'*' 运行结果更改tf=20运行结果Tf=40 11所以永远追不上！四、上机心得体会高等数学是工科学生的主干科目，它应用于生产生活的方方面面，通过建模，计算可以求出实际问题的最优化问题！因此我们需要掌握建模和利用专业软件处理实际问题的能力！ 12。
微分方程模型matlab
微分方程模型matlab微分方程是自然科学及工程领域中常用的数学工具，它描述了系统中的变化率与其当前状态之间的关系。
而Matlab是一个广泛使用的科学计算软件，它具有强大的数值计算能力和可视化功能。
在本文中，我们将讨论如何使用Matlab来求解微分方程模型，并给出一个实例。
一、Matlab中的微分方程求解函数Matlab中有几个可用于求解微分方程的函数，包括ode45、ode23、ode113、ode15s等。
其中前两者是比较常用的。
这些函数都采用了一些数值方法来近似求解微分方程，得到数值解。
这些求解函数都采用输入函数来表示微分方程。
二、小例子：一个单自由度振动模型作为一个简单的微分方程模型，我们来考虑一个单自由度振动模型。
它可以用下面的微分方程来表示：$$m\ddot{x}+c\dot{x}+kx=F(t)$$其中，$m$是质量，$c$是阻尼系数，$k$是弹簧刚度，$F(t)$是作用在系统上的外力，$x(t)$是位移，$\dot{x}(t)$是速度，$\ddot{x}(t)$是加速度。
我们假设质量$m$为1kg，阻尼系数$c$为0.2N-s/m，弹簧刚度$k$为100N/m，外力$F(t)$为3N。
首先，我们将这些参数定义为变量。
```matlabm = 1;c = 0.2;k = 100;F = 3;```然后，我们将微分方程表示为Matlab可以识别的形式。
我们定义一个名为odefun的函数，在此函数中计算微分方程的右侧。
```matlabfunction dxdt = odefun(t, x)dxdt = [x(2); (F - c*x(2) - k*x(1))/m];end```接着，我们要设置初始状态。
我们初始化$x_0$为0，$\dot{x}_0$为0.5m/s。
```matlabx0 = [0; 0.5];```我们使用ode45函数来解微分方程。
ode45函数需要我们提供微分方程函数和初始条件，然后它会自动计算微分方程的解。
使用Matlab进行微分方程求解的方法
使用Matlab进行微分方程求解的方法引言微分方程是数学中非常重要的一部分，广泛应用于物理、经济、工程等领域。
对于大部分微分方程的解析解往往难以求得，而数值解法则成为了一种常用的解决手段。
Matlab作为一种强大的科学计算软件，也提供了丰富的工具和函数用于求解微分方程，本文将介绍一些常见的使用Matlab进行微分方程求解的方法。
一、数值求解方法1. 欧拉方法欧拉方法是最简单的一种数值求解微分方程的方法，它将微分方程的微分项用差分的方式进行近似。
具体的公式为：y(n+1) = y(n) + hf(x(n), y(n))其中，y(n)表示近似解在第n个点的值，h为步长，f(x, y)为微分方程的右端项。
在Matlab中使用欧拉方法进行求解可以使用ode113函数，通过设定不同的步长，可以得到不同精度的数值解。
2. 中点法中点法是较为精确的一种数值求解微分方程的方法，它的计算公式为：k1 = hf(x(n), y(n))k2 = hf(x(n) + h/2, y(n) + k1/2)y(n+1) = y(n) + k2中点法通过计算两个斜率的平均值来得到下一个点的值，相较于欧拉方法，中点法能提供更精确的数值解。
3. 4阶龙格库塔法龙格库塔法是一类高阶数值求解微分方程的方法，其中4阶龙格库塔法是最常用的一种。
它的计算公式为：k1 = hf(x(n), y(n))k2 = hf(x(n) + h/2, y(n) + k1/2)k3 = hf(x(n) + h/2, y(n) + k2/2)k4 = hf(x(n) + h, y(n) + k3)y(n+1) = y(n) + (k1 + 2k2 + 2k3 + k4)/64阶龙格库塔法通过计算多个斜率的加权平均值来得到下一个点的值，相较于欧拉方法和中点法，它的精度更高。
二、Matlab函数和工具除了可以使用以上的数值方法进行微分方程求解之外，Matlab还提供了一些相关的函数和工具，方便用户进行微分方程的建模和求解。
实验二_基于Matlab的微分方程数值解法
实验二微分方程数值解法一．实验原理及实验内容：对微分方程描述的控制系统，利用欧拉法、二阶龙格-库塔法、四阶龙格-库塔法分别编写M文件，进行数值计算和作图。
1．分别用欧拉法、二阶龙格-库塔法、四阶龙格-库塔法求下面系统的输出响应y(t)在0≤t≤1上，h=0.1时的数值解。
'2,(0)1=-=y y y要求保留4位小数，并将三种方法的结果与真解2=进行比较。
()ty t e-2.若为如2y y y==何编程计算？',(0)1二．实验仪器：计算机Matlab软件三．实验数据记录：程序一：disp('欧拉算法');y=1;h=0.1;for i=0:0.1:1disp(y);y=y+h*(-2*y);enddisp('欧拉算法');ydisp('精确解');yy=exp(-2*t)h=0.1;disp('函数的2阶数值解为');disp('y=');y=1;for t=0:h:1;disp(y);k1=-2*y;k2=-2*(y+k1*h);y(i+1)=y(i)+(k1+k2)*h*1/2;endh=0.1;disp('函数的4阶数值解为');disp('y=');y=1;for t=0:h:1;disp(y);k1=-2*y;k2=-2*(y+k1*h*1/2);k3=-2*(y+k2*h*1/2);k4=-2*(y+k3*h);y=y+h*1/6*(k1+2*k2+2*k3+k4); end>>程序2：t=0:0.1:1;n=length(t);y(1)=1;h=0.1;for i=1:n-1y(i+1)=y(i)+h*(y(i)*y(i)); enddisp('欧拉算法');ydisp('精确解');yy=exp(-2*t)h=0.1;disp('函数的2阶数值解为');disp('y=');y=1;for t=0:h:1;disp(y);k1=y*y;k2=(y+k1*h)^2;y=y+(k1+k2)*h*1/2;endh=0.1;disp('函数的4阶数值解为');disp('y=');y=1;disp(y);k1=y*y;k2=(y+k1*h*1/2)^2;k3=(y+k2*h*1/2)^2;k4=(y+k3*h)^2;y=y+h*1/6*(k1+2*k2+2*k3+k4); end。
数学建模实验二：微分方程与差分方程模型Matlab求解
实验二：微分方程与差分方程模型Matlab 求解专业年级： 2014级信息与计算科学1班 姓名： 黄志锐 学号：201430120110一、实验目的[1] 掌握解析、数值解法，并学会用图形观察解的形态和进行解的定性分析； [2] 熟悉MATLAB 软件关于微分方程求解的各种命令；[3] 通过范例学习建立微分方程方面的数学模型以及求解全过程； [4] 熟悉离散 Logistic 模型的求解与混沌的产生过程。
二、实验内容1．求微分方程的解析解, 并画出图形，解：使用MATLAB 编程计算得出上述微分方程的解析解为：y =3e x −2x −2其解析解图示如图1所示：图1 解析解图示=+2,(0)1,01y y x y x '=<<MATLAB代码运行结果截图如下所示：2．求微分方程的数值解, 并画出图形，解：使用MATLAB编程计算上述微分方程的数值解，并作出其数值解的图示如图2所示：图2 数值解图示cos0,(0)1,(0)0y y x y y'''+===3．两种相似的群体之间为了争夺有限的同一种食物来源和生活空间而进行生存竞争时，往往是竞争力较弱的种群灭亡，而竞争力较强的种群达到环境容许的最大数量。
假设有甲、乙两个生物种群，当它们各自生存于一个自然环境中，均服从 Logistic 规律。
（1）是两个种群的数量； （2）是它们的固有增长率； （3）是它们的最大容量；（4）为种群乙(甲)占据甲(乙)的位置的数量，并且建立模型如下：)(),(21t x t x 21,r r 21,n n )(12m m .1122;x m x m βα==计算, 画出图形及相轨线图。
解释其解变化过程。
2）121,r r ==设12100,n n ==102010x x ==,=1.5，=0.7，计算, 画出图形及相轨线图。
解释其解变化过程。
解：（1）使用MATLAB 编程实现上述模型，并输入相关参数后，计算得出, 并画出图形及相轨线图如下所示：图3 数值解图示)(),(21t x t x αβ)(),(21t x t x )(),(21t x tx图4 相轨线图详细MATLAB代码如下：由图3、图4可以看出，甲、乙两个生物种群以几乎一致的趋势不断增长，直到达到一个相对稳定的数量。
matlab微分方程模型
matlab微分方程模型Matlab微分方程模型是一种基于Matlab软件的数学建模方法，用于解决微分方程相关的问题。
微分方程是描述物理、工程和数学问题的重要工具，通过建立微分方程模型，可以对各种现象进行定量分析和预测。
在Matlab中，可以使用ode45函数求解常微分方程（ODE）或者ode15s函数求解刚性ODE。
这些函数可以通过数值方法近似求解微分方程的解析解，从而得到问题的数值解。
具体来说，可以通过在Matlab中定义微分方程的右侧函数，然后使用相应的ode函数进行求解。
例如，考虑一个简单的一阶线性微分方程模型：dy/dx = -ky，其中k为常数。
我们可以通过在Matlab中定义这个微分方程的右侧函数，并使用ode45函数求解。
具体步骤如下：1. 在Matlab中定义微分方程的右侧函数：function dydx = myODE(x,y)k = 0.1; % 设定常数k的值dydx = -k*y;end2. 使用ode45函数求解微分方程：xspan = [0 10]; % 设定求解区间y0 = 1; % 设定初始条件[x,y] = ode45(@myODE, xspan, y0);3. 绘制得到的数值解：plot(x,y);xlabel('x');ylabel('y');title('Solution of dy/dx = -ky');通过以上步骤，我们可以得到微分方程dy/dx = -ky的数值解，并绘制出解的图像。
这个简单的例子展示了如何使用Matlab微分方程模型求解微分方程。
除了一阶线性微分方程，Matlab微分方程模型还可以用于解决更复杂的微分方程问题，包括高阶线性微分方程、非线性微分方程、偏微分方程等。
通过定义相应的微分方程函数和合适的求解方法，可以在Matlab中进行数值求解。
此外，Matlab还提供了丰富的绘图和分析工具，可以对微分方程的解进行可视化和进一步分析。
实验二MATLAB数值计算常微分方程（组）的求解
实验⼆MATLAB数值计算常微分⽅程（组）的求解实验⼆ MATLAB 数值计算：常微分⽅程(组)的求解⼀、实验⽬的在物理学和⼯程技术上，很多问题都可以⽤⼀个或⼀组常微分⽅程来描述，因此要解决相应的实际问题往往需要⾸先求解对应的微分⽅程。
在⼤多数情况下这些微分⽅程通常是⾮线性的或者是超越⽅程(⽐如范德堡⽅程，波导本征值⽅程等)，因此往往需要使⽤计算机数值求解。
MATLAB 作为⼀种强⼤的科学计算语⾔，其在数值计算和数据的可视化⽅⾯具有⽆以伦⽐的优势。
在解决常微分⽅程问题上，MATLAB 就提供了多种可适⽤于不同场合(如刚性和⾮刚性问题)下的求解器(Solver)，例如ode45，ode15s ，ode23，ode23s 等等。
本次实验将以范德堡⽅程的计算和地球卫星的运⾏轨道的仿真为例，练习使⽤MATLAB 的常微分⽅程求解器，以期达到如下⼏个⽬的：1. 熟悉常微分⽅程的求解⽅法，了解状态⽅程的概念；2. 能熟练使⽤dsolve 函数解析求解常微分⽅程；3. 能熟练运⽤ode45、ode15s 求解器分别数值求解⾮刚性和刚性常微分⽅程；4. 学习⽤求解器来绘制相图的⽅法。
⼆、实验的预备知识1．微分⽅程的概念未知的函数以及它的某些阶的导数连同⾃变量都由⼀已知⽅程联系在⼀起的⽅程称为微分⽅程。
如果未知函数是⼀元函数，称为常微分⽅程（Ordinary differential equations ，简称odes ）。
n 阶常微分⽅程的⼀般形式（隐式）为：0),,",',,()(=n y y y y t F （1）其中t 为⾃变量。
如果未知函数是多元函数，成为偏微分⽅程。
联系⼀些未知函数的⼀组微分⽅程组称为微分⽅程组。
微分⽅程中出现的未知函数的导数的最⾼阶解数称为微分⽅程的阶。
若⽅程中未知函数及其各阶导数都是⼀次的，称为线性常微分⽅程，⼀般表⽰为)()(')()(1)1(1)(t b y t a y t a y t a y n n n n =++++--若上式中的系数a i (t), i =1,2,…,n 均与t ⽆关，称之为常系数。
用MATLAB求解微分方程及微分方程组
dx dt X x ， 0 即： Y y dy dt w dx ( X x) dt 2 2 ( X x ) (Y y ) 消去λ 得： dy w (Y y ) 2 2 dt ( X x ) (Y y )
1 0.8 0.6 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -1 0 2 4 6 8 10 12
解
1、建立m-文件rigid.m如下： function dy=rigid(t,y) dy=zeros(3,1); dy(1)=y(2)*y(3); dy(2)=-y(1)*y(3); dy(3)=-0.51*y(1)*y(2);
用MATLAB求解微分方程
1. 微分方程的解析解
求微分方程（组）的解析解命令: dsolve(‘方程1’, ‘方程2’,…‘方程n’, ‘初始条件’, ‘自变量’)
记号: 在表达微分方程时，用字母 D 表示求微分，D2、D3 等 表示求高阶微分.任何 D 后所跟的字母为因变量，自变量可以指 定或由系统规则选定为确省. 例如，微分方程
任取k1、k2的一组初始值：k0=[2,1];
输入命令： k=lsqcurvefit('curvefun1',k0,t,c) 运行结果为： k =[ 1.3240 作图表示求解结果： t1=0:0.1:18; f=curvefun1(k,t1); plot(t,c,'ko',t1,f,'r-')
90 80 70 60 50 40 30 20 10 0
（2）
（3）
3．因乙舰以速度v0沿直线x=1运动，设v0=1，则w=5，X=1，Y=t
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
实验二： 微分方程模型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 。
（1）通解:[f,g]=dsolve('Df=3*f+4*g','Dg=-4*f+3*g')f =exp(3*t)*(C1*sin(4*t)+C2*cos(4*t)) g =exp(3*t)*(C1*cos(4*t)-C2*sin(4*t))特解：[f,g]=dsolve('Df=3*f+4*g','Dg=-4*f+3*g','f(0)=0,g(0)=1')f =exp(3*t)*sin(4*t) g =exp(3*t)*cos(4*t)数值解在微分方程(组)难以获得解析解的情况下，可以用Matlab 方便地求出数值解。
格式为:[t,y] = ode23('F',ts,y0,options)注意：➢ 微分方程的形式：y ' = F (t , y )，t 为自变量，y 为因变量（可以是多个，如微分方程组）；➢ [t, y]为输出矩阵，分别表示自变量和因变量的取值；➢ F 代表一阶微分方程组的函数名（m 文件，必须返回一个列向量，每个元素对应每个方程的右端）；➢ ts 的取法有几种，（1）ts=[t0, tf] 表示自变量的取值范围，（2）ts=[t0,t1,t2,…,tf],则输出在指定时刻t0,t1,t2,…,tf 处给出，（3）ts=t0:k:tf,则输出在区间[t0,tf]的等分点给出； ➢ y0为初值条件；➢ options 用于设定误差限（缺省是设定相对误差是10^(-3)，绝对误差是10^(-6)）； ode23是微分方程组数值解的低阶方法，ode45为中阶方法，与ode23类似。
例4 求解一个经典的范得波（Van Der pol ）微分方程：0)0(',1)0(0')1(''2===+-+u u u u u u ，解 形式转化：令)(');(21t u y t u y ==。
则以上方程转化一阶微分方程组：1221221)1(;y y y y y y --='='。
编写M 文件如下，必须是M 文件表示微分方程组，并保存，一般地，M 文件的名字与函数名相同，保存位置可以为默认的work 子目录，也可以保存在自定义文件夹，这时注意要增加搜索路径（File\Set Path\Add Folder ）function dot1=vdpol(t,y); dot1=[y(2); (1-y(1)^2)*y(2)-y(1)]; 在命令窗口写如下命令:[t,y]=ode23('vdpol',[0,20],[1,0]); y1=y(:,1);y2=y(:,2);plot(t,y1,t,y2,'--');title('Van Der Pol Solution ');xlabel('Time,Second');ylabel('y(1)andy(2)')注：Van der Pol方程描述具有一个非线性振动项的振动子的运动过程。
最初,由于它在非线性电路上的应用而引起广泛兴趣。
一般形式为-+u+uu。
u''2=')1(图形解无论是解析解还是数值解，都不如图形解直观明了。
即使是在得到了解析解或数值解的情况下，作出解的图形，仍然是一件深受欢迎的事。
这些都可以用Matlab方便地进行。
(1)图示解析解如果微分方程（组）的解析解为：y=f (x)，则可以用Matlab函数fplot作出其图形：fplot('fun',lims)其中：fun给出函数表达式；lims=[xmin xmax ymin ymax]限定坐标轴的大小。
例如fplot('sin(1/x)', [0.01 0.1 -11])(2)图示数值解设想已经得到微分方程（组）的数值解(x,y)。
可以用Matlab函数plot(x,y)直接作出图形。
其中x和y为向量（或矩阵）。
2、Volterra模型（食饵捕食者模型）意大利生物学家Ancona曾致力于鱼类种群相互制约关系的研究，他从第一次世界大战期间,地中海各港口捕获的几种鱼类捕获量百分比的资料中，发现鲨鱼的比例有明显增加（见下表）。
战争为什么使鲨鱼数量增加？是什么原因？因为战争使捕鱼量下降，食用鱼增加，显然鲨鱼也随之增加。
但为何鲨鱼的比例大幅增加呢？生物学家Ancona无法解释这个现象，于是求助于著名的意大利数学家V.Volterra，希望建立一个食饵—捕食者系统的数学模型，定量地回答这个问题。
1、符号说明：①x1(t), x2(t)分别是食饵、捕食者（鲨鱼）在t时刻的数量；②r1, r2是食饵、捕食者的固有增长率；③λ1是捕食者掠取食饵的能力，λ2是食饵对捕食者的供养能力；2、基本假设：① 捕食者的存在使食饵的增长率降低，假设降低的程度与捕食者数量成正比，即)(21111x r x dtdx λ-=② 食饵对捕食者的数量x 2起到增长的作用， 其程度与食饵数量x 1成正比，即)(12222x r x dtdx λ+-=综合以上①和②，得到如下模型：模型一：不考虑人工捕获的情况该模型反映了在没有人工捕获的自然环境中食饵与捕食者之间的制约关系，没有考虑食饵和捕食者自身的阻滞作用，是Volterra 提出的最简单的模型。
给定一组具体数据，用matlab 软件求解。
食饵： r 1= 1, λ1= 0.1, x 10= 25; 捕食者(鲨鱼)：r 2=0.5, λ2=0.02, x 20= 2;编制程序如下1、建立m-文件shier.m 如下：function dx=shier(t,x) dx=zeros(2,1); %初始化dx(1)=x(1)*(1-0.1*x(2));dx(2)=x(2)*(-0.5+0.02*x(1)); 2、在命令窗口执行如下程序：[t,x]=ode45('shier',0:0.1:15,[25 2]); plot(t,x(:,1),'-',t,x(:,2),'*'),grid⎪⎪⎩⎪⎪⎨⎧+-=-=)()(1222221111x r x dtdx x r x dt dx λλ2)0(,25)0()02.05.0()1.01(21122211==⎪⎪⎩⎪⎪⎨⎧+-=-=x x x x dtdx x x dt dx从图中可以看出它们的数量都呈现出周期性，而且鲨鱼数量的高峰期稍滞后于食饵数量的高峰期。
画出相轨迹图：模型二 考虑人工捕获的情况假设人工捕获能力系数为e ，相当于食饵的自然增长率由r 1 降为r 1-e ，捕食者的死亡率由r 2 增为 r 2+e ，因此模型一修改为：设战前捕获能力系数e=0.3, 战争中降为e=0.1, 其它参数与模型(一)的参⎪⎩⎪⎨⎧++-=--=])([])[(1222221111x e r x dtdxx e r x dt dx λλ数相同。
观察结果会如何变化？1）当e=0.3时： 2）当e=0.1时：分别求出两种情况下鲨鱼在鱼类中所占的比例。
即计算画曲线：plot(t,p1(t),t,p2(t),'*') MATLAB 编程实现 建立两个M 文件function dx=shier1(t,x) dx=zeros(2,1);dx(1)=x(1)*(0.7-0.1*x(2)); dx(2)=x(2)*(-0.8+0.02*x(1)); function dy=shier2(t,y) dy=zeros(2,1);dy(1)=y(1)*(0.9-0.1*y(2)); dy(2)=y(2)*(-0.6+0.02*y(1));运行以下程序：[t1,x]=ode45('shier1',[0 15],[25 2]); [t2,y]=ode45('shier2',[0 15],[25 2]); x1=x(:,1);x2=x(:,2); p1=x2./(x1+x2);y1=y(:,1);y2=y(:,2); p2=y2./(y1+y2);11222112(0.70.1)(0.80.02)(0)25,(0)2dx x x dt dx x x dtx x ⎧=-⎪⎪⎪=-+⎨⎪==⎪⎪⎩⎪⎪⎪⎩⎪⎪⎪⎨⎧==+-=-=2)0(,25)0()02.06.0()1.09.0(21122211y y y y dtdy y y dt dy )()()()(;)()()()(21222121t y t y t y t p t x t x t x t p +=+=plot(t1,p1,'-',t2,p2,'*')0510150.10.20.30.40.50.60.70.8图中‘*’曲线为战争中鲨鱼所占比例。