实验二微分方程与差分方程模型Matlab求解.docx
《实验二微分方程与差分方程模型Matlab求解.docx》由会员分享，可在线阅读，更多相关《实验二微分方程与差分方程模型Matlab求解.docx（20页珍藏版）》请在冰点文库上搜索。
实验二微分方程与差分方程模型Matlab求解
实验二：
微分方程与差分方程模型Matlab求解
一、实验目的
[1]掌握解析、数值解法，并学会用图形观察解的形态和进行解的定性分析；
[2]熟悉MATLAB软件关于微分方程求解的各种命令；
[3]通过范例学习建立微分方程方面的数学模型以及求解全过程；
[4]熟悉离散Logistic模型的求解与混沌的产生过程。
二、实验原理
1.微分方程模型与MATLAB求解
解析解
用MATLAB命令dsolve（‘eqn1’,’eqn2’,...）求常微分方程（组）的解析解。
其中‘eqni'表示第i个微分方程，Dny表示y的n阶导数,默认的自变量为t。
（1）微分方程
例1求解一阶微分方程
（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求解二阶微分方程
原方程两边都除以
，得
输入:
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求解 df/dx=3f+4g; dg/dx=-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求解一个经典的范得波（VanDerpol）微分方程：
解形式转化：
令
。
则以上方程转化一阶微分方程组：
。
编写M文件如下，必须是M文件表示微分方程组，并保存，一般地，M文件的名字与函数名相同，保存位置可以为默认的work子目录，也可以保存在自定义文件夹，这时注意要增加搜索路径（File\SetPath\AddFolder）
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（'VanDerPolSolution'）;
xlabel（'Time,Second'）;ylabel（'y
（1）andy
（2）'）
执行：
注：
VanderPol方程描述具有一个非线性振动项的振动子的运动过程。
最初,由于它在非线性电路上的应用而引起广泛兴趣。
一般形式为
。
图形解
无论是解析解还是数值解，都不如图形解直观明了。
即使是在得到了解析解或数值解的情况下，作出解的图形，仍然是一件深受欢迎的事。
这些都可以用Matlab方便地进行。
（1）图示解析解
如果微分方程（组）的解析解为：
y=f（x），则可以用Matlab函数fplot作出其图形：
fplot（'fun',lims）
其中：
fun给出函数表达式；lims=[xminxmaxyminymax]限定坐标轴的大小。
例如
fplot（'sin（1/x）',[-11]）
（2）图示数值解
设想已经得到微分方程（组）的数值解（x,y）。
可以用Matlab函数plot（x,y）直接作出图形。
其中x和y为向量（或矩阵）。
2、Volterra模型（食饵捕食者模型）
意大利生物学家Ancona曾致力于鱼类种群相互制约关系的研究，他从第一次世界大战期间,地中海各港口捕获的几种鱼类捕获量百分比的资料中，发现鲨鱼的比例有明显增加（见下表）。
年代
1914
1915
1916
1917
1918
百分比
年代
1919
1920
1921
1922
1923
百分比
战争为什么使鲨鱼数量增加？
是什么原因？
因为战争使捕鱼量下降，食用鱼增加，显然鲨鱼也随之增加。
但为何鲨鱼的比例大幅增加呢？
生物学家Ancona无法解释这个现象，于是求助于著名的意大利数学家，希望建立一个食饵—捕食者系统的数学模型，定量地回答这个问题。
1、符号说明：
① x1（t）,x2（t）分别是食饵、捕食者（鲨鱼）在t时刻的数量；
② r1,r2是食饵、捕食者的固有增长率；
③λ1是捕食者掠取食饵的能力，
λ2是食饵对捕食者的供养能力；
2、基本假设：
①捕食者的存在使食饵的增长率降低，假设降低的程度与捕食者数量成正比，即
② 食饵对捕食者的数量x2起到增长的作用，其程度与食饵数量x1成正比，即
综合以上①和②，得到如下模型：
模型一：
不考虑人工捕获的情况
该模型反映了在没有人工捕获的自然环境中食饵与捕食者之间的制约关系，没有考虑食饵和捕食者自身的阻滞作用，是Volterra提出的最简单的模型。
给定一组具体数据，用matlab软件求解。
食饵：
r1=1,λ1=,x10=25;
捕食者（鲨鱼）：
r2=,λ2=,x20=2;
编制程序如下
1、建立m-文件如下：
functiondx=shier（t,x）
dx=zeros（2,1）;%初始化
dx
（1）=x
（1）**x
（2））;
dx
（2）=x
（2）*+*x
（1））;
2、在命令窗口执行如下程序：
[t,x]=ode45（'shier',0:
:
15,[252]）;
plot（t,x（:
1）,'-',t,x（:
2）,'*'）,grid
图中，蓝色曲线和绿色曲线分别是食饵和鲨鱼数量随时间的变化情况，从图中可以看出它们的数量都呈现出周期性，而且鲨鱼数量的高峰期稍滞后于食饵数量的高峰期。
画出相轨迹图：
plot（x（:
1）,x（:
2））
模型二考虑人工捕获的情况
假设人工捕获能力系数为e，相当于食饵的自然增长率由r1降为r1-e，捕食者的死亡率由r2增为r2+e，因此模型一修改为：
设战前捕获能力系数e=,战争中降为e=,其它参数与模型
（一）的参数相同。
观察结果会如何变化？
1）当e=时：
2）当e=时：
分别求出两种情况下鲨鱼在鱼类中所占的比例。
即计算
画曲线：
plot（t,p1（t）,t,p2（t）,'*'）
MATLAB编程实现
建立两个M文件
functiondx=shier1（t,x）
dx=zeros（2,1）;
dx
（1）=x
（1）*
dx
（2）=x
（2）*+*x
（1））;
functiondy=shier2（t,y）
dy=zeros（2,1）;
dy
（1）=y
（1）*
dy
（2）=y
（2）*+*y
（1））;
运行以下程序：
[t1,x]=ode45（'shier1',[015],[252]）;
[t2,y]=ode45（'shier2',[015],[252]）;
x1=x（:
1）;x2=x（:
2）;
p1=x2./（x1+x2）;
y1=y（:
1）;y2=y（:
2）;
p2=y2./（y1+y2）;
plot（t1,p1,'-',t2,p2,'*'）
图中‘*’曲线为战争中鲨鱼所占比例。
结论：
战争中鲨鱼的比例比战前高。
3、Logistic映射
logistic映射---通向混沌的道路
混沌系统，由于其行为的复杂性，往往认为其动态特性（运动方程）也一定非常复杂，事实并非如此，一个参量很少、动态特性非常简单的系统有时也能够产生混沌现象，以一维虫口模型为例，假设某一区域内的现有虫口数为yn，昆虫的繁殖率为r，且第n代昆虫不能存活于第n+1代，既无世代交叠，则第n+1代虫口数为
，r＞1时，虫口会无限制地增长；r＜1时，虫口最终会趋于消亡，因此需要对模型进行修正。
由于环境的制约和食物有限，因争夺生存空间发生相互咬斗事件的最大次数为
，即制约虫口数的因素与
成正比，设咬斗事件的战死率为β则对虫口的修正项为
，则有：
.
令
，则
（1）
取最大虫口数为1，且虫口数不能为负，则
；当
=时，方程有极大值，
而
又必须小于1，因而r＜4，则参量r的取值范围为1到4，这就得到一个抽象的标准虫口方程
（1）。
记映射
为
（2）
方程
（1）可写为
（3）
这一迭代关系通常称为logistic映射。
从[0,1]内点x0出发，由Logistic映射的迭代形成
xn=fn（x0）,n=0,1,2,…
序列{xn}称为x0的轨道。
一个看似简单的系统，随着参量的不同会表现出截然不同的行为，当r的取值范围在1~3时，方程
（1）有定态解
即方程通过多次迭代后趋于一个稳定的不动点，此时系统是稳定的。
为方程
的解，称为周期2点。
当r在3~范围内取值时，经过多次迭代，方程
（1）出现周期2点
和
，即
和
是方程
的解，满足
是使解有意义的r最小值。
随着r的增大，r=；；…依次出现周期4、周期8、周期16…的振荡解，r的极限值约为。
这种行为称为倍周期分岔，直到r＞时，系统进入了混沌状态，如下图所示，此时系统的状态不再具有规律性，而是发生随机的波动，使图d的右侧的大部分区域被涂黑了，仔细观察发现，混沌区域并非一片涂斑，而是有粗粗细细的白色“竖线”，称为周期窗口，随着参量r的增大（如
）时，混沌突然消失，系统出现周期三的稳定状态，
Logistic映射分岔图
接着倍周期分岔以更快的速度进行，再次进入混沌状态。
如果将周期窗口放大，发现其结构与分岔图的整体结构具有相似性，而且是一种无限嵌套的自相似结构。
可以看出，通过改变系统参量，使系统进入混沌的第一种模式是倍周期分岔，即由不动点→周期二→周期四→…无限倍周期→进入混沌状态。
当然通向混沌的道路不只于此，第二种通向的道路是：
从平衡态到周期运动，再到拟周期运动，直到进入混沌状态。
第三种通向混沌的方式是阵发（或间歇）道路，即系统在近似周期运动的过程中，通过改变参量，系统会出现阵发性混沌过程，随着参量的调整，阵发性混沌越来越频繁，近似的周期运动越来越少，最后进入混沌。
Matlab程序
下面程序绘制r=2，x0=的轨道
clearall;clf;
x=;
r=2;
n=input（'Pleaseinputanumber:
'）;
fori=1:
n
x1=r*x*（1-x）;
x=x1;
plot（i,x1,'.'）;
holdon;
end
下面程序绘制logistic映射r=3,5的轨道，观察是否有周期点
clearall;clf;
x=;r=;
fori=1:
100
x1=r*x*（1-x）;
x=x1;
plot（i,x1,'.'）;
holdon;
end
下面程序绘制logistic映射r=4的轨道，观察其混沌
clearall;clf;
x=;
fori=1:
100
x1=4*x*（1-x）;
x=x1;
plot（i,x1,'.'）;
holdon;
end
下面程序绘制在混沌状态下不同初值x0=和x0=的轨道的差（对初值的敏感性）
clearall;clf;
x1=;x11=;
fori=1:
1000
x2=4*x1*（1-x1）;
x1=x2;
x22=4*x11*（1-x11）;
x11=x22;
plot（i,x11-x1）;
holdon;
end
下面程序绘制logistic映射的分岔图
clearall;clf;
x=;
forr=1:
:
4
fori=1:
18
x1=r*x*（1-x）;
x=x1;
ifi>9
plot（r,x）;
holdon;
end
end
end
蛛网迭代
clearall;clf;
x=;y=0;
fori=1:
200
XY=[x,y];
y=4*x*（1-x）;
XY1=[x,y];
line（XY,XY1）;holdon;
x=y;
XY2=[x,x];
line（XY1,XY2）;holdon;
end
三、实验内容
1．求微分方程的解析解,并画出图形，
2．求微分方程的数值解,并画出图形，
3．两种相似的群体之间为了争夺有限的同一种食物来源和生活空间而进行生存竞争时，往往是竞争力较弱的种群灭亡，而竞争力较强的种群达到环境容许的最大数量。
假设有甲、乙两个生物种群，当它们各自生存于一个自然环境中，均服从Logistic规律。
（1）
是两个种群的数量；
（2）
是它们的固有增长率；
（3）
是它们的最大容量；
（4）
为种群乙（甲）占据甲（乙）的位置的数量，并且
计算
画出图形及相轨迹图。
解释其解变化过程。
2）
=，
=，计算
画出图形及相轨迹图。
解释其解变化过程。
4.绘制Logistic映射的轨道图,分岔图和蛛网迭代图。
四、实验心得
升级会员 