数学建模第三讲幻灯片.ppt

《数学建模第三讲幻灯片.ppt》由会员分享，可在线阅读，更多相关《数学建模第三讲幻灯片.ppt（75页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、数学建模第三讲第1页，共75页，编辑于2022年，星期六第三讲第三讲 微分方程模型微分方程模型动态动态模型模型 描述对象特征随时间描述对象特征随时间(空间空间)的演变过程的演变过程 分析对象特征的变化规律分析对象特征的变化规律 预报对象特征的未来性态预报对象特征的未来性态 研究控制对象特征的手段研究控制对象特征的手段 根据函数及其变化率之间的关系确定函数根据函数及其变化率之间的关系确定函数微分微分方程方程建模建模 根据建模目的和问题分析作出简化假设根据建模目的和问题分析作出简化假设 按照内在规律或用类比法建立微分方程按照内在规律或用类比法建立微分方程第2页，共75页，编辑于2022年，星期六主

2、要内容主要内容生物单种群增长模型生物单种群增长模型 3.1 人口增长模型人口增长模型 3.2 传染病模型传染病模型生物多种群增长模型生物多种群增长模型 3.3 正规战与游击战正规战与游击战 3.4 捕食系统的捕食系统的Volterra方程方程 第3页，共75页，编辑于2022年，星期六 为了保持自然资料的合理开发与利用，人类必须保持并控制为了保持自然资料的合理开发与利用，人类必须保持并控制生态平衡，甚至必须控制人类自身的增长。生态平衡，甚至必须控制人类自身的增长。本节将建立几个简单的单种群增长模型，以简略分析一下本节将建立几个简单的单种群增长模型，以简略分析一下这方面的问题。一般生态系统的分析

3、可以通过一些简单模型的这方面的问题。一般生态系统的分析可以通过一些简单模型的复合来研究，大家若有兴趣可以根据生态系统的特征自行建立复合来研究，大家若有兴趣可以根据生态系统的特征自行建立相应的模型。相应的模型。美丽的大自然 种群的数量本应取离散值，但由于种群数量一种群的数量本应取离散值，但由于种群数量一般较大，为建立微分方程模型，可将种群数量看般较大，为建立微分方程模型，可将种群数量看作连续变量，甚至允许它为可微变量，由此引起作连续变量，甚至允许它为可微变量，由此引起的误差将是十分微小的。的误差将是十分微小的。离散化为连续，方便研离散化为连续，方便研究究3.1 如何预报人口的增长如何预报人口的增

4、长 -MalthusMalthus模型与模型与LogisticLogistic模型模型模型模型第4页，共75页，编辑于2022年，星期六背景背景 年年 1625 1830 1930 1960 1974 1987 1999人口人口(亿亿)5 10 20 30 40 50 60世界人口增长概况世界人口增长概况中国人口增长概况中国人口增长概况 年年 1908 1933 1953 1964 1982 1990 1995 2000人口人口(亿亿)3.0 4.7 6.0 7.2 10.3 11.3 12.0 13.0研究人口变化规律研究人口变化规律控制人口过快增长控制人口过快增长第5页，共75页，编辑于2

5、022年，星期六指数增长模型指数增长模型马尔萨斯提出马尔萨斯提出(1798)常用的计算公式常用的计算公式x(t):时刻时刻t的的人口人口基本假设基本假设:人口人口(相对相对)增长率增长率 r 是常数是常数,不考虑移民不考虑移民今年人口今年人口 x0,年增长率年增长率 rk年后人口年后人口随着时间增加，人口按指数规律无限增长随着时间增加，人口按指数规律无限增长第6页，共75页，编辑于2022年，星期六模型检验模型检验 比较历年的人口统计资料，可发现人口增长的实际情况比较历年的人口统计资料，可发现人口增长的实际情况与马尔萨斯模型的预报结果基本相符，例如，与马尔萨斯模型的预报结果基本相符，例如，19

6、61年世界人年世界人口数为口数为30.6（即（即3.06109），人口增长率约为），人口增长率约为2%，人口数大，人口数大约每约每35年增加一倍。检查年增加一倍。检查1700年至年至1961的的260年人口实际数年人口实际数量，发现两者几乎完全一致，且按马氏模型计算，人口数量量，发现两者几乎完全一致，且按马氏模型计算，人口数量每每34.6年增加一倍，两者也几乎相同。年增加一倍，两者也几乎相同。模型预测模型预测 假如人口数真能保持每假如人口数真能保持每34.6年增加一倍，那么人口数将年增加一倍，那么人口数将以几何级数的方式增长。例如，到以几何级数的方式增长。例如，到2510年，人口达年，人口达2

