如何调试Python中求解微分方程的四阶龙格-库塔法代码

API 探路者

2026-6-30

四阶龙格-库塔法求解微分方程组的代码问题排查

我尝试用Python的四阶龙格-库塔法求解如下微分方程组：
$$ \frac{du_l}{dx} = u,
\frac{du}{dx} = (V-w^2)u_l $$
但代码未能得到正确结果，希望排查问题所在。

代码说明

Rd：二阶导数计算函数，负责计算 $(V-w^2)u_l$
x0、x1：x坐标的边界值
V：字典类型，V[l]为随l变化的矩阵
ul0：解的初始值
n：决定一阶导数初始值的变量，取值为1或-1
h：计算步长
l：计算范围

代码实现

import math
import numpy as np

# 假设w是全局变量，xV函数已实现
def Rd(ul,V):
    return (V-w**2)*ul

def sol(Rd,x0,x1,V,ul0,n,h,l):   
    ulR_final={}
    n1=int((x1-x0)/h)  
    xVs1=xV(x0,x1,V,h,l)  
    xVs2=xV(x0,x1,V,h/2,l)

    for j in range(1,l,1):
        xV1=xVs1[j]
        xV2=xVs2[j]
        x_values=xV1[:,0] 
        ul_values=[ul0]
        u0=math.sqrt(abs(w**2-xV1[0,1]))*ul0*1j*n
        u_values=[u0]
        ulR_values=np.zeros((n1,4), dtype=complex) 
        for i in range(n1-1):
            x=xV1[i,0]
            V=xV1[i,1]
            x1=xV1[i+1,0]
            V_1=xV1[i+1,1] 
            x2=xV2[(2*i)+1,0]
            V_2=xV2[(2*i)+1,1]

            ul=ul_values[-1]
            u=u_values[-1]

            k1ul=h*u
            k2ul=h*(u+(k1ul/2))
            k3ul=h*(u+(k2ul/2))
            k4ul=h*(u+(k3ul/2))

            k1u=h*Rd(ul,V)
            k2u=h*Rd(ul+(k1u/2),V_2)
            k3u=h*Rd(ul+(k2u/2),V_2)
            k4u=h*Rd(ul+(k3u),V_1)

            ul_new=ul+((k1ul+(2*k2ul)+(2*k3ul)+k4ul)/6)
            u_new=u+((k1u+(2*k2u)+(2*k3u)+k4u)/6)

            ul_values.append(ul_new)
            u_values.append(u_new)

            ulR_values[i,0]=x
            ulR_values[i,1]=np.round(ul_new,decimals=5) 
            ulR_values[i,2]=np.round(u_new,decimals=5)  
            ulR_values[i,3]=j
            

        ulR_final[j]=ulR_values[:-1]
        
    return ulR_final

问题排查与修正建议

1. 四阶龙格-库塔法核心公式应用错误

原方程组是一阶耦合方程组，状态向量为 $[u_l, u]$，导数函数为 $f([u_l, u], x) = [u, (V-w^2)u_l]$。代码混淆了两个分量的增量计算，导致每一步近似值错误：

错误点1：k2ul、k3ul、k4ul的计算
当前用k1ul/2、k2ul/2更新u的近似值，但k1ul是$u_l$的增量，和u的增量无关。正确逻辑是：
- k2ul = h * (u + k1u/2)（基于k1u的中点u值）
- k3ul = h * (u + k2u/2)（基于k2u的中点u值）
- k4ul = h * (u + k3u)（基于k3u的终点u值）
错误点2：k2u、k3u、k4u的计算
当前用k1u/2、k2u/2更新u_l的近似值，但k1u是u的增量，和u_l的增量无关。正确逻辑是：
- k2u = h * Rd(ul + k1ul/2, V_2)（基于k1ul的中点u_l值）
- k3u = h * Rd(ul + k2ul/2, V_2)（基于k2ul的中点u_l值）
- k4u = h * Rd(ul + k3ul, V_1)（基于k3ul的终点u_l值）

修正后的RK4计算代码：

# 修正后的RK4步骤
k1ul = h * u
k1u = h * Rd(ul, V)

k2ul = h * (u + k1u / 2)
k2u = h * Rd(ul + k1ul / 2, V_2)

k3ul = h * (u + k2u / 2)
k3u = h * Rd(ul + k2ul / 2, V_2)

k4ul = h * (u + k3u)
k4u = h * Rd(ul + k3ul, V_1)

ul_new = ul + (k1ul + 2*k2ul + 2*k3ul + k4ul) / 6
u_new = u + (k1u + 2*k2u + 2*k3u + k4u) / 6

2. 初始值u0的计算不够严谨

手动处理复数的方式易出错，建议改用numpy的sqrt函数直接处理负数并返回复数结果：

u0 = np.sqrt(w**2 - xV1[0,1]) * ul0 * n

3. 结果存储的x值与解不匹配

当前ulR_values[i,0] = x存储的是步长起点，但ul_new是步长终点的解，两者不对应，应改为存储终点坐标：

ulR_values[i,0] = x1  # x1为当前步的终点xV1[i+1,0]

4. 潜在的步长计算误差

n1=int((x1-x0)/h)可能因浮点数精度问题导致最后一步x超过x1，建议改用np.linspace生成x点，或用n1 = round((x1-x0)/h)后验证终点x值。

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