使用四阶Runge-Kutta法求解微分方程遇OverflowError问题求助

问题分析 #

你遇到的OverflowError由两个核心错误导致：

• F函数的电容项计算错误：原微分方程中的电容项应为y/(L*C)，但你写成了1/L*C（等价于C/L），数值差了近10个数量级，直接导致方程失去稳定性，数值急剧膨胀。
• Runge-Kutta步骤实现错误：混淆了斜率与增量的计算逻辑，步长h的应用位置错误，进一步加速了数值溢出。

错误代码中的关键问题 #

1. F函数错误：
   原代码中：

   return ((Vo/L -(R0/L)*u -(R1/L)*u**3 - y*(1/L*C)))
   

   其中1/L*C应为1/(L*C)，否则电容项的数值和物理意义完全错误。

2. RK步骤错误：
   对于方程组：
   $$\frac{dy}{dt} = u$$
   $$\frac{du}{dt} = F(y, u, t)$$
   四阶RK的正确增量计算应将步长h与斜率相乘得到增量，而你在计算m1等变量时未正确应用h，导致更新逻辑混乱。

修正后的代码 #

import numpy as np
from matplotlib.pyplot import plot, show, legend

# 参数
R0 = 200
R1 = 250
L = 15
h = 0.002
Vo = 1000
C = 4.2 * 10**(-6)
t_end = 0.93

def F(y, u, x):
    # 修正电容项的计算：1/(L*C) 而非 1/L*C
    return (Vo/L - (R0/L)*u - (R1/L)*u**3 - y/(L*C))

t_points = np.arange(0, t_end, h)
y_points = []
u_points = []

# 初始条件
y = 0.0
u = Vo/L  # du/dt初始值为Vo/L，对应t=0时的电感电压突变

for t in t_points:
    y_points.append(y)
    u_points.append(u)
    
    # 四阶RK正确步骤：计算y和u的增量项
    k1_y = h * u
    k1_u = h * F(y, u, t)
    
    k2_y = h * (u + k1_u / 2)
    k2_u = h * F(y + k1_y / 2, u + k1_u / 2, t + h/2)
    
    k3_y = h * (u + k2_u / 2)
    k3_u = h * F(y + k2_y / 2, u + k2_u / 2, t + h/2)
    
    k4_y = h * (u + k3_u)
    k4_u = h * F(y + k3_y, u + k3_u, t + h)
    
    # 更新y和u
    y += (k1_y + 2*k2_y + 2*k3_y + k4_y) / 6
    u += (k1_u + 2*k2_u + 2*k3_u + k4_u) / 6

# 绘图
plot(t_points, u_points, label='u(t)')
plot(t_points, y_points, label='y(t)')
legend()
show()

说明 #

• 修正后的F函数恢复了正确的电容项计算，保证方程的物理稳定性。
• RK步骤严格按照四阶方法的标准流程实现，正确计算每个阶段的增量，避免数值发散。
• 合并了绘图代码，添加图例方便区分两条曲线。

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