如何在 Python 中调试微分方程的四阶 Runge-Kutta

Question

我一直在尝试在 Python 上使用 Runge Kutta 方法求解这个微分方程。我无法使用我编写的代码获得正确的结果。这里可能出了什么问题？等式为：

$$ rac{du_l}{dx}=u, rac{du}{dx}=(V-w^2)u_l $$

我写的代码是：（这里，Rd是二阶导数，x0和x1是边界x坐标值，V是字典，V[l]=some矩阵（V随l变化），ul0是解的初始值，n是决定一阶导数（1，-1）初始值的变量，h是步长，l是我需要计算值的某个范围）。 [这里的缩进不对，但在Python上工作正常]。

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]
    #x1data=pd.DataFrame(xV1[:,0])
    #x2data=pd.DataFrame(xV2[:,0])
    #data=pd.concat([x1data,x2data],axis=1)
    #print(data.head(50))
    #print(len(xV1),len(xV2))
    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]

        #print(x,'x+h/2',x2,x+(h/2),'x+h',x1,x+h)
        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

Answer 1

我认为您错误地实施了 RK 程序：

• 首先，您应该对两个变量一起执行 RK 的每一步，而不是一次只考虑一个变量。因此，在您的代码中，您希望成对计算 (k1ul, k1u) -> (k2ul, k2u)。由于您的系统是耦合的，因此您需要计算 (k2ul, k2u) 对 (k1ul, k1u) 的结果。尝试再次查看 RK 方案，您不能将其单独应用于变量。

• 您的 RK 计算中的各个步骤是错误的。以第二行为例：k2ul = h * (u + k1ul/2)。这不是微分方程所说的。 k1ul 是 ul 的导数，而不是 u 的导数。所以加法 u + k1ul/2 就像写 $u + .5