常微分方程解的几何解释-积分曲线族

字数统计: 2.2k字 | 阅读时长: 10分

2023-04-23

摘要: 常微分方程解的几何意义，方向场

【对数据分析、人工智能、金融科技、风控服务感兴趣的同学，欢迎关注我哈，阅读更多原创文章】
我的网站：潮汐朝夕的生活实验室
我的公众号：潮汐朝夕
我的知乎：潮汐朝夕
我的github：FennelDumplings
我的leetcode：FennelDumplings

在文章入门级常微分方程小结中，我们给出了常微分方程及其解的定义，并以物理问题为例说明了初值问题和边值问题这两种定解条件，最后罗列了可以通过初等积分法解决的方程类型。

本文我们以一阶微分方程为例，研究常微分方程其几何意义，并用 Python 进行一些绘图实践。

积分曲线、方向场

考虑一阶方程：

$\frac{\mathrm{d}y}{\mathrm{d}x} = f(x, y)$

其中 $f(x, y)$ 为平面区域 G 内的连续函数，假设 $y = \phi(x), x\in I$ 是方程的解（$I$ 是解的存在区间），则 $y = \phi(x)$ 在 $(x, y)$ 平面上的图形是一条光滑的曲线 $\Gamma$，称其为常微分方程的积分曲线。

任取一点 $P_{0}\in\Gamma$，坐标为 $(x_{0}, y_{0})$，则 $y_{0} = \phi(x_{0})$，由于 $y = \phi(x)$ 满足常微分方程，因此由微分的几何意义，积分曲线 $\Gamma$ 在 $P_{0}$ 的切线斜率为：

$\phi'(x_{0}) = f(x_{0}, \phi(x_{0}))$

因此，即使我们不知道积分曲线 $\Gamma: y = \phi(x)$ 是什么，依然可以知道积分曲线 $\Gamma$ 在$P_{0}$ 点的切线方程如下：

$y = y_{0} + f(x_{0}, y_{0})(x - x_{0})$

于是在 $G$ 内每一点 $P(x_{0}, y_{0})$，可以做一个以 $f(P)$ 为斜率的小线段 $l(P)$，表明积分曲线（如果存在的话）在该点的切线方向，称 $l(P)$ 为微分方程在点 $P$ 的线素。区域 $G$ 上全体线素称为微分方程的方向场。

常微分方程的任何积分曲线与该方程的方向场吻合。

求解初值问题 $\frac{\mathrm{d}y}{\mathrm{d}x} = f(x, y), y(x_{0}) = y_{0}$ 就是求经过点 $(x_{0}, y_{0})$ 并与方向场吻合的一条光滑曲线。

等斜线

在构造微分方程 $\frac{\mathrm{d}y}{\mathrm{d}x} = f(x, y)$ 的方向场时，通常由 $f(x, y) = k$ 确定曲线 $L_{k}$ 来辅助，称其为等斜线。

如果方向场的线素取的够密，方向场就会显示出积分曲线的草图，可以近似地描绘初值问题的积分曲线。这在无法求精确解时，通过方向场可以得到近似解。

奇异点

将一阶方程转换为恰当方程的形式：

$P(x, y)\mathrm{d}x + Q(x, y)\mathrm{d}y = 0$

如果在某个点 $(x_{0}, y_{0})$，有 $P(x_{0}, y_{0}) = Q(x_{0}, y_{0}) = 0$ 时，方向场在 $(x_{0}, y_{0})$ 没有意义，称其为微分方程的奇异点。

用 `quiver` 画箭头图

matplotlib 中的 quiver 方法可用于绘制箭头，前 5 个参数如下：

1	quiver(X, Y, U, V, C, **kwargs)

其中 X 和 Y 决定箭头尾部位置，U 和 V 决定箭头方向，C 是箭头对应的值。

对于二元函数 $z = f(x, y)$，当 $f(x, y) = C$ 时表示一个曲线族（等高线）。对曲线族 $f(x, y) = C$ 取全微分：

$\mathrm{d}z = \frac{\partial{f}}{\partial{x}}\mathrm{d}x + \frac{\partial{f}}{\partial{y}}\mathrm{d}y = 0$

也就是 $\nabla{f} = (\frac{\partial{f}}{\partial{x}}, \frac{\partial{f}}{\partial{y}})$ 与 $(\mathrm{d}x, \mathrm{d}y)$ 正交：

$(\frac{\partial{f}}{\partial{x}}, \frac{\partial{f}}{\partial{y}})\cdot(\mathrm{d}x, \mathrm{d}y) = 0$

