动力学实验1

数学建模与实验第一章 微分方程建模1.１ 生物动力学模型 1.１.1 单种群模型 1 、Malthus 模型群体的数量变化一个复杂的问题，但是如果把影响变化的许多次要因素忽略掉或简单化，可以用微分方程来近似描述群体数量规律，从而预测某一种群未来数量。

现假设某种群体处于除了本身的出生和死亡之外，既无迁出也无迁入的状态，且其增长过程是平稳的。

设在t 时间内这种群体的数量为N ，并设他们的出生率与死亡率分别为n 与m 。

假设他们的出生数与死亡数都和t 时的群体数及时间成正比。

现在讨论群体数N 与时间t 之间的函数关系。

设],[dt t t +时间间隔内群体数量的增量为dN ，由题意，在dt 时间内这类群体的出生量与死亡量分别为nNdt 与mNdt 。

根据增量=出生量-死亡量容易得到 mNdt nNdt dN −= 即如果处始条件为00N Nt ==，解上变量可分离方程，得 则或写成)(0m n t eN N −= (1.01) 更一般地取t=n-m,则上式可表为：0kt N N e = (1.02)从上两式看出，如果m n >或t>0,则种群数量将无限增加,此称为马尔萨斯规律，它反映了群体按几何级数增长的危险性；如果n m >或t>0，则种群数量将逐步减少，趋于灭亡。

这样的结论是非常粗糙的，马尔萨斯没有考虑环境的因素．事实上在自然界，种群不可能长期地指数增长．决不会如此简单地无限增加或趋于灭亡。

为此生物学家及数学家根据统计数据对m n ,作了修正，使结果能更符合事实。

2 、Logistic 模型()dN dt n m N =−00()t N N dN dt n m N=−∫∫1Nt n m N =−现对(1.01)式作如下修正，设,bN a n −=qN p m +=，式中q p b a ,,,均为正常数。

上两式说明出生率与死亡率已不再是常数。

而是N 的线性函数，前者均匀随减小，后者均匀随N 增加。

这时方程（1.01）化为令 则上式化为()dN kN l N dt =−即 ()dN kdt N l N =−积分上式（注意到11(11+=）得 或其中0N 是 0=t 时的动物数，不论初值0N 多少，当 时，N 的极限总为 l ，即环境的总容纳量为l 。

可以用实验的方法对不同的问题，像人口的增长、传染病的发生率等来确定（1）式的图形。

这个图形称为（Logistic ）曲线。

3 、药物代谢动力学中的房室模型药物代谢动力学是研究药物在体内运转和转化过程中浓度随时间变化的规律与药物输入方式，机体状况等的联系，对于认识药物的特点，指导临床用药都有较大的实用价值．假设１ 生命系统可以设想为一系列房室构成，体内每个器观为一房室．当药物通过静脉注射、肌注或口服进入动物体内很快均匀地分布到各房室之中，V 表示假想的该室的容积，称表现分布容积；D(t)和C(t)分别表示t 时刻室内的药量和浓度；V=D(t)/C(t)．假设2 为反映药物沙拉沙星的吸收过程，在分布室前增加一个吸收室，其中C a (t)， k a分别表示分布室的药物浓度和吸收速率常数．假设3 室内物质外流速度正比于该室内的药量． 根据上假设，可得一室模型：()()()dN a bN Ndt p qN Ndt a p dN b q N N dtb q =−−+−=+ +,a p k b q l b q−=+=+0()tNN dN kdt N l N =−∫∫001ln N Nt lk l N l N =−−0000()11kltklt lN l N N l N el e N −==+−+−t→+∞00()(0)dC kC t dt D C V D =− =其中,表示注射药物总量,k >0为清除速率,C (t )为浓度.00()()ktktD C t C e C t V−−==不难解出或其中V 与k 为要估计的模型参数.1.1.2 简单动力学模型求解的MATLAB 命令 1．解微分方程的的MA TLAB 命令常微分方程的初值问题的标准数学表述为：'(,),()0y f t y a t by a y =≤≤= 其中t 为自变量，y 为因变量（变量y 可以是向量，例如微分方程组），yo 为初始值。

任何高阶常微分方程都可以用替换法化为上式所示的一阶形式。

（１）符号解析解 Dsolve(‘eqn1’, ‘eqn2’,…,’x’)其中,’eqni’表示第i 个方程与初始条件等式,’X’代表微分方程(组)的自变量，默认自变量为t.方程中导数用Ｄ表示，２阶导数用Ｄ２表示（２）Matlab 利用数值方法来求解常微分方程的解，其思路如下：把求解的时间区间划分成有限步，对应于每一步将计算出一个解，如果求得的解不满足误差限制，则减少步长，再求解。

如此重复，直到满足误差限为止。

○a 刚性问题（stiff ）:方程组的解不同分量的数量级差别较大，对于数值求解是一大困难。

Matlab 既能解决非刚性问题，也能解决刚性问题。

○b 三个解决非刚性问题的函数：ode45,ode23,ode113 ○c 两个解刚性问题的函数：ode15s 和ode23s在MA TLAB 中，用2阶（3阶）龙格-库塔公式和4阶（5阶）龙格-库塔公式的程序分别为[t ，y]=ode23(‘F’,ts,y0,options) [t ，y]=ode45(‘F’,ts,y0,options)其中F 是由微分方程（组）写成的M 文件名．输入ts 的取法有几种，当ts=[t0，tf]，t0，tf 分别表示自变量的初值和终值；若ts=[t0，t1，t2，…，tf]，则输出在指定时刻t0，t1，t2，…，tf 处给出；对于等步长时用ts=t0：k ：tf ，则输出在区间[t0，tf]的等分点给出． y0为函数的初值，options 用于设定误差限．程序为options=odeset(‘reltol’,rt,’abstol’,at)[t ，y]为输出矩阵，分别表示自变量t 和因变量y 的取值．Ode23是微分方程组数值解的低阶方法，ode45为较高阶方法，与ode23类似．另外还有一些其他方法，如求解非刚性微分方程组的可变阶方法ode123．２．利用MA TLAB 的解微分方程的步骤 （1）：化方程为标准形式。

例如：y’’’-3y’’-y’y=0,y(0)=0,y’(0)=1,y’’(0)=-1.把微分方程的高阶导数写为低阶导数的算式，即：y’’’=3y’’+y’y,设：y1=y,y2=y’,y3=y’’,则原方程化为下列等价的方程组：+===1233'33'22'1y y y y y y y y 满足初值条件：−===1)0(31)0(20)0(1y y y 已把该方程化成了标准形式。

其中：y’ ->（y1’,y2’,y3’）,a->(0,0,0),y0->(0,1,-1),f(t,y)->(y2,y3,3y3+y2y1).(２)、把微分方程组编成m 函数文件(命令文件和函数文件)MA TLAB 提供了一个内置的具有编辑和调试功能的程序编辑器.从“File ”菜单中选择“New ”下的“M －file ”命令，即可进入程序编辑器（MA TLAB Edior/Debug ），,输入MA TLAB 程序，然后保存，即建立了一个M 文件. 该Ｍ文件都可被别的M 文件调用.如图1.4：图1.4M 文件的编写函数M 文件能够象库函数一样方便地调用，从而扩展MATLAB 的功能，其第一行比较特殊，其形式必须为： F u n c t i o n f ＝函数名（t ,y ） f=表达式；f=f(:);如上例中：function dy=F(t,y)dy=[y(2);y(3);3*y(3)+y(2)*y(1)];注意：A:在函数文件里，虽然写微分方程时并不同时包含参数t 和y ，但第一行必须包含这两个输入变量。

B:向量dy 必须为列向量。

或（含待定参数） f u n c t i o n [因变量列表]＝函数名（B e t a ,自变量列表） a =B e t a (1); b =B e t a (2); ………… c =B e t a (m );因变量表＝含自变量的表达式；当输出变量多于一个时，应该用方括号，自变量多于一个时应该用逗号隔开，编写完以后，必须以函数名存盘，否则不能被调用，函数M 文件不能访问工作区中的变量.（3）：调用一个微分方程的求解函数求解。

[t,Y]=Ｄsolver(‘F’,tspan,y0); 其中：solver:求解函数名；F:包含微分方程的m 文件;tspan 为积分的数据范围，其格式为：[t0,tfinal]; y0为t0时刻的初值列向量。

输出参数T 和Y 为列向量 t 为时刻向量。

Y 表是不同时刻的函数值。

(4)把函数M 文件的调入数值解作图形模拟例．求解方程y’’-3(1-y^2)y’+y=0在初值y’(0)=3,y(0)=2的解，写出其求解常微分方程初值问题的完整过程。

○1化成标准形式： 设y1=y,y2=y’,则：−−==12)2^11(3'22'1y y y y y y 初值为： ==3)0(22)0(1y y○2编写函数文件ode.m,类容为： function dy=ode(t,y)dy=[y(2);3*(1-y(1)^2)*y(2)-y(1)];○3调用函数ode45求解，时间区间为[0,20]: [T,Y]=ode45(‘ode’,[0,20],[2;3]);输出结果[T,Y]中T 为时间点组成的向量。

Y 为对应于T 中时间点的y(1)和y(2)的值。

○4绘制解的曲线，结果如图。

plot(T,Y(:,1),’-’,T,Y(:,2),’--’)title(‘Solution of ODE Equation’); xlabel(‘time T’) ylabel(‘solution Y’); legend(‘Y1’,’Y2’)Matlab 利用数值方法来求（4.1）．把函数M 文件的调入数值解作图形模拟 clear; close;y0=初始值; h=步长值; for i=1:循环数[t,y]=ode45('函数名',[t0,tf],y0); plot(t,y); hold on; y0=y0+h; end;利用MA TLAB 的非线性回归对微分方程的解析解作的参数估计（5）．把函数M 文件的调入非线性回归估计在MATLAB 命令窗口输入 Ｘ＝［自变量数据．．．］； Ｙ＝［因自变量数据．．．］； Beta0=[初值]；Beta=nlinfit(x,y,’函数名＇Beta0)[Beta,R,J]=nlinfit(x,y,’函数名＇Beta0) (6)模型的应用（求极值，预测等） 三、求函数的极值1．求解一元函数的最小值调用格式如下：x=fminbnd('fun',a ,b)函数fminbnd 来求一元函数)(x f y =在指定区间[a ,b ]上的函数局部极小值，该函数返回函数在极小值点时自变量x 的值，fun 是MATLAB 内置的M 文件函数2．求解多元函数的最小值函数fminsearch用于求多元函数在向量x0附近的最小值.它指定一个开始的向量（x0），并非指定一个区间.此函数返回一个向量为此多元函数局部最小函数值对应的自变量的取值，调用格式如下：x=fminsearch('fun',x0)例8 把一个3个自变量的函数创建在一个M文件里%three.mfunction b=three(v)x=v(1);y=v(2);z=v(3);b=x*x+2.5*sin(y)-z*z*x*y*y;求这个函数在[1,-1,0]点附近的最小值可以得到：>> v=[1 -1 0];>> fminsearch('three',v)ans =-0.0000 -1.5708 0.0008也可简单输入如下：f='x(1)^2+2.5*sin(x(2))-x(3)^2*x(1)*x(2)^2';>> x=fminsearch(f,[1 -1 0]),f=eval(f)x =-0.0000 -1.5708 0.0008 ％函数的最小值点f =-2.5000 ％函数的最小值例六、七种类型曲线。