自在学

我们与你共同进步

  • 分类课程
  • 文章
  • 工作台
  • 订阅

  • 关于我们
  • 隐私政策
  • 使用条款

探索

  • 分类课程
  • 文章
  • 工作台
  • 订阅

网站信息

  • 关于我们
  • 隐私政策
  • 使用条款

加入社区

自在学学习社区微信二维码

微信扫码,交流学习

株洲市自在学教育科技有限公司© 2025 - 2026 版权所有

© 2025 - 2026 株洲市自在学教育科技有限公司 版权所有

湘公网安备43020302000292号|湘ICP备2025148919号-1
分类课程工作台文章订阅
分类课程工作台文章价格

常微分方程 I

  1. 01常微分方程研究什么
  2. 02建模、初值问题与解的几何图像
  3. 03可分离变量方程
  4. 04一阶线性方程与积分因子
  5. 05精确方程、替换法与一阶方程工具箱
  6. 06存在唯一性、自治方程与相线
  7. 07一阶模型:从真实情境到方程
  8. 08数值方法:沿着变化率,一步步算出解
  9. 09二阶线性方程的结构
  10. 10常系数齐次方程与自由振动
  11. 11非齐次方程、待定系数法与共振
  12. 12变参数法与解的构造思想
  13. 13Laplace 变换:把初值和输入一起带进代数方程
  14. 14Laplace 变换 II:阶跃、冲击与卷积
  15. 15幂级数解法:把未知函数一项一项算出来
  16. 16一阶线性系统与矩阵方法
  17. 17平面线性系统、相图与稳定性
  18. 18非线性系统入门与综合建模
正在加载课程章节内容
课程数学常微分方程 I常微分方程研究什么

常微分方程研究什么

前面学方程时,我们常常在找一个数。比如 x2=4x^2=4x2=4,找到 x=2x=2x=2 和 x=−2x=-2x=−2,问题就解决了。到了常微分方程,这件熟悉的事会发生一次变化:要找的未知数,变成了一整个函数。

先想象一个培养皿。我们暂时不知道里面的细菌数量每天会是多少,却可以提出一个变化规律:细菌越多,同样一小段时间里新增的细菌也越多。如果把数量记成 P(t)P(t)P(t),这句话就把“现在有多少”和“现在增长多快”联系起来了。微分方程要做的,正是从这样的联系出发,找出数量随时间变化的函数。

咖啡冷却也有类似的关系,只是决定降温快慢的变成了温差;弹簧振动更复杂一些,当前位置和速度会一起影响加速度。我们会反复在现象、方程和解之间来回走。先把现象说清楚,方程才有来处;把解放回现象中检查,计算才有落点。


从求一个数,到找一整条曲线

不妨先离开具体单位,看看这个最简单的例子:

y′=2y.y'=2y.y′=2y.

这里的 yyy 是未知函数 y(t)y(t)y(t) 的简写,y′y'y′ 是它对 ttt 的导数。把式子读成一句话,就是:“在每一个时刻,函数的变化率都等于它当前数值的两倍。”

如果有人回答“y=3y=3y=3”,我们需要追问:是说某个时刻函数值为 333,还是说函数永远等于 333?前一种回答只给了曲线上的一个点,远远不够;后一种回答指的是常函数 y(t)=3y(t)=3y(t)=3,它的导数为 000,而右边 2y=62y=62y=6,并不满足方程。

那么,什么函数有希望?从微积分里,我们知道指数函数求导后仍然保持指数函数的样子。试一试 y(t)=e2ty(t)=e^{2t}y(t)=e2t:

y′(t)=2e2t=2y(t).y'(t)=2e^{2t}=2y(t).y′(t)=2e2t=2y(t).

等式对每个实数 ttt 都成立,所以它确实是一个解。再试 y(t)=3e2ty(t)=3e^{2t}y(t)=3e2t,求导得到 6e2t6e^{2t}6e2t,恰好也等于自身的两倍。实际上,对任意实常数 CCC,

y(t)=Ce2ty(t)=Ce^{2t}y(t)=Ce2t

都满足这个方程。C=0C=0C=0 给出一直贴在横轴上的零函数;C>0C>0C>0 给出向上增长的曲线;C<0C<0C<0 给出数值越来越负的曲线。数学方程允许这些解,若它具体描述人口,我们还要根据人口不能为负这一点限制初值。

请停下来比较一下两种求解。代数方程 x2=4x^2=4x2=4 要求某个数满足关系;这里要求函数和它的导数在一整段时间里始终配合。所以微分方程的答案通常要写成 y(t)=⋯y(t)=\cdotsy(t)=⋯,而不能只交出一个数。

导数已经告诉了我们什么

导数的含义并没有因为课程换了名字而改变。若一个量在时刻 ttt 的值是 y(t)y(t)y(t),经过一小段时间 hhh,可用

y(t+h)≈y(t)+y′(t)hy(t+h)\approx y(t)+y'(t)hy(t+h)≈y(t)+y′(t)h

估计新值。这里的近似来自可微性,只适合足够小的时间步长。对于 y′=2yy'=2yy′=2y,如果当前 y=3y=3y=3,当前斜率就是 666;取一个很小的正步长,曲线会从这个位置向上走。

注意这并不意味着以后一直以斜率 666 增长。函数值改变后,右边的 2y2y2y 也会改变,下一刻的斜率随之更新。这种“状态影响变化率,变化率又改变状态”的关系,是许多微分方程模型的共同结构。

