《实验二微分方程及差分方程模型Matlab求解》由会员分享，可在线阅读，更多相关《实验二微分方程及差分方程模型Matlab求解（10页珍藏版）》请在金锄头文库上搜索。
1、.实验二： 微分方程与差分方程模型Matlab求解一、实验目的1掌握解析、数值解法,并学会用图形观察解的形态和进行解的定性分析；2 熟悉MATLAB软件关于微分方程求解的各种命令；3通过范例学习建立微分方程方面的数学模型以及求解全过程；4 熟悉离散 Logistic模型的求解与混沌的产生过程。 二、实验原理1. 微分方程模型与MATLAB求解解析解用MATLAB命令dsolve 求常微分方程组的解析解。其中eqni表示第i个微分方程,Dny表示y的n阶导数,默认的自变量为t。1 微分方程例1 求解一阶微分方程 求通解输入：dsolve输出：ans =tan 2求特解输入：dsolveDy=1+
2、y2,y=1,x指定初值为1,自变量为x输出：ans =tan 例2求解二阶微分方程原方程两边都除以,得输入:dsolveD2y+*Dy+*y=0,y=2,Dy=-2/pi,xans =- exp*i/x + exp*exp*2*i/pi*x试试能不用用simplify函数化简输入: simplifyans =2*pi/x*sin 2微分方程组例3 求解 df/dx=3f+4g; dg/dx=-4f+3g。1通解:f,g=dsolvef =exp*C1*sin+C2*cosg =exp*C1*cos-C2*sin 特解：f,g=dsolveDf=3*f+4*g,Dg=-4*f+3*g,f=0,
3、g=1f =exp*sing =exp*cos 数值解在微分方程难以获得解析解的情况下,可以用Matlab方便地求出数值解。格式为:t,y = ode23注意： 微分方程的形式：y = F,t为自变量,y为因变量可以是多个,如微分方程组； t, y为输出矩阵,分别表示自变量和因变量的取值； F代表一阶微分方程组的函数名m文件,必须返回一个列向量,每个元素对应每个方程的右端； ts的取法有几种,1ts=t0, tf 表示自变量的取值范围,2ts=t0,t1,t2,tf,则输出在指定时刻t0,t1,t2,tf处给出,3ts=t0:k:tf,则输出在区间t0,tf的等分点给出； y0为初值条件； o
4、ptions用于设定误差限缺省是设定相对误差是10,绝对误差是10；ode23是微分方程组数值解的低阶方法,ode45为中阶方法,与ode23类似。例4 求解一个经典的范得波Van Der pol微分方程：解 形式转化：令。则以上方程转化一阶微分方程组：。编写M文件如下,必须是M文件表示微分方程组,并保存,一般地,M文件的名字与函数名相同,保存位置可以为默认的work子目录,也可以保存在自定义文件夹,这时注意要增加搜索路径FileSet PathAdd Folder function dot1=vdpol; dot1=y; 1-y2*y-y;在命令窗口写如下命令:t,y=ode23;y1=y;
5、y2=y;plot;title;xlabel;ylabelyandy执行：注：Van der Pol方程描述具有一个非线性振动项的振动子的运动过程。最初,由于它在非线性电路上的应用而引起广泛兴趣。一般形式为。图形解无论是解析解还是数值解,都不如图形解直观明了。即使是在得到了解析解或数值解的情况下,作出解的图形,仍然是一件深受欢迎的事。这些都可以用Matlab方便地进行。图示解析解如果微分方程组的解析解为：y=f ,则可以用Matlab函数fplot作出其图形：fplot其中：fun给出函数表达式；lims=xmin xmax ymin ymax限定坐标轴的大小。例如fplotsin, 0.01
6、 0.1 -11图示数值解设想已经得到微分方程组的数值解。可以用Matlab函数plot直接作出图形。其中x和y为向量或矩阵。2、Volterra模型食饵捕食者模型意大利生物学家Ancona曾致力于鱼类种群相互制约关系的研究,他从第一次世界大战期间,地中海各港口捕获的几种鱼类捕获量百分比的资料中,发现鲨鱼的比例有明显增加见下表。年代19141915191619171918百分比11.921.422.121.236.4年代19191920192119221923百分比27.316.015.914.819.7战争为什么使鲨鱼数量增加？是什么原因？因为战争使捕鱼量下降,食用鱼增加,显然鲨鱼也随之增加
7、。 但为何鲨鱼的比例大幅增加呢？生物学家Ancona无法解释这个现象,于是求助于著名的意大利数学家V.Volterra,希望建立一个食饵捕食者系统的数学模型,定量地回答这个问题。 1、符号说明：x1, x2分别是食饵、捕食者鲨鱼在t时刻的数量； r1, r2是食饵、捕食者的固有增长率；1是捕食者掠取食饵的能力, 2是食饵对捕食者的供养能力；2、基本假设：捕食者的存在使食饵的增长率降低,假设降低的程度与捕食者数量成正比,即食饵对捕食者的数量x2起到增长的作用,其程度与食饵数量x1成正比,即综合以上和,得到如下模型：模型一：不考虑人工捕获的情况 该模型反映了在没有人工捕获的自然环境中食饵与捕食者之
8、间的制约关系,没有考虑食饵和捕食者自身的阻滞作用,是Volterra提出的最简单的模型。给定一组具体数据,用matlab软件求解。 食饵：r1= 1, 1= 0.1, x10= 25; 捕食者：r2=0.5, 2=0.02, x20= 2;编制程序如下1、建立m-文件shier.m如下： function dx=shier dx=zeros; %初始化 dx=x*1-0.1*x; dx=x*-0.5+0.02*x;2、在命令窗口执行如下程序：t,x=ode45; plott,x,-,t,x,*,grid图中,蓝色曲线和绿色曲线分别是食饵和鲨鱼数量随时间的变化情况,从图中可以看出它们的数量都呈现
9、出周期性,而且鲨鱼数量的高峰期稍滞后于食饵数量的高峰期。画出相轨迹图：plotx,x模型二 考虑人工捕获的情况假设人工捕获能力系数为e,相当于食饵的自然增长率由r1 降为r1-e,捕食者的死亡率由r2 增为 r2+e,因此模型一修改为：设战前捕获能力系数e=0.3, 战争中降为e=0.1, 其它参数与模型的参数相同。观察结果会如何变化？1当e=0.3时：2当e=0.1时：分别求出两种情况下鲨鱼在鱼类中所占的比例。即计算画曲线：plott,p1,t,p2,*MATLAB编程实现建立两个M文件function dx=shier1 dx=zeros; dx=x*0.7-0.1*x; dx=x*-0.
10、8+0.02*x;function dy=shier2 dy=zeros; dy=y*0.9-0.1*y; dy=y*-0.6+0.02*y;运行以下程序：t1,x=ode45; t2,y=ode45; x1=x;x2=x; p1=x2./; y1=y;y2=y; p2=y2./; plot 图中*曲线为战争中鲨鱼所占比例。结论：战争中鲨鱼的比例比战前高。 3、 Logistic映射logistic映射-通向混沌的道路混沌系统,由于其行为的复杂性,往往认为其动态特性运动方程也一定非常复杂,事实并非如此,一个参量很少、动态特性非常简单的系统有时也能够产生混沌现象,以一维虫口模型为例,假设某一区域
11、内的现有虫口数为yn,昆虫的繁殖率为r,且第n代昆虫不能存活于第n+1代,既无世代交叠,则第n+1代虫口数为,r1时,虫口会无限制地增长；r1时,虫口最终会趋于消亡,因此需要对模型进行修正。由于环境的制约和食物有限,因争夺生存空间发生相互咬斗事件的最大次数为,即制约虫口数的因素与 成正比,设咬斗事件的战死率为则对虫口的修正项为 ,则有： .令,则 1取最大虫口数为1,且虫口数不能为负,则 ；当 =0.5时,方程有极大值,而 又必须小于1,因而r4,则参量r的取值范围为1到4,这就得到一个抽象的标准虫口方程1。记映射为 2方程1可写为 3 这一迭代关系通常称为logistic映射。从0,1内点x
12、0出发,由Logistic映射的迭代形成xn= f n, n = 0,1,2,序列xn称为x0的轨道。一个看似简单的系统,随着参量的不同会表现出截然不同的行为,当r的取值范围在13时,方程1有定态解即方程通过多次迭代后趋于一个稳定的不动点,此时系统是稳定的。为方程的解,称为周期2点。当r在33.448范围内取值时,经过多次迭代,方程1出现周期2点和,即和是方程的解,满足是使解有意义的r最小值。随着r的增大,r=3.449；3.544；3.564依次出现周期4、周期8、周期16的振荡解,r的极限值约为3.569。这种行为称为倍周期分岔,直到r3.5699时,系统进入了混沌状态,如下图所示,此时系统的状态不再具有规律性,而是发生随机的波动,使图d的右侧的大部分区域被涂黑
