Why doesn't the convolution of the impulse response match the system's output?

Question

If you define an LTI system sys in scipy, you may conveniently feed an input x to it to get its output y as follows:

from scipy import signal
import numpy as np
import matplotlib.pyplot as plt

# let's define as mass-damper-spring system:
# x'' + beta x' + w0**2 x = F/m
m = 10
w0 = 20*2*np.pi
beta = 1.5 

fs = 500 # sampling frequency (Hz)
ts = np.arange(0,100,1./fs) # time interval of simulation
sys = signal.lti(1./m,[1,beta, w0**2]) 

# simulate the response to an input
x = 10+np.random.normal(0,1,size=ts.shape) # input
_, y, _ = sys.output(U = x, T=ts, ) # output

However, it is also possible convolve the impulse response (which is also accessible via lti module) with the input and find the output to the input:

_, kern = sys.impulse(T=ts) # impulse response 
y_conv = signal.convolve(kern/fs,x)[:len(ts)] #1/fs=dt of integration

I realized that the normalized y_conv and normalized y differ even though they should be the same theoretically (look below):

fig, axs= plt.subplots(3,1, figsize=(10,5), sharex=True)
axs[0].plot(ts, y,label='output y') 
axs[1].scatter(ts, y-y_conv,marker='.', label='difference') 
axs[2].plot(ts, kern, label='kernel (impulse response)') 
plt.xlabel('time')
for ax in axs:
    ax.legend()

Indeed, scipy team didn't implement the output of the linear system using convolution.

I did a similar experiment with an even simpler system (sys = signal.lti(1,[1,1])) and realized that this error is indeed a function of the sampling frequency fs, as the following images indicate (note the ranges of middle panels). The higher the fs the lower the difference.

These pictures suggest that the difference vanishes asymptotically, as also commented below. However, my questions that concern finite sampling rate regimes are:

1. How does the sampling frequency cause such differences mathematically? Is there any expression for the error bound?
2. What is the role of bilinear mapping in this difference? Is it possible to tailor this mapping to mitigate such numerical discrepancy?
3. Which method is more accurate to estimate the response to a general input, and why (implementation-wise)?

Answer 1

How does the sampling frequency cause such differences?

The differences are caused by the fact that the discrete-time convolution between two discrete signals is not equal to the discrete signal of continuous-convolution between two continuous signals.

signal.convolve gives you the discrete-time convolution result, which refers to convolution sum, while sys.output returns the continuous-time convolution result, which is also known as convolution integral.

Although we can't implement a convolution integral using computer, the documentation of scipy.signal.lsim indicates that it interpolates the input signal first, and then apply a convolution sum to simulate the output of a continuous-time linear system.

Which method is more accurate to estimate the response to a general input, and why (implementation-wise)?

Since you are simulating a mass-damper-spring system which is a continuous-time system, sys.output is more accurate.

Answer 2

I investigated this further and turned out the main issue is not really the sampling rate, but the integration approximation. In the code above, I approximated the convolution $kern(t) * x(t)$ with the sum(kappa * u) *dt, corresponding to left Riemann sum of the integral. A better way to approximate is the trapezoidal rule which can be implemented by sum(kappa * u) /2 *dt. As the code below shows, this one matches the response of Scipy's impulse much better:

from scipy import signal
import numpy as np
import matplotlib.pyplot as plt

# let's define as mass-damper-spring system:
# x'' + beta x' + w0**2 x = F/m
m = 10
w0 = 10*2*np.pi
beta = 15 

fs = 2000 # sampling frequency (Hz)
T = 2
ts = t = np.linspace(0, T, int(T*fs), endpoint=False) # time interval of simulation
sys = signal.lti(1./m,[1,beta, w0**2]) 
_, kern = sys.impulse(T=ts) # impulse response 

# Generate an input: I substituted the previous input with a noise-free 
# impulse-like input for a better comparison
# x = 10 + np.random.normal(0, 1, size=ts.shape) # old input
x = np.zeros_like(ts)
x[0] = 1.

_, y, _ = sys.output(U=x, T=ts, ) # output
y_conv_rec = signal.convolve(x, kern)[:len(ts)] * 1./fs     # rectangular elements
y_conv_trap= signal.convolve(x, kern)[:len(ts)] * 1./fs / 2 # trapezoidal rule


fig, axs= plt.subplots(2,1, figsize=(8,5), sharex=True)
axs[0].plot(ts, y, label='y (LTI output)', linewidth=2) 
axs[0].plot(ts, y_conv_rec, '--', label='left Riemann sum',)
axs[0].plot(ts, y_conv_trap, ':', label='trapezoidal',)

axs[1].plot(ts, y-y_conv_rec, label='diff left sum')
axs[1].plot(ts, y-y_conv_trap, label='diff trap')

axs[0].set_ylabel('Response')
axs[1].set_ylabel('diff (relateive to LTI output)')
plt.xlabel('time')
for ax in axs:
    ax.legend()

Clearly, the trapezoidal rule isn't the perfect scheme and Scipy uses a more elaborate one. Yet, it's good enough to show that the discrepancy I observed above was indeed due to integration and not, e.g. the kernel estimation.