How to perform a $4^{th}$ Order Runge-Kutta Algorithm without knowing the explicit function.

Question

I was trying to simulate the dynamics of a coin with numerical integration. Using Lagrangian mechanics, I found the EOMs and now I can numerically integrate the functions to find a solution.

I've tried to use Euler's Method, but to I've been having numerical precision errors. I've tried to upgrade the integration method for a Runge-Kutta algorithm, but it seems the definition requires the derivatives to be evaluated at a later instant of time. For example, calculating $k_2$ (following wikipedia's notation): $$ k_2 = f \left(t_n + \frac{h}{2}, y_n + \frac{k_1}{2} \right) $$ where: $$ f(t, y) = \frac{dy}{dt} $$ To calculate the integral, I should know the derivative at a later instant in time. The problem is, the only way I can find the derivative at a later instant in time is when I first calculate the integral, so it is basically it is a recursive cycle.

As you can see from this line of code:

  d2theta = ((m * R * R) * dpsi * dphi * sin(theta) - (m * R * R + Ipsi - Itheta) * dphi * dphi * cos(theta) * sin(theta) - (m*g*R) * cos(theta)) / (m * R * R + Itheta);
  dtheta += d2theta * dt;
  theta += dtheta * dt;

I know the variables $\ddot \theta$ depends on, but I do not have an explicit form of the function. After this, I calculate $\dot \theta$ and $\theta$.

It is possible to implement a Runge-Kutta (or a similar method) in this case? If yes, how would I do it?

EDIT: for further calification here are the differential equations I'm trying to solve: $$ \ddot{\theta} = - \frac{m R^2 + I_{\psi} - I_{\theta}}{mR^2 + I_\theta} \dot{\varphi}^2 \cos \theta \sin \theta + \frac{m R^2}{mR^2 + I_\theta} \dot{\psi} \dot{\varphi} \sin \theta - \frac{m g R}{mR^2 + I_\theta} $$ $$ \dot{\varphi} = \frac{L_\varphi + m R^2 \dot{\psi} \cos \theta }{(m R^2 + I_\psi - I_\theta) \cos^2 \theta + I_\theta} $$ $$ \dot{\psi} = \frac{L_\psi + m R^2 \dot{\varphi} \cos \theta }{m R^2 \cos^2 \theta + I_\varphi} $$ where $\theta, \varphi$ and $\psi$ are variables. It turns out that $\theta$ is the most problematic variable due to being described by a second order differential equation.

Answer 1

You have to take two large steps on the way forward. The first is to conceptually separate the evaluation of the EOM from the specific numerical method. It is quite tempting to meddle around with Euler-like steps or other ideas from finite differences approaches on a component basis, but what you regularly will get is an order 1 method, and rarely an order 2 method. To get to higher orders, one needs to follow the unintuitive developments of the theory in code and data structures.

One part of that is to treat ODE systems as systems $\dot y(t)=f(t,y(t))$, to not look at the single trees but the forest. The state vector $y(t)$ is a vector, the interpretation of its components is left for specialized evaluation/display procedures. The numerical solver only needs vector operations, like copying and forming linear combinations. Some extended remarks and links on that in How to implement a Runge Kutta method (RK4) for a second order differential equation?

The second large step is to recognize implicit equations and the need to apply a dedicated solver to the resulting implicit system. Here it is only a $2\times 2$ system for $\dot\phi,\dot\psi$ which can be easily solved manually. For larger instances of mechanical systems the Euler-Lagrange equations will result in larger linear systems that should then be solved with a linear solver procedure. This could then look like Simplistic simulation of double pendulum doesn't work. Is this approach valid?

One could use symbolic data types to manipulate the Euler-Lagrange equations, see https://stackoverflow.com/a/61078311/3088138, the resulting condensed (stripped of all context) EOM procedure can then look like in https://stackoverflow.com/a/70236831/3088138

Answer 2

You are misunderstanding how the algorithm works. To solve for a discrete set of values of $y$, starting from $y_0 = y(t_0)$, you use the equation $y' = f(t,y)$, and you are going to approximate the successive values of $y$ by quadrature. The simplest version of this is Euler's method, where the integral equation

$$y(t_{i+1}) - y(t_i) = \int_{t_i}^{t_{i+1}} f(s, y(s))ds$$

is replaced by the crude approximation

$$y(t_{i+1}) - y(t_i) \approx f(t_i, y(t_i)) \Delta t_i$$

where $\Delta t_i = t_{i+1} -t_i$. The idea is that if $\Delta t_i$ is small enough, the area is well approximated by the rectangle resulting from a constant value of $f(s,y(s))$. You can use more refined numerical integration to get a more refined estimation.

The point of Runge--Kutta methods is to make the quadrature much more precise, but you never need to know $y'$.