How to debug fourth-order Runge-Kutta for a differential equation in Python

Question

I have been trying to solve this differential equation using Runge Kutta method on Python. I haven't been able to get the right result with the code I wrote. What could be wrong here? The equation is:

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

The code that I wrote is: (Here, Rd is the second derivative, x0 and x1 are the boundary x coordinate values, V is a dictionary, V[l]=some matrix (V varies with l), ul0 is the initial value for the solution, n is a variable which decides the initial value of the first derivative(1,-1). h is the step size, l is some range over which I need to calculate values). [The indentation here is not right, but it works fine on 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

I think you implemented the RK procedure wrongly:

• First, you should do each step of RK together on both variables instead of only considering one variable at a time. So in your code, you want to compute in pairs (k1ul, k1u) -> (k2ul, k2u). Because your system is coupled, the result of the pair (k1ul, k1u) you need to compute (k2ul, k2u). Try to take look at the RK scheme again, you can't apply it to the variables separately.

• The invidual steps in your RK computation are wrong. Take your second line for instance: k2ul = h * (u + k1ul/2). This is not what the differential equation says. k1ul is the derivative of ul, not of u. So the addition u + k1ul/2 is like writing $u + .5 \cdot \frac{du_l}{dx}$, which is wrong.