Runge Kutta (RK4) to solve coupled harmonic oscillators [duplicate]

Question

I have the same problem as in this question.

But can someone elaborate on the answer? The poster says that:

Solving this system of 4 ODE's with rk4 will solve for all your state variables simultaneously. position and velocity do not need to be handled seperately (and shouldn't be).

How is this actually implemented? There is still no expression that I can see for $v$, so I do not see how we can obtain the RK coefficients, and I cannot see how to handle position and velocity ‘not seperately’.

For a system of 2 coupled oscillators we have:

$m_1\ddot{x}_1=-k_1x_1+k_2(-x_1+x_2)$

$m_2\ddot{x}_2=-k_3x_2+k_2(x_1-x_2)$

where $x_1$ and $x_2$ are the displacements of masses 1 and 2 from their equilibrium positions. If we then substitute $v=\dot{x}$, we get:

$\dot{v}_1m_1=-k_1x_1+k_2(-x_1+x_2)$

and similar for $\dot{v}_2$. The answer to the question I linked suggested we combine the system to give (they have used a slightly different functional form but the underlying mechanics is the same):

$$ \left(\begin{array}{c} \dot{x}_1 \\ \dot{x}_2 \\ \dot{v}_1 \\ \dot{v}_2 \end{array}\right) = \left(\begin{alignat}{1} & v_1 \\ & v_2 \\ -\omega_1^2(x_1-R_1) &+ \omega_2^2(x_2-x_1-w-R_2)\\ &-\omega_2^2(x_2-x_1-w-R_2) \end{alignat}\right) $$

However I'm unsure how to actually proceed from here. RK4 algorithm needs to evaluate the RK coefficients at different timesteps, and my issue is $\dot{v}=f(x_1,x_2)$ instead of $\dot{v}=f(v)$ in normal (uncoupled) RK4. I can't find $v$ at a new timestep as I need $x$ to do that, and I can't find $x$ without $v$. Is it a matter of writing some operator matrix equation and solving that?

Answer 1

The Runge-Kutta approximates the solution of the initial value problem $$ \dot{y}(t) = F(t, y(t)), \quad y(0) = y_0. $$ The method expects $y_0$ and the function $F$ as input. If we set $$ y(t) = \begin{pmatrix} y_1(t) \\ y_2(t) \\ y_3(t) \\ y_4(t) \end{pmatrix} = \begin{pmatrix} x_1(t) \\ x_2(t) \\ v_1(t) \\ v_2(t) \end{pmatrix} $$ and define $F: [0, \infty) \times \mathbb{R}^4 \to \mathbb{R}^4$ by $$ F(t, y) = \begin{pmatrix} y_3 \\ y_4 \\ -\omega_1^2(y_1-R_1) + \omega_2^2(y_2-y_1-w-R_2) \\ -\omega_2^2(y_2-y_1-w-R_2) \end{pmatrix} \,, $$ then $\dot{y} = F(t, y(t))$ is your ODE. To use the Runge-Kutta method you just have to implement the function $F$ and hand it, together with the initial conditions, to the method. (Note that the right hand side of the function $F$ does not explicitly depend on $t$.)

Answer 2

Differential equations as solved by Runge–Kutta methods generally have the form¹:

$$\dot{y} = \frac{\mathrm{d}y(t)}{\mathrm{d}t} = f(y(t)) \quad \text{with} \quad y:ℝ→ℝ^n,~ f:ℝ^n→ℝ^n$$

For solving these equations, it doesn’t matter:

• what the physical interpretation of the components is, e.g., whether they represent velocity or position,
• whether any components form a subsystem, here, form one of the original oscillators,
• which components of the input of $f$ influence which components of the output, i.e., how the components interact with each other.

For a single oscillator ($\dot{x} = v$, $\dot{v} = -kx$), you have $n=2$ and have to substitute $x$ and $v$ with components of $y$ in the above, e.g., $y_1=x$ and $y_2=v$. With this, the function $f$ for your Runge–Kutta method is:

$$f(y) = \pmatrix{y_2\\-ky_1}.$$

For two oscillators, you can do the same, just that you now have $n=4$. If you do substitute $y_1=x_1$, $y_2=x_2$, $y_3=v_1$, and $y_4=v_2$ in the equations you give, you get:

$$ f(y) = \left(\begin{alignat}{1} & y_3 \\ & y_4 \\ -\omega_1^2(y_1-R_1) &+ \omega_2^2(y_2-y_1-w-R_2)\\ &- \omega_2^2(y_2-y_1-w-R_2) \end{alignat}\right). $$

---

^{¹ To keep things simple, I here ignore the possibility that $f$ depends on $t$ explicitly, as it is not relevant to this problem.}