Is there a generalized formula for solving systems of ODEs using 4th order Runge-Kutta?

Question

As far as I know, we numerically solve any system by reducing it to ODEs and somehow we manage to wire the new system into the RK algorithm. I have seen solutions where there are formulas to use the basic algorithm on a $2$- or $3$-dimensional problem, but so far I don't really understand how we do it.

I mean here basic differential equations so no partial differential equations, no delay, etc.

Answer 1

I hope you mean systems of differential equations of the ordinary type, no partial derivatives, no delay, function and derivative all evaluated at the same state vector.

These get indeed transformed into an explicit system of ODE of order 1 $$ \vec x'(t)=\vec F(t,\vec x(t)). $$ These you can now insert into the Runge-Kutta method of your choice. For the step from demonstrative but simple methods to serious applications the educational methods used are the stereotypical order 4 Heun-Kutta RK4 method with fixed step-size on the simple side and the embedded order 4(5) Fehlberg method with step-size control on the serious side.

_{(RKF45 is dated, the same can be said for the popular DoPri45. There are also simple methods like Fehlberg23 that demonstrate the variable-step principle with a simple but non-trivial step computation.)}

While as an intermediate step a component-wise implementation of the method step is often seen, the most clean and reusable implementation never leaves the vector formulation using vector or array data types. In some programming languages (with vector types with overloaded arithmetic operations) the same method step implementation can be used for the scalar and the vector case. The surrounding time loop has to make a distinction if the trajectory data is collected in a list.

Answer 2

I think you are asking for a general algorithm to implement an Rk4 scheme with some programming language. The implementation needs the following definitions:

1. An array/vector holding the state variables (all the DOF of the system).

   $$ \boldsymbol{Y} = \pmatrix{y_1 \\ y_2 \\ \vdots \\ y_n} $$

   Note that for a dynamical system half the state variables are going to represent positions and half are going to represent velocities (or momentum). In the example of $N$ connected spring-masses you have

   $$ \boldsymbol{Y} = \begin{bmatrix} \boldsymbol{X} \\\boldsymbol{\dot X} \end{bmatrix} = \pmatrix{x_1 \\ \vdots \\ x_N \\ \hline \dot{x}_1 \\ \vdots \\ \dot{x}_N} $$ with $n=2 N$ state variables.

2. A model of the ODE in terms of the rate of change of the state variable as a vector function of time and the current state

   $$ \boldsymbol{\dot Y} = \mathbf{f}(\,t,\;\boldsymbol{Y}) $$

   In the example of $N$ connected spring-masses the model is something like

   $$ \frac{\rm d}{{\rm d}t} \begin{bmatrix} \boldsymbol{X} \\\boldsymbol{\dot X} \end{bmatrix} = \begin{bmatrix} 0 \\ \mathbf{M}^{-1} \boldsymbol{F(t)} \end{bmatrix} + \begin{vmatrix}0 & \mathbf{1} \\ -\mathbf{M}^{-1} \mathbf{K} & -\mathbf{M}^{-1} \mathbf{C} \end{vmatrix} \begin{bmatrix} \boldsymbol{X} \\\boldsymbol{\dot X} \end{bmatrix}$$

   where $\mathbf{M}$ is the mass matrix, $\mathbf{K}$ is the stiffness matrix, $\mathbf{C}$ is the damping matrix and $\boldsymbol{F(t)}$ is the forcing vector.

3. The initial conditions for the ODE, given a starting state vector $\boldsymbol{Y}_0$ and time $t_0=0$.

4. To run the simulation until time $t_{\rm end}$ using $r$ steps, each tiome step is $h = t_{\rm end}/r$. The RK4 algorithm would look like this

   t = 0
   Y = Y_0
   ! store step [t,Y]
   do while(t < t_end)
       K_0 = f(t, Y)
       Y_1 = Y + h/2*K_0
       K_1 = f(t+h/2, Y_1)
       Y_2 = Y + h/2*K_1
       K_2 = f(t+h/2, Y_2)
       Y_3 = Y + h*K_2
       K_3 = f(t+h, Y_3)
       t = t + h
       Y = Y + h/6*(K_0 + 2*K_1 + 2*K_2 + K_3)
       ! store step [t,Y]
   end do
   

   I hope you can read the above Fortran code as the math operations are the easiest to understand IMHO. Assumes basic linear algebra is available where scalar-vector product $\lambda Y$ is defined, and vector-vector addition $Y_1 + Y_2$ is also defined.