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.