含奇点的一阶微分方程解曲线Python绘图方法

问题已知条件 #

• 待求解微分方程：$\displaystyle \frac{u}{v}\cdot\frac{dv}{du} = \frac{1-u}{u+v-3}$
• 边界条件：$u=3$时$v=0$，点$(3,0)$处导数值$\frac{dv}{du}=-\frac{5}{3}$
• 理论解特征：解曲线为环绕奇点$(1,2)$的螺旋线

实现逻辑 #

原方程在$(3,0)$处为0/0型奇点，直接从该点启动数值积分会触发除零错误，因此利用已知导数值取距离该点极小偏移的解上点作为积分初值，沿u轴正负两个方向积分，即可得到完整的螺旋线解。

可运行Python代码 #

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

# 整理为显式微分方程 dv/du = f(u,v)
def dv_du(u, v):
    return (v * (1 - u)) / (u * (u + v - 3))

# 取初值近点，避开(3,0)处0/0奇点
eps = 1e-6
u_start = 3 - eps
v_start = (-5/3) * (u_start - 3)  # 由已知点导数计算近点v值
# 校验近点导数值
print(f"初值近点导数：{dv_du(u_start, v_start):.4f}，预期值：{-5/3:.4f}")

# 双向积分覆盖完整曲线
sol_neg = solve_ivp(dv_du, [u_start, -1.5], [v_start], max_step=0.005, rtol=1e-9, atol=1e-12)
sol_pos = solve_ivp(dv_du, [u_start, 5.5], [v_start], max_step=0.005, rtol=1e-9, atol=1e-12)

# 拼接积分结果
u_total = np.concatenate([np.flip(sol_neg.t), sol_pos.t[1:]])
v_total = np.concatenate([np.flip(sol_neg.y[0]), sol_pos.y[0][1:]])

# 绘制曲线
plt.figure(figsize=(8,8))
plt.plot(u_total, v_total, color='#1f77b4', linewidth=1.2, label='u-v解曲线')
plt.scatter(1, 2, c='red', s=70, zorder=10, label='奇点(1,2)')
plt.scatter(3, 0, c='green', s=70, zorder=10, label='边界点(3,0)')
plt.xlabel('u', fontsize=12)
plt.ylabel('v', fontsize=12)
plt.axis('equal')
plt.grid(alpha=0.3)
plt.legend()
plt.title('一阶微分方程环绕奇点螺旋解曲线')
plt.show()

运行代码后会输出初值近点的导数校验值，与给定的-5/3一致，生成的图像中蓝色曲线即为环绕红色奇点的螺旋解。

内容的提问来源于stack exchange，提问作者Dibya Bharati