7、1014个，个，即使海洋全部变成陆地，每人也只有即使海洋全部变成陆地，每人也只有9.3平方英尺的活动范围，平方英尺的活动范围，而到而到2670年，人口达年，人口达361015个，只好一个人站在另一人的肩个，只好一个人站在另一人的肩上排成二层了。上排成二层了。故故马尔萨斯模型是不完善的。马尔萨斯模型是不完善的。几何级数的增长MalthusMalthus模型实际上只有在群体总数不模型实际上只有在群体总数不太大时才合理，到总数增大时，生物太大时才合理，到总数增大时，生物群体的各成员之间由于有限的生存空群体的各成员之间由于有限的生存空间，有限的自然资源及食物等原因，间，有限的自然资源及食物等原因，就可

8、能发生生存竞争等现象。就可能发生生存竞争等现象。所以所以MalthusMalthus模型假设的人口净增模型假设的人口净增长率不可能始终保持常数，它应当长率不可能始终保持常数，它应当与人口数量有关。与人口数量有关。第7页，共75页，编辑于2022年，星期六指数增长模型的应用及局限性指数增长模型的应用及局限性 与与19世纪以前欧洲一些地区人口统计数据吻合世纪以前欧洲一些地区人口统计数据吻合 适用于适用于19世纪后迁往加拿大的欧洲移民后代世纪后迁往加拿大的欧洲移民后代 可用于短期人口增长预测可用于短期人口增长预测 不符合不符合19世纪后多数地区人口增长规律世纪后多数地区人口增长规律 不能预测较长期的

9、人口增长过程不能预测较长期的人口增长过程1919世纪后人口数据世纪后人口数据人口增长率人口增长率r r不是常数不是常数(逐渐下降逐渐下降)第8页，共75页，编辑于2022年，星期六阻滞增长模型阻滞增长模型(Logistic模型模型)人口增长到一定数量后，增长率下降的原因：人口增长到一定数量后，增长率下降的原因：资源、环境等因素对人口增长的阻滞作用资源、环境等因素对人口增长的阻滞作用且阻滞作用随人口数量增加而变大且阻滞作用随人口数量增加而变大假设假设r固有增长率固有增长率(x很小时很小时)xm人口容量（资源、环境能容纳的最大数量）人口容量（资源、环境能容纳的最大数量）r是是x的减函数的减函数第9

10、页，共75页，编辑于2022年，星期六dx/dtx0 xmxm/2xmtx0 x(t)S形曲线形曲线,x增加先快后慢增加先快后慢x0 xm/2阻滞增长模型阻滞增长模型(Logistic模型模型)第10页，共75页，编辑于2022年，星期六参数估计参数估计用指数增长模型或阻滞增长模型作人口用指数增长模型或阻滞增长模型作人口预报，必须先估计模型参数预报，必须先估计模型参数 r 或或 r,xm 利用统计数据用最小二乘法作拟合利用统计数据用最小二乘法作拟合例：美国人口数据（单位例：美国人口数据（单位百万）百万）1860 1870 1880 1960 1970 1980 1990 31.4 38.6 5

11、0.2 179.3 204.0 226.5 251.4专家估计专家估计阻滞增长模型阻滞增长模型(Logistic模型模型)r=0.2557,xm=392.1继续继续第11页，共75页，编辑于2022年，星期六最小二乘法最小二乘法设经实际测量已得设经实际测量已得 到到n组数据（组数据（xi,yi），），i=1,n。将数据画在。将数据画在平面直角坐标系中，见平面直角坐标系中，见 图。如果建模者判断图。如果建模者判断 这这n个点很象是分个点很象是分布在某条直线附近，令布在某条直线附近，令 该直线方程该直线方程 为为y=ax+b，进而利用数据来，进而利用数据来求参求参 数数a和和b。由于该直线只是数据

12、近似满足的关系式，故。由于该直线只是数据近似满足的关系式，故 yi-(axi+b)=0一般不成立，但我们希望一般不成立，但我们希望 最小最小此式对此式对a和和b的偏导数均的偏导数均 为为0，解，解相应方程组，求得：相应方程组，求得：y=ax+byO(xi,yi)x其中其中 和和 分别为分别为xi和和yi的平均值的平均值 如果建模者判断变量间的关系并非线性关系而是其他类型的函数，如果建模者判断变量间的关系并非线性关系而是其他类型的函数，则可作则可作 变量替换变量替换使之转化为线性关系或用类似方使之转化为线性关系或用类似方 法法拟合拟合。第12页，共75页，编辑于2022年，星期六用用MATLAB

