如何选择 first-order 方法

first-order ODE 含有 first derivative $y'$，但不含 $y''$ 这样的更高阶 derivative。通常写成

$$y'=\frac{dy}{dt}=f(t,y),$$

其中 $t$ 是 independent variable，$y=y(t)$ 是我们要找的函数。

先别急着做 integration，先看方程的结构：能不能把变量分到两边？能不能把左边凑成一个乘积的 derivative？或者，能不能找到一个沿解曲线保持不变的函数？没有一种方法能解所有 first-order equations，这里主要用下面三种形式：

同一个方程可能同时符合不止一种形式。例如，一个方程既是 separable 又是 linear，两种方法都能用，选计算更简单的即可。

Separable equations

如果能把 first-order equation 整理成一边只含 $y$、另一边只含 $t$ 的形式，就称它为 separable。

常见形式是

$$\frac{dy}{dt}=f_1(t)f_2(y),$$

其中 $f_1$ 只依赖 $t$，$f_2$ 只依赖 $y$。

例如，

$$y'=e^{t+y}=e^t e^y$$

是 separable，因为右边可以拆成只含 $t$ 和只含 $y$ 的两个因子的乘积。而方程

$$y'=t^2+y^2$$

在这个形式下不是 separable，因为右边是两个项相加，而不是这样的乘积。

先检查 equilibrium solutions

分离变量时，需要除以 $f_2(y)$。这个因子为零时不能做除法，但对应的函数可能完全符合原方程，所以要先检查。

在除法之前，先解

$$f_2(y)=0.$$

如果 $y=y_*$ 使得 $f_2(y_*)=0$，那么 constant function $y(t)=y_*$ 满足

$$y'(t)=0=f_1(t)f_2(y_*).$$

这个 constant function 就是 equilibrium solution。符号 $y_*$ 读作 “$y$ star”。

不过，不能因此断定所有到达这个值的解都会一直停在那里。要排除“到达后又离开”的解，还需要 uniqueness theorem。

分离变量并做 integration

在 $f_2(y)\ne0$ 的地方，两边除以它，然后做 integration：

$$\boxed{ \int\frac{1}{f_2(y)}\,dy = \int f_1(t)\,dt+C. }$$

最后得到的式子里，$y$ 和 $t$ 可能仍然混在一起。这叫 implicit solution：它描述函数满足的关系，但没有把函数显式解出来。同一条曲线可能有多个 branch；必须选取可以求 derivative 的 branch，并限制在原 ODE 和刚才的除法都有意义的 interval 上。

Separation 为什么成立

写成

$$\frac{dy}{f_2(y)}=f_1(t)\,dt$$

看起来像是把 $dy/dt$ 当作普通分数来处理。这种写法方便记步骤，但真正支持它的是 chain rule。

给两边的 antiderivatives 分别命名为 $F$ 和 $G$，使它们满足

$$F'(y)=\frac{1}{f_2(y)}, \qquad G'(t)=f_1(t).$$

由于 $y$ 依赖 $t$，chain rule 给出

