如何用Python求解含f''(t)=r(f(t))的非线性微分方程组？

嗨，我来帮你搞定这个非线性微分方程组的求解问题！其实核心思路很简单——咱们把二阶方程拆成一阶方程组，就能用你熟悉的scipy工具来解整个系统了。

不管是单个二阶方程还是方程组，scipy.integrate.odeint或者更现代的solve_ivp都只能处理一阶微分方程组。所以咱们要给每个二阶方程引入新变量，把它拆成两个一阶方程，再把所有方程整合起来。

具体变量替换步骤 #

针对你的方程组：

1. 对于方程 f''(t) = r(f(t))：
   定义两个新变量：

   • y1 = f(t) （原函数）
   • y2 = f'(t) （原函数的一阶导数）
     转化为两个一阶方程：

   y1' = y2
   y2' = r(y1)
   

2. 对于方程 g''(t) = s(f(t), g(t))：
   同样定义两个新变量：

   • y3 = g(t) （原函数）
   • y4 = g'(t) （原函数的一阶导数）
     转化为：

   y3' = y4
   y4' = s(y1, y3)  # 因为y1就是f(t)，直接代入即可
   

现在整个系统的状态向量就是 y = [y1, y2, y3, y4]，对应的一阶方程组就完整了，接下来就能用scipy的积分工具求解。

示例代码（用odeint） #

咱们拿一个具体的例子来演示：假设r(f) = -f（类似单摆的恢复力），s(f,g) = f*g，初始条件设为：

• f(0) = 1，f'(0) = 0
• g(0) = 0，g'(0) = 1

代码如下：

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

# 定义整个系统的导数函数
def system(y, t):
    y1, y2, y3, y4 = y
    # 对应f的一阶方程
    dy1_dt = y2
    dy2_dt = -y1  # 这里r(f) = -f
    # 对应g的一阶方程
    dy3_dt = y4
    dy4_dt = y1 * y3  # 这里s(f,g) = f*g
    return [dy1_dt, dy2_dt, dy3_dt, dy4_dt]

# 设置初始条件：[f(0), f'(0), g(0), g'(0)]
initial_conditions = [1, 0, 0, 1]

# 设置求解的时间点
t = np.linspace(0, 20, 1000)

# 求解方程组
solution = odeint(system, initial_conditions, t)

# 提取结果：f(t), f'(t), g(t), g'(t)
f = solution[:, 0]
f_prime = solution[:, 1]
g = solution[:, 2]
g_prime = solution[:, 3]

# 绘图展示
plt.figure(figsize=(12, 6))
plt.subplot(211)
plt.plot(t, f, label='f(t)')
plt.plot(t, g, label='g(t)')
plt.xlabel('t')
plt.ylabel('函数值')
plt.legend()
plt.title('f(t)和g(t)的变化曲线')

plt.subplot(212)
plt.plot(t, f_prime, label="f'(t)")
plt.plot(t, g_prime, label="g'(t)")
plt.xlabel('t')
plt.ylabel('导数值')
plt.legend()
plt.tight_layout()
plt.show()

用solve_ivp的版本（推荐） #

如果你想用更现代、接口更灵活的solve_ivp（支持更多求解器，适合刚性方程组），代码可以这样调整：

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

def system(t, y):
    y1, y2, y3, y4 = y
    dy1_dt = y2
    dy2_dt = -y1
    dy3_dt = y4
    dy4_dt = y1 * y3
    return [dy1_dt, dy2_dt, dy3_dt, dy4_dt]

initial_conditions = [1, 0, 0, 1]
t_span = (0, 20)  # 求解的时间区间
t_eval = np.linspace(0, 20, 1000)  # 需要输出结果的时间点

solution = solve_ivp(system, t_span, initial_conditions, t_eval=t_eval)

# 提取结果
f = solution.y[0]
g = solution.y[2]

# 绘图
plt.plot(solution.t, f, label='f(t)')
plt.plot(solution.t, g, label='g(t)')
plt.xlabel('t')
plt.ylabel('函数值')
plt.legend()
plt.show()

关键注意事项 #

• 初始条件要完整：必须提供每个一阶变量的初始值，也就是f(0), f'(0), g(0), g'(0)，缺一不可。
• 自定义r和s函数：如果你的s(f,g)是cos(f)*sin(g)，直接把dy4_dt改成np.cos(y1)*np.sin(y3)就行，记得用numpy的函数来处理向量运算。
• 刚性方程组处理：如果你的方程组变化速度差异极大（刚性强），可以在solve_ivp里指定刚性求解器，比如method='Radau'或者'BDF'。

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