13、作线性最小二乘拟合作线性最小二乘拟合1.1.作多项式作多项式f(x)=a1xm+amx+am+1拟合拟合,可利用已有程序可利用已有程序:a=polyfit(x,y,m)输出拟合多项式系数输出拟合多项式系数a=a1,am,am+1(数组数组)）输入同长度输入同长度的数组的数组X，Y拟合多项拟合多项式次数式次数第13页，共75页，编辑于2022年，星期六1.lsqcurvefit1.lsqcurvefit已知数据点数据点：xdata=xdata=（xdata1，xdata2，xdataxdatan n），），ydata=ydata=（ydataydata1 1，ydataydata2 2，ydat

14、aydatan n）用用MATLAB作非线性最小二乘拟合作非线性最小二乘拟合 Matlab Matlab的提供了两个求非线性最小二乘拟合的函数：的提供了两个求非线性最小二乘拟合的函数：lsqcurvefitlsqcurvefit和lsqnonlinlsqnonlin。两个命令都要先建立。两个命令都要先建立M-M-文件文件fun.mfun.m，在其，在其中定义函数中定义函数f(x)f(x)，但两者定义，但两者定义f(x)f(x)的方式是不同的的方式是不同的,可参考例题可参考例题.lsqcurvefitlsqcurvefit用以求含参量用以求含参量x x（向量）的向量值函数（向量）的向量值函数F(

15、x,xdata)=F(x,xdata)=（F F（x x，xdataxdata1 1），），F F（x x，xdataxdatan n）T T中的参变量中的参变量x(x(向量向量),),使得使得 第14页，共75页，编辑于2022年，星期六 输入格式为输入格式为:(1)x=lsqcurvefit(fun,x0,xdata,ydata);(2)x=lsqcurvefit(fun,x0,xdata,ydata,options);(3)x=lsqcurvefit(fun,x0,xdata,ydata,options,grad);(4)x,options=lsqcurvefit(fun,x0,xdat

16、a,ydata,);(5)x,options,funval=lsqcurvefit(fun,x0,xdata,ydata,);(6)x,options,funval,Jacob=lsqcurvefit(fun,x0,xdata,ydata,);fun是一个事先建立的是一个事先建立的定义函数定义函数F(x,xdata)的的M-文件文件,自变量为自变量为x和和xdata说明：x=lsqcurvefit(fun,x0,xdata,ydata,options);迭代初值迭代初值已知数据点已知数据点选项见无选项见无约束优化约束优化第15页，共75页，编辑于2022年，星期六 lsqnonlin用以求含参

