实验二微分方程与差分方程模型Matlab求解文档格式.docx

实验二微分方程与差分方程模型Matlab求解文档格式.docx

《实验二微分方程与差分方程模型Matlab求解文档格式.docx》由会员分享，可在线阅读，更多相关《实验二微分方程与差分方程模型Matlab求解文档格式.docx（17页珍藏版）》请在冰豆网上搜索。

（2）求特解

'

y（0）=1'

x'

指定初值为1，自变量为x

tan（x+1/4*pi）

例2求解二阶微分方程

原方程两边都除以，得

输入:

D2y+（1/x）*Dy+（1-1/4/x^2）*y=0'

y（pi/2）=2,Dy（pi/2）=-2/pi'

-（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）

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（0）=0,g（0）=1'

exp（3*t）*sin（4*t）

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]限定坐标轴的大小。

例如

sin（1/x）'

[0.010.1-11]）

（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如下：

functiondx=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,[252]）;

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文件

functiondx=shier1（t,x）

dx

（1）=x

（1）*（0.7-0.1*x

（2））;

dx

（2）=x

（2）*（-0.8+0.02*x

（1））;

functiondy=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'

[015],[252]）;

[t2,y]=ode45（'

shier2'

x1=x（:

x2=x（:

p1=x2./（x1+x2）;

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，且虫口数不能为负，则；

当=0.5时，方程有极大值，而又必须小于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~3.448范围内取值时，经过多次迭代，方程

（1）出现周期2点和，即和是方程的解，满足

是使解有意义的r最小值。

随着r的增大，r=3.449；

3.544；

3.564…依次出现周期4、周期8、周期16…的振荡解，r的极限值约为3.569。

这种行为称为倍周期分岔，直到r＞3.5699时，系统进入了混沌状态，如下图所示，此时系统的状态不再具有规律性，而是发生随机的波动，使图d的右侧的大部分区域被涂黑了，仔细观察发现，混沌区域并非一片涂斑，而是有粗粗细细的白色“