2015年数学建模讲座4月20日.ppt

2015年数学建模讲座,微分方程模型讨论：基本思想、方法与数值模拟分析贺天兰,微分方程=自然定律,常微分方程：联系着自变量、未知函数、及未知函数的某些导数之间的关系式总是与变化率相联系：如速度、加速度、曲率、增长率、出生率、死亡率等,本讲座分两个部分,模型举例模型求解的基本理论,一、模型举例,体重变化模型(一阶线性常微分方程)减肥计划模型（差分方程）国家综合实力模型（一阶二维非线性微分方程组）,模型一、体重变化模型,某人食量10467J/Day,假设其中5038J/Day因新陈代谢而自动消耗，其它活动，如健身等所消耗的热量是69J/(kg.Day)与他体重的乘积。假设与脂肪形式存储的热量100%有效，且1kg脂肪含热量41868J.试研究此人体重随时间变化的规律。,1、问题分析：问题涉及的原则、理论或方法做适当分析,要研究此人体重随时间变化的规律，就是要找出体重随时间变化的函数关系式。注意各变量与它们的微小增量之间的关系。,2、问题假设：提出有利条件，对实际问题做理想化近似,(1)设t时刻某人的体重为W(t)，一天开始时此人的体重为W(0);(2)W(t)是t的光滑函数；(3)体重在某时间段内的变化等于输入与输出之差，其中，其中输入是指扣除基本新陈代谢之后的净食量吸收，输出是指如健身时的消耗；(4)只考虑一天的情况，且不会出现异常现象。,3.模型建立:实际问题-数学语言:函数方程方程组等,此人每天的体重变化=输入-输出输入热量=10467-5038=5429J输出热量=69W=69WJ,4、模型求解,对数学翻译进行求解，既避免了对原问题的直接处理，又利用了现有的数学知识。,体重变化曲线,ClearW0=60;t=0:3000;W=5429/69-(5429-69*w0)/69*exp(-23*t/12956);Plot(t,w,r),grid,legend(体重变化曲线)；,5、结果结论：将数学结果返回到实际中做相应的解释,给出明确可信的结论！,6、敏感性分析：讨论参数的变化对结果的影响,讨论：增加或减少热量的摄入，体重将如何变化？增加或减少健身活动，体重又将如何变化？那个参数的变化对体重比较敏感？,定性分析,模型二、减肥计划-节食与运动,背景,多数减肥食品达不到减肥目标，或不能维持,通过控制饮食和适当的运动，在不伤害身体的前提下，达到减轻体重并维持下去的目标,分析,体重变化由体内能量守恒破坏引起,饮食（吸收热量）引起体重增加,代谢和运动（消耗热量）引起体重减少,体重指数BMI=w(kg)/l2(m2).18.525超重;BMI30肥胖.,模型假设,1）体重增加正比于吸收的热量每8000千卡增加体重1千克；,2）代谢引起的体重减少正比于体重每周每公斤体重消耗200千卡320千卡(因人而异),相当于70千克的人每天消耗2000千卡3200千卡；,3）运动引起的体重减少正比于体重，且与运动形式有关；,4）为了安全与健康，每周体重减少不宜超过1.5千克，每周吸收热量不要小于10000千卡。,某甲体重100千克，目前每周吸收20000千卡热量，体重维持不变。现欲减肥至75千克。,第一阶段：每周减肥1千克，每周吸收热量逐渐减少，直至达到下限（10000千卡）；,第二阶段：每周吸收热量保持下限，减肥达到目标,2）若要加快进程，第二阶段增加运动，试安排计划。,1）在不运动的情况下安排一个两阶段计划。,减肥计划,3）给出达到目标后维持体重的方案。,确定某甲的代谢消耗系数,即每周每千克体重消耗20000/100=200千卡,基本模型,w(k)第k周(末)体重,c(k)第k周吸收热量,代谢消耗系数(因人而异),1）不运动情况的两阶段减肥计划,每周吸收20000千卡w=100千克不变,第一阶段:w(k)每周减1千克,c(k)减至下限10000千卡,第一阶段10周,每周减1千克，第10周末体重90千克,吸收热量为,1）不运动情况的两阶段减肥计划,第二阶段：每周c(k)保持Cm,w(k)减至75千克,1）不运动情况的两阶段减肥计划,基本模型,第二阶段：每周c(k)保持Cm,w(k)减至75千克,第二阶段19周,每周吸收热量保持10000千卡,体重按减少至75千克。,运动t=24(每周跳舞8小时或自行车10小时),14周即可。,2）第二阶段增加运动的减肥计划,t每周运动时间(小时),3）达到目标体重75千克后维持不变的方案,每周吸收热量c(k)保持某常数C，使体重w不变,不运动,运动(内容同前),模型三,综合实力的微分方程模型讨论,一个地区物质文明和精神文明的综合称为综合实力,物质文明的发展速度与现有水平及发展潜力之积成正比，还与精神文明的水平成比例，精神文明的发展速度与现有水平成比例，且与现有的物质文明相关，有如下微分方程模型：,综合实力的微分方程模型：,对综合实力的微分方程模型进行尺度变换得：,在以下三种参数条件下分析讨论模型奇点的类型并对比奇点处线性化前后相图：,（1）,（2）,（3）,（1）,模型为：,奇点为：（0,0）稳定焦点和（0.8,0.16）鞍点,（1）,相图,奇点(0,0)和(0.8,0.16)处线性化后模型相图,我们可以观察到，线性化之前模型奇点附近的渐近行为与对该模型在奇点处线性化后的两个模型在零点附近的渐近行为分别相同！,（2）,模型为：,奇点为：（0,0）和（0.875,0.2187）,（2）,相图,奇点(0,0)和（0.875,0.2187）处线性化后模型相图,我们可以观察到，线性化之前模型奇点附近的渐近行为与对该模型在奇点处线性化后的两个模型在零点附近的渐近行为分别相似！,（3）,模型为：,奇点为：（0,0）和（0.8,0.2）,（3）,相图,奇点(0,0)和(0.8,0.2)处线性化后模型相图,我们可以观察到，线性化之前模型奇点附近的渐近行为与对该模型在奇点处线性化后的两个模型在零点附近的渐近行为分别相同！,二、相关的理论知识,微分方程模型讨论：基本思想、方法与数值模拟分析贺天兰,基本概念：1、微分方程的解、通解,例：y=sin(arcsinx+c),(-1,1),2、初值问题：柯西问题与特解,3、常微分方程的阶、线性、非线性常微分方程,4、问题：哪些常微分方程可以用初等方法求解？,自治的常微分方程与非自治的常微分方程,可以用初等方法求解的方程（组）,可分离变量的一阶方程及求解,一阶线性微分方程求解：常数变易法,非齐次,一阶线性齐次微分方程,（1）,（2）,Firstorderlinearnon-homogeneousdifferentialequation；Variationofconstants,全微分方程(Total/complete/filly)DE,M(x,y)dx+N(x,y)dy=0是全微分方程的充要条件是,一阶二维线性方程组的一般理论,二阶常系数线性微分方程,二阶线性常系数非齐次微分方程,二阶线性常系数非齐次微分方程解的结构,齐次方程的通解,非齐次方程的一个特解,5、黎卡提方程可以求解吗？,6、一阶方程解的存在性、唯一性,f(x,y)在R0上连续,7、定性方法与数值方法,1、一阶微分方程的几何意义；2、斜率场两种特例；解析方法与定性方法相结合；举例3、数值方法：欧拉方法,一阶的自治与非自治方程,举例,相线：自治方程简化的斜率场,利用相线画解的图像简图,3,0,-/2,/2,/2,3,例画出相线，并研究解的趋势,平衡点的分类：汇、源、结点,具有Allee效应的Logistic模型,分歧,自然条件下海洋中鱼的总量：,单参数微分方程族,8、一阶二维非线性方程组模型的定性、数值讨论,模型举例,Lotka-Volterra竞争模型,解？,处理非线性自治常微分方程模型的定性方法,在平衡点处线性化，借助线性模型而得到结论,平衡解、线性化定理,（x0,y0）是平衡点,Jacobimatrix,原系统与线性化系统在平衡点处,特征值决定,线性系统与非线性系统在平衡点附近差不多,functiondy=ODEfun6(t,y)dy=zeros(2,1);dy(1)=-3*y(1)+y(2);dy(2)=-y(2);returncleart1,Y1=ode45(ODEfun6,-8:0.1:8,1;0);t2,Y2=ode45(ODEfun6,-8:0.1:8,1;3);t3,Y3=ode45(ODEfun6,-8:0.1:8,-1;0);subplot(3,1,1),plot(t1,Y1,-g);subplot(3,1,2),plot(t2,Y2,-b);subplot(3,1,3),plot(t3,Y3,-r);figure;subplot(3,1,1),plot(Y1(:,1),Y1(:,2),-g);subplot(3,1,2),plot(Y2(:,1),Y2(:,2),-b);subplot(3,1,3),plot(Y3(:,1),Y3(:,2),-r);figure;plot(Y1(:,1),Y1(:,2),-g);holdon;plot(Y2(:,1),Y2(:,2),-b);plot(Y3(:,1),Y3(:,2),-r);,非线性模型,functiondy=ODEfun7(t,y)dy=zeros(2,1);dy(1)=-3*y(1)+y(2);dy(2)=-y(2)+y(1)2;returncleart1,Y1=ode45(ODEfun7,-8:0.1:8,1;0);t2,Y2=ode45(ODEfun7,-8:0.1:8,1;3);t3,Y3=ode45(ODEfun7,-8:0.1:8,-1;0);numsol7=t1,Y1,t2,Y2,t3,Y3;plot(Y1(:,1),Y1(:,2),-g);holdon;plot(Y2(:,1),Y2(:,2),-b);plot(Y3(:,1),Y3(:,2),-r);,Anexampleofalimitcycle:VandelPolequation,程序,M-functiondy=ODEfun5(t,y)mu=0.1;dy=zeros(2,1);dy(1)=y(2)-mu*(y(1)3/3-y(1);dy(2)=-y(1);returncleart1,Y1=ode45(ODEfun5,-8:0.08:16,-2;-2);t2,Y2=ode45(ODEfun5,-8:0.08:16,0.7;-0.7);t3,Y3=ode45(ODEfun5,-8:0.08:16,2;0);numsol5=t1,Y1,t2,Y2,t3,Y3subplot(3,1,1),plot(t1,Y1,-g);subplot(3,1,2),plot(t2,Y2,-b);subplot(3,1,3),plot(t3,Y3,-r);figure;subplot(3,1,1),plot(Y1(:,1),Y1(:,2),-g);subplot(3,1,2),plot(Y2(:,1),Y2(:,2),-b);subplot(3,1,3),plot(Y3(:,1),Y3(:,2),-r);figure;plot(Y1(:,1),Y1(:,2),-g);holdon;plot(Y2(:,1),Y2(:,2),-b);plot(Y3(:,1),Y3(:,2),-r);,Lorenzequations,chengxu,functiondx=lorenzf(t,x,flag,sigma,r,b);dx=-sigma*(x(1)-x(2);r*x(1)-x(2)-x(1)*x(3);-b*x(3)+x(1)*x(2);returnclearx0=1,0,0;sigma=9;r=27;b=3;t,x=ode45(lorenzf,0,60,x0,sigma,r,b);subplot(221);plot(t,x(:,1);L(1)=xlabel(itt);L(2)=ylabel(itx);subplot(222);plot(t,x(:,2);L(3)=xlabel(itt);L(4)=ylabel(ity);subplot(223);plot(t,x(:,3);L(1)=xlabel(itt);L(2)=ylabel(itz);subplot(223);plot(t,x(:,3);L(5)=xlabel(itt);L(6)=ylabel(itz);t(1:1500)=;x(1:1500,:)=;subplot(224);plot(x(:,1),x(:,2),x(:,3);boxon;?Errorusing=plotNotenoughinputarguments.subplot(224);plot3(x(:,1),x(:,2),x(:,3);boxon;subplot(224);comet3(x(:,1),x(:,2),x(:,3);boxon;subplot(224);plot3(x(:,1),x(:,2),x(:,3);boxon;,几个特殊的轨道,周期轨道、同宿轨道、异宿轨道平衡点,注意：,数学模型中，等式的两端必须有相同的量纲-物理量的属性,差分方程的模型.求解.解的性质简介,差