生物数学-微分方程的数值解

1、生物数学生物数学 微分方程的数值解 学习目的 主要内容 学习目的 主要内容 MATLAB 2、学会用、学会用Matlab求微分方程的数值解求微分方程的数值解. 数学软件数学软件 1、学会用、学会用Matlab求简单微分方程的解析解求简单微分方程的解析解. 1、求简单微分方程的解析解1、求简单微分方程的解析解. 4、实验题目.4、实验题目. 2、求微分方程的数值解、求微分方程的数值解. 3、 数学建模实例、 数学建模实例 微分方程的数值解微分方程的数值解 （一）常微分方程数值解的定义（一）常微分方程数值解的定义 （二）建立数值解法的一些途径（二）建立数值解法的一些途径 （三）用（三）用Matla

2、b软件求常微分方程的数值解软件求常微分方程的数值解 1、目标跟踪问题一：追踪问题、目标跟踪问题一：追踪问题 2、目标跟踪问题二：慢跑者与狗、目标跟踪问题二：慢跑者与狗 3、地中海鲨鱼问题、地中海鲨鱼问题 数学建模实例数学建模实例 微分方程的解析解微分方程的解析解 求微分方程（组）的解析解命令： dsolve(方程方程1, 方程方程2,方程方程n, 初始条件初始条件,自变量自变量) 记号: 在表达微分方程时，用字母 D 表示求导数，D2、D3 等 表示求高阶微分.任何 D 后所跟的字母为因变量，自变量可以指 定或由系统规则选定为确省. 例如，微分方程 0 2 2 = dx d y 应表达为：D2

3、y=0. 例例 1 求 2 1u dt du = + 的通解. 解解 输入命令：dsolve(Du=1+u2,t) 结果：u = tg(t-c) 例例 2 求微分方程的特解. = += (0)0, (0)15 4290 2 2 yy y dx dy dx d y 解解 输入命令： y=dsolve(D2y+4*Dy+29*y=0,y(0)=0,Dy(0)=15,x) 结果为 ： y =3e-2xsin(5x) 例例 3 求微分方程组的通解. =+ =+ =+ xyz dt dz xyz dt dy xyz dt dx 442 453 233 解解 输入命令 ： x,y,z=dsolve(Dx=

4、2*x-3*y+3*z,Dy=4*x-5*y+3*z,Dz=4*x-4*y+2*z, t)； x=simple(x) % 将x化简 y=simple(y) z=simple(z) 结果为：x = (c1-c2+c3+c2e -3t-c3e-3t)e2t y = -c1e-4t+c2e-4t+c2e-3t-c3e-3t+c1-c2+c3)e2t z = (-c1e-4t+c2e-4t+c1-c2+c3)e2t 微分方程的数值解微分方程的数值解 （一）常微分方程数值解的定义（一）常微分方程数值解的定义 在生产和科研中所处理的微分方程往往很复杂且大多 得不出一般解。而在实际上对初值问题，一般是要求得

5、到 解在若干个点上满足规定精确度的近似值，或者得到一个 满足精确度要求的便于计算的表达式。 因此，研究常微分方程的数值解法是十分必要的因此，研究常微分方程的数值解法是十分必要的。 0 00 0121 212 ( , ) () (), (), () , n nn yf x y x y xy xxxxxy x y xy xy yy = = 对常微分方程 ：，其数值解是指由初始点开始 的若干离散的 值处，即对， 求出准确值 的相应近似值。 （二）建立数值解法的一些途径（二）建立数值解法的一些途径 = = = + 00 i 1 , 0,1,2,1, x y(x )y yf(x,y) xh in i 解

