How to implement a Runge Kutta method (RK4) for a second order differential equation?

Question

problem

I have the following system of differential equations; $$ \frac{dx}{dt} = y \tag{1}\label{1}\\ $$ $$ \frac{dy}{dt} = -\lambda^2 x -A \tag{2} \label{2} $$ where $\lambda$ and $A$ values are known values, we also further have the following initial condition values; $$ x_0=x(0) = 0.5 \\ y_0=y(0) =\left. \frac{dx}{dt} \right\vert_{x=0} = 0.25 $$ 1)How do I solve this system using RK4 method?

My attempt

According to RK4 method we get \begin{aligned} x_{n+1}&=x_{n}+{\frac {1}{6}}h\left(k_{1}+2k_{2}+2k_{3}+k_{4}\right),\\ t_{n+1}&=t_{n}+h \end{aligned} where \begin{aligned}k_{1}&=\ y_{n},\\k_{2}&=\ y_{n}+h{\frac {k_{1}}{2}},\\k_{3}&=\ y_{n}+h{\frac {k_{2}}{2}},\\k_{4}&=\ y_{n}+hk_{3}.\end{aligned}

2)After this, should I take $x_n$ or $x_{n+1}$ value for calculating $\eqref{2}$?

3)How do I use matrices to solve this system?

Answer 1

This question gets often asked. I'd like to speculate that there are 3 stages to understanding numerical ODE methods of the Runge-Kutta variety:

1. low-order methods applied to the scalar case,
2. transition to population models or mechanical second-order equations with 2 or 3 components,
3. the final insight that all methods for scalar first-order equations (except Kutta's 5th order method) apply without restriction to first-order systems, and that all ODE systems can be transformed to such first-order systems.

There appears to be a mental block to the idea to put all components of the system into one flat vector and then treat that vector like a scalar (in this context), also separating the numerical method from the dynamical model. I think that to overcome this block seems to be a short period in the learning curve. I believe that this is the reason that the example code in compendiums like wikipedia, rosetta-code etc. is mostly geared to the scalar case, phase 1 will not yet reliably understand the multi-component case and phase 3 knows that, up to some trivial overhead, the same code structure can be used, and the scalar-case code is a little more compact to write down.

Some links that contain the treatment of second order DE or systems of DE (like mechanical assemblages, celestial systems, ...) using the RK4 method

Applying the vectorized RK4 method to a first-order system constructed from a higher-order equation

• Example: Solve fourth order ODE using fourth order Runge-Kutta method, apply RK4 in vectorized version like in https://stackoverflow.com/questions/61333687/how-do-i-solve-a-2nd-order-differential-equation-for-projectile-motion-with-air

Applying RK4 to the first-order system component-wise

• Implementation example: https://stackoverflow.com/questions/65224923/i-want-to-have-the-pendulum-blob-in-my-double-pendulum, this should demonstrate that this approach is not very suitable if the number of components increases

Using the second-order structure of the equation or system

• Theory: https://scicomp.stackexchange.com/questions/34257/solving-coupled-odes-using-runge-kutta-method
• Implementation: https://math.stackexchange.com/questions/2615672/solve-fourth-order-ode-using-fourth-order-runge-kutta-method, 2-body gravity simulation with partitioned RK4 code