17、量用以求含参量x x（向量）的向量值函数（向量）的向量值函数 f(x)=(f1(x),f2(x),fn(x)T 中的参量中的参量x x，使得，使得 最小。最小。其中其中 fi（x）=f（x，xdatai，ydatai）=F(x,xdatai)-ydatai 2.lsqnonlin已知数据点：已知数据点：xdata=（xdata1，xdata2，xdatan）ydata=（ydata1，ydata2，ydatan）第16页，共75页，编辑于2022年，星期六输入格式为：输入格式为：1）x=lsqnonlin（fun，x0）；）；2）x=lsqnonlin（fun，x0，options）；）；3）

18、x=lsqnonlin（fun，x0，options，grad）；）；4）x，options=lsqnonlin（fun，x0，）；）；5）x，options，funval=lsqnonlin（fun，x0，）；说明：x=lsqnonlinlsqnonlin（fun，x0，options）；）；fun是一个事先建立的定是一个事先建立的定义函数义函数f(x)的的M-文件，文件，自变量为自变量为x迭代初值迭代初值选项有无选项有无约束优化约束优化第17页，共75页，编辑于2022年，星期六 例例 用下面一组数据拟合用下面一组数据拟合 中的参数中的参数a，b，k该问题即解最优化问题：该问题即解最优化问

19、题：第18页，共75页，编辑于2022年，星期六 1 1）编写）编写M-M-文件文件 curvefun1.m curvefun1.m function f=curvefun1(x,tdata)f=x(1)+x(2)*exp(-0.02*x(3)*tdata)%其中其中 x(1)=a;x(2)=b；x(3)=k;2）输入命令）输入命令tdata=100:100:1000tdata=100:100:1000cdata=cdata=1e-03*4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,*4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.3

20、9,6.50,6.59;6.50,6.59;x0=0.2,0.05,0.05;x0=0.2,0.05,0.05;x=lsqcurvefit(curvefun1,x0,tdata,cdata)x=lsqcurvefit(curvefun1,x0,tdata,cdata)f=f=curvefun1(x,tdata)F(xF(x，tdata)=tdata)=，x=(ax=(a，b b，k)k)解法解法1 1.用命令用命令lsqcurvefitlsqcurvefit第19页，共75页，编辑于2022年，星期六3）运算结果为）运算结果为：f=0.0043 0.0051 0.0056 0.0059 0.0

21、061 0.0062 0.0062 0.0063 0.0063 0.0063 x=0.0063 -0.0034 0.25424）结论）结论：a=0.0063,b=-0.0034,k=0.2542返回第20页，共75页，编辑于2022年，星期六模型检验模型检验用模型计算用模型计算2000年美国人口，与实际数据比较年美国人口，与实际数据比较实际为实际为281.4(百万百万)模型应用模型应用预报美国预报美国2010年的人口年的人口加入加入2000年人口数据后重新估计模型参数年人口数据后重新估计模型参数Logistic 模型在经济领域中的应用模型在经济领域中的应用(如耐用消费品的售量如耐用消费品的售量

22、)阻滞增长模型阻滞增长模型(Logistic模型模型)r=0.2490,xm=434.0 x(2010)=306.0第21页，共75页，编辑于2022年，星期六模型检验模型检验 用用LogisticLogistic模型来描述种群增长的规律效果如何呢？模型来描述种群增长的规律效果如何呢？19451945年克朗皮克（年克朗皮克（CrombicCrombic）做了一个人工饲养小谷虫的实验，数）做了一个人工饲养小谷虫的实验，数学生物学家高斯（学生物学家高斯（EFGaussEFGauss）也做了一个原生物草履虫实）也做了一个原生物草履虫实验，实验结果都和验，实验结果都和LogisticLogistic曲

23、线十分吻合。曲线十分吻合。大量实验资料表明用大量实验资料表明用LogisticLogistic模型来描述种群的增长，效模型来描述种群的增长，效果还是相当不错的。例如，高斯果还是相当不错的。例如，高斯把把5只草履虫放进一个盛有只草履虫放进一个盛有0.5cm3营养液的小试管，他发现，开始时草履虫以每天营养液的小试管，他发现，开始时草履虫以每天230.9%的速率增长，此后增长速度不断减慢，到第五天达到最大量的速率增长，此后增长速度不断减慢，到第五天达到最大量375个，实验数据与个，实验数据与r=2.309，a=0.006157，N(0)=5的的LogisticLogistic曲线：曲线：几乎完全吻合

24、，见图几乎完全吻合，见图3.6。图图3-6第22页，共75页，编辑于2022年，星期六MalthusMalthus模型和模型和LogisticLogistic模型的总结模型的总结 MalthusMalthus模型和模型和LogisticLogistic模型模型均为对微分方程（均为对微分方程（3.7）所作的模拟近似方程。前一模型假设了种群增长率所作的模拟近似方程。前一模型假设了种群增长率r为一常为一常数，（数，（r被称为该种群的内禀增长率）。后一模型则假设环被称为该种群的内禀增长率）。后一模型则假设环境只能供养一定数量的种群，从而引入了一个竞争项。境只能供养一定数量的种群，从而引入了一个竞争项。

25、用模拟近似法建立微分方程来研究实际问题时必须对用模拟近似法建立微分方程来研究实际问题时必须对求得的解进行检验，看其是否与实际情况相符或基本相符。求得的解进行检验，看其是否与实际情况相符或基本相符。相符性越好则模拟得越好，否则就得找出不相符的主要原相符性越好则模拟得越好，否则就得找出不相符的主要原因，对模型进行修改。因，对模型进行修改。Malthus Malthus模型与模型与LogisticLogistic模型虽然都是为了研究种群数量的增模型虽然都是为了研究种群数量的增长情况而建立的，但它们也可用来研究其他实际问题，只要这些实长情况而建立的，但它们也可用来研究其他实际问题，只要这些实际问题的数

26、学模型有相同的微分方程即可，下面我们来看两个较为际问题的数学模型有相同的微分方程即可，下面我们来看两个较为有趣的实例。有趣的实例。第23页，共75页，编辑于2022年，星期六 年龄分布对于人口预测的重要性年龄分布对于人口预测的重要性 只考虑自然出生与死亡，不计迁移只考虑自然出生与死亡，不计迁移人口发展方程人口发展方程第24页，共75页，编辑于2022年，星期六人口发展方程人口发展方程一阶偏微分方程一阶偏微分方程第25页，共75页，编辑于2022年，星期六例例2:传染病模型传染病模型问题问题 描述传染病的传播过程描述传染病的传播过程 分析受感染人数的变化规律分析受感染人数的变化规律 预报传染病高

27、潮到来的时刻预报传染病高潮到来的时刻 预防传染病蔓延的手段预防传染病蔓延的手段 按照传播过程的一般规律，用机按照传播过程的一般规律，用机理分析方法建立模型理分析方法建立模型第26页，共75页，编辑于2022年，星期六 已感染人数已感染人数(病人病人)i(t)每个病人每天有效接触每个病人每天有效接触(足以使人致病足以使人致病)人数为人数为 模型模型1 1假设假设若有效接触的是病人，则若有效接触的是病人，则不能使病人数增加不能使病人数增加必须区分已感染者必须区分已感染者(病病人人)和未感染者和未感染者(健康人健康人)建模建模?第27页，共75页，编辑于2022年，星期六模型模型2 2区分已感染者区

28、分已感染者(病人病人)和未感染者和未感染者(健康人健康人)假设假设1）总人数）总人数N不变，病人和健康不变，病人和健康 人的人的 比例分别为比例分别为 2）每个病人每天有效接触人数为）每个病人每天有效接触人数为,且且使接触的健康人致病使接触的健康人致病建模建模 日日接触率接触率SI 模型模型第28页，共75页，编辑于2022年，星期六模型模型21/2tmii010ttm传染病高潮到来时刻传染病高潮到来时刻 (日接触率日接触率)tm Logistic 模型病人可以治愈！病人可以治愈！?t=tm,di/dt 最大最大第29页，共75页，编辑于2022年，星期六模型模型3传染病无免疫性传染病无免疫性

29、病人治愈成为健病人治愈成为健康人，健康人可再次被感染康人，健康人可再次被感染增加假设增加假设SIS 模型模型3）病人每天治愈的比例为）病人每天治愈的比例为 日日治愈率治愈率建模建模 日接触率日接触率1/感染期感染期 一个感染期内一个感染期内每个病人的有每个病人的有效接触人数，称为效接触人数，称为接触数接触数。第30页，共75页，编辑于2022年，星期六1-1/idi/dt01 1i0i00ti 11-1/i0i0t 1di/dt 1/i(t)先升后降至先升后降至0P2:s01/i(t)单调降至单调降至01/阈值阈值P3P4P2S0第36页，共75页，编辑于2022年，星期六模型模型4SIR模型

30、模型预防传染病蔓延的手段预防传染病蔓延的手段 (日接触率日接触率)卫生水平卫生水平 (日日治愈率治愈率)医疗水平医疗水平 传染病不蔓延的条件传染病不蔓延的条件s01/降低降低 s0提高提高 r0 提高阈值提高阈值 1/降低降低 (=/),群体免疫群体免疫第37页，共75页，编辑于2022年，星期六 的估计的估计第38页，共75页，编辑于2022年，星期六模型模型4SIR模型模型被传染人数的估计被传染人数的估计记被传染人数比例记被传染人数比例xs0i0P1i0 0,s0 1 小小,s0 1提高阈值提高阈值1/降低降低被传被传染人数比例染人数比例 xs0-1/=第39页，共75页，编辑于2022年

31、，星期六例例3：正规战与游击战正规战与游击战战争分类：正规战争，游击战争，混合战争战争分类：正规战争，游击战争，混合战争只考虑双方兵力多少和战斗力强弱只考虑双方兵力多少和战斗力强弱兵力因战斗及非战斗减员而减少，因增援而增加兵力因战斗及非战斗减员而减少，因增援而增加战斗力与射击次数及命中率有关战斗力与射击次数及命中率有关建模思路和方法为用数学模型讨论社会领域的实际问题建模思路和方法为用数学模型讨论社会领域的实际问题提供了可借鉴的示例提供了可借鉴的示例第一次世界大战第一次世界大战Lanchester提出预测战役结局的模型提出预测战役结局的模型第40页，共75页，编辑于2022年，星期六一般模型一般

32、模型 每方战斗减员率取决于双方的兵力和战斗力每方战斗减员率取决于双方的兵力和战斗力 每方非战斗减员率与本方兵力成正比每方非战斗减员率与本方兵力成正比 甲乙双方的增援率为甲乙双方的增援率为u(t),v(t)f,g 取决于战争类型取决于战争类型x(t)甲方兵力，甲方兵力，y(t)乙方兵力乙方兵力模型模型假设假设模型模型第41页，共75页，编辑于2022年，星期六正规战争模型正规战争模型 甲方战斗减员率只取决于乙方的兵力和战斗力甲方战斗减员率只取决于乙方的兵力和战斗力双方均以正规部队作战双方均以正规部队作战 忽略非战斗减员忽略非战斗减员 假设没有增援假设没有增援f(x,y)=ay,a 乙方每个士兵的

33、杀伤率乙方每个士兵的杀伤率a=ry py,ry 射击率，射击率，py 命中率命中率第42页，共75页，编辑于2022年，星期六0正规战争模型正规战争模型为判断战争的结局，不求为判断战争的结局，不求x(t),y(t)而而在相平面上讨论在相平面上讨论 x 与与 y 的关系的关系平方律平方律 模型模型乙方胜乙方胜第43页，共75页，编辑于2022年，星期六游击战争模型游击战争模型双方都用游击部队作战双方都用游击部队作战 甲方战斗减员率还随着甲方兵力的增加而增加甲方战斗减员率还随着甲方兵力的增加而增加 忽略非战斗减员忽略非战斗减员 假设没有增援假设没有增援f(x,y)=cxy,c 乙方每个士兵的杀伤率

34、乙方每个士兵的杀伤率c=ry pyry射击率射击率py 命中率命中率py=sry/sxsx 甲方活动面积甲方活动面积sry 乙方射击有效面积乙方射击有效面积第44页，共75页，编辑于2022年，星期六0游击战争模型游击战争模型线性律线性律 模型模型第45页，共75页，编辑于2022年，星期六0混合战争模型混合战争模型甲方为游击部队，乙方为正规部队甲方为游击部队，乙方为正规部队乙方必须乙方必须10倍于甲方的兵力倍于甲方的兵力设设 x0=100,rx/ry=1/2,px=0.1,sx=1(km2),sry=1(m2)第46页，共75页，编辑于2022年，星期六微分方程的解析解微分方程的解析解 求微

35、分方程（组）的解析解命令求微分方程（组）的解析解命令:dsolve(方程方程1,方程方程2,方程方程n,初始条件初始条件,自变量自变量)结果：u=tg(t-c)第47页，共75页，编辑于2022年，星期六 解解 输入命令输入命令:y=dsolve(D2y+4*Dy+29*y=0,y(0)=0,Dy(0)=15,x)结结 果果 为为:y=3e-2xsin（5x）第48页，共75页，编辑于2022年，星期六解解 输入命令：x,y,z=dsolve(Dx=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

36、(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 第49页，共75页，编辑于2022年，星期六微分方程的数值解微分方程的数值解（一）常微分方程数值解的定义（一）常微分方程数值解的定义 在生产和科研中所处理的微分方程往往很复杂且大多得不出一在生产和科研中所处理的微分方程往往很复杂且大多得不出一般解。而在实际上对初值问题，一般是要求得到解在若干个点上满般解。而在实际上对初值问题，一般是要求得到解在若

37、干个点上满足规定精确度的近似值，或者得到一个满足精确度要求的便于计算足规定精确度的近似值，或者得到一个满足精确度要求的便于计算的表达式的表达式。因此，研究常微分方程的数值解法是十分必要的因此，研究常微分方程的数值解法是十分必要的。第50页，共75页，编辑于2022年，星期六（二）建立数值解法的一些途径1、用差商代替导数、用差商代替导数 若步长若步长h较小，则有较小，则有故有公式：此即欧拉法欧拉法。第51页，共75页，编辑于2022年，星期六2、使用数值积分、使用数值积分对方程y=f(x,y),两边由xi到xi+1积分，并利用梯形公式，有：实际应用时，与欧拉公式结合使用：实际应用时，与欧拉公式结

38、合使用：此即改进的欧拉法改进的欧拉法。故有公式：第52页，共75页，编辑于2022年，星期六3、使用泰勒公式、使用泰勒公式 以此方法为基础，有龙格以此方法为基础，有龙格-库塔法、线性多步法等方法。库塔法、线性多步法等方法。4、数值公式的精度、数值公式的精度 当一个数值公式的截断误差可表示为当一个数值公式的截断误差可表示为O（hk+1）时（）时（k为为正整数，正整数，h为步长），称它是一个为步长），称它是一个k阶公式。阶公式。k越大，则数值公式的精度越高。越大，则数值公式的精度越高。欧拉法是一阶公式，改进的欧拉法是二阶公式。欧拉法是一阶公式，改进的欧拉法是二阶公式。龙格龙格-库塔法有二阶公式和四

39、阶公式。库塔法有二阶公式和四阶公式。线性多步法有四阶阿达姆斯外插公式和内插公式。线性多步法有四阶阿达姆斯外插公式和内插公式。第53页，共75页，编辑于2022年，星期六（三）用（三）用Matlab软件求常微分方程的数值解软件求常微分方程的数值解t，x=solver（f,ts,x0,options）ode45 ode23 ode113ode15sode23s由待解方程写成的m-文件名ts=t0，tf，t0、tf为自变量的初值和终值函数的初值ode23：组合的组合的2/3阶龙格阶龙格-库塔库塔-芬尔格算法芬尔格算法ode45：运用组合的运用组合的4/5阶龙格阶龙格-库塔库塔-芬尔格算芬尔格算法法自

40、变量值函数值用于设定误差限用于设定误差限(缺省时设定相对误差缺省时设定相对误差10-3,绝对误差绝对误差10-6),命令为：命令为：options=odeset（reltol,rt,abstol,at）,rt，at：分别为设定的相对误差和绝对误差：分别为设定的相对误差和绝对误差.第54页，共75页，编辑于2022年，星期六 1、在解、在解n个未知函数的方程组时，个未知函数的方程组时，x0和和x均为均为n维向量，维向量，m-文件文件中的待解方程组应以中的待解方程组应以x的分量形式写成的分量形式写成.2、使用、使用Matlab软件求数值解时，高阶微分方程必须等价地软件求数值解时，高阶微分方程必须等

41、价地变换成一阶微分方程组变换成一阶微分方程组.注意注意:第55页，共75页，编辑于2022年，星期六解解:令 y1=x，y2=y11、建立建立m-文件文件vdp1000.m如下如下：function dy=vdp1000(t,y)dy=zeros(2,1);dy(1)=y(2);dy(2)=1000*(1-y(1)2)*y(2)-y(1);2、取t0=0，tf=3000，输入命令：T,Y=ode15s(vdp1000,0 3000,2 0);plot(T,Y(:,1),-)3、结果如图第56页，共75页，编辑于2022年，星期六解解 1、建立、建立m-文件文件rigid.m如下：如下：func

42、tion dy=rigid(t,y)dy=zeros(3,1);dy(1)=y(2)*y(3);dy(2)=-y(1)*y(3);dy(3)=-0.51*y(1)*y(2);2、取取t0=0，tf=12，输入命令：，输入命令：T,Y=ode45(rigid,0 12,0 1 1);plot(T,Y(:,1),-,T,Y(:,2),*,T,Y(:,3),+)3、结果如图、结果如图图中，图中，y1的图形为实线，的图形为实线，y2的图形为的图形为“*”线，线，y3的图形为的图形为“+”线线.第57页，共75页，编辑于2022年，星期六稳定性问题稳定性问题 在研究许多实际问题时，人们最为关心的也许并非

43、系统与时在研究许多实际问题时，人们最为关心的也许并非系统与时间有关的变化状态，而是系统最终的发展趋势。例如，在研究某间有关的变化状态，而是系统最终的发展趋势。例如，在研究某频危种群时，虽然我们也想了解它当前或今后的数量，但我们更频危种群时，虽然我们也想了解它当前或今后的数量，但我们更为关心的却是它最终是否会绝灭，用什么办法可以拯救这一种群，为关心的却是它最终是否会绝灭，用什么办法可以拯救这一种群，使之免于绝种等等问题。要解决这类问题，需要用到微分方程或使之免于绝种等等问题。要解决这类问题，需要用到微分方程或微分方程组的稳定性理论。在下两节，我们将研究几个与稳定性微分方程组的稳定性理论。在下两节

44、，我们将研究几个与稳定性有关的问题。有关的问题。第58页，共75页，编辑于2022年，星期六一般的微分方程或微分方程组可以写成：一般的微分方程或微分方程组可以写成：定义定义 称微分方程或微分方程组称微分方程或微分方程组 为自治系统或动力系统。为自治系统或动力系统。（3.28）若方程或方程组若方程或方程组f f(x x)=0)=0有解有解X Xo o，X=XX=Xo o显然满足（显然满足（3.28）。）。称点称点X Xo o为微分为微分方程或微分方程组（方程或微分方程组（3.28)3.28)的平衡点或奇点。的平衡点或奇点。第59页，共75页，编辑于2022年，星期六例例7 7 Logistic

