Recently I was working on a small project which effectively reduced to finding an approximation of the "steady-state" solutions to a simple linear differential equation (DE) whose fundamental solution in the frequency domain was previously approximated. If the differential operator is linear and time invariant, then the general approach is simple: find the fundamental solution, and then convolve it with the input function.

Admittedly, it took me longer to get the expected behavior. The issue turned out to be scaling and time shifting that was being done incorrectly. Hopefully this document can be helpful in case someone else runs into the same mistake I was making.

Quick Refresher

Fundamental solution

Consider the differential equation

$$ L y = f(t) $$

where $L=\sum{a_k D^{k}}$ is a translationaly invariant linear differential operator and $f(t)$ is a nice-enough function. If $f(t) = \delta(t)$, then the solution to the DE, i.e.

$$ L y_{\delta} = \delta(t), $$

is called the fundumental solution.

If $y_\delta$ is known, then solution of the DE subject to any $f(t)$ is the convolution of the fundamental solution and $f(t)$:

$$ \begin{align} L (f(t) * y_{\delta}) &= \left(\sum{a_k \frac{d^k}{dt^k}} \right) \left( \int_{-\infty}^{\infty}{ f(\tau) y_{\delta}(t-\tau) \, d\tau } \right) \\ &= \int_{-\infty}^{\infty}{ f(\tau) \left[ \left(\sum_k{a_k \frac{d^k}{dt^k}} \right) \right] y_{\delta}(t-\tau) \, d\tau } \\ &= (f * L y_{\delta}) \\ &= (f * \delta) \\ &= f(t) \end{align} $$

If the original assumption on the linearity and translational invariance of the differential operator does not hold, then this approach falls short.

Convolution theorem

The convolution $y_{\delta} * f(t)$ can be evaluated as point-wise product of their transforms:

$$ \mathscr{F}\{ y_{\delta} * f(t) \} = \mathscr{F} \{y_{\delta(t)} \} \cdot \mathscr{F} \{f(t) \} $$

CTFT, DTFT, and DFT

If $X_c(\omega)$ is the CTFT of $x(t)$ and $X(\omega)$ is the DTFT of $x(t)$ sampled at $\Omega_s$ or period $T$, then

$$ X(\omega) = \frac{1}{T} \sum_{k=-\infty}^{\infty}{ X_c(i(\frac{\omega}{T} - k \frac{2\pi}{T})) } $$

The DFS of a periodic (discrete) $\tilde{x}[n]$ signal coincides with the DTF of a finite discrete signal $x[n]$. I.e. the DFT may be thought of as the DFS of periodic extention of $x[n]$. With this, an approximation of the CTFT may be $T \cdot DTFT(x)$

Summary

Knowing the fundamental solution, the 'steady state' response of a linear system (with constant coefficients) to any (nice) input can be quickly determined in the frequency domain (and back in the time domain).

In practice, fundamental solutions to certain class of problems may readily be available. Or discrete approximations to the fundamental solution in the frequency domain may be known. For example, some software may provide solutions to specific DEs in the frequency domain. In that case, the solution of to a unit function $f(\omega) = 1$ will be an approximation of the fundamental solution in the DE in the frequency domain.

Approximation using DFT

  1. Sample $y_{\delta}(t)$ and $f(t)$ at frequency $f_s$. $$ \begin{align} y_{\delta}(t) \rightarrow \boxed{\mbox{sample}} \rightarrow &y_{\delta}[n], \hspace{20px} n \in (0, N_1] \\ f(t) \rightarrow \boxed{\mbox{sample}} \rightarrow &f[n], \hspace{20px} n \in (0, N_2], N_2>N_1 \end{align} $$ Depending on its behavior, $y_{\delta}(t)$ needs to be sampled for ±t, and $f(t)$ has to be delayed by half the length of $y_{\delta}[n]$
  2. Evaluate the discrete time (linear) convolution of $f[n]$ and $y_{\delta}[n]$ while taking into account for scaling.
  3. Discard the first and last $N$ points of convolution.
  4. Translate and sampled $t$ values, and discard the first and last $N$ points.

The example below includes a code snippet

Example

A simple example is $L = D^2 - 4$ or

$$ y''(t) - 4y(t) = f(t) $$

For $f_1(t) = 1$, the general solution to the DE is

$$ y_1(t) = e^{2t} \left( c_1 + c_2 e^{-4t} \right) - \frac{1}{4} $$

Continuous case

It can be verified that

$$ y_\delta(t) = -1/4 e^{-2|t|} $$

is the fundamental solution to the previously mentioned DE (Use ($e^{-2|t|} \delta (t) = \delta(t) $)). Or alternatively, working in the frequnecy domain, the fundamental solution ($\mathscr{F}\{ \delta(t)\} = 1 $) is

$$ \mathscr{F}\{y_\delta(t)\} = \hat{y}_\delta(\xi) = - \frac{1}{(2 \pi \xi)^2 + 4} $$

Here, $ \mathscr{F}\{f(t)\} := \int_{-\infty}^{\infty}{f(t) e^{-i 2 \pi \xi t}} \, dt $. If the solution is expressed in terms of the angular frequency $\omega$ i.e. $ y_\delta(\omega) = -1/(\omega^2 + 4)$ then there will be appropriate scaling of the frequency variable.