由于 $(\mathrm{d}x, \mathrm{d}y)$ 是曲线族 $f(x, y) = C$ 的切向量（方向向量），因此 $(\frac{\partial{f}}{\partial{x}}, \frac{\partial{f}}{\partial{y}})$ 是曲线 $f(x, y) = 0$ 的法向量。

也就是说由隐函数确定的平面曲线族 $f(x, y) = C$ （等高线）的法向量相当于二元函数 $z = f(x, y)$ 的梯度向量。

因此 quiver 方法可以绘制以下内容：

平面曲线族 $f(x, y) = C$ 的方向场
平面曲线族 $f(x, y) = C$ 的法向量场（相当于二元函数 $z = f(x, y)$ 的梯度场）

三元函数也有类似的结论，由隐函数确定的空间曲面族 $f(x, y, z) = C$ （等势面）的法向量相当于三元函数 $w = f(x, y, z)$ 的梯度向量。

二元函数的梯度场

我们以二元函数 $f(x, y) = xe^{x^{2} - y^{2}}$ 的梯度场为例，看看这个方法怎么用。

首先定义函数 $f$:

1 2	def f(x, y): return x * np.exp(- x2 - y2)

然后定义 $f$ 的梯度 $(\frac{\partial{f}}{\partial{x}}, \frac{\partial{f}}{\partial{y}})$：

def grad_field(f, x, y, d=1e-6):
    partial_x = (f(x + d, y) - f(x, y)) / d
    partial_y = (f(x, y + d) - f(x, y)) / d
    return partial_x, partial_y

指定画图的区域：

1	X, Y = np.meshgrid(np.linspace(x_start, x_end, 30), np.linspace(y_start, y_end, 30))

计算强度：

1	C = f(X, Y)

计算方向：

1	U, V = grad_field(f, X, Y)

画图：

1
2
3

plt.quiver(X, Y, U, V, C)
plt.colorbar()
plt.show()

完整代码：

import numpy as np
import matplotlib.pyplot as plt

def f(x, y):
    return x * np.exp(- x**2 - y**2)

def grad_field(f, x, y, d=1e-6):
    partial_x = (f(x + d, y) - f(x, y)) / d
    partial_y = (f(x, y + d) - f(x, y)) / d
    return partial_x, partial_y

def draw(x_start, x_end, y_start, y_end):
    X, Y = np.meshgrid(np.linspace(x_start, x_end, 30), np.linspace(y_start, y_end, 30))
    C = f(X, Y)
    partial_x, partial_y = grad_field(f, X, Y)

    plt.quiver(X, Y, partial_x, partial_y, C)
    plt.colorbar()
    plt.show()

draw(-2, 2, -2, 2)

结果如下：

$二元函数 $f(x, y) = xe^{x^{2} - y^{2}}$ 的梯度场$

上面我们画的是二元函数 $f(x, y) = xe^{x^{2} - y^{2}}$ 的梯度场，相当于平面曲线族（等高线） $f(x, y) = C$ 的法向量场。

这里我们对应地画一下平面曲线族 $f(x, y) = C$ 的方向场，对比看一下。

由于 $\frac{\partial{f}}{\partial{x}}\mathrm{d}x + \frac{\partial{f}}{\partial{y}}\mathrm{d}y = 0$，因此：

$\frac{\mathrm{d}y}{\mathrm{d}x} = -\frac{\frac{\partial{f}}{\partial{x}}}{\frac{\partial{f}}{\partial{y}}}$

这正是一阶常微分方程的形式。也就是说，平面曲线族（等高线） $f(x, y) = C$ 的方向场就是以上一阶常微分方程的积分曲线族的方向场。代码如下：

import numpy as np
import matplotlib.pyplot as plt

def f(x, y):
    return x * np.exp(- x**2 - y**2)

def grad_field(f, x, y, d=1e-6):
    partial_x = (f(x + d, y) - f(x, y)) / d
    partial_y = (f(x, y + d) - f(x, y)) / d
    return partial_x, partial_y

def draw(x_start, x_end, y_start, y_end):
    X, Y = np.meshgrid(np.linspace(x_start, x_end, 30), np.linspace(y_start, y_end, 30))
    C = f(X, Y)
    partial_x, partial_y = grad_field(f, X, Y)
    dx = partial_y
    dy = - partial_x

    plt.quiver(X, Y, dx, dy, C)
    plt.colorbar()
    plt.show()

draw(-2, 2, -2, 2)

$平面曲线族 $f(x, y) = xe^{x^{2} - y^{2}} = C$ 的方向场$

我们也可以通过 contour 直接画出平面曲线族，代码如下：

import numpy as np
import matplotlib.pyplot as plt

def f(x, y):
    return x * np.exp(- x**2 - y**2)