$$\frac{d}{dt}F(y(t)) =F'(y)y' =\frac{y'}{f_2(y)}.$$

原 ODE 告诉我们，这等于 $f_1(t)=G'(t)$。因此，

$$\frac{d}{dt}\bigl(F(y(t))-G(t)\bigr)=0.$$

一个量的 derivative 为零，就说明它是常数，所以

$$\boxed{F(y)-G(t)=C}.$$

这正是分离变量再做 integration 得到的 implicit equation。

Implicit solutions 与 level curves

刚才找到的表达式沿着每条解曲线保持不变。把它记作

$$\psi(t,y)=F(y)-G(t),$$

$\psi$ 读作 “psi”，这里称它为分离后方程的 potential function。ODE 告诉我们曲线的 slope，而这个函数告诉我们：沿着一条解曲线，哪个量保持不变。

$$\psi(t,y)=C.$$

如果一条曲线上每一点的 $\psi$ 值都相同，这条曲线就叫 level curve。每个解沿着其中一条曲线变化；选择不同的 $C$，就是选择不同的 level。

Level curves and the selected solution

例：圆形的 implicit solutions

考虑

$$\frac{dy}{dx}=-\frac{x}{y}.$$

分离变量并做 integration：

$$\begin{aligned} y\,dy&=-x\,dx,\\ \frac{y^2}{2}&=-\frac{x^2}{2}+C. \end{aligned}$$

把平方项移到同一边：

$$\boxed{\frac{x^2}{2}+\frac{y^2}{2}=C}.$$

当常数为正时，这个方程描述一个以原点为圆心的圆。但整个圆并不是横坐标的函数，所以还要选 branch。假设给定

$$y(0)=1.$$

代入 $x=0$ 和 $y=1$，得到 $C=1/2$，所以

$$x^2+y^2=1.$$

解出 $y$，得到两个 branches：

$$y=\pm\sqrt{1-x^2}.$$

我们需要 $y(0)=1$，因此选上半圆：

$$\boxed{y(x)=\sqrt{1-x^2}}.$$

这个解能延伸多远？maximal solution interval 是包含初始点、且这个公式始终满足原 ODE 的最大 interval。这里是

$$-1<x<1.$$

在 $x=\pm1$ 处，square root 仍然有实数值，但此时 $y=0$，而原 ODE 要除以 $y$。因此这两个端点不能包含在 solution interval 里。

例：logistic growth

设 $N(t)$ 是时刻 $t$ 的种群数量，$N'(t)$ 表示总数量变化的快慢。

一般的种群模型可以写成

$$\frac{dN}{dt}=R(N)N.$$

这里的 $R(N)$ 是平均每个个体对应的增长率，叫作 per-capita growth rate。乘以个体数量 $N$，就得到总数量的变化率。

Exponential growth 假设 $R(N)$ 不变。在食物和空间充足时，这可能合理，但资源不足时就不合适了。为了描述资源的限制，引入两个常数：

• $r>0$，表示种群数量很小时的 per-capita growth rate；
• $K>0$，称为 carrying capacity，表示环境能长期维持的种群数量。

因此，$R(N)$ 应满足

$$R(0)=r>0, \qquad R(K)=0.$$

对正的种群数量，$N<K$ 时应该增长，$N>K$ 时应该减少。让 per-capita rate 沿一条直线从 $r$ 降到 $0$，就能得到这种变化：

$$R(N)=r\left(1-\frac{N}{K}\right) =\frac{r}{K}(K-N).$$

检查 $K-N$ 的正负，就能判断种群数量如何变化：

• 若 $0<N<K$，则 $K-N>0$，种群数量增加。
• 若 $N=K$，则 $K-N=0$，种群数量不变。
• 若 $N>K$，则 $K-N<0$，种群数量减少。

carrying capacity 是一个 equilibrium，不是数量绝不能超过的上限。种群可以从高于 $K$ 的数量开始；这时模型预测它会逐渐降向 $K$。

把这个 $R(N)$ 代回种群模型，得到 logistic equation：

$$\boxed{\frac{dN}{dt}=\frac{r}{K}(K-N)N}.$$

这个方程是 first order，也是 separable。它是 nonlinear 的，因为 $(K-N)$ 与 $N$ 相乘后会出现 $N^2$ 项。

Equilibria 与变化方向

在 equilibrium solution 上，变化率为零。令 logistic equation 的右边等于零：

$$\frac{r}{K}(K-N)N=0.$$

得到两个 equilibrium populations：

$$N=0 \qquad\text{and}\qquad N=K.$$

从任何一个 equilibrium 开始，数量都会保持不变。从其他任意正值开始，$N'$ 的符号都表明数量会朝 $K$ 变化。附近的种群数量会靠近而不是远离这个 equilibrium，因此称它为 stable equilibrium。

Per-capita growth rate

Growth is positive below the carrying capacity K, zero at K, and negative above K.

Population trajectories

K is a stable equilibrium: every positive trajectory shown approaches the carrying capacity.

求解

现在已经知道种群会朝哪里变化，接下来求任意时刻的数量。设 $N(0)=N_0$，且 $N_0\notin\{0,K\}$，其中 $N_0$ 是初始种群数量。

方程右边是 polynomial，因而是 smooth 的，可以应用 uniqueness theorem。这保证了这些解在存在期间不会穿过任何一个 equilibrium。因此 $N\ne0$ 且 $N\ne K$，可以除过去并分离变量：

$$\frac{K}{N(K-N)}\,dN=r\,dt.$$

把左边拆成 partial fractions，会更容易做 integration：

$$\frac{K}{N(K-N)} =\frac1N+\frac1{K-N}.$$

这两项的 antiderivatives 都涉及 logarithms。注意第二项的符号：$K-N$ 对 $N$ 的 derivative 是 $-1$，所以对应的 integral 前面要有负号。

$$\begin{aligned} \int\left(\frac1N+\frac1{K-N}\right)dN &=\ln|N|-\ln|K-N|\\ &=\ln\left|\frac{N}{K-N}\right|. \end{aligned}$$

右边也做 integration，把两边合在一起：

$$\ln\left|\frac{N}{K-N}\right|=rt+C.$$

两边取 exponential 来消去 logarithm。把正负号和 $e^C$ 合并到一个非零常数 $C_1$ 中：

$$\frac{N}{K-N}=C_1e^{rt}.$$

为了解出 $N$，先把两边乘以 $K-N$：

$$N=C_1e^{rt}(K-N).$$

展开右边，再把含 $N$ 的项收集到左边：

$$N(1+C_1e^{rt})=KC_1e^{rt}.$$

用 $N(0)=N_0$ 求 $C_1$：

$$C_1=\frac{N_0}{K-N_0}.$$

代回 $C_1$，解出种群数量：

$$\boxed{ N(t)=\frac{KN_0}{N_0+(K-N_0)e^{-rt}} }.$$

推导时除以了 $N(K-N)$，所以没有包含两个 equilibria。答案里还要保留它们：

$$N(t)\equiv0 \qquad\text{and}\qquad N(t)\equiv K.$$

再把公式与之前根据符号作出的判断比较：

• 若 $0<N_0<K$，种群数量增加并趋向 $K$。
• 若 $N_0=K$，数量一直为 $K$。
• 若 $N_0>K$，种群数量减少并趋向 $K$。

Linear equations

在 linear equation 中，未知函数及其 first derivative 都只以一次幂出现。它们的 coefficients 可以依赖 independent variable，但不能依赖未知函数。题目可能给出这样的形式：

$$a_0(t)y'(t)+a_1(t)y(t)=h(t).$$

先把每一项除以 $a_0(t)$，并限制在 $a_0(t)\ne0$ 的 interval 上，保证除法合法。得到

$$y'(t)+P(t)y(t)=g(t),$$

其中

$$P(t)=\frac{a_1(t)}{a_0(t)}, \qquad g(t)=\frac{h(t)}{a_0(t)}.$$

这就是 standard form。找 integrating factor 之前，必须先让 $y'$ 的 coefficient 变成 $1$。

Integrating factor 的作用

standard form 的左边是

$$y'+P(t)y.$$

第二项里的函数还不知道，不能把它当作已知函数直接做 integration。我们要先把整个左边变成熟悉的 derivative。

把每一项都乘以一个函数 $\mu(t)$，读作 “mu of $t$”。目标是让

$$\mu y'+\mu P y$$

变成 $\mu y$ 的 derivative。

根据 product rule，

$$\frac{d}{dt}(\mu y)=\mu y'+\mu' y.$$

第一项已经对应上了。要让第二项也对应上，需要

$$\mu'=P\mu.$$

这样就得到一个关于 multiplier 本身的方程，而且它是 separable。除以非零的 multiplier，做 integration，再取 exponential，就能选取

$$\boxed{\mu(t)=e^{\int P(t)\,dt}}.$$

这就是 integrating factor。exponent 里不必写 $+C$：它只会让 $\mu$ 多乘一个非零常数，最后解出 $y$ 时会约掉。

选好这个 $\mu$ 后，product rule 就把乘过 multiplier 的方程变成

$$\frac{d}{dt}(\mu y)=\mu g.$$

现在两边做 integration，再除以 $\mu$：

$$\boxed{ y(t)=\frac{1}{\mu(t)} \left(\int \mu(t)g(t)\,dt+C\right). }$$

最后这个公式不必硬背。按“乘上 integrating factor、用 product rule、做 integration、再除回去”的顺序，就能重新推出来。

有些 WeBWorK 题目用 $u(t)$ 表示 integrating factor 的 reciprocal：

$$u(t)=\frac{1}{\mu(t)}.$$

这不是两种不同的方法。乘方程时仍然用 $\mu$；最后解出 $y$ 时，才会出现 reciprocal $u$。填写答案前，注意题目究竟要哪一个。

求解步骤

对 linear equation，可以按下面的顺序做：

1. 先除过去，使 $y'$ 的 coefficient 为 $1$。
2. 对照 $y'+P(t)y=g(t)$，确定 $P(t)$ 和 $g(t)$。
3. 求出 $\mu(t)=e^{\int P(t)\,dt}$。
4. 每一项都乘以 $\mu(t)$。
5. 把左边写成 $(\mu y)'$。
6. 两边做 integration，在这一步加入常数 $C$。
7. 除以 $\mu(t)$。
8. 如果给了 initial condition，用它确定 $C$。

检查答案时，先求 derivative，再代回原 ODE。

例：variable coefficient

考虑

$$y'+\frac{2t}{4+t^2}y=\frac{4t}{4+t^2}.$$

方程已经是 standard form，因此

$$P(t)=\frac{2t}{4+t^2}.$$

integrating factor 为

$$\begin{aligned} \mu(t) &=e^{\int \frac{2t}{4+t^2}\,dt}\\ &=e^{\ln(4+t^2)}\\ &=4+t^2. \end{aligned}$$

$4+t^2$ 始终为正，不会出现需要避开的零点，也不需要分 interval 讨论。每一项乘以 $4+t^2$：

$$(4+t^2)y'+2ty=4t.$$

左边正好成为我们需要的 derivative：

$$\frac{d}{dt}\bigl((4+t^2)y\bigr)=4t.$$

做 integration 得到

$$(4+t^2)y=2t^2+C,$$

因此 general solution 是

$$\boxed{y(t)=\frac{2t^2+C}{4+t^2}}.$$

不同的 $C$ 对应不同的解，initial condition 会从中选出一个。

例如，若 $y(0)=0$，则

$$0=\frac{C}{4},$$

所以 $C=0$，得到

$$y(t)=\frac{2t^2}{4+t^2}.$$

若 $y(0)=1$，则 $C=4$，得到

$$y(t)=\frac{2t^2+4}{4+t^2}.$$

例：exponential input

考虑

$$y'+\frac12y=\frac12e^{t/3}.$$

这里 $P(t)=1/2$，因此

$$\mu(t)=e^{\int 1/2\,dt}=e^{t/2}.$$

整个方程乘以 $e^{t/2}$：

$$\frac{d}{dt}\left(e^{t/2}y\right) =\frac12e^{t/2}e^{t/3} =\frac12e^{5t/6}.$$

做 integration：

$$\begin{aligned} e^{t/2}y &=\frac12\int e^{5t/6}\,dt+C\\ &=\frac35e^{5t/6}+C. \end{aligned}$$

除以 $e^{t/2}$，得到

$$\boxed{y(t)=\frac35e^{t/3}+Ce^{-t/2}}.$$

第一项是对应这个 input 的一个解，$Ce^{-t/2}$ 则让我们能调整初值而不改变 input。它满足把右边换成零后得到的方程，也就是 homogeneous equation：

$$y'+\frac12y=0.$$

例：下落物体

假设一个物体在重力作用下下落，air resistance 阻碍它的运动。取向下为正方向，设：

• $v(t)$ 是时刻 $t$ 的向下速度；
• $m>0$ 是物体质量；
• $g>0$ 是 gravitational acceleration；
• $\gamma>0$ 读作 “gamma”，是 linear drag coefficient。

重力贡献 $+mg$。这里把 drag 建模为 $-\gamma v$，所以它始终与 velocity 的方向相反。Newton's second law 给出

$$mv'(t)=mg-\gamma v(t).$$

两边除以 $m$，得到 acceleration：

$$\boxed{v'(t)=g-\frac{\gamma}{m}v(t)}.$$

物体向下运动时，可以比较两个力的大小。若重力大于 drag，则 $v'>0$，物体加速；若 drag 大于重力，则 $v'<0$，物体减速。

Terminal velocity

acceleration 为零（$v'=0$）时，velocity 不再改变。令右边为零：

$$0=g-\frac{\gamma}{m}v,$$

得到

$$\boxed{v_T=\frac{mg}{\gamma}},$$

这里的 $v_T$ 就是 terminal velocity。

把 $\frac{\gamma}{m}$ 提出来，用 $v_T$ 改写方程：

$$v'=\frac{\gamma}{m}(v_T-v).$$

这样就可以直接判断变化方向：

• 若 $v<v_T$，则 $v'>0$，velocity 增大。
• 若 $v>v_T$，则 $v'<0$，velocity 减小。
• 若 $v=v_T$，则 $v'=0$，velocity 保持不变。

无论哪种情况，velocity 都趋向或保持在 $v_T$。

求解

为了求 $v(t)$，先把方程写成 standard first-order linear form：

$$v'+\frac{\gamma}{m}v=g.$$

integrating factor 是

$$\mu(t)=e^{\gamma t/m}.$$

每一项乘以 $\mu(t)$，左边就成为乘积的 derivative：

$$\frac{d}{dt}\left(e^{\gamma t/m}v(t)\right) =ge^{\gamma t/m}.$$

两边做 integration：

$$e^{\gamma t/m}v(t) =\frac{mg}{\gamma}e^{\gamma t/m}+C.$$

设初速度为 $v(0)=v_0$。代入 $t=0$，求出常数：

$$C=v_0-\frac{mg}{\gamma}.$$

再乘以 $e^{-\gamma t/m}$，并代入 $C$：

$$\boxed{ v(t)=\frac{mg}{\gamma} +\left(v_0-\frac{mg}{\gamma}\right)e^{-\gamma t/m} }.$$

可以把这个答案看成 equilibrium 值，再加上随时间变化的差值：

• $mg/\gamma$ 是稳定的 terminal velocity；
• 另一项是初速度与 terminal velocity 的差，乘以一个不断衰减的 exponential。

因为 $m>0$ 且 $\gamma>0$，exponential 项会衰减：

$$e^{-\gamma t/m}\longrightarrow0 \qquad\text{as }t\to\infty.$$

差值趋向零时，velocity 就趋向 terminal velocity：

$$\boxed{ \lim_{t\to\infty}v(t)=\frac{mg}{\gamma} }.$$

如果物体从静止释放，令 $v_0=0$，答案简化为

$$v(t)=\frac{mg}{\gamma} \left(1-e^{-\gamma t/m}\right).$$

velocity 从零开始平滑上升，趋向 terminal velocity，但不会超过它。

Velocity under linear drag

例：mixing tank

一个水箱起初装有 $100$ gallons 充分搅匀的溶液，其中含盐 $Q_0$ pounds。液体不断流入和流出，同时带进和带走盐。设：

• $Q(t)$ 是经过 $t$ minutes 后的盐量，单位为 pounds；
• $r>0$ 是流入和流出的 flow rate，两者相同，单位为 gallons per minute；
• 流入溶液的盐浓度为 $\frac14$ lb/gal。

流入与流出的液体量相同，因此水箱始终有 $100$ gallons。盐量的变化率等于流入的盐量速率减去流出的盐量速率：

$$\text{rate of change}=\text{rate in}-\text{rate out}.$$

流入的盐量速率等于浓度乘以 flow rate：

$$\left(\frac14\frac{\mathrm{lb}}{\mathrm{gal}}\right) \left(r\frac{\mathrm{gal}}{\mathrm{min}}\right) =\frac r4\frac{\mathrm{lb}}{\mathrm{min}}.$$

由于充分搅拌，盐均匀分布在溶液中。因此，流出液体的浓度与水箱中液体的浓度相同：

$$\frac{Q(t)}{100}\frac{\mathrm{lb}}{\mathrm{gal}}.$$

所以，流出的盐量速率是

$$\left(\frac{Q(t)}{100}\frac{\mathrm{lb}}{\mathrm{gal}}\right) \left(r\frac{\mathrm{gal}}{\mathrm{min}}\right) =\frac r{100}Q(t)\frac{\mathrm{lb}}{\mathrm{min}}.$$

两者相减，得到 initial-value problem：

$$\boxed{ Q'+\frac r{100}Q=\frac r4, \qquad Q(0)=Q_0 }.$$

检查单位：每一项都应是 pounds per minute。这里描述的是盐量变化的快慢，不是浓度变化的快慢。

求解

这已经是 linear standard form：derivative 的 coefficient 为一，盐量只以一次幂出现。integrating factor 是

$$\mu(t)=e^{\int r/100\,dt}=e^{rt/100}.$$

把方程乘以 $\mu(t)$，左边变成乘积的 derivative：

$$\frac d{dt}\left(e^{rt/100}Q(t)\right) =\frac r4e^{rt/100}.$$

两边做 integration：

$$\begin{aligned} e^{rt/100}Q(t) &=\frac r4\int e^{rt/100}\,dt+C\\ &=25e^{rt/100}+C. \end{aligned}$$

除以 integrating factor：

$$Q(t)=25+Ce^{-rt/100}.$$

再用 initial condition 确定常数：

$$Q_0=Q(0)=25+C, \qquad C=Q_0-25.$$

于是

$$\boxed{ Q(t)=25+(Q_0-25)e^{-rt/100} }.$$

这个 equilibrium 也可以从浓度来理解。当水箱与流入液体的浓度相同时，盐流入和流出的速率相等。此时箱中的盐量应为

$$100\,\mathrm{gal}\left(\frac14\frac{\mathrm{lb}}{\mathrm{gal}}\right) =25\,\mathrm{lb}.$$

因为 $r>0$，exponential 因子趋向零。若 $Q_0<25$，盐量增加并趋向 $25$；若 $Q_0>25$，盐量减少并趋向 $25$；若 $Q_0=25$，盐量就保持不变。

不追踪总盐量，而是追踪它与 equilibrium 的差，也能看出同样的变化。定义 $E(t)=Q(t)-25$，模型变成

$$E'=-\frac{r}{100}E, \qquad E(t)=E(0)e^{-rt/100}.$$

如果起初不在 equilibrium，$E(t)/E(0)=e^{-rt/100}$ 就表示初始差值还剩下多少比例。如果起初就在 equilibrium，差值会一直为零。

Stable salt equilibrium at 25 pounds

• Q(0) = 50 lb
• Q(0) = 25 lb
• Q(0) = 10 lb

达到目标需要多久

取 $Q_0=50$ lb、$r=3$ gal/min，则

$$Q(t)=25+25e^{-3t/100}.$$

与 $25$-lb equilibrium 相差不超过 $2\%$，意味着

$$|Q(t)-25|\le0.02(25)=0.5.$$

盐量从上方逐渐降向 equilibrium。第一次达到要求时，$Q(T)=25.5$，其中 $T$ 就是所求时间。解得：

$$\begin{aligned} 25+25e^{-3T/100}&=25.5\\ e^{-3T/100}&=\frac1{50}\\ T&=\frac{100}{3}\ln50\\ &\approx130.4\ \mathrm{min}. \end{aligned}$$

如果要求在 $45$ minutes 内达到这个范围呢？这次时间已知，要找的是 flow rate。恰好够用的最小 flow rate 会使盐量在最后一刻达到目标：

$$25+25e^{-45r/100}=25.5.$$

因此，

$$\boxed{ r_{\min}=\frac{100}{45}\ln50 \approx8.69\ \mathrm{gal/min} }.$$

任何 $r\ge r_{\min}$ 都能在不晚于 $45$ minutes 时达到允许的误差范围。

例：cooling

热物体与环境的温差越大，降温越快。随着温度接近室温，温差缩小，降温也会变慢。设 $T(t)$ 为物体温度，$T_{\mathrm{room}}$ 为保持不变的环境温度。Newton's law of cooling 假设温度变化率与这个温差成正比：

$$\boxed{ \frac{dT}{dt}=-k\left(T-T_{\mathrm{room}}\right), \qquad k>0 }.$$

负号保证温度会朝室温变化：

• 若 $T>T_{\mathrm{room}}$，则 $T'<0$，物体降温。
• 若 $T<T_{\mathrm{room}}$，则 $T'>0$，物体升温。
• 若 $T=T_{\mathrm{room}}$，则 $T'=0$，温度处于 equilibrium。

这里可以用 separation，也可以用 linear equation 的方法。采用后者时，先写成 standard form：

$$T'+kT=kT_{\mathrm{room}}.$$

integrating factor 是

$$\mu(t)=e^{kt}.$$

把方程乘以 $e^{kt}$，左边变成乘积的 derivative：

$$\frac{d}{dt}\left(e^{kt}T\right) =kT_{\mathrm{room}}e^{kt}.$$

做 integration 后再除以 $e^{kt}$，得到

$$T(t)=T_{\mathrm{room}}+Ce^{-kt}.$$

若 $T(0)=T_0$，则 $C=T_0-T_{\mathrm{room}}$，因此

$$\boxed{ T(t)=T_{\mathrm{room}} +\left(T_0-T_{\mathrm{room}}\right)e^{-kt} }.$$

作为比较，separation 会给出

$$\int\frac{dT}{T-T_{\mathrm{room}}} =\int-k\,dt,$$

所以

$$\ln|T-T_{\mathrm{room}}|=-kt+C.$$

这里除以了温差，因此没有包含物体一开始就处于室温的情况。要另外保留这个 constant solution。对于其他解，两边取 exponential，再把正负号合并进任意常数：

$$T(t)=T_{\mathrm{room}}+C_1e^{-kt}.$$

按 exponential 衰减的是温差 $T(t)-T_{\mathrm{room}}$。在这个模型中，$T(t)$ 会逐渐接近 $T_{\mathrm{room}}$，但不会越过它。

咖啡降温

一杯刚倒出的咖啡温度为 $90^\circ\mathrm C$，室温为 $20^\circ\mathrm C$。一分钟后，咖啡降至 $85^\circ\mathrm C$。模型为

$$T(t)=20+70e^{-kt}.$$

利用 $T(1)=85$，有

$$85=20+70e^{-k}, \qquad e^{-k}=\frac{65}{70},$$

所以

$$k=-\ln\left(\frac{65}{70}\right) \approx0.0741\ \mathrm{min}^{-1}.$$

要找咖啡何时降到 $60^\circ\mathrm C$，解方程

$$60=20+70e^{-kt}.$$

得到

$$\boxed{ t=-\frac1k\ln\left(\frac{40}{70}\right) \approx7.55\ \mathrm{min} }.$$

下落物体、mixing tank 和 cooling 的解有一个共同形式：

$$\text{current value} =\text{equilibrium} +\text{initial difference}\times e^{-\lambda t}.$$

其中 $\lambda>0$ 描述差值衰减的快慢；三个模型中分别是 $\gamma/m$、$r/100$ 和 $k$。这是这些 constant-coefficient models 的特点，并不是所有 linear solutions 都有这种形式。

Exact and non-exact equations

对于 exact equation，我们要找一个沿解曲线保持不变的函数。虽然还不知道这个函数，但方程给出了它的两个 partial derivatives，可以据此把它找回来。先明确 partial derivative 的含义。

Partial derivatives

像 $\psi(x,y)$ 这样的函数依赖两个变量。记号

$$\psi_x=\frac{\partial\psi}{\partial x}$$

表示对 $x$ 求 derivative，同时把 $y$ 当作常数。同样地，

$$\psi_y=\frac{\partial\psi}{\partial y}$$

表示对 $y$ 求 derivative，同时把 $x$ 当作常数。

例如，设

$$\psi(x,y)=3x^2y+y^3.$$

固定 $y$，得到

$$\psi_x=6xy,$$

固定 $x$，则得到

$$\psi_y=3x^2+3y^2.$$

求 partial derivative 时，我们固定一个输入，只改变另一个。沿着一条曲线变化时，两个输入都可能改变，因此 chain rule 需要同时用到两个 partial derivatives。

Exact equations

把 first-order equation 写成

$$M(x,y)+N(x,y)\frac{dy}{dx}=0.$$

如果存在函数 $\psi(x,y)$，满足

$$M=\psi_x, \qquad N=\psi_y.$$

这个方程就称为 exact。这样做的用处在于：沿着解 $y=y(x)$，chain rule 给出

$$\frac{d}{dx}\psi(x,y(x)) =\psi_x+\psi_y\frac{dy}{dx}.$$

若 $M=\psi_x$ 且 $N=\psi_y$，那么 differential equation 就是在说

$$\frac{d}{dx}\psi(x,y(x))=0.$$

因此，$\psi$ 沿每条解曲线保持不变：

$$\boxed{\psi(x,y)=C}.$$

从一个已知函数出发，可以看清这个形式是怎么来的。设

$$\psi(x,y)=x^2y+3x.$$

它的 partial derivatives 是

$$\psi_x=2xy+3, \qquad \psi_y=x^2.$$

令 $d\psi/dx=0$，就得到

$$(2xy+3)+x^2\frac{dy}{dx}=0.$$

这个方程是对一个已知函数求 derivative 得到的，所以我们已经知道它来自哪个函数。它的解满足

$$x^2y+3x=C.$$

如何检查 exactness

通常我们知道的是 $M$ 和 $N$，不知道 $\psi$。在尝试求它之前，先检查必要条件。若 $M=\psi_x$ 且 $N=\psi_y$，再分别求 derivative：

$$M_y=(\psi_x)_y=\psi_{xy}, \qquad N_x=(\psi_y)_x=\psi_{yx}.$$

当函数足够 smooth 时，mixed partial derivatives 相等，因此得到检验条件

$$\boxed{M_y=N_x}.$$

这不是要求 $M$ 与 $N$ 相等，而是要求它们的 cross-derivatives 相等。

假设 $M$、$N$、$M_y$、$N_x$ 在一个 rectangular region 上连续。那么，方程在这个区域上是 exact，当且仅当整个区域内都满足 $M_y=N_x$。

求解 exact equation

通过检验后，就可以反过来找 $\psi$。对一个 coefficient 做 integration，先找回函数的一部分，再用另一个 coefficient 确定缺少的部分：

1. 将 $M=\psi_x$ 对 $x$ 做 integration。
2. 加上一个待定函数 $h(y)$，而不只是一个常数。因为对 $x$ 求 derivative 时，任何只依赖 $y$ 的函数都会被当作常数。
3. 将结果对 $y$ 求 derivative。
4. 令这个 derivative 等于 $N=\psi_y$，求出 $h'(y)$。
5. 对 $h'(y)$ 做 integration，得到 $h(y)$，最后把完整的函数 $\psi(x,y)$ 设为 $C$。

例如，考虑

$$2x+y^2+2xy\frac{dy}{dx}=0.$$

这里

$$M=2x+y^2, \qquad N=2xy.$$

cross-derivatives 相等：

$$M_y=2y=N_x.$$

将 $M$ 对 $x$ 做 integration：

$$\psi(x,y)=x^2+xy^2+h(y).$$

再对 $y$ 求 derivative：

$$\psi_y=2xy+h'(y).$$

它必须等于 $N=2xy$，因此 $h'(y)=0$。缺少的函数只是一个常数，可以合并进 $C$。得到

$$\boxed{x^2+xy^2=C}.$$

也可以从另一个 coefficient 开始：将 $N$ 对 $y$ 做 integration，加上一个只依赖 $x$ 的待定函数，再用 $M$ 确定它。哪边的 integral 更容易算，就从哪边开始。

例：含 trigonometric 与 exponential 项的方程

考虑这个 exact differential equation：

$$\left(y\cos x+2xe^y\right) +\left(\sin x+x^2e^y-1\right)\frac{dy}{dx}=0.$$

这里

$$M_y=\cos x+2xe^y=N_x,$$

所以方程是 exact。将 $M$ 对 $x$ 做 integration：

$$\psi(x,y)=y\sin x+x^2e^y+h(y).$$

对 $y$ 求 derivative，并与 $N$ 对照：

$$\psi_y=\sin x+x^2e^y+h'(y) =\sin x+x^2e^y-1.$$

于是 $h'(y)=-1$，可取 $h(y)=-y$。implicit solution 是

$$\boxed{y\sin x+x^2e^y-y=C}.$$

Non-exact equations

如果 $M_y\ne N_x$，方程就不是 exact。有时可以给整个方程乘上一个 integrating factor $\mu$，让新的 coefficients $\mu M$ 和 $\mu N$ 满足 exactness test。

对一般的因子 $\mu(x,y)$，exactness 要求

$$\frac{\partial(\mu M)}{\partial y} = \frac{\partial(\mu N)}{\partial x}.$$

如果 multiplier 同时依赖两个变量，我们就要处理两个新的 partial derivatives。因此，先尝试只依赖一个变量的 multiplier。

假设 $\mu=\mu(x)$，我们需要

$$(\mu M)_y=(\mu N)_x.$$

因为 $\mu$ 只依赖 $x$，它对 $y$ 的 derivative 为零。应用 product rule：

$$\mu M_y=\mu' N+\mu N_x.$$

把 $\mu N_x$ 移到左边，并在 $\mu N$ 非零的地方除以它：

$$\frac{\mu'}{\mu} = \frac{M_y-N_x}{N}.$$

左边只依赖选定的那个变量。如果右边化简后仍然依赖两个变量，这种 multiplier 就行不通。

如果

$$\frac{M_y-N_x}{N}$$

只依赖 $x$，把它记为 $F(x)$，就有

$$\boxed{\mu(x)=e^{\int F(x)\,dx}}.$$

若想找只依赖 $y$ 的因子，令 $\mu=\mu(y)$。product rule 给出

$$(\mu M)_y=\mu' M+\mu M_y, \qquad (\mu N)_x=\mu N_x.$$

令两式相等，并在 $\mu M$ 的两个因子都非零的地方除以它：

$$\frac{\mu'}{\mu} = \frac{N_x-M_y}{M}.$$

因此，如果

$$\frac{N_x-M_y}{M}$$

只依赖 $y$，把它记为 $G(y)$，就有

$$\boxed{\mu(y)=e^{\int G(y)\,dy}}.$$

这里减法的顺序反了过来，是因为 product rule 多出来的那一项换到了另一边。如果忘记符号，不必猜；重新写出 exactness condition，再展开 product rule 就能确定。

例：把 non-exact equation 变成 exact

考虑这个方程：

$$(3xy+y^2)+(x^2+xy)\frac{dy}{dx}=0.$$

先确定

$$M=3xy+y^2, \qquad N=x^2+xy.$$

cross-derivatives 是

$$M_y=3x+2y, \qquad N_x=2x+y.$$

它们不相等，所以方程不是 exact。试着找只依赖 $x$ 的 integrating factor：

$$\begin{aligned} \frac{M_y-N_x}{N} &=\frac{(3x+2y)-(2x+y)}{x^2+xy}\\ &=\frac{x+y}{x(x+y)}\\ &=\frac1x. \end{aligned}$$

在被约去的因子和 $x$ 都非零的区域上，结果只依赖 $x$。因此

$$\mu(x)=e^{\int 1/x\,dx}=e^{\ln|x|}.$$

integrating factor 乘以非零常数，不会改变变换后方程的解。所以，在不跨过 $x=0$ 的 interval 上，可以选

$$\boxed{\mu(x)=x}.$$

把原方程的每一项乘以 $x$：

$$(3x^2y+xy^2)+(x^3+x^2y)\frac{dy}{dx}=0.$$

新的 coefficients 满足

$$(3x^2y+xy^2)_y =3x^2+2xy =(x^3+x^2y)_x,$$

因此变换后的方程是 exact。

在 $x+y=0$ 上，刚才约去的因子为零，那个比值没有定义。但直接检查新的 coefficients，可以看到 exactness 在这条线上仍然成立。因此 integrating factor $\mu=x$ 可以延伸到 $x+y=0$ 上。不过，如果要把两个方程当作等价方程，就仍然不能跨过 $x=0$，因为 multiplier 本身在那里为零。

现在将新的 $M$ 对 $x$ 做 integration：

$$\psi(x,y)=x^3y+\frac12x^2y^2+h(y).$$

于是

$$\psi_y=x^3+x^2y+h'(y).$$

与新的 $N=x^3+x^2y$ 对照，可得 $h'(y)=0$。implicit solution 为

$$\boxed{x^3y+\frac12x^2y^2=C}.$$

选择方法

1. 右边能否拆成只含 $t$ 和只含 $y$ 的两个因子的乘积？能的话，尝试 separation。
2. 能否写成 $y'+P(t)y=g(t)$？能的话，使用 linear integrating factor。
3. 能否写成 $M+Ny'=0$？检查是否满足 $M_y=N_x$。
4. 如果不是 exact，检查两个 single-variable integrating-factor ratios，看看是否有一个适用。
5. 保留因除法而漏掉的 equilibrium solutions，并说明最终公式在哪个 interval 上满足原方程。

Solution intervals 与 uniqueness

求出公式还没有完成 initial-value problem。我们还要知道：经过初始点的解是否存在，是否只有一个，以及它能在时间上延伸多远。

Picard's local existence and uniqueness theorem

从下面的 IVP 出发：

$$y'=f(t,y), \qquad y(0)=y_0.$$

取有限的数 $a,b>0$，以初始点为中心画一个 open rectangle：

$$R=(-a,a)\times(y_0-b,y_0+b).$$

$a$ 是沿 $t$ 方向的水平半宽，$b$ 是沿 $y$ 方向的竖直半高。

假设 $f$ 及其 partial derivative

$$f_y=\frac{\partial f}{\partial y}$$

在整个 $R$ 内都连续。Picard's theorem（也称 Picard–Lindelöf theorem）说明，存在一个数 $h$，满足

$$0<h\le a$$

使得这个 IVP 在 $(-h,h)$ 上恰好有一个解。在这个 local interval 上，解的图像保持在 rectangle 内。

What the local theorem guarantees

$h$ 可以小于或等于 $a$，并不是只要方程 nonlinear，就必须严格小于。还要区分：$(-h,h)$ 是 theorem 保证有解的 interval，不一定是解能存在的最大 interval。

两个假设分别起着不同的作用：

• $f$ 连续，已经足以保证至少有一个 local solution。
• $f_y$ 连续，是保证 local uniqueness 的一个方便的 sufficient condition，因为它限制了 $f$ 在 $y$ 方向上变化的剧烈程度。

这些是 sufficient conditions，不是双向判据。某个条件不满足时，只能说这个 theorem 无法给出相应结论；IVP 仍可能只有一个解，也可能有多个解。

例：选择 Picard rectangle

考虑

$$y'=\frac{y^2}{1+t}+\cos t, \qquad y(0)=y_0.$$

这里

$$f(t,y)=\frac{y^2}{1+t}+\cos t, \qquad f_y(t,y)=\frac{2y}{1+t}.$$

只要 $t\ne-1$，这两个函数都连续，对 $y$ 没有额外限制。由于初始时刻为 $0$，可以选

$$R=(-1,1)\times(y_0-b,y_0+b)$$

其中 $b>0$ 可以是任意有限值。水平方向的边界标为 $-a$ 和 $a$，此处 $a=1$；$b$ 表示从 $y_0$ 向上或向下的距离。singular time $t=-1$ 落在边界上，不属于这个 open rectangle。

事实上，两个函数在整个 unbounded strip $(-1,1)\times\mathbb R$ 上都连续，所以 $b$ 没有自然的上限。但 theorem 中的 rectangle 仍应取任意有限的 $b$，而不是写成 $b=\infty$。

Choosing a rectangle that avoids the singularity

现在 Picard's theorem 保证，在某个 $(-h,h)$ 上存在唯一解，其中 $0<h\le1$。它不会告诉我们 $h$ 的具体值，也不会给出解的公式或 maximal solution interval。

为什么需要 uniqueness 的假设

考虑

$$y'=y^{2/3}, \qquad y(0)=0.$$

对实数 $y$，这里的 fractional power 要用 real cube root 来解释：

$$y^{2/3}:=\left(\sqrt[3]{y}\right)^2=|y|^{2/3}.$$

当 $y<0$ 时，这个约定尤其重要。按这个约定，constant function

$$y_1(t)=0$$

是一个解，因为 $y_1'(t)=0=(y_1(t))^{2/3}$，且 $y_1(0)=0$。与此同时，

$$\boxed{y_2(t)=\frac{t^3}{27}}$$

是另一个解。直接把第二个函数代回原方程检查：

$$\begin{aligned} y_2'(t)&=\frac{t^2}{9},\\ \left(\sqrt[3]{y_2(t)}\right)^2 &=\left(\sqrt[3]{\frac{t^3}{27}}\right)^2 =\left(\frac t3\right)^2 =\frac{t^2}{9},\\ y_2(0)&=0. \end{aligned}$$

两条不同的解曲线经过同一个初始点，因此 uniqueness 不成立。

这套计算对负的 $t$ 也成立：$t^3/27$ 的 real cube root 是 $t/3$，虽然为负，但平方后仍然是 $t^2/9$。

用 separation 解这个方程时，要除以 $y^{2/3}$，因此在 $y=0$ 处这一步不合法，也会漏掉零解。这样找到的非零候选解，都必须像刚才那样代回原 ODE 验证。甚至还存在先在零处停留一段时间、然后沿着平移后的 cubic 变化的解，所以一个初值并不能只选出一条曲线。

右边 $f(t,y)=y^{2/3}$ 是连续的，但当 $y\ne0$ 时，

$$f_y(t,y)=\frac{2}{3\sqrt[3]{y}},$$

这个 partial derivative 在 $y=0$ 处没有定义。因此，在任何围绕 $(0,0)$ 的 rectangle 上，$f_y$ 都不连续。不过，单凭这一点还不能证明 nonuniqueness；真正的证明是我们已经验证了两个不同的解。它们并不违背 Picard's theorem，因为 theorem 所需的 uniqueness 假设没有满足。

Local intervals、maximal intervals 与 finite-time blowup

local theorem 给出的 interval，与 maximal solution interval 是两个概念：

• local interval，例如 $(-h,h)$，是保证 existence 和 uniqueness 的一个 neighborhood。
• maximal interval 则是包含初始时刻、且这个特定的解能始终满足 ODE 的最大 connected interval。

maximal interval 不一定关于初始时刻对称，下面的例子可以说明这一点。

考虑

$$y'=y^2, \qquad y(0)=y_0>0.$$

这里 $f(t,y)=y^2$ 和 $f_y(t,y)=2y$ 处处连续，所以在每个点附近都满足 Picard 的 local existence and uniqueness 假设。对这个正的初值，沿解曲线可以合法地分离变量：

$$\frac{dy}{y^2}=dt.$$

做 integration：

$$-\frac1y=t+C.$$

在 $t=0$ 处，用 initial condition 得到 $C=-1/y_0$。代回去并解出 $y$：

$$\boxed{ y(t)=\frac{1}{\frac1{y_0}-t} =\frac{y_0}{1-y_0t}. }$$

分母在 $t=1/y_0$ 时变成零，并且

$$\lim_{t\to(1/y_0)^-}y(t)=+\infty.$$

所以，包含初始时刻的 maximal interval 是

$$\boxed{\left(-\infty,\frac1{y_0}\right)}.$$

在这个 domain 内，以 $0$ 为中心的最大 symmetric interval 只是

$$\left(-\frac1{y_0},\frac1{y_0}\right).$$

这个 symmetric interval 仍然不是 maximal：解可以一直向左延伸。因此，local theorem 中对称的 $(-h,h)$ 不能与完整的 solution interval 混为一谈。

A smooth equation can still blow up in finite time

这种情况叫 finite-time blowup。它不违背 Picard's theorem：theorem 保证的是 local 的唯一解，并没有说 nonlinear solution 在所有时刻都必须保持有限。相比之下，$y'=y$ 的解在每个有限时刻都保持有限，只有当 $t\to\infty$ 时才会无界增长。

Linear equations 有更强的 interval 结论

考虑 linear IVP：

$$y'+p(t)y=g(t), \qquad y(t_0)=y_0,$$

假设 $p$ 和 $g$ 在 open interval $I$ 上连续，且 $t_0\in I$。那么，在包含 $t_0$ 的这个完整 interval $I$ 上，恰好有一个解，而不只是保证在一个可能更小的 neighborhood 上有唯一解。

integrating-factor formula 可以解释原因。定义

$$P(t)=\int_{t_0}^t p(s)\,ds.$$

那么唯一解为

$$\boxed{ y(t)=e^{-P(t)} \left[ y_0+\int_{t_0}^t e^{P(s)}g(s)\,ds \right]. }$$

连续性保证这些 integrals 对每个 $t\in I$ 都有意义，因此解不会在 $I$ 的内部发生 blowup。同一个 interval 并不意味着整个 $\mathbb R$：某个 coefficient 不再连续时，这个保证就到此为止。具体问题的解可能还能延伸得更远，但仅靠这个 theorem 不能断定。

例：不求解，先确定 interval

考虑

$$y'+\frac{1}{t-2}y=\cos\left(\frac1t\right), \qquad y(1)=4.$$

coefficient 和 forcing term 分别是

$$p(t)=\frac{1}{t-2}, \qquad g(t)=\cos\left(\frac1t\right).$$

$p$ 在 $t=2$ 处不连续，$g$ 在 $t=0$ 处没有定义。它们共同的 continuity set 分成三个 intervals：

$$(-\infty,0), \qquad (0,2), \qquad (2,\infty).$$

初始时刻 $t_0=1$ 只落在其中一个 interval 里，因此 theorem 保证唯一解存在于

$$\boxed{I=(0,2)}.$$

Intersect the coefficient intervals, then keep the piece containing the initial time

不必求 explicit solution。这道题只需要找到包含初始时刻、且两个函数都连续的 interval。