第3讲 微分方程模型

hold on和hold off
hold on 作用:当前轴及图形保持而不被刷新，准备 接受此后的绘制 例如，在作图时，想在一张图上同时显示 多组数据(便于观察等等原因)的时候用hold on (就是等一下,我还想在本图上画个东西) hold off 自然就是取消这个功能
指数增长模型的应用及局限性
[2] 建立模型
当 x xm 时，增长率应为0
r ( xm ) 0
r s xm
代入
r x r sx
得：
x r x r 1 x m
(3)
将（3）式代入（1）得：
|
i 1
m
2 i
| ( ( x i ) y
i 1
m
i)
2
最小,此即称为最小二乘原理
用Matlab作非线性最小二乘拟和
Matlab的最优化工具箱中提供了求非线性最小二乘拟合的函数： lsqcurvefit，调用格式为： x = lsqcurvefit (‘fun’, x0, xdata, ydata,…); fun是一个事先建立的 定义函数F(x,xdata) 的 M-文件, 自变量为x和 xdata
再编写一个M文件

clear all; close all; t=0:1:20; x=[3.9 5.3 7.2 9.6 12.9 17.1 23.2 31.4 38.6 50.2 62.9 76.0 92.0 106.5 123.2 131.7 150.7 179.3 204.0 226.5 251.4]; plot(t,x,'o'); hold on; c(1)=3.9; c(2)=mean(diff(x)./diff(t)./x(1:20)); k=lsqcurvefit(@fun,c,t,x) tt=[21 22]; xx=fun(k,tt); plot(tt,xx,'r*'); tt=0:0.1:22; xx=fun(k,tt); plot(tt,xx,'r'); hold off
x0初始值自然应取t=0时的x的值3.9 r初始值取增长率的平均值 mean(diff(x)./diff(t)./x(1:20))
Matlab求解过程
首先编写M文件fun.m，其中参数c(1)表示 x0,c(2)表示r。 function x=fun(c,t) x=c(1)*exp(c(2)*t);
杂的，建模时在不同的假设下去模拟实际的现
象，建立能近似反映问题的微分方程，然后从 数学上求解或分析所建方程及其解的性质，再 去同实际情况对比，检验此模型能否刻画、模 拟某些实际现象。
3.1 年代鉴定问题
在巴基斯坦一个洞穴里，发现了具有古代某种动 物特征的骨碎片，科学家把它带到实验室，作碳 14年代测定，分析表明，
y0 e
t 8000
y ke
.
当y 0.0624y0 时
求得 t 8000ln 0.0624 22194
此即所求死亡年数。
用Matlab求解对数
log(x) x的自然对数 log2(x) x的以2为底的对数 log10(x) x的以10为底的对数 以a为底x的对数，根据换底公式可表示为： log(x)/log(a)
建立微分方程模型的方法
（1）根据规律列方程
利用数学、力学、物理、化学等学科中的定理或 经过实验检验的规律等来建立微分方程模型。 （2）微元分析法 利用已知的定理与规律寻找微元之间的关系 式，与第一种方法不同的是对微元而不是直 接对函数及其导数应用规律。
（3）模拟近似法
在生物、经济等学科的实际问题中，许多现象 的规律性不很清楚，即使有所了解也是极其复
c j 103 4.54 4.99 5.35 5.65 5.90 6.10 6.26 6.39 6.50 6.59
该问题即解最优化问题：
min F (a, b, k ) [a be
j 1 10 0.02 kt j
c j ]2
解法
1）编写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:1000 cdata=1e-03*[4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39, 6.50,6.59]; x0=[0.2,0.05,0.05]; x=lsqcurvefit ('curvefun1',x0,tdata,cdata)
x 0 x0
[3] 模型求解
解微分方程得 x t x ert 0 t

(2)
x0 (1 r )
(2)被称为指数增长模型——马尔萨斯提出 (1798)，表示 随着时间增加，人口按指数规律无限增长
x0 (1 r )t 被称为等比数列模型
[4] 模型的参数估计
第3讲 微分方程建模方法及案例
在研究实际问题时，常常会联系到某些变 量的变化率或导数，这样所得到变量之间 的关系式就是微分方程模型。 微分方程模型反映的是变量之间的间接关 系，因此，要得到直接关系，就得求微分 方程。
求解微分方程的方法
求解微分方程有三种方法 1）求精确解； 2）求数值解（近似解）； 3）定性理论方法。
间的年代：
真正的年代＝
c
14
年 1.4 900
3.2
人口演化模型
预备知识:最小二乘原理
当由实验提供了大量数据时,不能要求拟合函数 ( x ) 在数据点( xi , yi ) 处的偏差,即 i ( xi ) y i (i=1,2,…,m) 严格为零,但为了使近似曲线尽量反映所给数据点 的变化趋势 ,需对偏差有所要求.通常要求偏差平方和
14
c
12
的数量成定比，
生物体死亡后，交换过程就停止了，放射性碳便以
设 t 为死后年数，
r (t ) xc14 (t ) xc12
14 12
则t 0时, r r0 , 即活体中 c 与 c 数量的比例 .
dxc14 dt
积分得