45、Logistic模型模型 共有两个平衡点：共有两个平衡点：N=0和和N=K，分别对应微分方程的两两个，分别对应微分方程的两两个特殊解。前者为特殊解。前者为No=0时的解而后者为时的解而后者为No=K时的解。时的解。当当NoK时，则位于时，则位于N=K的上方。从图的上方。从图3-17中不难看出，若中不难看出，若No0，积分曲线在，积分曲线在N轴上的投轴上的投影曲线（称为轨线）将趋于影曲线（称为轨线）将趋于K。这说明，平衡点。这说明，平衡点N=0和和N=K有着极大的区有着极大的区别。别。图3-17 定义定义1 1 自治系统自治系统 的相空间是指以（的相空间是指以（x1,xn）为坐标）为坐标 的空间

46、的空间Rn。特别，当特别，当n=2时，称相空间为相平面。时，称相空间为相平面。空间空间Rn的点集的点集(x1,xn)|xi=xi(t)满足满足(3.28)，i=1,n称为系统的轨线，所有轨线在相称为系统的轨线，所有轨线在相空间的分布图称为相图。空间的分布图称为相图。第60页，共75页，编辑于2022年，星期六定义定义2 2 设设x x0 0是（是（3.283.28）的平衡点，称：）的平衡点，称：（1 1）x x0 0是稳定的，如果对于任意的是稳定的，如果对于任意的00，存在一个，存在一个00，只要，只要|x x(0)-(0)-x x0 0|，就有，就有|x x(t t)-)-x x0 0|对所

