数学建模实验二：微分方程模型Matlab求解与分析.doc

完整word)数学建模实验二：微分方程模型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*f3＊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）微分方程：

解 形式转化：令。则以上方程转化一阶微分方程组： 。
编写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方程描述具有一个非线性振动项的振动子的运动过程。最初,由于它在非线性电路上的应用而引起广泛兴趣。一般形式为。
图形解
无论是解析解还是数值解，都不如图形解直观明了。即使是在得到了解析解或数值解的情况下，作出解的图形，仍然是一件深受欢迎的事。这些都可以用Matlab方便地进行。
(1)图示解析解
如果微分方程(组）的解析解为：y=f (x），则可以用Matlab函数fplot作出其图形：
fplot('fun',lims）
其中:fun给出函数表达式；lims=[xmin xmax ymin ymax]限定坐标轴的大小。例如
fplot('sin（1/x)’, [0.01 0。1 —1 1])

(2)图示数值解
设想已经得到微分方程(组）的数值解(x,y)。可以用Matlab函数plot(x,y）直接作出图形。其中x和y为向量（或矩阵）。
2、Volterra模型（食饵捕食者模型）
意大利生物学家Ancona曾致力于鱼类种群相互制约关系的研究,他从第一次世界大战期间,地中海各港口捕获的几种鱼类捕获量百分比的资料中,发现鲨鱼的比例有明显增加（见下表）.
年代
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)分别是食饵、捕食者(鲨鱼)在t时刻的数量；
②　r1, r2是食饵、捕食者的固有增长率；
③λ1是捕食者掠取食饵的能力，
λ2是食饵对捕食者的供养能力;
2、基本假设：
① 捕食者的存在使食饵的增长率降低,假设降低的程度与捕食者数量成正比，即

②　食饵对捕食者的数量x2起到增长的作用， 其程度与食饵数量x1成正比,即

综合以上①和②,得到如下模型：
模型一：不考虑人工捕获的情况

该模型反映了在没有人工捕获的自然环境中食饵与捕食者之间的制约关系,没有考虑食饵和捕食者自身的阻滞作用，是Volterra提出的最简单的模型。
给定一组具体数据,用matlab软件求解。
食饵: r1= 1, λ1= 0。1， x10= 25；
捕食者（鲨鱼）：r2=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

图中，蓝色曲线和绿色曲线分别是食饵和鲨鱼数量随时间的变化情况,从图中可以看出它们的数量都呈现出周期性,而且鲨鱼数量的高峰期稍滞后于食饵数量的高峰期。

画出相轨迹图：
plot(x（：,1）,x(：，2）)

模型二 考虑人工捕获的情况
假设人工捕获能力系数为e，相当于食饵的自然增长率由r1 降为r1—e，捕食者的死亡率由r2 增为 r2+e，因此模型一修改为:
设战前捕获能力系数e=0.3, 战争中降为e=0.1, 其它参数与模型（一）的参数相同。观察结果会如何变化？
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）
plot(t1，p1，'—'，t2,p2,'*'）

图中‘*’曲线为战争中鲨鱼所占比例。
结论：战争中鲨鱼的比例比战前高。 
三、实验内容
1．求微分方程的解析解, 并画出图形，

解：
求解代码：
>> y=dsolve('Dy=y+2*x', ’y（0)=1’， ’x’)
输出：

画图代码：
〉> x=[0 1]；
〉> fplot（@（x)3.*exp(x)-2。*x—2，x）
输出：

2．求微分方程的数值解, 并画出图形，

解：
函数M代码:
function y=fun（t，x）
y= [x(2）；—x(1）*cos(t）］;
求解代码:
t0=0；tf=10；
［t，y]=ode23（'fun’, [t0，tf］， [1, 0］）
输出：

画图代码：
y1=y(:，1);
y2=y（:,2）;
plot（t，y1，t，y2，':'）;
输出：

3．两种相似的群体之间为了争夺有限的同一种食物来源和生活空间而进行生存竞争时，往往是竞争力较弱的种群灭亡，而竞争力较强的种群达到环境容许的最大数量.
假设有甲、乙两个生物种群,当它们各自生存于一个自然环境中，均服从 Logistic 规律。
（1）是两个种群的数量；
(2是它们的固有增长率；
(3)是它们的最大容量;
（4）为种群乙(甲）占据甲(乙)的位置的数量，并且

1）计算， 画出图形及相轨迹图。解释其解变化过程。
2），=1.5,=0。7，计算, 画出图形及相轨迹图。解释其解变化过程。
（1）
函数M代码：
function xdot=fun(～，x）
r=［1,1］；
n=［100，100］;
m=[0.5,2]；
xdot=zeros(2，1)；
xdot（1）=r（1）*x（1)＊(1-（x（1)+m（2）)/n（1)）;
xdot（2）=r（2）*x（2)＊(1—（x（2)+m(1)）/n（2））；
画图代码：
［t,x］ = ode45（'fun’, 0：0.1:15， [10 10］);
subplot（1,2，1）；plot（t， x(:,1）， t， x(:，2)， ’：’）；
legend('x1’， 'x2’）;title（’ÊýÖµ½âͼ’)；
subplot(1,2,2)；plot（x(:，1）, x(:,2）);
title(’Ïà¹ìÏßͼ'）；
输出：

分析：甲、乙两个种群一开始增长率较快，在达到某个临界值时,增长率迅速减缓,然后达到一个相对稳定的水平，此时就是环境容许的最大数量，此外，两个种群的增长率在全部时间内几乎完全一致，说明两个竞争种群在同一环境下的增长率几乎一致.
（2）
函数M代码：
function xdot=fun(~， x)
r=[1，1]；
n=［100,100］；
a=1。5；
b=0。7;
m=[a＊x(2）, b*x（1）］;
xdot=zeros（2， 1）;
xdot（1） = r（1)*x(1）＊(1—(x(1）+m（2)）/n（1））;
xdot(2） = r(2)＊x（2)＊（1-（x（2）+m(1）)/n(2));
画图代码：
[t，x］ = ode45(’fun’， 0:0。1:15, [10 10]）;
subplot(1,2，1）;plot(t， x（:,1), t, x(：,2), ’：’）;
legend（'x1'， 'x2’）；title（’ÊýÖµ½âͼ’）;
subplot(1,2，2);plot(x(：,1）, x(:，2));
title(’Ïà¹ìÏßͼ’);
输出：

分析：同样,甲、乙两个种群一开始增长率较快，在达到某个临界值时，增长率迅速减缓，然后达到一个相对稳定的水平,不同的是，竞争力较强的甲种群数量远大于乙种群的，且乙种群达到稳定水平的时间比甲早。

四、实验心得
本次实验主要学习了利用MATLAB求解常微分方程的符号解和数值解，学习了利用MATLAB画常微分方程的解的图像，此外,还学习了logistics模型的求解方法。
在这次实验过程中，本身MATLAB求解微分方程非常方便，主要难点在于对于高阶的常微分方曾，需要将它先转化为一阶常微分方程，即状态方程，才可进行求解。此外在画常微分方程的解图像和相轨图时，要注意输入的参数值和一般的函数画图不同。