《实验二微分方程与差分方程模型Matlab求解》由会员分享，可在线阅读，更多相关《实验二微分方程与差分方程模型Matlab求解（15页珍藏版）》请在金锄头文库上搜索。
1、实验二： 微分方程与差分方程模型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、2) 输出：ans =tan(t+C1) （2）求特解输入：dsolve(Dy=1+y2,y(0)=1,x) 指定初值为1，自变量为x输出：ans =tan(x+1/4*pi) 例2 求解二阶微分方程 原方程两边都除以，得输入: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/2)*2*i)/(pi*x(1/2)试试能不用用simplify函数化简输入: simplify(a
3、ns)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) 数值解在微分方程(组)难以获得解
4、析解的情况下，可以用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为初值条件； optio
5、ns用于设定误差限（缺省是设定相对误差是10(-3)，绝对误差是10(-6)）；ode23是微分方程组数值解的低阶方法，ode45为中阶方法，与ode23类似。例4 求解一个经典的范得波（Van Der pol）微分方程：解 形式转化：令。则以上方程转化一阶微分方程组： 。编写M文件如下，必须是M文件表示微分方程组，并保存，一般地，M文件的名字与函数名相同，保存位置可以为默认的work子目录，也可以保存在自定义文件夹，这时注意要增加搜索路径（FileSet PathAdd Folder） function dot1=vdpol(t,y); dot1=y(2); (1-y(1)2)*y(2)-y
6、(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方便地进行。(
7、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曾致力于鱼类种群相互制约关系的研究，他从第一次世界大战期间,地中海各港口捕获的几种鱼类捕获量
8、百分比的资料中，发现鲨鱼的比例有明显增加（见下表）。年代19141915191619171918百分比11.921.422.121.236.4年代19191920192119221923百分比27.316.015.914.819.7战争为什么使鲨鱼数量增加？是什么原因？因为战争使捕鱼量下降，食用鱼增加，显然鲨鱼也随之增加。 但为何鲨鱼的比例大幅增加呢？生物学家Ancona无法解释这个现象，于是求助于著名的意大利数学家V.Volterra，希望建立一个食饵捕食者系统的数学模型，定量地回答这个问题。 1、符号说明：x1(t), x2(t)分别是食饵、捕食者（鲨鱼）在t时刻的数量； r1, r2是食
9、饵、捕食者的固有增长率；1是捕食者掠取食饵的能力， 2是食饵对捕食者的供养能力；2、基本假设： 捕食者的存在使食饵的增长率降低，假设降低的程度与捕食者数量成正比，即食饵对捕食者的数量x2起到增长的作用， 其程度与食饵数量x1成正比，即综合以上和，得到如下模型：模型一：不考虑人工捕获的情况 该模型反映了在没有人工捕获的自然环境中食饵与捕食者之间的制约关系，没有考虑食饵和捕食者自身的阻滞作用，是Volterra提出的最简单的模型。给定一组具体数据，用matlab软件求解。 食饵： r1= 1, 1= 0.1, x10= 25; 捕食者(鲨鱼)：r2=0.5, 2=0.02, x20= 2;编制程序
10、如下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) 模型二 考
11、虑人工捕获的情况假设人工捕获能力系数为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*
12、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,*) 图中*曲线为战争中鲨鱼所占比例。结论：战争中鲨鱼的比例比战前高。 3、 Log
13、istic映射logistic映射-通向混沌的道路 混沌系统，由于其行为的复杂性，往往认为其动态特性（运动方程）也一定非常复杂，事实并非如此，一个参量很少、动态特性非常简单的系统有时也能够产生混沌现象，以一维虫口模型为例，假设某一区域内的现有虫口数为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内点x0出发，由Logistic映射的迭代形成xn= f n(x0), n = 0,1,2,序列xn称为x0的轨道。一个看似简单的系统，随着参量的不同会表现出截然不同的行为，当r的取值范围在13时，方程（1）有定态解 即方