47、有的对所有的t t都成立。都成立。（2 2）x x0 0是渐近稳定的，如果它是稳定的且是渐近稳定的，如果它是稳定的且 。微分方程平衡点的稳定性除了几何方法，还可以通过解微分方程平衡点的稳定性除了几何方法，还可以通过解析方法来讨论，所用工具为以下一些定理。析方法来讨论，所用工具为以下一些定理。（3 3）x x0 0是不稳定的，如果（是不稳定的，如果（1 1）不成立。）不成立。根据这一定义，根据这一定义，LogisticLogistic方程方程的平衡点的平衡点N=KN=K是稳定的且为渐近是稳定的且为渐近稳定的，而平衡点稳定的，而平衡点N=0N=0则是不稳则是不稳定的。定的。第61页，共75页，编辑

48、于2022年，星期六解析方法解析方法定理定理1 1 设设x xo o是微分方程是微分方程 的平衡点：的平衡点：若若 ,则则x xo o是渐近稳定的是渐近稳定的若若 ,则则x xo o是渐近不稳定的是渐近不稳定的证证 由泰勒公式，当由泰勒公式，当x x与与x xo o充分接近时，有：充分接近时，有：由于由于x xo o是平衡点，故是平衡点，故f f(x xo o)=0)=0。若。若 ，则当，则当x x 0)0，从而，从而x x单增；当单增；当x x x xo o时，又有时，又有f f(x x)0)00，可能出现以下情形：，可能出现以下情形：若若q0，120。当当p0时，时，零点不稳定；零点不稳定

