常微分方程研究什么
前面学方程时,我们常常在找一个数。比如 x 2 = 4 x^2=4 x 2 = 4 ,找到 x = 2 x=2 x = 2 和 x = − 2 x=-2 x = − 2 ,问题就解决了。到了常微分方程,这件熟悉的事会发生一次变化:要找的未知数,变成了一整个函数。
先想象一个培养皿。我们暂时不知道里面的细菌数量每天会是多少,却可以提出一个变化规律:细菌越多,同样一小段时间里新增的细菌也越多。如果把数量记成 P ( t ) P(t) P ( t ) ,这句话就把“现在有多少”和“现在增长多快”联系起来了。微分方程要做的,正是从这样的联系出发,找出数量随时间变化的函数。
咖啡冷却也有类似的关系,只是决定降温快慢的变成了温差;弹簧振动更复杂一些,当前位置和速度会一起影响加速度。我们会反复在现象、方程和解之间来回走。先把现象说清楚,方程才有来处;把解放回现象中检查,计算才有落点。
从求一个数,到找一整条曲线
不妨先离开具体单位,看看这个最简单的例子:
y ′ = 2 y . y'=2y. y ′ = 2 y .
这里的 y y y 是未知函数 y ( t ) y(t) y ( t ) 的简写,y ′ y' y ′ 是它对 t t t 的导数。把式子读成一句话,就是:“在每一个时刻,函数的变化率都等于它当前数值的两倍。”
如果有人回答“y = 3 y=3 y = 3 ”,我们需要追问:是说某个时刻函数值为 3 3 3 ,还是说函数永远等于 3 3 3 ?前一种回答只给了曲线上的一个点,远远不够;后一种回答指的是常函数 y ( t ) = 3 y(t)=3 y ( t ) = 3 ,它的导数为 0 0 0 ,而右边 2 y = 6 2y=6 2 y = 6 ,并不满足方程。
那么,什么函数有希望?从微积分里,我们知道指数函数求导后仍然保持指数函数的样子。试一试 y ( t ) = e 2 t y(t)=e^{2t} y ( t ) = e 2 t :
y ′ ( t ) = 2 e 2 t = 2 y ( t ) . y'(t)=2e^{2t}=2y(t). y ′ ( t ) = 2 e 2 t = 2 y ( t ) .
等式对每个实数 t t t 都成立,所以它确实是一个解。再试 y ( t ) = 3 e 2 t y(t)=3e^{2t} y ( t ) = 3 e 2 t ,求导得到 6 e 2 t 6e^{2t} 6 e 2 t ,恰好也等于自身的两倍。实际上,对任意实常数 C C C ,
y ( t ) = C e 2 t y(t)=Ce^{2t} y ( t ) = C e 2 t
都满足这个方程。C = 0 C=0 C = 0 给出一直贴在横轴上的零函数;C > 0 C>0 C > 0 给出向上增长的曲线;C < 0 C<0 C < 0 给出数值越来越负的曲线。数学方程允许这些解,若它具体描述人口,我们还要根据人口不能为负这一点限制初值。
请停下来比较一下两种求解。代数方程 x 2 = 4 x^2=4 x 2 = 4 要求某个数满足关系;这里要求函数和它的导数在一整段时间里始终配合 。所以微分方程的答案通常要写成 y ( t ) = ⋯ y(t)=\cdots y ( t ) = ⋯ ,而不能只交出一个数。
导数已经告诉了我们什么
导数的含义并没有因为课程换了名字而改变。若一个量在时刻 t t t 的值是 y ( t ) y(t) y ( t ) ,经过一小段时间 h h h ,可用
y ( t + h ) ≈ y ( t ) + y ′ ( t ) h y(t+h)\approx y(t)+y'(t)h y ( t + h ) ≈ y ( t ) + y ′ ( t ) h
估计新值。这里的近似来自可微性,只适合足够小的时间步长。对于 y ′ = 2 y y'=2y y ′ = 2 y ,如果当前 y = 3 y=3 y = 3 ,当前斜率就是 6 6 6 ;取一个很小的正步长,曲线会从这个位置向上走。
注意这并不意味着以后一直以斜率 6 6 6 增长。函数值改变后,右边的 2 y 2y 2 y 也会改变,下一刻的斜率随之更新。这种“状态影响变化率,变化率又改变状态”的关系,是许多微分方程模型的共同结构。
更一般的一阶方程常写成
d y d t = f ( t , y ) . \frac{dy}{dt}=f(t,y). d t d y = f ( t , y ) .
右边的 f f f 是已知规则。它可以只看时间,如 y ′ = cos t y'=\cos t y ′ = cos t ;可以只看当前状态,如 y ′ = 2 y y'=2y y ′ = 2 y ;也可以两者都看,如 y ′ = t − y y'=t-y y ′ = t − y 。求解时未知的是沿着时间走下去的函数 y ( t ) y(t) y ( t ) ,并不是右边这个规则。
下面的交互把几种变化放在一起。移动时间探针时,可以先看当前函数值,再看切线斜率:数值高,不一定意味着变化快;变化快慢最终要由各自的方程决定。
判断一个函数是不是解,要把它放回去
现在即使还不会系统求解,也已经能做一件很有用的事:检查别人给出的答案。方法只有一个核心动作——求出需要的导数,代回原方程,看等式是否在所讨论的区间内处处成立。
一个算完的验证例子
假设要检查函数
y ( t ) = 2 + 5 e − 3 t y(t)=2+5e^{-3t} y ( t ) = 2 + 5 e − 3 t
是否满足方程
y ′ + 3 y = 6. y'+3y=6. y ′ + 3 y = 6.
先求导。常数 2 2 2 的导数是 0 0 0 ,指数项求导要乘上指数中 t t t 的系数 − 3 -3 − 3 ,所以
y ′ = − 15 e − 3 t . y'=-15e^{-3t}. y ′ = − 15 e − 3 t .
再把函数和导数一起代入左边:
y ′ + 3 y = − 15 e − 3 t + 3 ( 2 + 5 e − 3 t ) = 6. y'+3y=-15e^{-3t}+3(2+5e^{-3t})=6. y ′ + 3 y = − 15 e − 3 t + 3 ( 2 + 5 e − 3 t ) = 6.
左边确实等于右边,而且没有限制 t t t 的分母、根号或对数,所以这个函数在整个实数轴上都是方程的解。它还有一个容易看懂的趋势:随着 t t t 增大,指数项逐渐变小,y ( t ) y(t) y ( t ) 靠近 2 2 2 。这个趋势也能从方程改写成的 y ′ = 3 ( 2 − y ) y'=3(2-y) y ′ = 3 ( 2 − y ) 中读出:在 y > 2 y>2 y > 2 的地方,导数为负。
如果候选函数换成 y = 2 + 5 e − 2 t y=2+5e^{-2t} y = 2 + 5 e − 2 t ,那么
y ′ + 3 y = − 10 e − 2 t + 6 + 15 e − 2 t = 6 + 5 e − 2 t , y'+3y=-10e^{-2t}+6+15e^{-2t}=6+5e^{-2t}, y ′ + 3 y = − 10 e − 2 t + 6 + 15 e − 2 t = 6 + 5 e − 2 t ,
结果不等于 6 6 6 。它虽然也靠近 2 2 2 ,图像看起来甚至很相似,却不是这个方程的解。趋势相符可以帮助检查,不能代替代回验证。
在一个点上碰巧成立,还不够
对于 y ′ = 2 y y'=2y y ′ = 2 y ,试着代入 y = t 2 y=t^2 y = t 2 ,会得到
2 t = 2 t 2 . 2t=2t^2. 2 t = 2 t 2 .
这个等式在 t = 0 t=0 t = 0 和 t = 1 t=1 t = 1 处成立,在别处通常不成立。因此 y = t 2 y=t^2 y = t 2 不是这个微分方程在任何非空开区间上的解。微分方程要求的是持续遵守变化规律,几个孤立时刻的巧合不算。
我们据此把“解”说得完整一些:在某个区间上,函数具有方程要求的各阶导数,代入后又在区间内每一点满足等式,这个函数就是方程在该区间上的解。它的图像叫解曲线。研究初值附近的解时,通常用包含初始时刻的开区间;实际模型则常只讨论其中从初始时刻向后的部分。
解的区间也是答案的一部分
考虑另一个候选函数:
y ( t ) = 1 2 − t . y(t)=\frac{1}{2-t}. y ( t ) = 2 − t 1 .
它的导数是
y ′ ( t ) = 1 ( 2 − t ) 2 = y ( t ) 2 , y'(t)=\frac{1}{(2-t)^2}=y(t)^2, y ′ ( t ) = ( 2 − t ) 2 1 = y ( t ) 2 ,
所以它满足 y ′ = y 2 y'=y^2 y ′ = y 2 。但在 t = 2 t=2 t = 2 处,函数根本没有定义。我们可以说它在 ( − ∞ , 2 ) (-\infty,2) ( − ∞ , 2 ) 上是解,也可以说同一个表达式在 ( 2 , ∞ ) (2,\infty) ( 2 , ∞ ) 上是解,却不能说它在整个实数轴上是一条解。
若另外要求 y ( 0 ) = 1 / 2 y(0)=1/2 y ( 0 ) = 1/2 ,所讨论的解区间就要包含 0 0 0 ,因此对应的最大区间是 ( − ∞ , 2 ) (-\infty,2) ( − ∞ , 2 ) 。虽然公式还能在 t > 2 t>2 t > 2 计算,那一段与经过初始点的这一段之间隔着无法跨越的发散点,不能把两段直接连成一个初值问题的解。
这个例子还提醒我们:方程右边 y 2 y^2 y 2 本身处处有意义,并不保证每个解都能一直存在下去。解可能在有限时间内越长越大,最终无法以有限函数值延伸。更系统的存在性与唯一性条件,我们会在后面单独讨论。
检查一个答案时,要同时检查函数、导数、原方程的定义域和初始条件。分母为零、对数无定义、函数不可微,都可能限制解的区间。仅仅写出一个代数表达式,还没有把解交代完整。
“常”说的是自变量,不是变化速度
“常微分方程”这个名字很容易让人误会成“变化率是常数的方程”。其实其中的“常”是在区分普通导数与偏导数,和变化率是否恒定没有关系。
如果把一杯充分搅拌的咖啡看成温度均匀的整体,用 T ( t ) T(t) T ( t ) 描述它,我们只让温度依赖一个自变量——时间。方程里出现 d T / d t dT/dt d T / d t ,这是常微分方程的情形。哪怕 T ′ T' T ′ 每时每刻都在变,仍然叫常微分方程。
如果研究一块金属板,板中心和边缘在同一时刻可能温度不同,就需要用 u ( x , y , t ) u(x,y,t) u ( x , y , t ) 记录不同位置、不同时间的温度。它有多个自变量,描述这些变化的方程可能同时出现 ∂ u / ∂ t \partial u/\partial t ∂ u / ∂ t 、∂ 2 u / ∂ x 2 \partial^2u/\partial x^2 ∂ 2 u / ∂ x 2 等偏导数,那就进入了偏微分方程的范围。
常微分方程关注一个自变量上的变化,偏微分方程关注多个自变量共同决定的变化。
这里还有一个区别:未知函数有几个,与自变量有几个,是两回事。 同时研究两杯互相传热的水,可以有 T 1 ( t ) T_1(t) T 1 ( t ) 、T 2 ( t ) T_2(t) T 2 ( t ) 两个未知函数;只要它们都只依赖时间,描述它们的联立方程仍然是常微分方程组。后面讲线性系统时,我们会把多个量一起放进状态向量里研究。
我们在本课通常用 t t t 表示时间,但常微分方程的自变量也可以是空间位置、距离等。判断依据是未知函数依赖的自变量结构,并不取决于字母选了什么。
先认清方程的结构,再选择解法
遇到一条新方程,先看它需要几阶导数、未知函数怎样出现,会让后面的求解更有方向。就像修理一个装置前先认清部件,分类是在识别结构。
阶数:看最高求到了几次导数
位置函数 x ( t ) x(t) x ( t ) 的一次导数是速度,二次导数是加速度。因此,描述“速度由位置决定”的模型和描述“加速度由位置、速度决定”的模型,在信息结构上不同。
方程中实际出现的最高阶导数,决定方程的阶数。例如
y ′ + 4 y = t y'+4y=t y ′ + 4 y = t
是一阶方程,而
x ′ ′ + 3 x ′ + 2 x = 0 x''+3x'+2x=0 x ′′ + 3 x ′ + 2 x = 0
是二阶方程。若最高出现的是 u ′ ′ ′ u''' u ′′′ ,则是三阶。
不要把“导数的次数”与“导数整体的幂次”混在一起。方程 ( y ′ ) 2 + y = 0 (y')^2+y=0 ( y ′ ) 2 + y = 0 只出现一阶导数,因而仍是一阶方程;平方使它具有另一种结构上的复杂性,却没有把一阶导数变成二阶导数。
还要先看清是否有项可以消去,以及最高阶导数的系数在讨论区间内是否非零。例如 t y ′ ′ + y ′ = 0 t y''+y'=0 t y ′′ + y ′ = 0 在 t ≠ 0 t\ne0 t = 0 的区间上是通常的二阶方程;在 t = 0 t=0 t = 0 处,最高阶项的系数消失,不能直接除以 t t t 后假装什么都没发生。后面把方程整理成标准形式时,这类点会影响我们能使用的结论和解的区间。
线性:看未知函数如何参与运算
看看下面两个方程:
y ′ + t 2 y = sin t , y ′ + y 2 = sin t . y'+t^2y=\sin t,\qquad y'+y^2=\sin t. y ′ + t 2 y = sin t , y ′ + y 2 = sin t .
第一个有 t 2 t^2 t 2 ,却是线性的;第二个把 t 2 y t^2y t 2 y 换成 y 2 y^2 y 2 ,就变成非线性的。原因是线性所限制的对象是未知函数及其导数 。自变量 t t t 是我们输入的数,t 2 t^2 t 2 可以当作已知系数;y y y 才是正在寻找的函数,把它平方会改变方程的结构。
一个 n n n 阶线性方程可以写成
a n ( t ) y ( n ) + a n − 1 ( t ) y ( n − 1 ) + ⋯ + a 1 ( t ) y ′ + a 0 ( t ) y = g ( t ) , a_n(t)y^{(n)}+a_{n-1}(t)y^{(n-1)}+\cdots+a_1(t)y'+a_0(t)y=g(t), a n ( t ) y ( n ) + a n − 1 ( t ) y ( n − 1 ) + ⋯ + a 1 ( t ) y ′ + a 0 ( t ) y = g ( t ) ,
其中各个 a j ( t ) a_j(t) a j ( t ) 和 g ( t ) g(t) g ( t ) 都是已知函数,所讨论的通常区间上要求 a n ( t ) ≠ 0 a_n(t)\ne0 a n ( t ) = 0 。未知函数及其导数只以一次方出现,彼此不相乘,也不被放进平方根、对数、正弦等函数里。
用几个例子把判断做实。y ′ ′ + ( cos t ) y ′ + e t y = t 2 y''+(\cos t)y'+e^t y=t^2 y ′′ + ( cos t ) y ′ + e t y = t 2 是二阶线性方程,系数复杂不妨碍线性。y ′ ′ + y y ′ = 0 y''+yy'=0 y ′′ + y y ′ = 0 是二阶非线性方程,因为 y y y 与 y ′ y' y ′ 相乘;y ′ ′ + sin y = 0 y''+\sin y=0 y ′′ + sin y = 0 也是二阶非线性方程,因为正弦作用在未知函数上。y ′ = y ( 1 − y ) y'=y(1-y) y ′ = y ( 1 − y ) 展开后是 y ′ = y − y 2 y'=y-y^2 y ′ = y − y 2 ,所以是一阶非线性方程。
“线性”也不意味着解曲线必须是直线。刚刚见到的 y ′ = 2 y y'=2y y ′ = 2 y 是线性方程,它的解却是指数曲线。线性说的是方程怎样使用未知函数,不是在描述图像长得直不直。
线性方程只看未知函数及其导数如何出现;平方、相乘或套函数都会破坏线性。
齐次:在线性结构里看有没有独立输入
在线性方程中,把含 y y y 及其导数的项放在左边,剩下的已知项放在右边。若右边是零函数,就叫齐次线性方程;若不是零函数,就叫非齐次线性方程。例如
y ′ + 2 y = 0 y'+2y=0 y ′ + 2 y = 0
是齐次的,而
y ′ + 2 y = cos t y'+2y=\cos t y ′ + 2 y = cos t
是非齐次的。后一条方程的右边在某些时刻会等于零,仍然不改变分类,因为判断的是它是否在整个区间上恒等于零。
在模型里,这个独立项常常对应外部输入,例如持续供热或电源驱动。但分类首先是一件数学上的事,要根据选定变量下的方程判断,不能只凭“现实里有没有外力”来猜。后面一阶方程中还会出现另一个术语“齐次型方程”,含义不同,届时我们会重新说明。
下面的练习器可以帮助你练习分类。每次都先说出依据:最高导数是哪一个,未知函数有没有相乘或被套入其他函数,再看线性结构右端是否恒为零。
通解是一族可能,初始条件说明从哪里开始
前面验证过,y = C e 2 t y=Ce^{2t} y = C e 2 t 对每个常数 C C C 都满足 y ′ = 2 y y'=2y y ′ = 2 y 。这一族覆盖了该方程所有解的表达式,叫作它的通解。任意选定一个常数,例如 C = 4 C=4 C = 4 ,就得到一个具体解,也常叫特解。
这里说“覆盖了所有解”不能只靠猜。对于任意一个满足 y ′ = 2 y y'=2y y ′ = 2 y 的可微函数,利用乘积求导法则,有
d d t ( e − 2 t y ( t ) ) = e − 2 t ( y ′ − 2 y ) = 0. \frac{d}{dt}\left(e^{-2t}y(t)\right)
=e^{-2t}(y'-2y)=0. d t d ( e − 2 t y ( t ) ) = e − 2 t ( y ′ − 2 y ) = 0.
因此 e − 2 t y ( t ) e^{-2t}y(t) e − 2 t y ( t ) 在所讨论的区间上是常数,记成 C C C ,就得到 y = C e 2 t y=Ce^{2t} y = C e 2 t 。这说明任何解都没有跑出这个函数族,而不只是说明这个函数族中的成员都是解。我们这里只借用了乘积求导,如何主动找到这样的乘子,会留到一阶线性方程中解释。
一个初始点怎样选定常数
假如我们要求这条曲线在 t = 0 t=0 t = 0 时取值为 4 4 4 ,就增加条件
y ( 0 ) = 4. y(0)=4. y ( 0 ) = 4.
代入通解得到 C e 0 = C = 4 Ce^0=C=4 C e 0 = C = 4 ,因此选出的函数是
y ( t ) = 4 e 2 t . y(t)=4e^{2t}. y ( t ) = 4 e 2 t .
把变化规律和起点合起来,写成
y ′ = 2 y , y ( 0 ) = 4 , y'=2y,\qquad y(0)=4, y ′ = 2 y , y ( 0 ) = 4 ,
就是一个初值问题。它要求解同时满足两件事:沿途的变化率服从方程,在指定时刻经过指定状态。验证答案时,两件事都要检查。
初始时刻也不一定是零。如果题目给出 y ( 1 ) = 6 y(1)=6 y ( 1 ) = 6 ,由 C e 2 = 6 Ce^2=6 C e 2 = 6 得到 C = 6 e − 2 C=6e^{-2} C = 6 e − 2 ,于是
y ( t ) = 6 e 2 ( t − 1 ) . y(t)=6e^{2(t-1)}. y ( t ) = 6 e 2 ( t − 1 ) .
写成这种形式,一眼就能看见它在 t = 1 t=1 t = 1 时等于 6 6 6 。所谓“初始”,说的是我们用来指定状态的那个时刻,并不要求坐标原点一定放在那里。
图中展示的是初值能唯一确定解的典型情形。对刚才的指数方程,我们已经通过常数 C C C 的计算说明了唯一性;对任意微分方程,不能只凭“给了一个点”就断定一定有且只有一个解。有些方程没有满足指定初值的解,也有些会有多条。后面学存在唯一性时,我们会解释需要检查什么条件。
同样,看到一个含任意常数的解族,也别立刻认定它已经包含全部解。某些求解步骤会丢掉常数解,某些非线性方程还可能有不属于所求函数族的解。现在先养成一个习惯:区分“我找到了一族解”和“我证明这就是所有解”。
为什么二阶问题通常还要给初始速度
想象一个小车,在 t = 0 t=0 t = 0 时位于轨道的同一个位置。它可能静止,也可能正在向右运动,还可能正在向左运动。只给位置,显然不能区分接下来会发生什么。
用一个可以直接积分的模型算一遍。假设小车在一段时间里保持加速度 2 m / s 2 2\,\mathrm{m/s^2} 2 m/ s 2 ,位置为 x ( t ) x(t) x ( t ) ,时间以秒计,那么
x ′ ′ = 2. x''=2. x ′′ = 2.
积分一次得到速度 x ′ = 2 t + C 1 x'=2t+C_1 x ′ = 2 t + C 1 ,再积分得到位置
x = t 2 + C 1 t + C 2 . x=t^2+C_1t+C_2. x = t 2 + C 1 t + C 2 .
如果已知 x ( 0 ) = 3 m x(0)=3\,\mathrm m x ( 0 ) = 3 m ,只能确定 C 2 = 3 C_2=3 C 2 = 3 ,C 1 C_1 C 1 仍然没有定。再给初始速度 x ′ ( 0 ) = − 1 m / s x'(0)=-1\,\mathrm{m/s} x ′ ( 0 ) = − 1 m/s ,才得到 C 1 = − 1 C_1=-1 C 1 = − 1 ,从而
x ( t ) = t 2 − t + 3. x(t)=t^2-t+3. x ( t ) = t 2 − t + 3.
核对一下:x ′ ′ = 2 x''=2 x ′′ = 2 ,x ( 0 ) = 3 x(0)=3 x ( 0 ) = 3 ,x ′ ( 0 ) = − 1 x'(0)=-1 x ′ ( 0 ) = − 1 ,三项要求都满足。由于初始速度是负的,小车开始时向左走,虽然它的加速度始终向右;到 t = 1 / 2 t=1/2 t = 1/2 秒,速度变成零,随后才向右走。这也说明导数的正负必须对应到正确的物理量,正加速度不等于位置一定在增加。
一般的 n n n 阶初值问题,会在同一时刻指定函数值和前 n − 1 n-1 n − 1 阶导数。在适当条件下,这些信息能确定唯一解。上述二阶例子中,每积分一次出现一个常数,是这个规则最容易看见的起点。
变化率始终等于函数值的两倍,这一条件先给出一族指数函数;再指定初始时刻 t=0 的函数值为 3,就能确定其中唯一一个。
四个模型,把现象一步步写成方程
现在我们已经会认方程、验答案,也知道初值的作用。接下来把这些语言放回真实情境中。这里的参数都是为了练习设置的假想数值;每个模型都建立在明确假设上,我们先练习把假设翻译准确,再讨论什么时候它还适用。
人口增长:每个个体平均贡献多少新增数量
培养皿里细菌越多,短时间内发生分裂的机会通常也越多。我们先假设营养充足、环境保持稳定,没有迁入迁出,并用连续可微的函数近似大量个体的总数量。个体数本来是整数,把它连续化,是为了描述整体趋势的一个近似。
设 P ( t ) P(t) P ( t ) 为数量,时间以小时计。若每个个体平均每小时贡献的净增长比例是常数 r r r ,那么在很短的 Δ t \Delta t Δ t 小时内,新增数量近似为 r P ( t ) Δ t rP(t)\Delta t r P ( t ) Δ t 。除以 Δ t \Delta t Δ t ,再用导数表示瞬时变化率,就得到
P ′ = r P . P'=rP. P ′ = r P .
P ′ P' P ′ 的单位是“个每小时”,因此 r r r 的单位是“每小时”。这句话说的是相对增长率 P ′ / P = r P'/P=r P ′ / P = r 保持不变,并不是每小时增加的个数保持不变。当前数量翻倍时,当前增长速度也翻倍。
设 r = 0.3 h − 1 r=0.3\,\mathrm{h}^{-1} r = 0.3 h − 1 ,初始数量 P ( 0 ) = 800 P(0)=800 P ( 0 ) = 800 。照着前面的指数函数验证思路,候选解是
P ( t ) = 800 e 0.3 t . P(t)=800e^{0.3t}. P ( t ) = 800 e 0.3 t .
求导得 P ′ = 240 e 0.3 t = 0.3 P P'=240e^{0.3t}=0.3P P ′ = 240 e 0.3 t = 0.3 P ,初值也对。在开始时,增长速度是 240 240 240 个每小时;模型预测数量翻倍的时刻满足 e 0.3 t = 2 e^{0.3t}=2 e 0.3 t = 2 ,所以
t = ln 2 0.3 ≈ 2.31 小时 . t=\frac{\ln2}{0.3}\approx2.31\text{ 小时}. t = 0.3 ln 2 ≈ 2.31 小时 .
这里的 r r r 是净增长率,可以理解为平均出生率减去平均死亡率。若 r < 0 r<0 r < 0 ,同一种方程描述衰减;若 r = 0 r=0 r = 0 ,数量不变。人口模型通常只考虑 P ≥ 0 P\ge0 P ≥ 0 。
恒定正增长率会使指数函数一直增大,真实培养皿却没有无限的营养和空间。若要描述资源压力,一种后续修正是让每个个体的净增长率随着数量增加而下降,例如引入承载量 K > 0 K>0 K > 0 ,写成
P ′ = r P ( 1 − P K ) , r > 0. P'=rP\left(1-\frac{P}{K}\right),\qquad r>0. P ′ = r P ( 1 − K P ) , r > 0.
现在先读懂这句话就够了:当 P P P 远小于 K K K 时,括号接近 1 1 1 ,仍然近似指数增长;当 P P P 接近 K K K 时,增长逐渐变慢;P = K P=K P = K 时增长率为零。这就是后面要学习的逻辑斯蒂模型。我们会进一步求解,但此刻最该看清的是:换一条现实假设,方程的结构也会跟着改变。
咖啡冷却:决定快慢的是温差
一杯热咖啡放在房间里,通常刚开始降温较快,后来越来越慢。为了把这个现象写成方程,我们把咖啡看成内部温度均匀的整体,假设室温保持不变,而且这段温度范围内的换热效果可以用一个常数表示。
设咖啡温度为 T ( t ) T(t) T ( t ) ,环境温度为 T a T_a T a ,冷却系数为 k > 0 k>0 k > 0 。先把“温差越大,变化越快”写成变化率的大小与 ∣ T − T a ∣ |T-T_a| ∣ T − T a ∣ 成正比,再把变化方向放进去,就得到
T ′ = − k ( T − T a ) . T'=-k(T-T_a). T ′ = − k ( T − T a ) .
负号的作用可以在两种情形下检查。咖啡比房间热时,T − T a > 0 T-T_a>0 T − T a > 0 ,于是 T ′ < 0 T'<0 T ′ < 0 ,温度下降;一杯冷饮比房间冷时,T − T a < 0 T-T_a<0 T − T a < 0 ,于是 T ′ > 0 T'>0 T ′ > 0 ,温度回升。同一个方程能描述冷却和回温,因为它让温度朝环境温度靠近。
设室温 T a = 24 ∘ C T_a=24^\circ\mathrm C T a = 2 4 ∘ C ,咖啡初温 84 ∘ C 84^\circ\mathrm C 8 4 ∘ C ,时间以分钟计,k = 0.1 m i n − 1 k=0.1\,\mathrm{min}^{-1} k = 0.1 min − 1 。初值问题为
T ′ = − 0.1 ( T − 24 ) , T ( 0 ) = 84. T'=-0.1(T-24),\qquad T(0)=84. T ′ = − 0.1 ( T − 24 ) , T ( 0 ) = 84.
开始时 T ′ ( 0 ) = − 0.1 ( 84 − 24 ) = − 6 T'(0)=-0.1(84-24)=-6 T ′ ( 0 ) = − 0.1 ( 84 − 24 ) = − 6 ,单位是摄氏度每分钟。它表示初始瞬间的降温速率,并不表示以后每分钟都降 6 6 6 度。
利用指数函数,我们可以检查候选解
T ( t ) = 24 + 60 e − 0.1 t . T(t)=24+60e^{-0.1t}. T ( t ) = 24 + 60 e − 0.1 t .
它在 t = 0 t=0 t = 0 时为 84 84 84 ,导数是 − 6 e − 0.1 t -6e^{-0.1t} − 6 e − 0.1 t ;右边 − 0.1 ( T − 24 ) -0.1(T-24) − 0.1 ( T − 24 ) 也等于 − 6 e − 0.1 t -6e^{-0.1t} − 6 e − 0.1 t ,所以验证通过。十分钟后模型给出 T ( 10 ) = 24 + 60 / e ≈ 46.07 ∘ C T(10)=24+60/e\approx46.07^\circ\mathrm C T ( 10 ) = 24 + 60/ e ≈ 46.0 7 ∘ C 。若问多久降到 54 ∘ C 54^\circ\mathrm C 5 4 ∘ C ,就是温差从 60 60 60 度变为 30 30 30 度,于是
60 e − 0.1 t = 30 , t = 10 ln 2 ≈ 6.93 分钟 . 60e^{-0.1t}=30,
\qquad t=10\ln2\approx6.93\text{ 分钟}. 60 e − 0.1 t = 30 , t = 10 ln 2 ≈ 6.93 分钟 .
在这个理想模型里,只要初始温度不等于室温,温差就会在每个有限时刻保持非零,逐渐趋近于零。实际测量有分辨率,我们会在某个时候读到“已经是室温”,两句话并不冲突。若环境也被显著加热,或者咖啡内部温度并不均匀,就要调整原来的建模假设。
RC 电路:储存电荷让响应需要时间
电源接上以后,电容上的电压通常不会立刻跳到电源电压。原因是电容需要逐渐积累电荷,而电流正是电荷积累的速度。这个故事只要写清楚两个元件的关系,就会自然出现导数。
考虑电阻 R > 0 R>0 R > 0 与电容 C > 0 C>0 C > 0 串联的理想电路,设输入电压为 V i n ( t ) V_{\mathrm{in}}(t) V in ( t ) ,电容电压为 V C ( t ) V_C(t) V C ( t ) 。电容上的电荷满足 Q = C V C Q=CV_C Q = C V C 。假设电容值不随时间改变,电流就满足
I = Q ′ = C V C ′ . I=Q'=CV_C'. I = Q ′ = C V C ′ .
电阻上的电压是 R I RI R I 。沿回路按一致方向计量电压,电阻电压与电容电压之和等于输入电压,所以
R I + V C = V i n ( t ) . RI+V_C=V_{\mathrm{in}}(t). R I + V C = V in ( t ) .
把 I = C V C ′ I=CV_C' I = C V C ′ 代进去,得到
R C V C ′ + V C = V i n ( t ) . RCV_C'+V_C=V_{\mathrm{in}}(t). R C V C ′ + V C = V in ( t ) .
这是一阶线性方程,未知函数是电容电压。若想得到电流,可以在得到 V C V_C V C 后再通过 I = C V C ′ I=CV_C' I = C V C ′ 计算。不要把“研究电路”自动等同于“未知量一定是电流”。
假设 R = 1000 Ω R=1000\,\Omega R = 1000 Ω ,C = 0.002 F C=0.002\,\mathrm F C = 0.002 F ,则 R C = 2 s RC=2\,\mathrm s R C = 2 s 。从 t = 0 t=0 t = 0 起接入恒定 10 V 10\,\mathrm V 10 V 电源,电容初始未充电,那么在接通后的时段,模型写为
2 V C ′ + V C = 10 , V C ( 0 ) = 0. 2V_C'+V_C=10,\qquad V_C(0)=0. 2 V C ′ + V C = 10 , V C ( 0 ) = 0.
候选解 V C ( t ) = 10 ( 1 − e − t / 2 ) V_C(t)=10(1-e^{-t/2}) V C ( t ) = 10 ( 1 − e − t /2 ) 在开始时确实是零,且
2 V C ′ + V C = 2 ⋅ 5 e − t / 2 + 10 − 10 e − t / 2 = 10. 2V_C'+V_C=2\cdot5e^{-t/2}+10-10e^{-t/2}=10. 2 V C ′ + V C = 2 ⋅ 5 e − t /2 + 10 − 10 e − t /2 = 10.
所以它满足模型。电流为 I = 0.002 ⋅ 5 e − t / 2 = 0.01 e − t / 2 A I=0.002\cdot5e^{-t/2}=0.01e^{-t/2}\,\mathrm A I = 0.002 ⋅ 5 e − t /2 = 0.01 e − t /2 A ,开始时最大,随后衰减。电容电压逐渐接近 10 V 10\,\mathrm V 10 V ,电阻两端的压差越来越小,电流自然也越来越小。
把方程整理成 V C ′ = ( 10 − V C ) / 2 V_C'=(10-V_C)/2 V C ′ = ( 10 − V C ) /2 ,你会发现它与冷却模型很相似:变化率由“目标值与当前值的差”决定。这个相似性会在后面的解法中得到统一解释。这里的 R C RC R C 具有时间单位,指数中的 t / ( R C ) t/(RC) t / ( R C ) 因而没有单位,这也是检查公式是否合理的一条线索。
弹簧振子:位置和速度一起决定加速度
把一个小物块连在水平弹簧上,拉开一点再松手。物块离平衡位置越远,弹簧往回拉的力通常越大;若还有阻尼,运动越快,阻碍运动的力也越大。两种力作用在一起,决定的是加速度,因此方程会出现二阶导数。
设 x ( t ) x(t) x ( t ) 是相对平衡位置的位移,向右为正。质量为 m > 0 m>0 m > 0 ,弹簧刚度为 k > 0 k>0 k > 0 ,阻尼系数为 c ≥ 0 c\ge0 c ≥ 0 。这里的 k k k 表示弹簧刚度,和前面冷却模型中的 k k k 是不同参数,只是各自在本模型中使用的常见记号。
在弹簧伸缩不大、回复力可以近似与位移成正比的范围内,弹簧力是 − k x -kx − k x 。若阻尼力近似与速度成正比且方向相反,则为 − c x ′ -cx' − c x ′ 。再设外部驱动力为 F ( t ) F(t) F ( t ) ,由“质量乘加速度等于合力”得到
m x ′ ′ = − c x ′ − k x + F ( t ) , mx''=-cx'-kx+F(t), m x ′′ = − c x ′ − k x + F ( t ) ,
也就是
m x ′ ′ + c x ′ + k x = F ( t ) . mx''+cx'+kx=F(t). m x ′′ + c x ′ + k x = F ( t ) .
每一项都具有力的单位:m m m 用千克,x x x 用米,t t t 用秒时,k k k 的单位是牛顿每米,c c c 的单位是牛顿秒每米。若用竖直弹簧并从重力作用下的平衡位置测位移,重力与静态伸长产生的弹簧力抵消,也会得到相同形式;选对位置原点,能让方程更清楚。
先取一个没有阻尼、也没有外力的例子。设 m = 1 k g m=1\,\mathrm{kg} m = 1 kg 、k = 4 N / m k=4\,\mathrm{N/m} k = 4 N/m ,物块初始向右偏离 0.1 m 0.1\,\mathrm m 0.1 m 后静止释放,则
x ′ ′ + 4 x = 0 , x ( 0 ) = 0.1 , x ′ ( 0 ) = 0. x''+4x=0,\qquad x(0)=0.1,\qquad x'(0)=0. x ′′ + 4 x = 0 , x ( 0 ) = 0.1 , x ′ ( 0 ) = 0.
我们不急着推导全部解法,先验证 x ( t ) = 0.1 cos ( 2 t ) x(t)=0.1\cos(2t) x ( t ) = 0.1 cos ( 2 t ) 。求导两次得到
x ′ ( t ) = − 0.2 sin ( 2 t ) , x ′ ′ ( t ) = − 0.4 cos ( 2 t ) . x'(t)=-0.2\sin(2t),\qquad x''(t)=-0.4\cos(2t). x ′ ( t ) = − 0.2 sin ( 2 t ) , x ′′ ( t ) = − 0.4 cos ( 2 t ) .
因此 x ′ ′ + 4 x = 0 x''+4x=0 x ′′ + 4 x = 0 ,初始位置和初始速度也都正确。这个解描述来回振动,周期为 π \pi π 秒。它一直保持相同振幅,是因为模型特意忽略了阻尼。加入阻尼后振幅怎样变化,加入周期外力后会不会越振越大,将是二阶方程部分要回答的问题。
从文字到答案,中间还有几次检查
四个模型的故事不同,写方程时的思考顺序却可以重复使用。先选定正在变化的量,说明它的输入、输出和单位;接着把影响变化的因素说清楚,再把这句话变成导数关系。然后加入起始状态,最后把方程与候选解放回现象中,检查符号、单位和适用范围。
从现象出发建立一阶微分方程模型的基本闭环。
比如一段文字说“数量每小时增加 20 20 20 个”,它对应 P ′ = 20 P'=20 P ′ = 20 ;说“每小时净增长率为当前数量的 20 % 20\% 20% ”,在连续变化率的解释下对应 P ′ = 0.2 P P'=0.2P P ′ = 0.2 P ,其中 0.2 0.2 0.2 的单位是每小时。前者每小时增加同样多,后者当前数量越大就增加越多。一句话中的“个”和“比例”,会让模型从直线变化转向指数变化。
再比如小车匀加速模型 x = t 2 − t + 3 x=t^2-t+3 x = t 2 − t + 3 在数学上对所有实数 t t t 有定义,但如果题目说小车只在前五秒保持恒定加速度,那么它对现实的预测范围就只能先放在 0 ≤ t ≤ 5 0\le t\le5 0 ≤ t ≤ 5 。函数的数学定义域、微分方程解的区间、现实假设的有效时间范围,彼此相关,却不能混为一谈。
我们也不必要求每个方程都立刻写出一个漂亮公式。一个模型建立以后,可以先从右端正负判断增加还是减少,从右端为零寻找可能不变的状态,再想办法求解或近似计算。下一章的方向场就是沿着这个想法,把许多位置上的斜率画出来,让局部规则逐渐显出整条曲线的走向。
练习:把定义用到具体问题上
验证函数,而不是只看外形
判断 y 1 ( t ) = 3 + 2 e − 4 t y_1(t)=3+2e^{-4t} y 1 ( t ) = 3 + 2 e − 4 t 和 y 2 ( t ) = 3 + 2 e − 2 t y_2(t)=3+2e^{-2t} y 2 ( t ) = 3 + 2 e − 2 t 是否满足 y ′ + 4 y = 12 y'+4y=12 y ′ + 4 y = 12 。若满足,再检查初值 y ( 0 ) = 5 y(0)=5 y ( 0 ) = 5 ,并说明解的区间。
显示推导答案 对第一个函数,y 1 ′ = − 8 e − 4 t y_1'=-8e^{-4t} y 1 ′ = − 8 e − 4 t ,因此
y 1 ′ + 4 y 1 = − 8 e − 4 t + 12 + 8 e − 4 t = 12. y_1'+4y_1=-8e^{-4t}+12+8e^{-4t}=12. y 1 ′ + 4 y 1 = − 8 e − 4 t + 12 + 8 e − 4 t = 12. 同时 y 1 ( 0 ) = 3 + 2 = 5 y_1(0)=3+2=5 y 1 ( 0 ) = 3 + 2 = 5 ,所以它是所给初值问题的解,最大解区间为 ( − ∞ , ∞ ) (-\infty,\infty) ( − ∞ , ∞ ) 。第二个函数虽然同样满足初值,但 y 2 ′ = − 4 e − 2 t y_2'=-4e^{-2t} y 2 ′ = − 4 e − 2 t ,代入后得到 y 2 ′ + 4 y 2 = 12 + 4 e − 2 t y_2'+4y_2=12+4e^{-2t} y 2 ′ + 4 y 2 = 12 + 4 e − 2 t ,并不恒等于 12 12 12 ,因此不是方程的解。满足初始点,不能代替满足沿途的变化规律。
分类时盯住未知函数
判断下面四个方程的阶数与线性;对线性方程再判断齐次性:y ′ ′ + t 3 y ′ = e t y''+t^3y'=e^t y ′′ + t 3 y ′ = e t ,( y ′ ) 3 + y = t (y')^3+y=t ( y ′ ) 3 + y = t ,y ′ ′ + sin ( t ) y = 0 y''+\sin(t)y=0 y ′′ + sin ( t ) y = 0 ,y ′ ′ + sin ( y ) = 0 y''+\sin(y)=0 y ′′ + sin ( y ) = 0 。
显示推导答案 第一个最高导数是 y ′ ′ y'' y ′′ ,是二阶线性非齐次方程;t 3 t^3 t 3 是已知系数,不破坏线性,右端 e t e^t e t 不恒为零。第二个只含一阶导数,所以是一阶方程,但 ( y ′ ) 3 (y')^3 ( y ′ ) 3 使它非线性,不能在此套用线性方程的齐次分类。第三个是二阶线性齐次方程,sin t \sin t sin t 只是 y y y 的已知系数。第四个是二阶非线性方程,正弦作用在未知函数上。两条方程只改变了正弦括号里的对象,线性结构就不同了。
初始时刻可以不在零点
求解可以直接积分的初值问题 y ′ = 6 t y'=6t y ′ = 6 t ,y ( 2 ) = 7 y(2)=7 y ( 2 ) = 7 。再说明它的图像与通解之间的关系。
显示推导答案 对 6 t 6t 6 t 积分,通解为 y = 3 t 2 + C y=3t^2+C y = 3 t 2 + C 。由 y ( 2 ) = 12 + C = 7 y(2)=12+C=7 y ( 2 ) = 12 + C = 7 得 C = − 5 C=-5 C = − 5 ,所以
y ( t ) = 3 t 2 − 5. y(t)=3t^2-5. y ( t ) = 3 t 2 − 5. 求导得到 y ′ = 6 t y'=6t y ′ = 6 t ,再代入 t = 2 t=2 t = 2 得到 7 7 7 ,两项都满足。通解是一族上下平移的抛物线,初始条件要求其中一条经过 ( 2 , 7 ) (2,7) ( 2 , 7 ) ,因而固定了平移量。所求解的最大区间为整个实数轴。
同一个表达式,初值决定讨论哪一段
验证 y ( t ) = 1 / ( 3 − t ) y(t)=1/(3-t) y ( t ) = 1/ ( 3 − t ) 满足 y ′ = y 2 y'=y^2 y ′ = y 2 。若给定 y ( 0 ) = 1 / 3 y(0)=1/3 y ( 0 ) = 1/3 ,最大解区间是什么?若改为 y ( 4 ) = − 1 y(4)=-1 y ( 4 ) = − 1 呢?
显示推导答案 求导得到 y ′ = 1 / ( 3 − t ) 2 y'=1/(3-t)^2 y ′ = 1/ ( 3 − t ) 2 ,与 y 2 y^2 y 2 相等,验证在 t ≠ 3 t\ne3 t = 3 时成立。第一个初值确实满足,包含 0 0 0 的最大开区间是 ( − ∞ , 3 ) (-\infty,3) ( − ∞ , 3 ) 。第二个初值也满足,但包含 4 4 4 的最大开区间是 ( 3 , ∞ ) (3,\infty) ( 3 , ∞ ) 。这两个区间不能跨过 t = 3 t=3 t = 3 拼接,因为函数在那里无定义,而且从两侧靠近时都没有有限极限。
从回温故事写出初值问题
一杯温度为 8 ∘ C 8^\circ\mathrm C 8 ∘ C 的饮料放进 26 ∘ C 26^\circ\mathrm C 2 6 ∘ C 的房间,假设饮料内部温度均匀、室温不变,温度变化率与温差成正比,系数为 0.05 m i n − 1 0.05\,\mathrm{min}^{-1} 0.05 min − 1 。写出初值问题,求初始变化率,并验证候选函数 T ( t ) = 26 − 18 e − 0.05 t T(t)=26-18e^{-0.05t} T ( t ) = 26 − 18 e − 0.05 t 。
显示推导答案 令 T ( t ) T(t) T ( t ) 为 t t t 分钟后的温度,初值问题为
T ′ = − 0.05 ( T − 26 ) , T ( 0 ) = 8. T'=-0.05(T-26),\qquad T(0)=8. T ′ = − 0.05 ( T − 26 ) , T ( 0 ) = 8. 初始导数 T ′ ( 0 ) = − 0.05 ( 8 − 26 ) = 0.9 T'(0)=-0.05(8-26)=0.9 T ′ ( 0 ) = − 0.05 ( 8 − 26 ) = 0.9 ,单位是摄氏度每分钟,正号表示回温。候选函数的导数为 0.9 e − 0.05 t 0.9e^{-0.05t} 0.9 e − 0.05 t ;代入方程右端得到 − 0.05 ( − 18 e − 0.05 t ) = 0.9 e − 0.05 t -0.05(-18e^{-0.05t})=0.9e^{-0.05t} − 0.05 ( − 18 e − 0.05 t ) = 0.9 e − 0.05 t ,并且 T ( 0 ) = 8 T(0)=8 T ( 0 ) = 8 ,所以通过验证。它在 t ≥ 0 t\ge0 t ≥ 0 时从下方靠近 26 ∘ C 26^\circ\mathrm C 2 6 ∘ C ,温度上升但上升速度逐渐减小。
两次积分与两条初始信息
假想一辆小车在研究时段内满足 x ′ ′ = − 4 x''=-4 x ′′ = − 4 ,时间以秒计,位置以米计,初始位置 x ( 0 ) = 2 x(0)=2 x ( 0 ) = 2 ,初始速度 x ′ ( 0 ) = 6 x'(0)=6 x ′ ( 0 ) = 6 。求位置函数,并求它第一次速度为零的时刻及此时的位置。
显示推导答案 第一次积分得到 x ′ = − 4 t + C 1 x'=-4t+C_1 x ′ = − 4 t + C 1 ,由初始速度得 C 1 = 6 C_1=6 C 1 = 6 。第二次积分得到 x = − 2 t 2 + 6 t + C 2 x=-2t^2+6t+C_2 x = − 2 t 2 + 6 t + C 2 ,由初始位置得 C 2 = 2 C_2=2 C 2 = 2 ,因此
x ( t ) = − 2 t 2 + 6 t + 2. x(t)=-2t^2+6t+2. x ( t ) = − 2 t 2 + 6 t + 2. 令速度 − 4 t + 6 = 0 -4t+6=0 − 4 t + 6 = 0 ,得 t = 1.5 t=1.5 t = 1.5 秒;此时位置为 − 2 ( 1.5 ) 2 + 6 ( 1.5 ) + 2 = 6.5 -2(1.5)^2+6(1.5)+2=6.5 − 2 ( 1.5 ) 2 + 6 ( 1.5 ) + 2 = 6.5 米。虽然加速度始终为负,但在 0 < t < 1.5 0<t<1.5 0 < t < 1.5 时速度仍为正,小车仍朝正方向运动,只是越来越慢。核对得到 x ′ ′ = − 4 x''=-4 x ′′ = − 4 且两个初值都满足,结果只用于题设恒加速度假设成立的时段。
做完这些练习,再回到一开始的 y ′ = 2 y y'=2y y ′ = 2 y :你现在能说明它在找什么,怎样检查一个答案,为什么会有任意常数,以及初始点怎样进入问题。接下来我们会把“每个状态给出一个斜率”真的画成图。即使一时求不出函数公式,也可以先看清曲线从哪里出发、向哪里走。