Runge-Kutta 4th order method to solve second-order ODES of the coupled RLC circuit

I am trying to solve the system of the second-order differential equations of the coupled oscillatory circuit using Runge-Kutta 4th order method in python.

The system

And the initial conditions are:

Condition	Value
L_1 =L_2	2,5E-04 H
C_1	2,5Е-11 F
C_2	5,0Е-11 F
C_св	1,0Е-09 F
I_10	0 A
I_20	0 A
q_10	2,5Е-10 C
q_20	0 C
t_0	0 s
t_max	3,0Е-04 s

My python code is below:


import numpy as np

L = 2.5e-04
C = 1.0e-09
C_1 = 2.5e-11
C_2 = 5.0e-11
i_1_0 = 0
i_2_0 = 0
q_1_0 = 2.5e-10
q_2_0 = 0
t_0 = 0
t_max = 3.0e-04
h = (t_max - t_0) / 100

t = np.arange(t_0, t_max, h)

def f1(q_1, q_2):
    return -q_1 / (L * C_1) + (q_2 - q_1) / (L * C)

def f2(q_1, q_2):
    return -q_2 / (L * C_2) - (q_2 - q_1) / (L * C)

def runge_kutta(q_1, i_1, q_2, i_2, t, h):
    q_1_list = [q_1]
    i_1_list = [i_1]
    q_2_list = [q_2]
    i_2_list = [i_2]

    for _ in t:
        m_1_1 = h * i_1
        k_1_1 = h * f1(q_1, q_2)
        m_2_1 = h * i_2
        k_2_1 = h * f2(q_1, q_2)

        m_1_2 = h * (i_1 + 0.5 * k_1_1)
        k_1_2 = h * f1(q_1 + 0.5 * m_1_1, q_2 + 0.5 * m_2_1)
        m_2_2 = h * (i_2 + 0.5 * k_2_1)
        k_2_2 = h * f2(q_1 + 0.5 * m_1_1, q_2 + 0.5 * m_2_1)

        m_1_3 = h * (i_1 + 0.5 * k_1_2)
        k_1_3 = h * f1(q_1 + 0.5 * m_1_2, q_2 + 0.5 * m_2_2)
        m_2_3 = h * (i_2 + 0.5 * k_2_2)
        k_2_3 = h * f2(q_1 + 0.5 * m_1_2, q_2 + 0.5 * m_2_2)

        m_1_4 = h * (i_1 + k_1_3)
        k_1_4 = h * f1(q_1 + m_1_3, q_2 + m_2_3)
        m_2_4 = h * (i_2 + k_2_3)
        k_2_4 = h * f2(q_1 + m_1_3, q_2 + m_2_3)

        q_1 += (1.0 / 6.0) * (m_1_1 + 2 * m_1_2 + 2 * m_1_3 + m_1_4)
        i_1 += (1.0 / 6.0) * (k_1_1 + 2 * k_1_2 + 2 * k_1_3 + k_1_4)
        q_2 += (1.0 / 6.0) * (m_2_1 + 2 * m_2_2 + 2 * m_2_3 + m_2_4)
        i_2 += (1.0 / 6.0) * (k_2_1 + 2 * k_2_2 + 2 * k_2_3 + k_2_4)

        q_1_list.append(q_1)
        i_1_list.append(i_1)
        q_2_list.append(q_2)
        i_2_list.append(i_2)

    return q_1_list, i_1_list, q_2_list, i_2_list


q_1_points, i_1_points, q_2_points, i_2_points = runge_kutta(q_1_0, i_1_0, q_2_0, i_2_0, t, h)

print(q_1_points)
print(i_1_points)
print(q_2_points)
print(i_2_points)

The code is running but the results I get look suspicious. All quantities quickly acquire huge powers and go to NaN.

Solution

With the given constants the Lipschitz constant of the system is in the region of Lip=1e+14. With an adapted norm that gives the i components a weight of 1e+7, this can be reduced to about Lip=1e+7.

For the step size in Runge-Kutta 4 you want Lip*h=0.5 or smaller, so h=5e-8 or smaller.

sqrt(L*C_k) is also about the wave length of the oscillations, and you want at least 6, better 10 or more points per period to get samples close to the maxima and minima. This amounts to about the same bounds for the step size.

Using solve_ivp with max_step=1e-8 gives smooth looking graphs, both quite periodic, q_1 a simple oscillation, q_2 indeed a superposition, but also periodic looking. No amplitude degradation was observed.

The actual frequency is lower, in q_1 only 2.1e+6, so only 630 oscillations in the time window. The interference pattern in q_2 has a frequency of 3e+5. Despite the lower frequency, using max_step=5e-8 for the sampling resolution frequently misses the tops of the waves, giving visual artifacts.

I distributed the factors a little between q and i, so my i is not your i, the q remain unchanged. This brings the scales of the components closer together, which gives better results from the step size controller.

t_max = 1.0e-05
h = 1e-8
u_0 = [ q_1_0, q_2_0, i_1_0*L, i_2_0*L ]

def f(t,u):
    q_1,q_2,i_1,i_2 = u
    f_1 = -q_1 / C_1 + (q_2 - q_1) / C
    f_2 = -q_2 / C_2 - (q_2 - q_1) / C
    return i_1/L,i_2/L,f_1,f_2

res = solve_ivp(f,(t_0,t_max), u_0, atol=1e-7, rtol=1-10, max_step=h )
print(res.message)

fig,ax=plt.subplots(4,1,figsize=(7,5))
for k in range(4): ax[k].plot(res.t,res.y[k]); ax[k].grid()
plt.tight_layout(); plt.show()