物理

[栖岸计划][数学四大工具之一-微分方程] 第四章：微分方程组与相平面分析

4.1 微分方程组的基本概念

一阶正规方程组（标准形式）：$$\begin{cases}\dfrac{dx_1}{dt} = f_1(t, x_1, x_2, \dots, x_n) \\ \dfrac{dx_2}{dt} = f_2(t, x_1, x_2, \dots, x_n) \\ \vdots \\ \dfrac{dx_n}{dt} = f_n(t, x_1, x_2, \dots, x_n)\end{cases}$$

向量形式：令 $ \mathbf{x} = (x_1, x_2, \dots, x_n)^T $，$ \mathbf{f} = (f_1, f_2, \dots, f_n)^T $，则$$\frac{d\mathbf{x}}{dt} = \mathbf{f}(t, \mathbf{x})$$

高阶方程化为一阶方程组：对于 $ n $ 阶方程 $ y^{(n)} = f(t, y, y', \dots, y^{(n-1)}) $，令$$ x_1 = y,\quad x_2 = y',\quad \dots,\quad x_n = y^{(n-1)} $$则化为：$$\begin{cases}x_1' = x_2 \x_2' = x_3 \vdots \x_{n-1}' = x_n \x_n' = f(t, x_1, x_2, \dots, x_n)\end{cases}$$

4.2 线性微分方程组

齐次线性方程组：$$\frac{d\mathbf{x}}{dt} = A(t) \mathbf{x}$$其中 $ A(t) $ 为 $ n \times n $ 函数矩阵。

非齐次线性方程组：$$\frac{d\mathbf{x}}{dt} = A(t) \mathbf{x} + \mathbf{f}(t)$$

（质心论坛你拉完了）

叠加原理：若 $ \mathbf{x}_1, \mathbf{x}_2 $ 是齐次方程的解，则 $ C_1 \mathbf{x}_1 + C_2 \mathbf{x}_2 $ 也是解。

基解矩阵：以 $ n $ 个线性无关解为列的矩阵 $ \Phi(t) $，满足 $ \Phi'(t) = A(t) \Phi(t) $。齐次通解为 $ \mathbf{x} = \Phi(t) \mathbf{C} $。

常数变易法（非齐次方程）：设 $ \mathbf{x} = \Phi(t) \mathbf{C}(t) $，代入得$$ \Phi(t) \mathbf{C}'(t) = \mathbf{f}(t) \quad \Rightarrow \quad \mathbf{C}(t) = \int \Phi^{-1}(t) \mathbf{f}(t) \, dt $$特解为 $ \mathbf{x}^* = \Phi(t) \int \Phi^{-1}(t) \mathbf{f}(t) \, dt $。

4.3 常系数线性微分方程组

方程形式：$$\frac{d\mathbf{x}}{dt} = A \mathbf{x}$$其中 $ A $ 为 $ n \times n $ 常数矩阵。

解法一：特征值法

设解为 $ \mathbf{x} = \mathbf{v} e^{\lambda t} $，代入得特征方程：$$ (A - \lambda I) \mathbf{v} = 0 $$$ \lambda $ 为特征值，$ \mathbf{v} $ 为对应特征向量。

根据特征根写出解：

解法二：矩阵指数法

定义矩阵指数：$$ e^{At} = \sum_{k=0}^{\infty} \frac{A^k t^k}{k!} $$性质：

• $ \frac{d}{dt} e^{At} = A e^{At} = e^{At} A $
• $ e^{A(t+s)} = e^{At} e^{As} $
• $ (e^{At})^{-1} = e^{-At} $

齐次方程通解：$ \mathbf{x}(t) = e^{At} \mathbf{x}(0) $，基解矩阵 $ \Phi(t) = e^{At} $。

计算 $ e^{At} $ 的方法：

1. Jordan 标准形法：若 $ A = P J P^{-1} $，则 $ e^{At} = P e^{Jt} P^{-1} $。
2. Cayley-Hamilton 法：$ e^{At} = \sum_{k=0}^{n-1} \alpha_k(t) A^k $，系数由特征值确定。
3. 拉普拉斯变换法：$ e^{At} = \mathcal{L}^{-1}[(sI - A)^{-1}] $。

非齐次方程 $ \frac{d\mathbf{x}}{dt} = A \mathbf{x} + \mathbf{f}(t) $ 的解：$$ \mathbf{x}(t) = e^{At} \mathbf{x}_0 + \int_0^t e^{A(t-\tau)} \mathbf{f}(\tau) \, d\tau $$

例：求解 $ \begin{cases} x_1' = 2x_1 + x_2 \\ x_2' = x_1 + 2x_2 \end{cases} $

特征方程 $ \begin{vmatrix} 2-\lambda & 1 \\ 1 & 2-\lambda \end{vmatrix} = (\lambda-1)(\lambda-3)=0 $，特征值 $ \lambda_1=1, \lambda_2=3 $。

$ \lambda=1 $：特征向量 $ (1,-1)^T $，解 $ \mathbf{x}^{(1)} = \begin{pmatrix}1 \\ -1\end{pmatrix} e^t $。

$ \lambda=3 $：特征向量 $ (1,1)^T $，解 $ \mathbf{x}^{(2)} = \begin{pmatrix}1 \\ 1\end{pmatrix} e^{3t} $。

通解：$ \mathbf{x} = C_1 \begin{pmatrix}1 \\ -1\end{pmatrix} e^t + C_2 \begin{pmatrix}1 \\ 1\end{pmatrix} e^{3t} $。

4.4 相平面分析

基本概念：对于二维自治系统$$\begin{cases}\dfrac{dx}{dt} = P(x, y) \\ \dfrac{dy}{dt} = Q(x, y)\end{cases}$$相平面是以 $ (x, y) $ 为坐标的平面，相轨线是解曲线在相平面上的投影，相图是轨线的全体。

平衡点：满足 $ P(x, y) = 0 $ 且 $ Q(x, y) = 0 $ 的点 $ (x_0, y_0) $。

线性化：在平衡点附近，令 $ u = x - x_0 $，$ v = y - y_0 $，则$$\begin{pmatrix} u' \\ v' \end{pmatrix} = J \begin{pmatrix} u \\ v \end{pmatrix} + \text{高阶项}$$其中 Jacobi 矩阵：$$ J = \begin{pmatrix} P_x & P_y \\ Q_x & Q_y \end{pmatrix}_{(x_0, y_0)} $$

平衡点分类（对于 $ \begin{pmatrix} u' \\ v' \end{pmatrix} = A \begin{pmatrix} u \\ v \end{pmatrix} $，$ A = \begin{pmatrix} a & b \\ c & d \end{pmatrix} $）：

记 $ p = \operatorname{tr}(A) = a + d $，$ q = \det(A) = ad - bc $，判别式 $ \Delta = p^2 - 4q $。

例：判断 $ \begin{cases} x' = y \\ y' = -x - y \end{cases} $ 的平衡点类型

平衡点 $ (0,0) $，$ A = \begin{pmatrix} 0 & 1 \\ -1 & -1 \end{pmatrix} $，$ p = -1 $，$ q = 1 $，$ \Delta = 1-4 = -3 < 0 $。

故 $ (0,0) $ 为稳定焦点。

4.5 非线性系统的近似分析

极限环：相平面上的孤立闭轨线。

保守系统：$ \begin{cases} x' = y \\ y' = -g(x) \end{cases} $，存在能量函数 $ E(x,y) = \frac{1}{2}y^2 + \int_0^x g(s) ds $ 守恒，轨线为等能曲线。

耗散系统：存在阻尼项，能量递减，趋于平衡点或极限环。

Liénard 系统：$$\frac{d^2x}{dt^2} + f(x)\frac{dx}{dt} + g(x) = 0$$或等价方程组$$\begin{cases}x' = y - F(x) \y' = -g(x)\end{cases}$$其中 $ F(x) = \int_0^x f(s) ds $。

Liénard 定理：若 $ f(x) $ 偶，$ g(x) $ 奇，$ xg(x) > 0 $（$ x \neq 0 $），且 $ F(x) \to \infty $（$ x \to \infty $），且 $ F(x) $ 在 $ x=0 $ 处有唯一零点，则系统存在唯一稳定极限环。

Van der Pol 方程：$$\frac{d^2x}{dt^2} + \varepsilon(x^2-1)\frac{dx}{dt} + x = 0$$当 $ \varepsilon > 0 $ 时，存在稳定极限环。

（我服了我python不认识汉字）

原码自取

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# Lotka-Volterra 相图绘制
def lotka_volterra(t, z, alpha=1, beta=0.1, delta=0.02, gamma=0.5):
    x, y = z
    dxdt = alpha * x - beta * x * y
    dydt = delta * x * y - gamma * y
    return [dxdt, dydt]

t_span = (0, 50)
y0_list = [[20, 5], [30, 4], [15, 8]]

plt.figure(figsize=(8, 6))
for y0 in y0_list:
    sol = solve_ivp(lotka_volterra, t_span, y0, t_eval=np.linspace(0, 50, 3000))
    plt.plot(sol.y[0], sol.y[1])

plt.xlabel('猎物数量 x')
plt.ylabel('捕食者数量 y')
plt.title('Lotka-Volterra 相图')
plt.grid()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# 定义 Van der Pol 系统
def vanderpol(t, y, eps=1.0):
    x, y_vel = y
    dxdt = y_vel
    dydt = eps * (1 - x**2) * y_vel - x
    return [dxdt, dydt]

# 数值求解
t_span = (0, 30)
y0_list = [[0.5, 0], [1, 0], [2, 0]]  # 不同初始条件

plt.figure(figsize=(8, 6))
for y0 in y0_list:
    sol = solve_ivp(vanderpol, t_span, y0, t_eval=np.linspace(0, 30, 2000))
    plt.plot(sol.y[0], sol.y[1], label=f'初值 ({y0[0]}, {y0[1]})')

plt.xlabel('x')
plt.ylabel('y')
plt.title('Van der Pol 振荡器相图 (ε=1)')
plt.legend()
plt.grid()
plt.axis('equal')
plt.show()