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

## 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?

• $x$ and $y$ are coupled variables. You have to update them both in each sub-step. The way you are doing it here, you are trying to update $x$ alone. This will give an incorrect answer. Jun 18, 2021 at 17:51
• Follow en.wikipedia.org/wiki/…. The independent variable is $t$.
– lhf
Jun 18, 2021 at 18:12
• See for instance math.stackexchange.com/questions/256393/… for a problem with a correctly implemented RK4 for a second order equation. Or scicomp.stackexchange.com/questions/34257/… for a "condensed" implementation for second order systems. Jun 19, 2021 at 10:48

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

Applying RK4 to the first-order system component-wise

Using the second-order structure of the equation or system