1

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.

yigal
  • 23
  • 3
  • You should try a smaller step size. The larger the step size, the larger your error. Since errors propagate forward in this method, they grow exponentially with each step, eventually breaking the math. Smaller step sizes are more accurate, but will take longer to compute. – the_strange Oct 24 '22 at 08:57
  • With your constants you get Lipschitz constants in the region `Lip=1e+14`. With an adapted norm that gives the `i` components a weight of `1e7`, this can be reduced to about `Lip=1e7`. 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. – Lutz Lehmann Oct 24 '22 at 09:12
  • @LutzLehmann I changed my step to h=5e-8 and my results became more natural. However I tried to build [q_1(t)](https://i.imgur.com/oYLYLFM.png) and [q_2(t)](https://imgur.com/JrjpPxm) plots. Are they correct or I have some mistakes in my calculations? – yigal Oct 26 '22 at 09:10
  • By the rough frequency/wavelength estimate above, you have `3e-4*1e7=3000` oscillations in the time interval. This are still multiple oscillations per pixel column, so the filled graphs are to be expected. As is some exponential or long-wave dynamic in the amplitudes as the two oscillations are coupled. The pattern in the second could be super-imposition of multiple frequencies, or more likely an aliasing issue, that is, purely a feature of the plot and not of the solution. – Lutz Lehmann Oct 26 '22 at 09:18
  • @LutzLehmann So as I understood my solution is correct. Is it possible to prettify my plots? Sorry if my questions are too dumb, I'm not experienced in physics and solving ODEs, as well as in using matplotlib. – yigal Oct 26 '22 at 09:32
  • There is indeed superposition and interference in the q_2 plot. Reducing the time interval to only plot a few dozen oscillation gives more informative plots. – Lutz Lehmann Oct 26 '22 at 10:13

1 Answers1

1

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()

enter image description here

Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51