Getting torsion and curvature out of ODE solution skeleton

Question

Suppose I have solved an ODE $v'(t) = f(t,x)$ via some adaptive stepper, such as RK4 or Dormand-Prince, generating a list of points $\{(t_i, v_i, v_i' = f(t_i, v_i))\}_{i=0}^{n-1}$.

I wish to use this data structure to compute the curvature and torsion of the ODE solution, however, computation of curvature requires access to second derivatives, and torsion requires access to third derivatives.

Currently, I differentiate an interpolant to compute these quantities; namely, for a given $t$, I identify the subinterval $[t_i, t_{i+1}]$ which contains it and constructs a cubic Hermite spline using $v_i, v_{i+1}, v_i'$ and $v_{i+1}'$. This gets nice interpolated values, and if I differentiate the interpolant, then I also get pretty nice derivatives.

But second derivatives are sketchy, and third derivatives are hopelessly inaccurate.

In there a better way to compute the torsion and curvature of an ODE solution?

Answer 1

DifferentialEquations.jl includes high order interpolations for many of its schemes. They are dependent on the method, so if you check for example in the ODE Solvers page you'll find things like "Vern9 - Verner's "Most Efficient" 9/8 Runge-Kutta method. (lazy 9th order interpolant)". That is a 9th order method with a 9th order interpolation. Sources can be found by using the help, i.e. ?Vern9 in the REPL. The interpolations include analytical solutions to the derivatives as well, so on the solution handling page it describes how you can use sol(t,Val{3}) to get the third derivative at time t. That would be more than sufficient for your use case.

Answer 2

So I've tried to compute successive time derivatives for the position $x_1$ of a mass from a coupled mass-spring system (adapted from a previous question). The system is linear: $d_t x = Ax$, hence the analytical $n$-th time derivative is $\frac{d^n x}{dt^n} = A^n x$, which permits an easy estimation of the error.

In a "try hard" approach which I am not necessarily proud of, I've tested different ways to compute these derivatives numerically:

• using a spline fit which can then be differentiated
• using finite differences (with gradient from the Python package Numpy) on the solver's output time vector
• using finite differences again, but on the solution interpolated on a finer time grid with the solver's continuous extension
• using finite differences again, but on the solution interpolated on a finer time grid using a spline interpolator

A quick side not on continuous extension / dense output for Runge-Kutta methods: Many RK methods are equipped with a dense output capability which uses the internal stage values to compute a high-order interpolant a the discrete solution, with an interpolation error of order lower than or equal to the order of the method. I'll see if I can find a reference.

Here's the code:

import numpy as np
import matplotlib.pyplot as plt
import numpy.linalg
import scipy.integrate

# System of two coupled masses linked by springs
k1=20;k2=300;m1=1;m2=5
A = np.array([[0,1,0,0],
              [-(k1+k2)/m1, 0, k2/m1, 0],
              [0,0,0,1],
              [k2/m2,0,-k2/m2,0]])
def f(t,x):
    return A.dot(x)
  
def jacfun(t,x):
  return A

def deriv_sol(t,x,order):
  """ Compute the n-th time derivative of the solution"""
  temp = np.copy(x)
  for i in range(order):
    temp = A.dot(temp)
  return temp

x0 = np.array([1, 0, 0, 0]) # initial condition
tf = 0.25 # physical time simulated
tol = 1e-8
# compute numerical solution with adaptive time stepping
sol = scipy.integrate.solve_ivp(fun=f, y0=x0, t_span=(0,tf), method='RK45',
                                atol=tol, rtol=tol, jac=jacfun,
                                dense_output=True)
x1,v1,x2,v2 = sol.y


#%% various ways to ompute the high-order time derivatives
t_test = sol.t
sol_interp = sol.sol(t_test) # high-order continuous extension of the solve_ivp method

## 1 - true solution using the fact that the system is linear
true_result = [deriv_sol(t_test, sol_interp, order=i)[0,:] for i in range(4)]

## 2 - using a spline interpolator
import scipy.interpolate as interp
tck = interp.splrep(sol.t, x1) # get the spline fit
spline_result = [interp.splev(t_test, tck, der=i) for i in range(4)]

## 3 - Finite difference on the solver time grid
grad_result = [sol_interp[0,:]]
for ider in range(1,4):
  grad_result.append( np.gradient(grad_result[-1], t_test))
  
## 3.5 - Finite difference of the solution interpolated on a finer grid (using continuous extension)
t_test_fine = np.linspace(0, tf, int(1e4))
sol_interp_fine= sol.sol(t_test_fine) # high-order continuous extension of the solve_ivp method
grad_result_fine = [sol_interp_fine[0,:]]
for ider in range(1,4):
  grad_result_fine.append( np.gradient(grad_result_fine[-1], t_test_fine))
# reinterp on the initial grid so that we can compare with the other solutions
for ider in range(4):
  tck = interp.splrep(t_test_fine, grad_result_fine[ider]) # get the spline fit
  grad_result_fine[ider] = interp.splev(t_test, tck)

