1、实验二： 微分方程与差分方程模型 Matlab 求解一、实验目的1 掌握解析、数值解法，并学会用图形观察解的形态和进行解的定性分析；2 熟悉 MATLAB 软件关于微分方程求解的各种命令；3 通过范例学习建立微分方程方面的数学模型以及求解全过程；4 熟悉离散 Logistic 模型的求解与混沌的产生过程。 二、实验原理1. 微分方程模型与 MATLAB 求解解析解用 MATLAB 命令 dsolve(eqn1,eqn2, .) 求常微分方程（组）的解析解。其中eqni表示第 i 个微分方程，Dny 表示 y 的 n 阶导数,默认的自变量为 t。（1） 微分方程例 1 求解一阶微分方程 21yd
2、x(1) 求通解输入：dsolve(Dy=1+y2) 输出：ans =tan(t+C1) （2）求特解输入：dsolve(Dy=1+y2,y(0)=1,x) 指定初值为 1，自变量为 x输出：ans =tan(x+1/4*pi) 例 2 求解二阶微分方程 221()04(/)/xyxy原方程两边都除以 ，得2x21()04yyx输入:dsolve(D2y+(1/x)*Dy+(1-1/4/x2)*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、(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 求解 df/dx=3f+4g; dg/dx=-4 f+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)
4、 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,则输出在指
5、定时刻 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)(,1)0()1(2 uuu，解 形式转化：令 。则以上方程转化一阶微分方程组： );21tyt221(;yy。编写 M 文件如下，必须是 M 文件表示微分方程组，并保存，一般地，M 文件的名字与函数名相同，保存位置可以
6、为默认的 work 子目录，也可以保存在自定义文件夹，这时注意要增加搜索路径（FileSet PathAdd 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) 执行：0 2 4 6 8 10 12 14 16 18 20-3-2-1012
7、3 Van Der Pol Solution Time,Secondy(1)andy(2)注：Van der Pol 方程描述具有一个非线性振动项的振动子的运动过程。最初,由于它在非线性电路上的应用而引起广泛兴趣。一般形式为。0)1(2uu图形解无论是解析解还是数值解，都不如图形解直观明了。即使是在得到了解析解或数值解的情况下，作出解的图形，仍然是一件深受欢迎的事。这些都可以用 Matlab 方便地进行。(1)图示解析解如果微分方程（组）的解析解为：y= f (x)，则可以用 Matlab 函数 fplot 作出其图形：fplot(fun,lims)其中：fun 给出函数表达式；lims=xm
8、in xmax ymin ymax限定坐标轴的大小。例如fplot(sin(1/x), 0.01 0.1 -1 1) 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1-1-0.8-0.6-0.4-0.200.20.40.60.81(2)图示数值解设想已经得到微分方程（组）的数值解(x,y)。可以用 Matlab 函数plot(x,y)直接作出图形。其中 x 和 y 为向量（或矩阵）。2、Volterra 模型（食饵捕食者模型）意大利生物学家 Ancona 曾致力于鱼类种群相互制约关系的研究，他从第一次世界大战期间,地中海各港口捕获的几种鱼类捕获量
9、百分比的资料中，发现鲨鱼的比例有明显增加（见下表） 。年代 1914 1915 1916 1917 1918百分比 11.9 21.4 22.1 21.2 36.4年代 1919 1920 1921 1922 1923百分比 27.3 16.0 15.9 14.8 19.7战争为什么使鲨鱼数量增加？是什么原因？因为战争使捕鱼量下降，食用鱼增加，显然鲨鱼也随之增加。 但为何鲨鱼的比例大幅增加呢？生物学家 Ancona 无法解释这个现象，于是求助于著名的意大利数学家 V.Volterra，希望建立一个食饵捕食者系统的数学模型，定量地回答这个问题。1、符号说明： x1(t), x2(t)分别是食饵、
10、捕食者（鲨鱼）在 t 时刻的数量； r 1, r2 是 食饵、捕食者的固有增长率； 1 是捕食者掠取食饵的能力， 2 是食饵对捕食者的供养能力；2、基本假设： 捕食者的存在使食饵的增长率降低，假设降低的程度与捕食者数量成正比，即 )(211xrdtx 食饵对捕食者的数量 x2 起到增长的作用， 其程度与食饵数量 x1 成正比，即 )(122rxdt综合以上和，得到如下模型：模型一：不考虑人工捕获的情况该模型反映了在没有人工捕获的自然环境中食饵与捕食者之间的制约关系，没有考虑食饵和捕食者自身的阻滞作用，是 Volterra 提出的最简单的模型。给定一组具体数据，用 matlab 软件求解。 食饵
11、： r1= 1, 1= 0.1, x10= 25; 捕食者(鲨鱼)：r 2=0.5, 2=0.02, x20= 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 )(12211xrxdt2)0(,25)0()02.5(. 1121 xxxdtt0
12、 5 10 150102030405060708090100图中，蓝色曲线和绿色曲线分别是食饵和鲨鱼数量随时间的变化情况，从图中可以看出它们的数量都呈现出周期性，而且鲨鱼数量的高峰期稍滞后于食饵数量的高峰期。画出相轨迹图：plot(x(:,1),x(:,2) 0 10 20 30 40 50 60 70 80 90 100051015202530模型二 考虑人工捕获的情况假设人工捕获能力系数为 e，相当于食饵的自然增长率由 r1 降为 r1-e，捕食者的死亡率由 r2 增为 r 2+e，因此模型一修改为：设战前捕获能力系数 e=0.3, 战争中降为 e=0.1, 其它参数与模型(一)的参)(1
13、222 2111 xerxdtt 数相同。观察结果会如何变化？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+
14、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); 1 22 112(0.7.)(.8.0)(0)5,)dxxtxdt2)0(,25)0( ).6.)1.09.(21 12 211yyydtytd )()()(;)()()( 212211 tyttptxttp plot(t1,p1,-,t2,p2,*) 0 5 10 1500.10.20.30.4
15、0.50.60.70.8图中*曲线为战争中鲨鱼所占比例。结论：战争中鲨鱼的比例比战前高。 3、 Logistic 映射logistic 映射-通向混沌的道路混沌系统，由于其行为的复杂性，往往认为其动态特性（运动方程）也一定非常复杂，事实并非如此，一个参量很少、动态特性非常简单的系统有时也能够产生混沌现象，以一维虫口模型为例，假设某一区域内的现有虫口数为 yn，昆虫的繁殖率为 r，且第 n 代昆虫不能存活于第 n+1 代，既无世代交叠，则第 n+1 代虫口数为 ，r1 时，虫口会1ny无限制地增长；r1 时，虫口最终会趋于消亡，因此需要对模型进行修正。由于环境的制约和食物有限，因争夺生存空间发生
16、相互咬斗事件的最大次数为 ，即制约虫()2n口数的因素与 成正比，设咬斗事件的战死率为 则对虫口的修正项为 ，则有：2ny 2ny.21nnyry令 ，则 nxyr（1） 1()nnxrx取最大虫口数为 1，且虫口数不能为负，则 ；当 =0.5 时，方程有极大值， 而 又必须小于 1，因而 r4，则参量 r 的取值范围为 1 到 4，这就得14nxr到一个抽象的标准虫口方程（1）。记映射 为f（2）()1)fxr方程（1）可写为（3） 1()nnxf这一迭代关系通常称为 logistic 映射。从0,1内点 x0出发，由 Logistic 映射的迭代形成xn= f n(x0), n = 0,1
17、,2,序列x n称为 x0 的轨道。一个看似简单的系统，随着参量的不同会表现出截然不同的行为，当 r 的取值范围在 13 时，方程（1）有定态解 1xr即方程通过多次迭代后趋于一个稳定的不动点，此时系统是稳定的。为方程 的解，称为周期 2 点。1xr()fx当 r 在 33.448 范围内取值时，经过多次迭代，方程（1）出现周期 2 点 和 ，1x2即 和 是方程 的解，满足1x22()fx12()rx1是使解有意义的 r 最小值。3r随着 r 的增大，r=3.449；3.544；3.564依次出现周期 4、周期 8、周期 16的振荡解，r 的极限值约为 3.569。这种行为称为倍周期分岔，直
18、到 r3.5699 时，系统进入了混沌状态，如下图所示，此时系统的状态不再具有规律性，而是发生随机的波动，使图 d的右侧的大部分区域被涂黑了，仔细观察发现，混沌区域并非一片涂斑，而是有粗粗细细的白色“竖线”，称为周期窗口，随着参量 r 的增大（如 ）时，混沌突然消失，1系统出现周期三的稳定状态， Logistic 映射分岔图接着倍周期分岔以更快的速度进行，再次进入混沌状态。如果将周期窗口放大，发现其结构与分岔图的整体结构具有相似性，而且是一种无限嵌套的自相似结构。 可以看出，通过改变系统参量，使系统进入混沌的第一种模式是倍周期分岔，即由不动点周期二周期四无限倍周期进入混沌状态。当然通向混沌的道
19、路不只于此，第二种通向的道路是：从平衡态到周期运动，再到拟周期运动，直到进入混沌状态。第三种通向混沌的方式是阵发（或间歇）道路，即系统在近似周期运动的过程中，通过改变参量，系统会出现阵发性混沌过程，随着参量的调整，阵发性混沌越来越频繁，近似的周期运动越来越少，最后进入混沌。Matlab 程序下面程序绘制 r=2，x 0=0.3 的轨道clear all;clf;x=0.3;r=2;n=input(Please input a number: ); for i=1:nx1=r*x*(1-x);x=x1; plot(i,x1,.);hold on;end下面程序绘制 logistic 映射 r=3
20、,5 的轨道，观察是否有周期点clear all;clf;x=0.3;r=3.5;for i=1:100x1=r*x*(1-x); 0 10 20 30 40 50 60 70 80 90 1000.40.50.60.70.80.910 10 20 30 40 50 60 70 80 90 1000.420.430.440.450.460.470.480.490.50.51x=x1;plot(i,x1,.);hold on;end下面程序绘制 logistic 映射 r=4 的轨道，观察其混沌clear all;clf;x=0.007;for i=1:100x1=4*x*(1-x); x=x1
21、;plot(i,x1,.);hold on;end下面程序绘制在混沌状态下不同初值 x0=0.100001 和 x0=0.1 的轨道的差（对初值的敏感性）clear all;clf;x1=0.100001;x11=0.1;for i=1:1000x2=4*x1*(1-x1); x1=x2;x22=4*x11*(1-x11);x11=x22;plot(i,x11-x1);hold on;end下面程序绘制 logistic 映射的分岔图clear all;clf;x=0.2;for r=1:0.001:4for i=1:18x1=r*x*(1-x); x=x1;if i9plot(r,x); h
22、old on;endendend0 10 20 30 40 50 60 70 80 90 10000.10.20.30.40.50.60.70.80.910 100 200 300 400 500 600 700 800 900 1000-1-0.8-0.6-0.4-0.200.20.40.60.811 1.5 2 2.5 3 3.5 400.10.20.30.40.50.60.70.80.91蛛网迭代clear all;clf;x=0.007;y=0;for i=1:200 XY=x,y;y=4*x*(1-x);XY1=x,y;line(XY,XY1);hold on;x=y;XY2=x,x
23、;line(XY1,XY2);hold on;end0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 100.10.20.30.40.50.60.70.80.91三、实验内容1求微分方程的解析解, 并画出图形，=+2,(0)1,yxx 2求微分方程的数值解, 并画出图形，cos,(),(0)yxy 3两种相似的群体之间为了争夺有限的同一种食物来源和生活空间而进行生存竞争时，往往是竞争力较弱的种群灭亡，而竞争力较强的种群达到环境容许的最大数量。假设有甲、乙两个生物种群，当它们各自生存于一个自然环境中，均服从 Logistic 规律。（1） )(,21tx是两个种群的数量；（2） ,r是它们的固有增长率；（3） 21,n是它们的最大容量；（4） )(m为种群乙(甲)占据甲(乙)的位置的数量，并且 .12;x计算 )(,21tx, 画出图形及相轨迹图。解释其解变化过程。2） ,=1.5， =0.7，计算 )(,21tx, ,r设 120,n120x画出图形及相轨迹图。解释其解变化过程。4. 绘制 Logistic 映射的轨道图,分岔图和蛛网迭代图。四、实验心得
