微分方程建模课件

微分方程建模,1,卢长娜changnalu,当我们描述实际对象的某些特性随时间（空间）而演变的过程、分析它的变化规律、预测它的未来形态、研究它的控制手段时。通常要建立对象的动态模型。,在许多实际问题中，当直接导出变量之间的函数关系较为困难，但导出包含未知函数的导数或微分的关系式较为容易时，可用建立微分方程模型的方法来研究该问题。,2,一、微分方程模型,涉及“改变”、“变化”、“增加”、“减少”、“衰变”、“边际”、“速度”、“运动”、“追赶”、“逃跑”、等等词语的确定性连续问题。,b、微分方程建模的基本手段微元法等,a、微分方程建模的对象,3,1、寻找改变量一般说来微分方程问题都遵循这样的文字等式变化率（微商）=单位增加量-单位减少量等式通常是利用已有的原则或定律。,c、微分方程建模的基本规则,2、对问题中的特征进行数学刻画3、用微元法建立微分方程；4、确定微分方程的定解条件（初边值条件）；5、求解或讨论方程（数值解或定性理论）6、模型和结果的讨论与分析。,4,论文格式及组成,题目摘要，关键词问题重述,模型假设分析与建立数学模型数学模型的求解模型检验(总结与评价)模型推广参考文献附录（若有）,正文,二、导弹跟踪问题,1、实验目的本实验主要涉及常微分方程。通过实验复习微分方程的建模和求解；介绍两种微分方程的数值方法：Euler法和改进的Euler法；还介绍了仿真方法。,6,3、数学模型微分方程建模的方法主要是依据守恒律来建立等量关系式。对于这个问题，寻求等量关系是比较简单的。设坐标系如图3.1所示，取导弹基地为原点O(0,0)，x轴指向正东方,y轴指向正北方。,2、实际问题某军的一导弹基地发现正北方向120km处海面上有敌舰一艘以90km/h的速度向正东方向行驶。该基地立即发射导弹跟踪追击敌艇，导弹速度为450km/h，自动导航系统使导弹在任一时刻都能对准敌艇。试问导弹在何时何处击中敌舰？,当t=0时，导弹位于点O，敌艇位于点A(0,H)，其中H=120（km）。,7,设导弹在t时刻的位置为P(x(t),y(t)，由题意,方程(3.1),(3.3)连同初值条件构成了一个关于时间变量t的一阶微分方程组的初值问题。,其中另外在t时刻,敌艇位置为,其中。由于导弹轨迹的切线方向必须指向敌艇，即直线PM的方向就是导弹轨迹上点P的切线方向,故有或写为,8,两边对t求导,为了寻求x与y的关系,要设法消去变量t,由式(3.2)得,即有把式(3.1)写为代入上式,就得到轨迹方程.这是一个二阶非线性微分方程,加上初值条件,则初值问题,9,就是导弹轨迹的数学模型。值得注意的是，前面导出的一阶微分方程组(3.1)，(3.3)和(3.4)实际上已经是一个数学模型了，不过多一个变量（或说参数）而已。,10,则（3.5）化为一阶可分离变量方程,即易得,由初值条件(3.7)即,4、解析方法方程（3.5）可以降阶.令,注意到上式可改写为于是有这样我们又得到一个可分离变量方程积分得利用于是导弹轨迹方程为,11,设导弹击中敌艇于B(L,H),以Y=H代入(3.9)式,得而导弹击中敌艇的时刻将数据代入(3.10),(3.11)式,得,12,5、数值方法将初值问题(3.5)(3.7)化为一阶微分方程组,取自变量y的步长为,于是得节点相应点上的x的值和p的值记为显然,有初值条件我们将介绍两种近似算法来进行数值处理.,13,设导弹到达处的时刻为那么得到计算的迭代格式,Euler方法Euler方法十分简单,就是用差商代替微商,即以代之以而代之以这样有,于是表3.1是取n=4时的计算结果,读者可以用来检验程序或应用软件的正确性.表3.1kykxkpk000013000.052601.50.123905.00.22412011.50.42此时,14,表3.2是对于不同的n值所对应的计算结果.显然,n越大(即h越小),结果就越精确.表3.2n4812244896120240L11.5215.9617.9720.5522.2523.3323.5824.15T0.1280.1770.2000.2280.2470.2590.2620.268,15,注意，由问题（3.1）,(3.3),(3.4)消去t推导出问题（3.5）(3.7)是较为巧妙和偶然的.一般而言，一个微分方程组未必能消去一些变元而减少方程的个数。那么，我们能否直接对初值问题(3.1),(3.3),(3.4)进行数值处理呢？答案是肯定的。实际上，只要由方程(3.1),(3.3)解出和的表达式，这样问题变为,取时间步长对应时导弹轨迹上点的坐标为则Euler格式为当计算到即停止,于是,16,表3.3和表3.4分别列出了取步长为0.1和0.05时的计算结果:表3.3ktkxkyk00.00.000000.0000010.10.0000045.0000020.25.3615489.6794630.322.67495131.21553此时取,17,表3.4ktkxkyk10.050.0000022.5000020.101.0373642.9760730.153.4120567.3504140.207.6461589.4484350.2514.86790110.7579660.3029.19480128.10702此时取表3.5是对应不同的,用Euler法所得相应的步长推进次数n和计算结果.表3.50.10.050.0050.001n3656278L22.6749529.1948025.6673125.04935T0.251940.324390.285190.27833,18,Euler方法较为简单,但也较为粗糙,对形式较复杂的微分方程更易有较大的误差.人们设计了不少更精确的近似算法,这里我们介绍其中的一种,进一步研究可看参考书.,19,.改进的Euler方法(预报-校正法)以一维情况为例,对问题Euler迭代格式是,由积分表达式的几何意义看,右边为下方的曲边梯形,从图3.2我们可以看出Euler法是用矩形来代替曲边梯形,而改进的Euler法则是用梯形来代替曲边梯形.对问题(3.18)(3.20),我们写出相应的改进Euler迭代格式,20,21,表3.6和表3.7分别列出了取步长为0.1和0.05时的计算结果:表3.6ktkxkyk00.00.000000.0000010.12.6807744.8397320.212.5752488.2867930.322.07242130.25569此时取表3.7ktkxkyk10.050.5186822.4880420.102.1059644.9219530.155.0982167.2021340.2010.1609689.0690650.2519.65646108.9789860.3024.24089130.99030,22,此时取表3.8是对应不同的,用改进的Euler法所得相应的步长推进次数n和计算结果.表3.80.10.050.0050.001n3656278L27.0724224.2408925.1355224.98112T0.300800.269340.279280.27757,图3.3画出了导弹轨迹由解析式所给出的精确曲线以及由Euler法和改进的Euler法进行数值计算所给出的近似曲线.,23,6、仿真方法如果建立微分方程很困难，或者微分方程很复杂而较难作出数值处理，常常可以用仿真方法。所谓仿真方法，顾名思义，指的是模仿真实行为和过程的方法。在这个具体问题中，就是一步步地模拟导弹追踪敌艇的实际过程。而计算机仿真，则是在计算机上通过相应的程序和软件来实现对事件运行的实际过程的模拟。设导弹和敌艇在初始时刻（即t=0时）分别位于P0(0,0)和M0(0,H)。此时，导弹指向M0。而在t=时，导弹的位置P1(x1,y1)，其中，敌艇的位置则为这时导弹沿P1M1方向飞行，P1M1的倾角为在t=2时，导弹的位置为P(x,y)，其中,24,此时敌艇位置为,导弹沿P2M2方向飞行(见图3.4).以此方式,一般地,设时,导弹位置为敌艇的位置则为导弹将沿PkMk方向飞行,那么,PkMk的倾角为,25,从而时,导弹位置为,其中而敌艇位置为仍然可以如前那样,当时,仿真停止;或者事先给定误差界,当时,仿真停止,这时对于我们用仿真迭代格式(3.33)(3.36)进行计算,结果与Euler迭代格式的结果完全一致(见表3.33.5).这两种迭代格式实际上确实是相同的,建议读者自己验证一下.值得注意的是,在仿真方法中,我们根本没有用到微分方程组(3.18)(3.20),却得到了它的一种离散形式,这是十分有意思的.,26,三、传染病模型,传染病是人类的大敌，通过疾病传播过程中若干重要因素之间的联系建立微分方程加以讨论，研究传染病流行的规律并找出控制疾病流行的方法显然是一件十分有意义的工作。在本节中，我们将主要用多房室系统的观点来看待传染病的流行，并建立起相应的多房室模型。,医生们发现，在一个民族或地区，当某种传染病流传时，波及到的总人数大体上保持为一个常数。即既非所有人都会得病也非毫无规律，两次流行（同种疾病）的波及人数不会相差太大。如何解释这一现象呢？试用建模方法来加以证明。,问题的提出：,设某地区共有n+1人，最初时刻共有io人得病，t时刻已感染（infective）的病人数为i(t)，假定每一已感染者在单位时间内将疾病传播给k个人（k称为该疾病的传染强度），且设此疾病既不导致死亡也不会康复,模型1,此模型即Malthus模型，它大体上反映了传染病流行初期的病人增长情况，在医学上有一定的参考价值，但随着时间的推移，将越来越偏离实际情况。,已感染者与尚未感染者之间存在着明显的区别，有必要将人群划分成已感染者与尚未感染的易感染，对每一类中的个体则不加任何区分，来建立两房室系统。,模型2,记t时刻的病人数与易感染人数（susceptible）分别为i(t)与s(t)，初始时刻的病人数为i。根据病人不死也不会康复的假设及（竞争项）统计筹算律，,其中：,统计结果显示，(3.17)预报结果比(3.15)更接近实际情况。医学上称曲线为传染病曲线，并称最大值时刻t1为此传染病的流行高峰。,模型2仍有不足之处，它无法解释医生们发现的现象，且当时间趋与无穷时，模型预测最终所有人都得病，与实际情况不符。,为了使模型更精确，有必要再将人群细分，建立多房室系统,（3.18）,求解过程如下：,对（3）式求导，由（1）、（2）得：,解得：,将人群划分为三类（见右图）：易感染者、已感染者和已恢复者（recovered）。分别记t时刻的三类人数为s(t)、i(t)和r(t)，则可建立下面的三房室模型：,模型3,由(1)式可得：,从而解得：,为揭示产生上述现象的原因（3.18）中的第（1）式改写成：,其中通常是一个与疾病种类有关的较大的常数。,下面对进行讨论，请参见右图,如果，则开始时，i(t)单增。但在i(t)增加的同时，伴随地有s(t)单减。当s(t)减少到小于等于时，i(t)开始减小，直至此疾病在该地区消失。,鉴于在本模型中的作用，被医生们称为此疾病在该地区的阀值。的引入解释了为什么此疾病没有波及到该地区的所有人。,图3-14,综上所述，模型3指出了传染病的以下特征：,（1）当人群中有人得了某种传染病时，此疾病并不一定流传，仅当易受感染的人数与超过阀值时，疾病才会流传起来。,（2）疾病并非因缺少易感染者而停止传播，相反，是因为缺少传播者才停止传播的，否则将导致所有人得病。,（3）种群不可能因为某种传染病而绝灭。,模型检验：,医疗机构一般依据r(t)来统计疾病的波及人数，从广义上理解，r(t)为t时刻已就医而被隔离的人数，是康复还是死亡对模型并无影响。,通常情况下，传染病波及的人数占总人数的百分比不会太大，故一般是小量。利用泰勒公式展开取前三项，有：,代入（3.20）得近似方程：,积分得：,其中：,这里双曲正切函数：,而：,图3-14（a）给出了（3.21）式曲线的图形，可用医疗单位每天实际登录数进行比较拟合得最优曲线。,图3-14（a）,图3-14（b）记录了1905年下半年至1906年上半年印度孟买瘟疫大流行期间每周死亡人数，不难看出两者有较好的一致性。,微分方程数值解,35,卢长娜changnalu,36,常微分方程(组)初值问题的数值解法,考虑一阶常微分方程的初值问题:,（1）（2）,要计算出解函数y(x)在一系列离散点a=x0x1xn=b处的近似值,节点间距为步长，通常采用等距节点，即取hi=h(常数)。,37,欧拉法的建立,差商方法,1欧拉方法,欧拉法的几何意义,隐式欧拉法,由于未知数yi+1同时出现在等式的两边，不能直接得到，故称为隐式欧拉公式，而前者称为显式欧拉公式。,39,局部截断误差：,欧拉法的局部截断误差：,欧拉法具有1阶精度。,隐式欧拉法的局部截断误差：,即隐式欧拉公式具有1阶精度。,40,改进的欧拉法,注：此法亦称为预测-校正法。可以证明该算法具有2阶精度，同时可以看到它是个单步递推格式，比隐式公式的迭代求解过程简单。后面将看到，它的稳定性高于显式欧拉法。,改进的欧拉法的“嵌套”表示形式,41,改进的欧拉法也可以表示为“平均化形式”：,xi,xi+1,P,Q,R,S,42,单步递推法的基本思想是从(xi,yi)点出发，以某一斜率沿直线达到(xi+1,yi+1)点。欧拉法及其各种变形所能达到的最高精度为2阶。,2龙格-库塔法,斜率一定取K1K2的平均值吗？步长一定是一个h吗？,43,一个比较简单而又重要的三阶龙格库塔公式，称为库塔公式，其为：,继续上述这个过程,可以进一步讨论四阶龙格库塔公式。,最常用的是4阶经典龙格-库塔法：,44,3微分方程组与高阶方程,一阶微分方程组,IVP的一般形式为：,前述所有公式皆适用于向量形式。,45,高阶微分方程,化作一阶微分方程组求解。,引入新变量,初值条件为：,46,2阶常微分方程边值问题,有限差分法,将求解区间a,b等分为N份，取节点xi=a+ih(i=0,N)，在每一个节点处将y和y离散化。,47,差分方法的基本思想“就是以差商代替微商”,考虑如下两个Taylor公式：,（1）,（2）,从（1）得到：,48,从（2）得到：,从（1）-（2）得到：,从（1）+（2）得到：,2020/6/7,NUIST,49,Programs,Example1:Solvetheinitialvalueproblem,2020/6/7,NUIST,50,2020/6/7,NUIST,51,2020/6/7,NUIST,52,2020/6/7,NUIST,53,Poincarmapsareimportanttoolsforstudyingthebehaviourofdynamicalsystems.Theygivesnapshotsofthestateofthedynamicalsystematpre-definedinstants.However,ratherthantakingsnapshotsatfixedintervalsoftime,valuesareoftensoughtforwhichalinearcombinationofsolutioncomponentsvanishes;e.g.,whenagivencomponentis0.,Weconsiderapairofcoupledharmonicoscillators,Where,Example2:,2020/6/7,NUIST,54,2020/6/7,NUIST,55,2020/6/7,NUIST,56,2020/6/7,NUIST,57,Theterm“butterflyeffect”itselfisrelatedtotheworkofEdwardLorenz,anditisbasedinchaostheoryandsensitivedependenceoninitialconditions.TheLorenzeequationisderivedfromthesimplifiedequationsofconvectionrollsa