如何用Python求解非线性微分方程4(y')³ - y' = 1/x²

这个方程属于关于导数y'的隐式一阶非线性ODE，无法直接写成y' = f(x,y)的显式形式，所以不能直接用odeint或solve_ivp的常规调用方式。我们可以分两步解决：先对每个x求解y'的实根，再积分得到y(x)。

步骤1：拆解方程，求解y'的表达式 #

令p = y'，方程转化为三次代数方程：

4p³ - p - 1/x² = 0

对于每个给定的x，我们可以用数值方法求解这个三次方程的实根p(x)。三次方程的实根数量取决于x的取值：

• 当|x| > 0.98左右时，方程仅有1个实根；
• 当|x| < 0.98时，方程有3个实根，对应方程的3个解分支。

步骤2：Python实现方案 #

方法1：先求p(x)再积分 #

用scipy.optimize.root_scalar对每个x求解p的实根，再用solve_ivp积分得到y(x)。示例代码如下（假设初始条件为x=1时y=0，选择其中一个实根分支）：

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

# 定义三次方程：4p³ - p - 1/x² = 0
def cubic_eq(p, x):
    return 4*p**3 - p - 1/(x**2)

# 定义y'=p(x)，对当前x求解三次方程的实根
def dy_dx(x, y):
    # 用牛顿法求解，初始猜测值决定获取哪个实根分支
    sol = root_scalar(cubic_eq, args=(x,), method='newton', x0=0.5)
    return sol.root

# 初始条件：x=1时y=0
x_range = (1, 5)
y_initial = [0]

# 用solve_ivp进行数值积分
sol = solve_ivp(dy_dx, x_range, y_initial, t_eval=np.linspace(1, 5, 100))

# 绘制结果
plt.plot(sol.t, sol.y[0])
plt.xlabel('x')
plt.ylabel('y(x)')
plt.title('Solution of 4(y\')³ - y\' = 1/x²')
plt.show()

注意事项 #

• 实根分支选择：root_scalar的初始猜测值x0决定了获取哪个实根，x较小时需要调整初始值来切换不同的解分支；
• 奇点处理：当x趋近于0时，1/x²趋近于无穷大，方程解会出现奇点，需避开x接近0的区间；
• 初始条件适配：如果给定的是y(x₀)而非y'(x₀)，需要先根据x₀和y(x₀)反推对应的p(x₀)，再启动积分。

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