如何用Scipy的solve_ivp数值求解含平方项的微分方程

首先从原方程解出适配scipy.integrate.solve_ivp的显式形式，这是核心步骤：
$$x'(t) = \pm \sqrt{b - a x^2(t)}$$
这对应两个解分支（正号、负号），你之前遇到的定义域错误，本质是根号内的表达式$b - a x^2(t)$在数值计算中出现负数——要么是初始条件不符合约束，要么是浮点误差导致。下面是具体解决方法和代码实现：

核心处理逻辑 #

1. 解的存在性约束

   • 若$a>0$：必须满足$b>0$，且初始值$x(t_0)$符合$a x(t_0)^2 \leq b$，否则无实数解
   • 若$a<0$：$b - a x^2 = b + |a|x^2$，只要$b\geq0$，根号内恒为正，不会出现定义域问题

2. 数值误差规避
   在计算平方根前，强制将根号内的表达式截断为非负数，避免math.sqrt()抛出定义域错误。

代码实现 #

情况1：$a>0$（需处理非负约束） #

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

# 正号分支的右端函数
def f_pos(t, y, a, b):
    sqrt_arg = max(b - a * y**2, 0.0)  # 截断负数，规避定义域错误
    return math.sqrt(sqrt_arg)

# 负号分支的右端函数
def f_neg(t, y, a, b):
    sqrt_arg = max(b - a * y**2, 0.0)
    return -math.sqrt(sqrt_arg)

# 示例参数与初始条件
a = 1.0
b = 4.0
t_span = [0, 10]
x0 = [1.0]  # 满足a*x0²=1 ≤4，符合约束

# 求解两个分支的解
sol_pos = solve_ivp(f_pos, t_span, x0, args=(a, b), dense_output=True)
sol_neg = solve_ivp(f_neg, t_span, x0, args=(a, b), dense_output=True)

# 可视化结果
t_eval = np.linspace(t_span[0], t_span[1], 100)
plt.plot(t_eval, sol_pos.sol(t_eval)[0], label='正号分支')
plt.plot(t_eval, sol_neg.sol(t_eval)[0], label='负号分支')
plt.xlabel('t')
plt.ylabel('x(t)')
plt.legend()
plt.show()

情况2：$a<0$（无定义域问题） #

def f_a_neg(t, y, a, b):
    sqrt_arg = b - a * y**2  # a为负，-a为正，根号内恒正
    return math.sqrt(sqrt_arg)

# 示例参数
a = -2.0
b = 3.0
t_span = [0, 5]
x0 = [0.0]

sol = solve_ivp(f_a_neg, t_span, x0, args=(a, b), dense_output=True)

t_eval = np.linspace(t_span[0], t_span[1], 100)
plt.plot(t_eval, sol.sol(t_eval)[0], label='a<0时的解')
plt.xlabel('t')
plt.ylabel('x(t)')
plt.legend()
plt.show()

关键说明 #

• 当$a>0$时，解是周期性的（类似简谐运动），正负分支对应不同的运动方向
• 使用max(b - a*y**2, 0.0)是为了处理浮点计算的微小误差，避免因精度问题导致根号内出现极小负数
• 若需要完整的周期解，可在解达到极值点（x’=0）时自动切换分支，但这需要额外的逻辑判断；仅需单方向解的话，上述代码即可满足需求

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