## 3.5.5 - Finite difference of the solution interpolated on a finer grid (using splines)
tck = interp.splrep(sol.t, x1) # get the spline fit
sol_interp_fine_spline = interp.splev(t_test_fine, tck)
grad_result_fine_spline = [ sol_interp_fine[0,:] ]
for ider in range(1,4):
  grad_result_fine_spline.append( np.gradient(grad_result_fine_spline[-1], t_test_fine))
# reinterp on the initial grid so that we can compare with the other solutions
for ider in range(4):
  tck = interp.splrep(t_test_fine, grad_result_fine_spline[ider]) # get the spline fit
  grad_result_fine_spline[ider] = interp.splev(t_test, tck)

## plot the values
various_solutions = ((spline_result, 'spline', '-'),
                     (grad_result, 'grad coarse', '-'),
                     (grad_result_fine, 'grad fine (cont ext)', '-'),
                     (grad_result_fine_spline, 'grad fine (spline)', '--'))

selected_derivatives = range(1,4)
fig, ax = plt.subplots(len(selected_derivatives),1,sharex=True, dpi=300)
for iax, ider in enumerate(selected_derivatives):
  cax = ax[iax]
  for val,name, linestyle in various_solutions:
      cax.plot(t_test, val[ider], label=name, linestyle=linestyle)
  cax.plot(t_test, true_result[ider], label='analytical')
  cax.set_ylim(1.1*np.min(true_result[ider]), 1.1*np.max(true_result[ider]))
  if ider==1:
    cax.legend()
  cax.grid()
  cax.set_ylabel(r'$\frac{{ d^{{{}}} x_1 }}{{ dt^{{{}}} }}$'.format(ider, ider),  rotation='horizontal')
fig.suptitle(r'Successive time derivatives of $x_1$')
plt.tight_layout()

## plot the error
fig, ax = plt.subplots(len(selected_derivatives),1,sharex=True, dpi=300)
for iax, ider in enumerate(selected_derivatives):
  cax = ax[iax]
  for val,name,linestyle in various_solutions:
      cax.semilogy(t_test, np.abs((val[ider]-true_result[ider])/true_result[ider]),
                        label=name, linestyle=linestyle)
  if ider==1:
    cax.legend()
  cax.grid()
  cax.set_ylabel(r'$\frac{{ d^{{{}}} x_1 }}{{ dt^{{{}}} }}$'.format(ider, ider),  rotation='horizontal')
  cax.yaxis.get_major_locator().numticks = 4 
fig.suptitle('Relative errors wrt to analytical values')
plt.tight_layout()

And here are the plots it produces when the solution is computed with integration tolerances set to $1e-8$:

And here for a looser integration tolerance ($1e-6$):

This shows that with these techniques, you need to have a good temporal resolution if you want to compute sensible high-order time derivatives.

Answer 3

Let's for concreteness interpret the ODE as advection of a particle in a given velocity field, $d_t \xi(t) =f(t,x)$. Then the solution of the ODE is the particle's position vs. time $\xi(t)$, and you can plug it in the right-hand side to have the identity $d_t \xi(t) = f(t,\xi(t))$. Then the second derivative of the position along the trajectory is

$ \frac{d^2}{dt^2} \xi(t) = f_t + f_x \xi'(t) = f_t + f_x f(t,x), $

evaluated at $x=\xi(t)$

For the third derivative $d^3 \xi(t)$ one can construct an expression in a similar way.

Answer 4

A good reference is Shampine's paper Interpolation for Runge-Kutta Methods. Shampine contends that there are no "natural" interpolations for Runge-Kutta methods, but a few desiderata constrain the choice of interpolant.

First, Shampine recommends use of a local interpolant (i.e., using just a few points near the desired abscissa to do interpolation; not the entire solution skeleton), because the error of Runge-Kutta methods is only controlled locally anyway. In addition, the complexity of evaluating a local interpolant is lower (perhaps only $\mathcal{O}(\log(N))$ rather than $\mathcal{O}(N)$).

Some of his suggestions:

it has occurred to many people to [use] Hermite interpolation to $y_{n+1}, y_{n+1}'; y_{n}, y_{n}'; y_{n-1}, y_{n-1}'$. Using both solution and derivative approximations is natural here because they are available . .

This is an attractive option if the rhs $f(t, x)$ is not available during interpolation and analysis of the solution.

If the ODE stepper has accuracy of order $p$, then it is reasonable to require the interpolant to also have this accuracy. So quartic interpolation is required for 4th order accurate steppers and quintic is required for 5th-order accurate steppers.

Another question to answer is whether the stages of the method are available at the time the interpolation is desired, or if they have been discarded. If the stages are available, equation (4.2) of the reference can be used.

The 5th order Hermite spline can then be differentiated in order to build the curvature and the torsion.

Another reference which provides some improvements over the suggestions of Shampine is D.J. Higham, Highly Continuous Runge-Kutta Interpolants.