6、微分方程：可用以下离散化方法求设 1、用差商代替导数、用差商代替导数 若步长h较小，则有 h y xhy x y x ()( ) ( ) + 故有公式： i0,1,2,n -1 () (,) 00 1 = = =+ + yy x yyhf xy iiii 此即欧拉法欧拉法。 2、使用数值积分、使用数值积分 对方程y=f(x,y), 两边由xi到xi+1积分，并利用梯形公式，有： 实际应用时，与欧拉公式结合使用： 此即改进的欧拉法改进的欧拉法。 故有公式： (1)(1 1111i 2 kkk iii i yyyy y + + + （2） 消去得： + = + = () ()() () ()()

7、22 22 Yy XxYy w dt dy Xx XxYy w dt dx （3） 3因乙舰以速度v0沿直线x=1运动，设v0=1，则w=5，X=1，Y=t 因此运动轨迹的参数方程为： = + = + = (0)0, (0)0 () (1)() 5 (1) (1)() 5 22 22 xy ty xty dt dy x xty dt dx 4. 解运动轨迹的参数方程 建立m-文件eq2.m如下： function dy=eq2(t,y) dy=zeros(2,1); dy(1)=5*(1-y(1)/sqrt(1-y(1)2+(t-y(2)2); dy(2)=5*(t-y(2)/sqrt(1-y

8、(1)2+(t-y(2)2); 取t0=0，tf=2，建立主程序chase2.m如下： t,y=ode45(eq2,0 2,0 0); Y=0:0.01:2; plot(1,Y,-)， hold on plot(y(:,1),y(:,2),*) 5. 结果见图1 大致在（1, 0.2）处击中乙舰，与前面的结论一致. 图1 图2 在chase2.m中，按二分法逐步修改tf，即分别取tf=1, 0.5, 0.25, , 直到tf=0.21时，得图2. 结论：时刻结论：时刻t=0.21时，在（时，在（1，0.21）处击中乙舰。）处击中乙舰。 慢跑者与狗慢跑者与狗 一个慢跑者在平面上沿椭圆以恒定的速率

9、v=1跑步,设 椭圆方程为: x=10+20cost, y=20+5sint. 突然有一只狗攻击他. 这 只狗从原点出发,以恒定速率w跑向慢跑者,狗的运动方向始终指向 慢跑者.分别求出w=20, w=5时狗的运动轨迹. 1. 模型建立 设时刻t慢跑者的坐标为(X(t)，Y(t)，狗的坐标为(x(t)，y(t). 则X=10+20cost, Y=20+15sint, 狗从(0,0)出发,与追踪问题类似，建立狗 的运动轨迹的参数方程: = + + = + + = (0)0, (0)0 (2015sin) (1020cos)(2015sin) (1020cos) (1020cos)(2015sin)

10、 22 22 yx ty tytx w dt dy tx tytx w dt dx 2. 模型求解 (1) w=20时时,建立m-文件eq3.m如下: function dy=eq3(t,y) dy=zeros(2,1); dy(1)=20*(10+20*cos(t)-y(1)/sqrt (10+20*cos(t)-y(1)2+(20+15*sin(t)-y(2)2); dy(2)=20*(20+15*sin(t)-y(2)/sqrt (10+20*cos(t)-y(1)2+(20+15*sin(t)-y(2)2); 取t0=0，tf=10，建立主程序chase3.m如下： t0=0;tf=1

11、0; t,y=ode45(eq3,t0 tf,0 0); T=0:0.1:2*pi; X=10+20*cos(T); Y=20+15*sin(T); plot(X,Y,-) hold on plot(y(:,1),y(:,2),*) 在chase3.m，不断修改tf的值,分别取tf=5, 2.5, 3.5,至3.15时, 狗刚好追上慢跑者. 建立m-文件eq4.m如下: function dy=eq4(t,y) dy=zeros(2,1); dy(1)=5*(10+20*cos(t)-y(1)/sqrt (10+20*cos(t)-y(1)2+(20+15*sin(t)-y(2)2); dy(

12、2)=5*(20+15*sin(t)-y(2)/sqrt (10+20*cos(t)-y(1)2+(20+15*sin(t)-y(2)2); 取t0=0，tf=10，建立主程序chase4.m如下： t0=0;tf=10; t,y=ode45(eq4,t0 tf,0 0); T=0:0.1:2*pi; X=10+20*cos(T); Y=20+15*sin(T); plot(X,Y,-) hold on plot(y(:,1),y(:,2),*) 在chase3.m，不断修改tf的值,分别取tf=20, 40, 80, 可以看出,狗永远追不上慢跑者. (2) w=5时时 地中海鲨鱼问题地中海鲨

13、鱼问题 意大利生物学家Ancona曾致力于鱼类种群相互制约关 系的研究，他从第一次世界大战期间,地中海各港口捕获的 几种鱼类捕获量百分比的资料中，发现鲨鱼等的比例有明 显增加（见下表），而供其捕食的食用鱼的百分比却明显 下降.显然战争使捕鱼量下降，食用鱼增加，鲨鱼等也随之 增加，但为何鲨鱼的比例大幅增加呢？ 他无法解释这个现象，于是求助于著名的意大利数学 家V.Volterra，希望建立一个食饵捕食系统的数学模型， 定量地回答这个问题. 年代19141915191619171918 百分比11.921.422.121.236.4 年代19191920192119221923 百分比27.316

14、.015.914.819.7 1符号说明：符号说明： ( ) 1 x t食饵在 t 时刻的数量； ( ) 2 x t捕食者在 t 时刻的数量； 1 r食饵独立生存时的增长率； 2 r捕食者独自存在时的死亡率； 1 捕食者掠取食饵的能力； 2 食饵对捕食者的供养能力. e捕获能力系数 2基本假设：基本假设： （1）食饵由于捕食者的存在使增长率降低，假设降低的程度与捕食者数量成正比； （2）捕食者由于食饵为它提供食物的作用使其死亡率降低或使之增长，假定增长 的程度与食饵数量成正比。 3模型建立与求解模型建立与求解 模型（一） 不考虑人工捕获 () 1112 1 x rx dt dx = () 22

15、21 2 xrx dt dx =+ 该 模型反映了在没有人工捕获的自然环境中食饵与捕食者之间的制约 关系，没有考虑食饵和捕食者自身的阻滞作用，是Volterra提出的最简单的 模型. 首先，建立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); 其次，建立主程序shark.m如下： t,x=ode45(shier,0 15,25 2); plot(t,x(:,1),-,t,x(:,2),*) plot(x(:,1),x(:,2) 针对一组

16、具体的数据用 Matlab 软件进行计算. 设食饵和捕食者的初始数量分别为 101 )0(xx=， 202 )0(xx= 对于数据1,0.1,0.5,0.02,25,2 20102211 =xxrr， t的终值经试验后确定为 15，即模型为： 112 221 12 (1 0.1) ( 0.50.02 ) (0)25,(0)2 xxx xxx xx = =+ = 相图(,) 12 x x为： 051015 0 10 20 30 40 50 60 70 80 90 100 020406080100 0 5 10 15 20 25 30 数值解如下图： ( ) 1 xt为实线，( ) 2 xt为“*

17、”线. 求解结果： 左图反映了x1(t)与x2(t)的关系。 可以猜测： x1(t)与x2(t)都是周期函数。 模型（二）： 考虑人工捕获 设表示捕获能力的系数为e，相当于食饵的自然增长率 由r1降为r1-e，捕食者的死亡率由r2增为 r2+e =+ = () () 2122 2 1211 1 xrex dt dx xrex dt dx = =+ = (0)25,(0)2 ( 0.80.02) (0.70.1) 21 12 2 21 1 xx xx dt dx xx dt dx = =+ = (0)25,(0)2 ( 0.60.02) (0.90.1) 21 12 2 21 1 xx xx dt dx xx dt dx 112212 1,0.1,0.5,0.02,025,02rrxx=仍取（ ）（ ） 设战前捕获能力系数e=0.3, 战争中降为e=0.1, 则战前与 战争中的模型分别为: 模型求解: 1、分别用m-文件shier1.m和shier2.m定义上述两个方程 2、建立主程序shark1.m, 求解两个方程，并画出两种情况下 鲨鱼数在鱼类总数中所占比例 x2(t)/x1(t)+x2(t) 051015 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 实线为战前的鲨 鱼比例，“*”线为 战争中的鲨鱼比例 结