A Two-Point Cubic Spline

In the specific case of interpolating between GADGET-2 galaxy simulation snapshots, it is possible to use a third-order interpolation scheme to interpolate between two snapshots. A third-order interpolation is possible in this case because the velocities of the particles are known at the end-points. In the simpler case of quadratic spline interpolation, if the acceleration of a particle changes sign even once during the interval, the interpolation is poor and the derivative of the spline at the end points is generally discontinuous. Using a third-order interpolation scheme (cubic spline) allows us to make smoother animations, since the velocity can be kept continuous between stapshot intervals. The only real disadvantage of the cubic spline is that it adds complexity, compared with the quadratic spline.

Cubic spline interpolation is quite popular. Derivations of the general method can be found in most textbooks on numerical methods. I’ll just discuss the specific result that’s useful for the application here—the cubic spline in terms of the end-point positions and velocities. To derive the result, we just have to slightly modify the derivation given in the previous post. We add a jerk term $\mathbf{j}(t)$ (third derivative of the position with respect to time) to the third and fourth Taylor expansions,

$$\mathbf{r}_0 \approx \mathbf{r}(t) + (t_0 - t) \mathbf{v}(t) + \frac{1}{2}(t_0 - t)^2 \mathbf{a}(t)$$

$$\mathbf{r}_1 \approx \mathbf{r}(t) + (t_1 - t) \mathbf{v}(t) + \frac{1}{2}(t_1 - t)^2 \mathbf{a}(t)$$

$$\mathbf{v}_0 \approx \mathbf{v}(t) + (t_0 - t) \mathbf{a}(t) + \frac{1}{2}(t_1 - t)^2 \mathbf{j}(t)$$

$$\mathbf{v}_1 \approx \mathbf{v}(t) + (t_1 - t) \mathbf{a}(t) + \frac{1}{2}(t_1 - t)^2 \mathbf{j}(t)$$

and solve the system of equations for $\mathbf{r}(t)$. There are, of course, several ways of doing this. I’ll skip the details and show the end result. I defined

$$\delta t \equiv t_1-t_0,\qquad \tau_0 \equiv t_0-t,\qquad \tau_1 \equiv t_1-t$$

In terms of these definitions,

$$\mathbf{r}(t)=\frac{\mathbf{r}_0\tau_1^2+\mathbf{r}_1\tau_0 ^2+\tau_0\tau_1(\mathbf{v}_1\tau_0 ^2-\mathbf{v}_0\tau_1^2)/\delta t}{\tau_1^2+\tau_0^2}$$

This can be written more compactly as

$$\mathbf{r}(t)=\frac{\mathbf{r}_0+\rho \mathbf{r}_1+(\rho \mathbf{v}_1-\mathbf{v}_0)\tau_1\tau_0/\delta t}{\rho+1}$$

where

$$\rho\equiv\frac{\tau_0^2}{\tau_1^2}=\frac{(t_0-t)^2}{(t_1-t)^2}$$

As an example, the cubic function generated by this technique when given values corresponding to $\sin(x)$ at $x=0$ and $x=2\pi$ is

$$ s(x)=\frac{(2\pi-x)(\pi-x)x}{x^2-2\pi x+2\pi^2}$$

In the plot below, $\sin(x)$ is the blue curve and $s(x)$ is the dotted black curve.