xc14 8000
t 8000
dy y dt 8000
迭代初值
已知数据点
注意: 1.如无合理初值，那就只能给出一个猜想初值 2.拟合结果是初值敏感的，因为找到的不一定是全局最优 而可能是初值附近的局部最优.
例 用下面一组数据拟合 中的参数a，b，k
tj
c(t ) a be0.02 kt
100 200 300 400 500 600 700 800 900 1000
从整体来说，人口的变化由两个因素决定：出生 和死亡。出生使得人口增加，死亡使得人口减少。 对于局部地区来说，除了出生和死亡外，影响人 口的变化还有两个因素：迁入和迁出。迁入使得 局部人口增加，迁出使得局部人口减少。 在迁入、迁出人口的差别不大时，人口的变化主 要由出生率和死亡率决定。 根据上面的分析，不难建立起人口演化模型。
下表是近两百年的美国人口统计数据，试依此建立美国人 口增长的数学模型，最后用它预报2000年、2010年美国 人口.
年（公元） 人口（百万） 年（公元） 人口（百万） 年（公元） 人口（百万）
1790 3.9 1860 31.4 1930 123.2
1800 5.3 1870 38.6 1940 131.7
wenku.baidu.com题的提出
人口的增长情况是当前世界引起普遍关注的问题， 早在18世纪人们就开始进行人口预报工作了。 20世纪90年代，我们经常可以在报刊上看见关于 人口增长的预报，说到20世纪末，或21世纪中叶， 全世界（或某地区）的人口将达到多少多少亿。 这些人口预报的数值是从哪里来的？准确不准确？ 你能不能对某地人口数目的演化进行一下估算？
• 与19世纪以前欧洲一些地区人口统计数据吻合 • 适用于19世纪后迁往加拿大的欧洲移民后代
• 可用于短期人口增长预测
• 不符合19世纪后多数地区人口增长规律 • 不能预测较长期的人口增长过程 19世纪后人口数据
人口增长率r不是常数(逐渐下降)
美国实际人口与按指数增长模型计算 的人口比较
年 实际人口 （百万） 指数增长模型 预测人口（百万） 误差（%） 1790 1800 1810 1820 1830 1840 1850 1860 1870 1880 1890 3.9 5.3 7.2 9.6 12.9 17.1 23.2 31.4 38.6 50.2 62.9 7.3 10.0 13.7 18.7 25.6 35.0 47.8 65.5 89.6 1.4 4.2 6.2 9.4 10.3 10.8 23.8 30.5 42.4
Matlab中的求解方法
-8000*log(0.0624) ans = 2.2194e+004
c
14
年代测定的修订
1966年，耶鲁实验室的Minze Stuiver和加利福尼亚 大学圣地亚哥分校的HansE.Suess在一份报告中指出：在 2500到10000年前这段时间中测得的结果有差异，其根本 原因在于那个年代，宇宙射线的放射性强度减弱了，偏差的 峰值发生在大约6000年以前。他们提出了一个很成功的误 差公式，用来校正根据碳测定出的2300年到6000年前这期
求解步骤
[1] 假设 [2] 建立模型 [3] 模型求解 [4] 模型的参数估计
[1] 假设
人口增长率r(即出生率b-死亡率d)是常数 （或单位时间内人口的增长量与当时的人 口成正比）.
[2] 建立模型
记时刻t=0时人口数为x0, 时刻t的人口为x(t),由于 量大,可视为连续、可微函数. t到t t时间内人口的增量为：xt t xt rxt t 于是x(t)满足微分方程： dx rx dt (1)
c
14
与
c
12
的比例仅仅
是活组织内的6.24%，能否判断此动物生活在多
少年前？
背景
c
14
年代测定：活体中的碳有一小部分是放射性
同位素
c
14
，这种放射性碳是由于宇宙射线在高层
大气中的撞击引起的，经过一系列交换过程进入活
组织内，直到在生物体内达到平衡浓度，这意味着
在活体中， c 的数量与稳定的 每年八千分之一的速度减少。
[1] 假设
0 最简单假定 rx r sx, r, s （线性函数）， r 为固有增长率(x很小时). 1838年,荷兰生物数学家韦尔侯斯特 (Verhulst)引入常数 ,用来表示自然环境 条件所能容许的最大人口数（一般说来,一 个国家工业化程度越高,它的生活空间就越 大,食物就越多,从而 x m 就越大）
要用模型的结果来预报人口，必须对其中 的参数r进行估计，这可以用表中所给的数 据通过拟合得到.
拟合的关键问题
1. t的变换
指数函数exp(t)当t很大时可能会溢出,为了减小 数据误差,首先将时间域变换至[0,20],所用的变 换为: t=1800+(t-1800)/10
2. x0和r初值的确定
1900
1910 1920
76.0
92.0 106.5
122.5
167.6 229.3
61.2
82.1 115.3
阻滞增长模型(Logistic模型)
马尔萨斯模型为什么不能预测未来的人口呢？ 这主要是地球上的各种资源只能供一定数量的人 生活,随着人口的增加,自然资源、环境条件等因素 对人口增长的阻滞作用越来越显著 如果当人口较少时,人口的自然增长率可以看作常 数的话,那么当人口增加到一定数量以后,这个增长 r是x的减函数 率就要随人口的增加而减小. 因此,应对马尔萨斯模型中关于净增长率为常数的 假设进行修改.
1810 7.2 1880 50.2 1950 150.7
1820 9.6 1890 62.9 1960 179.3
1830 12.9 1900 76.0 1970 204.0
1840 17.1 1910 92.0 1980 226.5
1850 23.2 1920 106.5 1990 251.4
问题的分析