# Getting torsion and curvature out of ODE solution skeleton

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?

• If you use tight integration tolerances (rtol=atol=1e-12 for example, if that's not too heavy), maybe you can get more temporal resolution and hence better derivatives from the interpolation. Does that help ? Do you have a precise example for us to test ? – Laurent90 Nov 26 '20 at 21:37
• Is it feasible to compute derivatives of $f$? For example, $v''(t) = \frac{d}{dt} f(t,v(t)) = f_v(t, v(t)) f(t, v(t)) + f_t(t, v(t))$. Also are you restrict to RK methods? High order linear multistep methods naturally have high order interpolants. – Steven Roberts Nov 26 '20 at 22:24
• @StevenRoberts: Good idea, but $f$ is a black box; moreover, it is unavailable when users analyze the data. Could you explain the natural high order interpolants idea? I can restrict indeed to RK45 and Dormand-Prince. – user14717 Nov 27 '20 at 2:33
• Looks like what @StevenRoberts proposed is essentially the same as my answer. One can just calculate and save the partial derivatives of $f(t,x)$ along the solution trajectory, that would be enough to reconstruct the curvature etc. in the post-processing. – Maxim Umansky Nov 27 '20 at 4:18
• @Laurent90: Sure, the Lorenz ODE would be fine. But literature references and explanation of the high-order interpolants would be preferred. I believe my problem was that I was unaware of the existence of "natural" interpolation for Runge-Kutta methods. – user14717 Nov 27 '20 at 17:43

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.

• Am I correct in thinking that the stages of the stepper need to be stored to use those interpolators? – user14717 Jan 20 at 13:14
• Yes that is true for high order interpolations from RK methods. Though if you do it on the fly via the integrator interface that is not true, since it can just locally interpolate without requiring the saving. – Chris Rackauckas Jan 20 at 14:29

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
for ider in range(1,4):

## 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
for ider in range(1,4):
# 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

## 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)
for ider in range(1,4):
# 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

## plot the values
various_solutions = ((spline_result, '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.

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.

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.