更一般的一阶方程常写成

dydt=f(t,y).\frac{dy}{dt}=f(t,y).dtdy​=f(t,y).

右边的 fff 是已知规则。它可以只看时间,如 y′=cos⁡ty'=\cos ty′=cost;可以只看当前状态,如 y′=2yy'=2yy′=2y;也可以两者都看,如 y′=t−yy'=t-yy′=t−y。求解时未知的是沿着时间走下去的函数 y(t)y(t)y(t),并不是右边这个规则。

下面的交互把几种变化放在一起。移动时间探针时,可以先看当前函数值,再看切线斜率:数值高,不一定意味着变化快;变化快慢最终要由各自的方程决定。


判断一个函数是不是解,要把它放回去

现在即使还不会系统求解,也已经能做一件很有用的事:检查别人给出的答案。方法只有一个核心动作——求出需要的导数,代回原方程,看等式是否在所讨论的区间内处处成立。

一个算完的验证例子

假设要检查函数

y(t)=2+5e−3ty(t)=2+5e^{-3t}y(t)=2+5e−3t

是否满足方程

y′+3y=6.y'+3y=6.y′+3y=6.

先求导。常数 222 的导数是 000,指数项求导要乘上指数中 ttt 的系数 −3-3−3,所以

y′=−15e−3t.y'=-15e^{-3t}.y′=−15e−3t.

再把函数和导数一起代入左边:

y′+3y=−15e−3t+3(2+5e−3t)=6.y'+3y=-15e^{-3t}+3(2+5e^{-3t})=6.y′+3y=−15e−3t+3(2+5e−3t)=6.

左边确实等于右边,而且没有限制 ttt 的分母、根号或对数,所以这个函数在整个实数轴上都是方程的解。它还有一个容易看懂的趋势:随着 ttt 增大,指数项逐渐变小,y(t)y(t)y(t) 靠近 222。这个趋势也能从方程改写成的 y′=3(2−y)y'=3(2-y)y′=3(2−y) 中读出:在 y>2y>2y>2 的地方,导数为负。

如果候选函数换成 y=2+5e−2ty=2+5e^{-2t}y=2+5e−2t,那么

y′+3y=−10e−2t+6+15e−2t=6+5e−2t,y'+3y=-10e^{-2t}+6+15e^{-2t}=6+5e^{-2t},y′+3y=−10e−2t+6+15e−2t=6+5e−2t,

结果不等于 666。它虽然也靠近 222,图像看起来甚至很相似,却不是这个方程的解。趋势相符可以帮助检查,不能代替代回验证。

在一个点上碰巧成立,还不够

对于 y′=2yy'=2yy′=2y,试着代入 y=t2y=t^2y=t2,会得到

2t=2t2.2t=2t^2.2t=2t2.

这个等式在 t=0t=0t=0 和 t=1t=1t=1 处成立,在别处通常不成立。因此 y=t2y=t^2y=t2 不是这个微分方程在任何非空开区间上的解。微分方程要求的是持续遵守变化规律,几个孤立时刻的巧合不算。

我们据此把“解”说得完整一些:在某个区间上,函数具有方程要求的各阶导数,代入后又在区间内每一点满足等式,这个函数就是方程在该区间上的解。它的图像叫解曲线。研究初值附近的解时,通常用包含初始时刻的开区间;实际模型则常只讨论其中从初始时刻向后的部分。

解的区间也是答案的一部分

考虑另一个候选函数:

y(t)=12−t.y(t)=\frac{1}{2-t}.y(t)=2−t1​.

它的导数是

y′(t)=1(2−t)2=y(t)2,y'(t)=\frac{1}{(2-t)^2}=y(t)^2,y′(t)=(2−t)21​=y(t)2,

所以它满足 y′=y2y'=y^2y′=y2。但在 t=2t=2t=2 处,函数根本没有定义。我们可以说它在 (−∞,2)(-\infty,2)(−∞,2) 上是解,也可以说同一个表达式在 (2,∞)(2,\infty)(2,∞) 上是解,却不能说它在整个实数轴上是一条解。

若另外要求 y(0)=1/2y(0)=1/2y(0)=1/2,所讨论的解区间就要包含 000,因此对应的最大区间是 (−∞,2)(-\infty,2)(−∞,2)。虽然公式还能在 t>2t>2t>2 计算,那一段与经过初始点的这一段之间隔着无法跨越的发散点,不能把两段直接连成一个初值问题的解。

这个例子还提醒我们:方程右边 y2y^2y2 本身处处有意义,并不保证每个解都能一直存在下去。解可能在有限时间内越长越大,最终无法以有限函数值延伸。更系统的存在性与唯一性条件,我们会在后面单独讨论。

检查一个答案时,要同时检查函数、导数、原方程的定义域和初始条件。分母为零、对数无定义、函数不可微,都可能限制解的区间。仅仅写出一个代数表达式,还没有把解交代完整。


“常”说的是自变量,不是变化速度

“常微分方程”这个名字很容易让人误会成“变化率是常数的方程”。其实其中的“常”是在区分普通导数与偏导数,和变化率是否恒定没有关系。

如果把一杯充分搅拌的咖啡看成温度均匀的整体,用 T(t)T(t)T(t) 描述它,我们只让温度依赖一个自变量——时间。方程里出现 dT/dtdT/dtdT/dt,这是常微分方程的情形。哪怕 T′T'T′ 每时每刻都在变,仍然叫常微分方程。