49、；当当p p00时，零点稳定时，零点稳定 若若q q00，1200时，零点不时，零点不 稳定稳定 当当p0时，零点稳定时，零点稳定第64页，共75页，编辑于2022年，星期六（2）0，此时此时若若a0,零点不稳定零点不稳定若若a=0，有零点为中心的周期解，有零点为中心的周期解 综上所述：仅当综上所述：仅当p0时，时，（3.30）零点才是渐近稳定的；当）零点才是渐近稳定的；当p=0且且q0时（时（3.30）有周期解，零点是稳定的中心（非渐近稳定）；在其他）有周期解，零点是稳定的中心（非渐近稳定）；在其他情况下，零点均为不稳定的情况下，零点均为不稳定的。非线性方程组（非线性方程组（3.29）平衡点

50、稳定性讨论可以证明有下面定）平衡点稳定性讨论可以证明有下面定理成立理成立:定理定理2 若（若（3.30）的零点是渐近稳定的，则（）的零点是渐近稳定的，则（3.29）的平衡点）的平衡点 也是渐近稳定的；若（也是渐近稳定的；若（3.30）的零点是不稳定的，则（）的零点是不稳定的，则（3.29）的平衡点也是不稳定的。的平衡点也是不稳定的。第65页，共75页，编辑于2022年，星期六捕食系统的捕食系统的VolterraVolterra方程方程 问题背景：问题背景：意大利生物学家意大利生物学家DAncona曾致力于鱼类种群相互制约关系的研曾致力于鱼类种群相互制约关系的研究，在研究过程中他无意中发现了一些