def draw(x_start, x_end, y_start, y_end):
    X, Y = np.meshgrid(np.linspace(x_start, x_end, 30), np.linspace(y_start, y_end, 30))
    Z = f(X, Y)

    plt.contour(X, Y, Z, levels=np.linspace(-1, 1, 20))
    plt.colorbar()
    plt.show()

draw(-2, 2, -2, 2)

$平面曲线族 $f(x, y) = xe^{x^{2} - y^{2}} = C$$

微分方程的积分曲线族

考虑微分方程：

$(x^{2} + 1)(y^{2} - 1)\mathrm{d}x + xy\mathrm{d}y = 0$

求解后，该方程有通解 $y^{2} = 1 + C\frac{e^{-x^{2}}}{x^{2}}$ 和特解 $x=0$。

从微分方程 $\frac{\mathrm{d}y}{\mathrm{d}x} = f(x, y) = -\frac{(x^{2}+1)(y^{2}-1)}{xy}$ 出发，通过线素的方式画积分曲线族大致图形。

import math
import numpy as np
import matplotlib.pyplot as plt

def f(x, y):
    return - (y**2 - 1) * (x**2 + 1) / (x * y)

def derivative(x, y):
    return f(x, y)

def draw(x_start, x_end, y_start, y_end):
    X, Y = np.meshgrid(np.linspace(x_start, x_end, 20), np.linspace(y_start, y_end, 20))
    Z = f(X, Y)
    U = 1.0
    V = derivative(X, Y)

    # 归一化箭头长度
    N = np.sqrt(U ** 2 + V ** 2)
    U = U / N
    V = V / N

    plt.quiver(X, Y, U, V, Z)
    plt.colorbar()

    plt.show()

draw(-3, 3, -3, 3)

$微分方程 $\frac{\mathrm{d}y}{\mathrm{d}x} = f(x, y) = -\frac{(x^{2}+1)(y^{2}-1)}{xy}$ 的积分曲线族的方向场$

从二元函数 $f(x, y) = (y^{2} - 1)\frac{x^{2}}{e^{-x^{2}}} = C$ 出发，通过算偏微分，画出积分曲线族（等高线）的方向场：

import math
import numpy as np
import matplotlib.pyplot as plt

def f(x, y):
    return (y**2 - 1) * x**2 / np.power(math.e, -x**2)

def grad_field(f, x, y, d=1e-6):
    partial_x = (f(x + d, y) - f(x, y)) / d
    partial_y = (f(x, y + d) - f(x, y)) / d
    return partial_x, partial_y

def draw(x_start, x_end, y_start, y_end):
    X, Y = np.meshgrid(np.linspace(x_start, x_end, 20), np.linspace(y_start, y_end, 20))
    Z = f(X, Y)

    partial_x, partial_y = grad_field(f, X, Y)
    dx = partial_y
    dy = - partial_x

    # 归一化箭头长度
    U = dx
    V = dy
    N = np.sqrt(U ** 2 + V ** 2)
    U = U / N
    V = V / N

    plt.quiver(X, Y, U, V, Z)
    plt.colorbar()
    plt.show()

draw(-3, 3, -3, 3)

$平面曲线族 $f(x, y) = (y^{2} - 1)\frac{x^{2}}{e^{-x^{2}}} = C$ 的方向场$

从二元函数 $f(x, y) = (y^{2} - 1)\frac{x^{2}}{e^{-x^{2}}} = C$ 出发，画出积分曲线族（等高线）大致图形：

import math
import numpy as np
import matplotlib.pyplot as plt

def f(x, y):
    return (y**2 - 1) * x**2 / np.power(math.e, -x**2)

def draw(x_start, x_end, y_start, y_end):
    X, Y = np.meshgrid(np.linspace(x_start, x_end, 20), np.linspace(y_start, y_end, 20))
    Z = f(X, Y)

    plt.contour(X, Y, Z, levels=np.linspace(-100, 100, 20))
    plt.colorbar()
    plt.show()

draw(-3, 3, -3, 3)

$平面曲线族 $f(x, y) = (y^{2} - 1)\frac{x^{2}}{e^{-x^{2}}} = C$$

潮汐朝夕的生活实验室 \Doge 陪伴一个算法工程师的职业生涯

数方程

常微分方程解的几何解释-积分曲线族

积分曲线、方向场

等斜线

奇异点

用 `quiver` 画箭头图

二元函数的梯度场

微分方程的积分曲线族

数方程

积分曲线、方向场

等斜线

奇异点

用 quiver 画箭头图

二元函数的梯度场

微分方程的积分曲线族

用 `quiver` 画箭头图