如果研究一块金属板,板中心和边缘在同一时刻可能温度不同,就需要用 u(x,y,t)u(x,y,t)u(x,y,t) 记录不同位置、不同时间的温度。它有多个自变量,描述这些变化的方程可能同时出现 ∂u/∂t\partial u/\partial t∂u/∂t、∂2u/∂x2\partial^2u/\partial x^2∂2u/∂x2 等偏导数,那就进入了偏微分方程的范围。

常微分方程与偏微分方程对照图,左侧为 T(t) 随时间变化曲线,右侧为薄板温度分布热图 u(x,y,t)。
常微分方程关注一个自变量上的变化,偏微分方程关注多个自变量共同决定的变化。

这里还有一个区别:未知函数有几个,与自变量有几个,是两回事。同时研究两杯互相传热的水,可以有 T1(t)T_1(t)T1​(t)、T2(t)T_2(t)T2​(t) 两个未知函数;只要它们都只依赖时间,描述它们的联立方程仍然是常微分方程组。后面讲线性系统时,我们会把多个量一起放进状态向量里研究。

我们在本课通常用 ttt 表示时间,但常微分方程的自变量也可以是空间位置、距离等。判断依据是未知函数依赖的自变量结构,并不取决于字母选了什么。


先认清方程的结构,再选择解法

遇到一条新方程,先看它需要几阶导数、未知函数怎样出现,会让后面的求解更有方向。就像修理一个装置前先认清部件,分类是在识别结构。

阶数:看最高求到了几次导数

位置函数 x(t)x(t)x(t) 的一次导数是速度,二次导数是加速度。因此,描述“速度由位置决定”的模型和描述“加速度由位置、速度决定”的模型,在信息结构上不同。

方程中实际出现的最高阶导数,决定方程的阶数。例如

y′+4y=ty'+4y=ty′+4y=t

是一阶方程,而

x′′+3x′+2x=0x''+3x'+2x=0x′′+3x′+2x=0