Now, $$ \begin{align} y(t) &= \mathscr{F}^{-1} \left\{ y_\delta(\xi) \cdot \mathscr{F} \left\{ f(t) \right\} \right\} \\ y(t) &= \mathscr{F}^{-1} \left\{ - \frac{1}{ (2\pi \xi)^2 + 4} \cdot \delta(\xi) \right\} \\ &= -\frac{1}{4} \end{align} $$

Here, the forward transform $\mathscr{F} \{ 1 \} = \delta(\xi)$ and the inverse transform $ \mathscr{F}^{-1} \{ \hat{f} (\xi) \} := \int_{-\infty}^{\infty}{ \hat{f} (\xi) e^{i 2 \pi \xi t} \, d\xi} $ are consistent with the earlier defintion.

Alternatively, of course, the convolution could have been carried out in the time-domain instead.

Using DFT

Now suppose that we have access to the fundamental solution in the frequency domain that is continuously sampled at $f$. This can be obtained for example through another software by inputing a constant forcing function for a given differential operator. Then the solution in the time-domain can be approximated as before using linear convolution. Here's an example: 

def convolve_diy(x, y, dt):
    ''' Approximate (continuous) linear convolution of x and y

    Convolution is carried out in the frequency domain.

    @param {np.ndarray} x vector of sampled fundamental sol in t-domain
    @param {np.ndarray} y vector of sampled input function in t-domain
    @param {float} dt sampling period
    @return {np.ndarray} approximation to the continuous convolution of x and y
    '''
    M = len(x) + len(y)

    # Approximation to CTFT of input and delta_fn
    x_dft = np.fft.fft(x, M) * dt
    y_dft = np.fft.fft(y, M) * dt

    # Linear convolution
    conv_f = x_dft * y_dft
    # Approximation to ICTFT -- i.e. the continuous-variable linear convolution in t
    return 1./dt * np.real(np.fft.ifft(conv_f))

def cos_example():
    # Fundamental solution in t-domain, y_delta
    deltaFn = lambda t : -1/4.0 * np.exp(-2*np.abs(t))
    # Input function in t-domain, f(t)
    inputFn = lambda t : t*np.cos(t)
    # Analytical solution to $y'' - 4y = f(t)$
    analyticalFn = lambda t : 2/25.*np.sin(t) - 1./5*t*np.cos(t)

    # duration of fundamental solution in t-domain
    W = 30.0
    # Number of samples of y_delta over W
    N = 2**10

    # Sampling points of the fundamental solution in t-domain
    t_delta = np.linspace(-W/2.0, W/2.0, N)

    # Sampling period
    T = abs(t_delta[1] - t_delta[0])

    # Sampling points for f(t)
    #
    # If we are interested in solution for t>=0, delay the sampling by
    # W/2
    #
    # Must be sampled at the same rate as the delta function
    t_fn = np.arange(-W/2.0, -W/2.0 + T*N*10, T)

    # Sampled fundamental solution in t-domain
    y_delta = deltaFn(t_delta)
    # Sampled input function in t-domain
    input_fn = inputFn(t_fn)

    # Approximate the linear convolution of y_delta and f(t) in
    # f-domain
    conv_t = convolve_diy(y_delta, input_fn, T)

    # Discard first and last N points to remove windowing effect
    conv_trunc = conv_t[N-1:-N]

    # `N-1` corresponds to the first point that includes the full
    # (discrete) convolution of the fundamental solution for t=0
    # i.e. where N/2 lies. So shift the `x` points by N/2*dx
    xp = np.arange(-W/2.0, -W/2.0 + len(conv_t) * T, T) \
        - 0.5 * T * (N-1)
    xp = xp[N-1 : -N]
    yp = conv_trunc

    # Plot
    fig, ax = plt.subplots()
    ax.grid(True)
    ax.plot(xp, yp, label='DFT')
    ax.plot(xp, analyticalFn(xp), '--r', linewidth=2, label='Analytical')
    ax.legend(loc='best', fancybox=True, framealpha=0.5)
    ax.set_title('$(D^2 - 4) y(t) = t \cdot \cos(t)$')
    plt.show()

if __name__ == '__main__':
    cos_example()

Graph below shows comparison between analytical and discrete solution where $f(t) = t \cos(t)$ is sampled using 10240 points with sampling period, $T$, of 0.0293.

Next graph investigates whether, as one might expect, some convergent behavior exists with increasing sampling frequency.

It may not be obvious, but at the very initial portion of the curve, there is an initial jump in error with increasing sampling frequency. This behavior is most likely attributed to the initial sampling frequency falling below the Nyquist frequency.

Interactive Demo

Here's a little interactive demo to show the solution computed for an input $f(t)$ for the DE given in the example above. The field accepts common functions, e.g.:

For the exponential function, use Math.pow(expression, power), e.g. Math.pow(t, 2) * cos(2*pi*t)

The interactive demo below may be a little sluggish since I implemented the forward and inverse dfts using a quick (not computationally!) and dirty naive DFT algorithm using brute force since the input discrete vectors are not radix 2.