是二阶方程。若最高出现的是 u′′′u'''u′′′,则是三阶。

不要把“导数的次数”与“导数整体的幂次”混在一起。方程 (y′)2+y=0(y')^2+y=0(y′)2+y=0 只出现一阶导数,因而仍是一阶方程;平方使它具有另一种结构上的复杂性,却没有把一阶导数变成二阶导数。

还要先看清是否有项可以消去,以及最高阶导数的系数在讨论区间内是否非零。例如 ty′′+y′=0t y''+y'=0ty′′+y′=0 在 t≠0t\ne0t=0 的区间上是通常的二阶方程;在 t=0t=0t=0 处,最高阶项的系数消失,不能直接除以 ttt 后假装什么都没发生。后面把方程整理成标准形式时,这类点会影响我们能使用的结论和解的区间。

线性:看未知函数如何参与运算

看看下面两个方程:

y′+t2y=sin⁡t,y′+y2=sin⁡t.y'+t^2y=\sin t,\qquad y'+y^2=\sin t.y′+t2y=sint,y′+y2=sint.

第一个有 t2t^2t2,却是线性的;第二个把 t2yt^2yt2y 换成 y2y^2y2,就变成非线性的。原因是线性所限制的对象是未知函数及其导数。自变量 ttt 是我们输入的数,t2t^2t2 可以当作已知系数;yyy 才是正在寻找的函数,把它平方会改变方程的结构。

一个 nnn 阶线性方程可以写成

an(t)y(n)+an−1(t)y(n−1)+⋯+a1(t)y′+a0(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),an​(t)y(n)+an−1​(t)y(n−1)+⋯+a1​(t)y′+a0​(t)y=g(t),

其中各个 aj(t)a_j(t)aj​(t) 和 g(t)g(t)g(t) 都是已知函数,所讨论的通常区间上要求 an(t)≠0a_n(t)\ne0an​(t)=0。未知函数及其导数只以一次方出现,彼此不相乘,也不被放进平方根、对数、正弦等函数里。

用几个例子把判断做实。y′′+(cos⁡t)y′+ety=t2y''+(\cos t)y'+e^t y=t^2y′′+(cost)y′+ety=t2 是二阶线性方程,系数复杂不妨碍线性。y′′+yy′=0y''+yy'=0y′′+yy′=0 是二阶非线性方程,因为 yyy 与 y′y'y′ 相乘;y′′+sin⁡y=0y''+\sin y=0y′′+siny=0 也是二阶非线性方程,因为正弦作用在未知函数上。y′=y(1−y)y'=y(1-y)y′=y(1−y) 展开后是 y′=y−y2y'=y-y^2y′=y−y2,所以是一阶非线性方程。

“线性”也不意味着解曲线必须是直线。刚刚见到的 y′=2yy'=2yy′=2y 是线性方程,它的解却是指数曲线。线性说的是方程怎样使用未知函数,不是在描述图像长得直不直。

线性与非线性微分方程对照图,左侧展示线性方程中 y 和 y' 只一次方且互不相乘,右侧展示 y²、sin y、yy' 等非线性形式。
线性方程只看未知函数及其导数如何出现;平方、相乘或套函数都会破坏线性。

齐次:在线性结构里看有没有独立输入

在线性方程中,把含 yyy 及其导数的项放在左边,剩下的已知项放在右边。若右边是零函数,就叫齐次线性方程;若不是零函数,就叫非齐次线性方程。例如

y′+2y=0y'+2y=0y′+2y=0

是齐次的,而

y′+2y=cos⁡ty'+2y=\cos ty′+2y=cost

是非齐次的。后一条方程的右边在某些时刻会等于零,仍然不改变分类,因为判断的是它是否在整个区间上恒等于零。

在模型里,这个独立项常常对应外部输入,例如持续供热或电源驱动。但分类首先是一件数学上的事,要根据选定变量下的方程判断,不能只凭“现实里有没有外力”来猜。后面一阶方程中还会出现另一个术语“齐次型方程”,含义不同,届时我们会重新说明。

下面的练习器可以帮助你练习分类。每次都先说出依据:最高导数是哪一个,未知函数有没有相乘或被套入其他函数,再看线性结构右端是否恒为零。


通解是一族可能,初始条件说明从哪里开始

前面验证过,y=Ce2ty=Ce^{2t}y=Ce2t 对每个常数 CCC 都满足 y′=2yy'=2yy′=2y。这一族覆盖了该方程所有解的表达式,叫作它的通解。任意选定一个常数,例如 C=4C=4C=4,就得到一个具体解,也常叫特解。

这里说“覆盖了所有解”不能只靠猜。对于任意一个满足 y′=2yy'=2yy′=2y 的可微函数,利用乘积求导法则,有

ddt(e−2ty(t))=e−2t(y′−2y)=0.\frac{d}{dt}\left(e^{-2t}y(t)\right) =e^{-2t}(y'-2y)=0.dtd​(e−2ty(t))=e−2t(y′−2y)=0.

因此 e−2ty(t)e^{-2t}y(t)e−2ty(t) 在所讨论的区间上是常数,记成 CCC,就得到 y=Ce2ty=Ce^{2t}y=Ce2t。这说明任何解都没有跑出这个函数族,而不只是说明这个函数族中的成员都是解。我们这里只借用了乘积求导,如何主动找到这样的乘子,会留到一阶线性方程中解释。

一个初始点怎样选定常数

假如我们要求这条曲线在 t=0t=0t=0 时取值为 444,就增加条件

y(0)=4.y(0)=4.y(0)=4.

代入通解得到 Ce0=C=4Ce^0=C=4Ce0=C=4,因此选出的函数是

y(t)=4e2t.y(t)=4e^{2t}.y(t)=4e2t.

把变化规律和起点合起来,写成

y′=2y,y(0)=4,y'=2y,\qquad y(0)=4,y′=2y,y(0)=4,

就是一个初值问题。它要求解同时满足两件事:沿途的变化率服从方程,在指定时刻经过指定状态。验证答案时,两件事都要检查。

初始时刻也不一定是零。如果题目给出 y(1)=6y(1)=6y(1)=6,由 Ce2=6Ce^2=6Ce2=6 得到 C=6e−2C=6e^{-2}C=6e−2,于是

y(t)=6e2(t−1).y(t)=6e^{2(t-1)}.y(t)=6e2(t−1).

写成这种形式,一眼就能看见它在 t=1t=1t=1 时等于 666。所谓“初始”,说的是我们用来指定状态的那个时刻,并不要求坐标原点一定放在那里。

图中展示的是初值能唯一确定解的典型情形。对刚才的指数方程,我们已经通过常数 CCC 的计算说明了唯一性;对任意微分方程,不能只凭“给了一个点”就断定一定有且只有一个解。有些方程没有满足指定初值的解,也有些会有多条。后面学存在唯一性时,我们会解释需要检查什么条件。

同样,看到一个含任意常数的解族,也别立刻认定它已经包含全部解。某些求解步骤会丢掉常数解,某些非线性方程还可能有不属于所求函数族的解。现在先养成一个习惯:区分“我找到了一族解”和“我证明这就是所有解”。

为什么二阶问题通常还要给初始速度

想象一个小车,在 t=0t=0t=0 时位于轨道的同一个位置。它可能静止,也可能正在向右运动,还可能正在向左运动。只给位置,显然不能区分接下来会发生什么。

用一个可以直接积分的模型算一遍。假设小车在一段时间里保持加速度 2 m/s22\,\mathrm{m/s^2}2m/s2,位置为 x(t)x(t)x(t),时间以秒计,那么

x′′=2.x''=2.x′′=2.

积分一次得到速度 x′=2t+C1x'=2t+C_1x′=2t+C1​,再积分得到位置

x=t2+C1t+C2.x=t^2+C_1t+C_2.x=t2+C1​t+C2​.

如果已知 x(0)=3 mx(0)=3\,\mathrm mx(0)=3m,只能确定 C2=3C_2=3C2​=3,C1C_1C1​ 仍然没有定。再给初始速度 x′(0)=−1 m/sx'(0)=-1\,\mathrm{m/s}x′(0)=−1m/s,才得到 C1=−1C_1=-1C1​=−1,从而

x(t)=t2−t+3.x(t)=t^2-t+3.x(t)=t2−t+3.

核对一下:x′′=2x''=2x′′=2,x(0)=3x(0)=3x(0)=3,x′(0)=−1x'(0)=-1x′(0)=−1,三项要求都满足。由于初始速度是负的,小车开始时向左走,虽然它的加速度始终向右;到 t=1/2t=1/2t=1/2 秒,速度变成零,随后才向右走。这也说明导数的正负必须对应到正确的物理量,正加速度不等于位置一定在增加。

一般的 nnn 阶初值问题,会在同一时刻指定函数值和前 n−1n-1n−1 阶导数。在适当条件下,这些信息能确定唯一解。上述二阶例子中,每积分一次出现一个常数,是这个规则最容易看见的起点。

未知数从数变成函数:代数方程得到两个数,微分方程得到指数函数族,初值选出其中一个函数。
变化率始终等于函数值的两倍,这一条件先给出一族指数函数;再指定初始时刻 t=0 的函数值为 3,就能确定其中唯一一个。

四个模型,把现象一步步写成方程

现在我们已经会认方程、验答案,也知道初值的作用。接下来把这些语言放回真实情境中。这里的参数都是为了练习设置的假想数值;每个模型都建立在明确假设上,我们先练习把假设翻译准确,再讨论什么时候它还适用。

人口增长:每个个体平均贡献多少新增数量

培养皿里细菌越多,短时间内发生分裂的机会通常也越多。我们先假设营养充足、环境保持稳定,没有迁入迁出,并用连续可微的函数近似大量个体的总数量。个体数本来是整数,把它连续化,是为了描述整体趋势的一个近似。

设 P(t)P(t)P(t) 为数量,时间以小时计。若每个个体平均每小时贡献的净增长比例是常数 rrr,那么在很短的 Δt\Delta tΔt 小时内,新增数量近似为 rP(t)ΔtrP(t)\Delta trP(t)Δt。除以 Δt\Delta tΔt,再用导数表示瞬时变化率,就得到

P′=rP.P'=rP.P′=rP.

P′P'P′ 的单位是“个每小时”,因此 rrr 的单位是“每小时”。这句话说的是相对增长率 P′/P=rP'/P=rP′/P=r 保持不变,并不是每小时增加的个数保持不变。当前数量翻倍时,当前增长速度也翻倍。

设 r=0.3 h−1r=0.3\,\mathrm{h}^{-1}r=0.3h−1,初始数量 P(0)=800P(0)=800P(0)=800。照着前面的指数函数验证思路,候选解是

P(t)=800e0.3t.P(t)=800e^{0.3t}.P(t)=800e0.3t.

求导得 P′=240e0.3t=0.3PP'=240e^{0.3t}=0.3PP′=240e0.3t=0.3P,初值也对。在开始时,增长速度是 240240240 个每小时;模型预测数量翻倍的时刻满足 e0.3t=2e^{0.3t}=2e0.3t=2,所以

t=ln⁡20.3≈2.31 小时.t=\frac{\ln2}{0.3}\approx2.31\text{ 小时}.t=0.3ln2​≈2.31 小时.

这里的 rrr 是净增长率,可以理解为平均出生率减去平均死亡率。若 r<0r<0r<0,同一种方程描述衰减;若 r=0r=0r=0,数量不变。人口模型通常只考虑 P≥0P\ge0P≥0。

恒定正增长率会使指数函数一直增大,真实培养皿却没有无限的营养和空间。若要描述资源压力,一种后续修正是让每个个体的净增长率随着数量增加而下降,例如引入承载量 K>0K>0K>0,写成

P′=rP(1−PK),r>0.P'=rP\left(1-\frac{P}{K}\right),\qquad r>0.P′=rP(1−KP​),r>0.

现在先读懂这句话就够了:当 PPP 远小于 KKK 时,括号接近 111,仍然近似指数增长;当 PPP 接近 KKK 时,增长逐渐变慢;P=KP=KP=K 时增长率为零。这就是后面要学习的逻辑斯蒂模型。我们会进一步求解,但此刻最该看清的是:换一条现实假设,方程的结构也会跟着改变。

咖啡冷却:决定快慢的是温差

一杯热咖啡放在房间里,通常刚开始降温较快,后来越来越慢。为了把这个现象写成方程,我们把咖啡看成内部温度均匀的整体,假设室温保持不变,而且这段温度范围内的换热效果可以用一个常数表示。

设咖啡温度为 T(t)T(t)T(t),环境温度为 TaT_aTa​,冷却系数为 k>0k>0k>0。先把“温差越大,变化越快”写成变化率的大小与 ∣T−Ta∣|T-T_a|∣T−Ta​∣ 成正比,再把变化方向放进去,就得到

T′=−k(T−Ta).T'=-k(T-T_a).T′=−k(T−Ta​).

负号的作用可以在两种情形下检查。咖啡比房间热时,T−Ta>0T-T_a>0T−Ta​>0,于是 T′<0T'<0T′<0,温度下降;一杯冷饮比房间冷时,T−Ta<0T-T_a<0T−Ta​<0,于是 T′>0T'>0T′>0,温度回升。同一个方程能描述冷却和回温,因为它让温度朝环境温度靠近。

设室温 Ta=24∘CT_a=24^\circ\mathrm CTa​=24∘C,咖啡初温 84∘C84^\circ\mathrm C84∘C,时间以分钟计,k=0.1 min−1k=0.1\,\mathrm{min}^{-1}k=0.1min−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)=−6T'(0)=-0.1(84-24)=-6T′(0)=−0.1(84−24)=−6,单位是摄氏度每分钟。它表示初始瞬间的降温速率,并不表示以后每分钟都降 666 度。

利用指数函数,我们可以检查候选解

T(t)=24+60e−0.1t.T(t)=24+60e^{-0.1t}.T(t)=24+60e−0.1t.

它在 t=0t=0t=0 时为 848484,导数是 −6e−0.1t-6e^{-0.1t}−6e−0.1t;右边 −0.1(T−24)-0.1(T-24)−0.1(T−24) 也等于 −6e−0.1t-6e^{-0.1t}−6e−0.1t,所以验证通过。十分钟后模型给出 T(10)=24+60/e≈46.07∘CT(10)=24+60/e\approx46.07^\circ\mathrm CT(10)=24+60/e≈46.07∘C。若问多久降到 54∘C54^\circ\mathrm C54∘C,就是温差从 606060 度变为 303030 度,于是

60e−0.1t=30,t=10ln⁡2≈6.93 分钟.60e^{-0.1t}=30, \qquad t=10\ln2\approx6.93\text{ 分钟}.60e−0.1t=30,t=10ln2≈6.93 分钟.

在这个理想模型里,只要初始温度不等于室温,温差就会在每个有限时刻保持非零,逐渐趋近于零。实际测量有分辨率,我们会在某个时候读到“已经是室温”,两句话并不冲突。若环境也被显著加热,或者咖啡内部温度并不均匀,就要调整原来的建模假设。

RC 电路:储存电荷让响应需要时间

电源接上以后,电容上的电压通常不会立刻跳到电源电压。原因是电容需要逐渐积累电荷,而电流正是电荷积累的速度。这个故事只要写清楚两个元件的关系,就会自然出现导数。

考虑电阻 R>0R>0R>0 与电容 C>0C>0C>0 串联的理想电路,设输入电压为 Vin(t)V_{\mathrm{in}}(t)Vin​(t),电容电压为 VC(t)V_C(t)VC​(t)。电容上的电荷满足 Q=CVCQ=CV_CQ=CVC​。假设电容值不随时间改变,电流就满足

I=Q′=CVC′.I=Q'=CV_C'.I=Q′=CVC′​.

电阻上的电压是 RIRIRI。沿回路按一致方向计量电压,电阻电压与电容电压之和等于输入电压,所以

RI+VC=Vin(t).RI+V_C=V_{\mathrm{in}}(t).RI+VC​=Vin​(t).

把 I=CVC′I=CV_C'I=CVC′​ 代进去,得到

RCVC′+VC=Vin(t).RCV_C'+V_C=V_{\mathrm{in}}(t).RCVC′​+VC​=Vin​(t).

这是一阶线性方程,未知函数是电容电压。若想得到电流,可以在得到 VCV_CVC​ 后再通过 I=CVC′I=CV_C'I=CVC′​ 计算。不要把“研究电路”自动等同于“未知量一定是电流”。

假设 R=1000 ΩR=1000\,\OmegaR=1000Ω,C=0.002 FC=0.002\,\mathrm FC=0.002F,则 RC=2 sRC=2\,\mathrm sRC=2s。从 t=0t=0t=0 起接入恒定 10 V10\,\mathrm V10V 电源,电容初始未充电,那么在接通后的时段,模型写为

2VC′+VC=10,VC(0)=0.2V_C'+V_C=10,\qquad V_C(0)=0.2VC′​+VC​=10,VC​(0)=0.

候选解 VC(t)=10(1−e−t/2)V_C(t)=10(1-e^{-t/2})VC​(t)=10(1−e−t/2) 在开始时确实是零,且

2VC′+VC=2⋅5e−t/2+10−10e−t/2=10.2V_C'+V_C=2\cdot5e^{-t/2}+10-10e^{-t/2}=10.2VC′​+VC​=2⋅5e−t/2+10−10e−t/2=10.

所以它满足模型。电流为 I=0.002⋅5e−t/2=0.01e−t/2 AI=0.002\cdot5e^{-t/2}=0.01e^{-t/2}\,\mathrm AI=0.002⋅5e−t/2=0.01e−t/2A,开始时最大,随后衰减。电容电压逐渐接近 10 V10\,\mathrm V10V,电阻两端的压差越来越小,电流自然也越来越小。

把方程整理成 VC′=(10−VC)/2V_C'=(10-V_C)/2VC′​=(10−VC​)/2,你会发现它与冷却模型很相似:变化率由“目标值与当前值的差”决定。这个相似性会在后面的解法中得到统一解释。这里的 RCRCRC 具有时间单位,指数中的 t/(RC)t/(RC)t/(RC) 因而没有单位,这也是检查公式是否合理的一条线索。

弹簧振子:位置和速度一起决定加速度

把一个小物块连在水平弹簧上,拉开一点再松手。物块离平衡位置越远,弹簧往回拉的力通常越大;若还有阻尼,运动越快,阻碍运动的力也越大。两种力作用在一起,决定的是加速度,因此方程会出现二阶导数。

设 x(t)x(t)x(t) 是相对平衡位置的位移,向右为正。质量为 m>0m>0m>0,弹簧刚度为 k>0k>0k>0,阻尼系数为 c≥0c\ge0c≥0。这里的 kkk 表示弹簧刚度,和前面冷却模型中的 kkk 是不同参数,只是各自在本模型中使用的常见记号。

在弹簧伸缩不大、回复力可以近似与位移成正比的范围内,弹簧力是 −kx-kx−kx。若阻尼力近似与速度成正比且方向相反,则为 −cx′-cx'−cx′。再设外部驱动力为 F(t)F(t)F(t),由“质量乘加速度等于合力”得到

mx′′=−cx′−kx+F(t),mx''=-cx'-kx+F(t),mx′′=−cx′−kx+F(t),

也就是

mx′′+cx′+kx=F(t).mx''+cx'+kx=F(t).mx′′+cx′+kx=F(t).

每一项都具有力的单位:mmm 用千克,xxx 用米,ttt 用秒时,kkk 的单位是牛顿每米,ccc 的单位是牛顿秒每米。若用竖直弹簧并从重力作用下的平衡位置测位移,重力与静态伸长产生的弹簧力抵消,也会得到相同形式;选对位置原点,能让方程更清楚。

先取一个没有阻尼、也没有外力的例子。设 m=1 kgm=1\,\mathrm{kg}m=1kg、k=4 N/mk=4\,\mathrm{N/m}k=4N/m,物块初始向右偏离 0.1 m0.1\,\mathrm m0.1m 后静止释放,则

x′′+4x=0,x(0)=0.1,x′(0)=0.x''+4x=0,\qquad x(0)=0.1,\qquad x'(0)=0.x′′+4x=0,x(0)=0.1,x′(0)=0.

我们不急着推导全部解法,先验证 x(t)=0.1cos⁡(2t)x(t)=0.1\cos(2t)x(t)=0.1cos(2t)。求导两次得到

x′(t)=−0.2sin⁡(2t),x′′(t)=−0.4cos⁡(2t).x'(t)=-0.2\sin(2t),\qquad x''(t)=-0.4\cos(2t).x′(t)=−0.2sin(2t),x′′(t)=−0.4cos(2t).

因此 x′′+4x=0x''+4x=0x′′+4x=0,初始位置和初始速度也都正确。这个解描述来回振动,周期为 π\piπ 秒。它一直保持相同振幅,是因为模型特意忽略了阻尼。加入阻尼后振幅怎样变化,加入周期外力后会不会越振越大,将是二阶方程部分要回答的问题。


从文字到答案,中间还有几次检查

四个模型的故事不同,写方程时的思考顺序却可以重复使用。先选定正在变化的量,说明它的输入、输出和单位;接着把影响变化的因素说清楚,再把这句话变成导数关系。然后加入起始状态,最后把方程与候选解放回现象中,检查符号、单位和适用范围。

ODE 建模流程闭环图,包含选变量和单位、写变化率假设、给初始条件、解释解的趋势,中间为 dy/dt=f(t,y)。
从现象出发建立一阶微分方程模型的基本闭环。

比如一段文字说“数量每小时增加 202020 个”,它对应 P′=20P'=20P′=20;说“每小时净增长率为当前数量的 20%20\%20%”,在连续变化率的解释下对应 P′=0.2PP'=0.2PP′=0.2P,其中 0.20.20.2 的单位是每小时。前者每小时增加同样多,后者当前数量越大就增加越多。一句话中的“个”和“比例”,会让模型从直线变化转向指数变化。

再比如小车匀加速模型 x=t2−t+3x=t^2-t+3x=t2−t+3 在数学上对所有实数 ttt 有定义,但如果题目说小车只在前五秒保持恒定加速度,那么它对现实的预测范围就只能先放在 0≤t≤50\le t\le50≤t≤5。函数的数学定义域、微分方程解的区间、现实假设的有效时间范围,彼此相关,却不能混为一谈。

我们也不必要求每个方程都立刻写出一个漂亮公式。一个模型建立以后,可以先从右端正负判断增加还是减少,从右端为零寻找可能不变的状态,再想办法求解或近似计算。下一章的方向场就是沿着这个想法,把许多位置上的斜率画出来,让局部规则逐渐显出整条曲线的走向。


练习:把定义用到具体问题上

验证函数,而不是只看外形

判断 y1(t)=3+2e−4ty_1(t)=3+2e^{-4t}y1​(t)=3+2e−4t 和 y2(t)=3+2e−2ty_2(t)=3+2e^{-2t}y2​(t)=3+2e−2t 是否满足 y′+4y=12y'+4y=12y′+4y=12。若满足,再检查初值 y(0)=5y(0)=5y(0)=5,并说明解的区间。

对第一个函数,y1′=−8e−4ty_1'=-8e^{-4t}y1′​=−8e−4t,因此

y1′+4y1=−8e−4t+12+8e−4t=12.y_1'+4y_1=-8e^{-4t}+12+8e^{-4t}=12.y1′​+4y1​=−8e−4t+12+8e−4t=12.

同时 y1(0)=3+2=5y_1(0)=3+2=5y1​(0)=3+2=5,所以它是所给初值问题的解,最大解区间为 (−∞,∞)(-\infty,\infty)(−∞,∞)。第二个函数虽然同样满足初值,但 y2′=−4e−2ty_2'=-4e^{-2t}y2′​=−4e−2t,代入后得到 y2′+4y2=12+4e−2ty_2'+4y_2=12+4e^{-2t}y2′​+4y2​=12+4e−2t,并不恒等于 121212,因此不是方程的解。满足初始点,不能代替满足沿途的变化规律。

分类时盯住未知函数

判断下面四个方程的阶数与线性;对线性方程再判断齐次性:y′′+t3y′=ety''+t^3y'=e^ty′′+t3y′=et,(y′)3+y=t(y')^3+y=t(y′)3+y=t,y′′+sin⁡(t)y=0y''+\sin(t)y=0y′′+sin(t)y=0,y′′+sin⁡(y)=0y''+\sin(y)=0y′′+sin(y)=0。

第一个最高导数是 y′′y''y′′,是二阶线性非齐次方程;t3t^3t3 是已知系数,不破坏线性,右端 ete^tet 不恒为零。第二个只含一阶导数,所以是一阶方程,但 (y′)3(y')^3(y′)3 使它非线性,不能在此套用线性方程的齐次分类。第三个是二阶线性齐次方程,sin⁡t\sin tsint 只是 yyy 的已知系数。第四个是二阶非线性方程,正弦作用在未知函数上。两条方程只改变了正弦括号里的对象,线性结构就不同了。

初始时刻可以不在零点

求解可以直接积分的初值问题 y′=6ty'=6ty′=6t,y(2)=7y(2)=7y(2)=7。再说明它的图像与通解之间的关系。

对 6t6t6t 积分,通解为 y=3t2+Cy=3t^2+Cy=3t2+C。由 y(2)=12+C=7y(2)=12+C=7y(2)=12+C=7 得 C=−5C=-5C=−5,所以

y(t)=3t2−5.y(t)=3t^2-5.y(t)=3t2−5.

求导得到 y′=6ty'=6ty′=6t,再代入 t=2t=2t=2 得到 777,两项都满足。通解是一族上下平移的抛物线,初始条件要求其中一条经过 (2,7)(2,7)(2,7),因而固定了平移量。所求解的最大区间为整个实数轴。

同一个表达式,初值决定讨论哪一段

验证 y(t)=1/(3−t)y(t)=1/(3-t)y(t)=1/(3−t) 满足 y′=y2y'=y^2y′=y2。若给定 y(0)=1/3y(0)=1/3y(0)=1/3,最大解区间是什么?若改为 y(4)=−1y(4)=-1y(4)=−1 呢?

求导得到 y′=1/(3−t)2y'=1/(3-t)^2y′=1/(3−t)2,与 y2y^2y2 相等,验证在 t≠3t\ne3t=3 时成立。第一个初值确实满足,包含 000 的最大开区间是 (−∞,3)(-\infty,3)(−∞,3)。第二个初值也满足,但包含 444 的最大开区间是 (3,∞)(3,\infty)(3,∞)。这两个区间不能跨过 t=3t=3t=3 拼接,因为函数在那里无定义,而且从两侧靠近时都没有有限极限。

从回温故事写出初值问题

一杯温度为 8∘C8^\circ\mathrm C8∘C 的饮料放进 26∘C26^\circ\mathrm C26∘C 的房间,假设饮料内部温度均匀、室温不变,温度变化率与温差成正比,系数为 0.05 min−10.05\,\mathrm{min}^{-1}0.05min−1。写出初值问题,求初始变化率,并验证候选函数 T(t)=26−18e−0.05tT(t)=26-18e^{-0.05t}T(t)=26−18e−0.05t。

令 T(t)T(t)T(t) 为 ttt 分钟后的温度,初值问题为

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.9T'(0)=-0.05(8-26)=0.9T′(0)=−0.05(8−26)=0.9,单位是摄氏度每分钟,正号表示回温。候选函数的导数为 0.9e−0.05t0.9e^{-0.05t}0.9e−0.05t;代入方程右端得到 −0.05(−18e−0.05t)=0.9e−0.05t-0.05(-18e^{-0.05t})=0.9e^{-0.05t}−0.05(−18e−0.05t)=0.9e−0.05t,并且 T(0)=8T(0)=8T(0)=8,所以通过验证。它在 t≥0t\ge0t≥0 时从下方靠近 26∘C26^\circ\mathrm C26∘C,温度上升但上升速度逐渐减小。

两次积分与两条初始信息

假想一辆小车在研究时段内满足 x′′=−4x''=-4x′′=−4,时间以秒计,位置以米计,初始位置 x(0)=2x(0)=2x(0)=2,初始速度 x′(0)=6x'(0)=6x′(0)=6。求位置函数,并求它第一次速度为零的时刻及此时的位置。

第一次积分得到 x′=−4t+C1x'=-4t+C_1x′=−4t+C1​,由初始速度得 C1=6C_1=6C1​=6。第二次积分得到 x=−2t2+6t+C2x=-2t^2+6t+C_2x=−2t2+6t+C2​,由初始位置得 C2=2C_2=2C2​=2,因此

x(t)=−2t2+6t+2.x(t)=-2t^2+6t+2.x(t)=−2t2+6t+2.

令速度 −4t+6=0-4t+6=0−4t+6=0,得 t=1.5t=1.5t=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.50<t<1.50<t<1.5 时速度仍为正,小车仍朝正方向运动,只是越来越慢。核对得到 x′′=−4x''=-4x′′=−4 且两个初值都满足,结果只用于题设恒加速度假设成立的时段。

做完这些练习,再回到一开始的 y′=2yy'=2yy′=2y:你现在能说明它在找什么,怎样检查一个答案,为什么会有任意常数,以及初始点怎样进入问题。接下来我们会把“每个状态给出一个斜率”真的画成图。即使一时求不出函数公式,也可以先看清曲线从哪里出发、向哪里走。

下一章建模、初值问题与解的几何图像