Algorithm to compute a convolution recursively

Question

Let $$ f(t) = \int_0^t k(t-s)g(s) \, ds. $$ Assume that $g$ is only given in a grid $t_j = j\delta_t$, and that we wish to compute similarly $f$ on the same grid.

What's an efficient algorithm to approximate $f(t_{j+1})$ given (an approximation of) $f(t_j)$? I tried to write \begin{align} f(t_{j+1})&=f(t_j)+\int_0^{t_{j+1}}k(t_{j+1}-s)g(s) \, ds - \int_0^{t_j}k(t_j-s)g(s) \, ds \\ &=f(t_j)+\int_0^{t_j} \left(k(t_{j+1}-s) - k(t_{j}-s)\right)g(s) \, ds + \int_{t_j}^{t_{j+1}}k(t_{j+1}-s)g(s) \, ds. \\ \end{align} The second integral can be approximated easily at each step with a trapezoidal rule. For the first integral, we would need to recompute looking in the past until 0 at each step?

Answer 1

Let $\displaystyle\{t_j\}_{j=0}^{N} \subset \mathbb{R}$ be a discrete grid with $\displaystyle t_j = j \delta_t$, where $\displaystyle \delta_t$ is the step size. Let $\displaystyle f : \mathbb{R} \to \mathbb{R}$ be:

$$ \displaystyle f(t) = \int_0^t k(t-s) g(s) \, ds, $$

where $\displaystyle k : \mathbb{R} \to \mathbb{R}$ and $\displaystyle g : \mathbb{R} \to \mathbb{R}$ are given functions. We have

$$ \displaystyle f(t_{j+1}) = \int_0^{t_{j+1}} k(t_{j+1} - s) g(s) \, ds = \int_0^{t_j} k(t_{j+1} - s) g(s) \, ds + \int_{t_j}^{t_{j+1}} k(t_{j+1} - s) g(s) \, ds. $$

Since

$$ \displaystyle f(t_j) = \int_0^{t_j} k(t_j - s) g(s) \, ds, $$

we have

$$ \displaystyle f(t_{j+1}) = f(t_j) + \int_0^{t_j} \left[ k(t_{j+1} - s) - k(t_j - s) \right] g(s) \, ds + \int_{t_j}^{t_{j+1}} k(t_{j+1} - s) g(s) \, ds. $$

Using the trapezoidal rule for 2nd integral:

$$ \displaystyle\int_{t_j}^{t_{j+1}} k(t_{j+1} - s) g(s) \, ds \approx \frac{\delta_t}{2} \left[ k(\delta_t) g(t_j) + k(0) g(t_{j+1}) \right]. $$

For the 1st integral, we assume $k(t_{j+1} - s) - k(t_j - s)$ is small and approximate it as:

$$ \displaystyle\int_0^{t_j} \left[ k(t_{j+1} - s) - k(t_j - s) \right] g(s) \, ds \approx \delta_t \sum_{m=0}^{j-1} \left[ k(t_{j+1} - t_m) - k(t_j - t_m) \right] g(t_m). $$

Hence, $$ \displaystyle f(t_{j+1}) \approx f(t_j) + \delta_t \sum_{m=0}^{j-1} \left[ k(t_{j+1} - t_m) - k(t_j - t_m) \right] g(t_m) + \frac{\delta_t}{2} \left[ k(\delta_t) g(t_j) + k(0) g(t_{j+1}) \right]. $$

Implementation.

If we use this method, for $g(t)=\sin{(t)}$ and $k(t)=\exp{(-t)}$ and noting:

$$\int_{0}^{t}e^{s-t}\sin(s)\:ds=\sin(t)-\int_{0}^{t}e^{s-t}\cos(s)\:ds \implies \int_{0}^{t}e^{s-t}\cos(s)\:ds=-\sin(t)+\frac{1}{2}(e^{-t}+\sin(t)-\cos(t))\implies\int_{0}^{t}e^{s-t}\sin(s)\:ds=\frac{1}{2}(e^{-t}+\sin(t)-\cos(t)). $$

#include <stdio.h>
#include <math.h>

#define DELTA_T 0.01
#define MAX_T 10
#define N ((int)(MAX_T / DELTA_T) + 1)

double k(double t) {
    return exp(-t);
}

double g(double t) {
    return sin(t);
}

double analytical_solution(double t) {
    return 0.5 * (exp(-t) - cos(t) + sin(t));
}

int main() {
    double tGrid[N];
    double fAnalyticalValues[N];
    double fNumerical[N];

    for (int i = 0; i < N; i++) {
        tGrid[i] = i * DELTA_T;
    }

    for (int i = 0; i < N; i++) {
        fAnalyticalValues[i] = analytical_solution(tGrid[i]);
    }

    fNumerical[0] = 0.0;

    for (int j = 1; j < N; j++) {
        double integral1 = 0.0;
        for (int i = 0; i < j - 1; i++) {
            integral1 += (k(tGrid[j] - tGrid[i]) - k(tGrid[j - 1] - tGrid[i])) * g(tGrid[i]) * DELTA_T;
        }
        double integral2 = 0.5 * DELTA_T * (k(DELTA_T) * g(tGrid[j - 1]) + k(0) * g(tGrid[j]));
        fNumerical[j] = fNumerical[j - 1] + integral1 + integral2;
    }

    printf("t\tAnalytical f(t)\tNumerical f(t)\n");
    for (int i = 0; i < N; i++) {
        printf("%12.6lf\t%20.12lf\t%20.12lf\n", tGrid[i], fAnalyticalValues[i], fNumerical[i]);
    }

    return 0;
}

which

$$ \left( \begin{array}{ccc} \text{t} & \text{Analytical f(t)} & \text{Numerical f(t)} \\ \hline 0.000000 & 0.000000000000 & 0.000000000000 \\ 0.010000 & 0.000049833333 & 0.000049999167 \\ 0.020000 & 0.000198666667 & 0.000199494167 \\ 0.030000 & 0.000445500001 & 0.000447484951 \\ 0.040000 & 0.000789333339 & 0.000792971424 \\ 0.050000 & 0.001229166688 & 0.001234953442 \\ 0.060000 & 0.001764000064 & 0.001772430824 \\ 0.070000 & 0.002392833495 & 0.002404403347 \\ \cdots & \cdots & \cdots \\ 9.930000 & 0.195556984477 & 0.204886234848 \\ 9.940000 & 0.188751735036 & 0.198056873899 \\ 9.950000 & 0.181927615399 & 0.191208209748 \\ 9.960000 & 0.175085307925 & 0.184340927209 \\ 9.970000 & 0.168225496791 & 0.177455712955 \\ 9.980000 & 0.161348867926 & 0.170553255456 \\ 9.990000 & 0.154456108941 & 0.163634244905 \\ 10.000000 & 0.147547909058 & 0.156699373152 \\ \end{array} \right). $$

More accurate integral methods includes (wiki).

Use

$$ f(t_i) \approx \sum_{j=0}^{i-1} \int_{t_j}^{t_{j+1}} \exp(-(t_i - s)) \sin(s) \, ds. $$

By Simpson's rule, the integral over each subinterval $[t_j, t_{j+1}]$ can be: $$ \int_{t_j}^{t_{j+1}} \exp(-(t_i - s)) \sin(s) \, ds \approx \frac{\Delta t}{6} \left[ \exp(-(t_i - t_j)) \sin(t_j) + 4 \exp(-(t_i - (t_j + \frac{\Delta t}{2}))) \sin(t_j + \frac{\Delta t}{2}) + \exp(-(t_i - t_{j+1})) \sin(t_{j+1}) \right]. $$

So $$ f(t_i) = \sum_{j=0, j+2 \leq i}^{i-2} \frac{\Delta t}{6} \left[ \exp(-(t_i - t_j)) \sin(t_j) + 4 \exp(-(t_i - (t_j + \frac{\Delta t}{2}))) \sin(t_j + \frac{\Delta t}{2}) + \exp(-(t_i - t_{j+1})) \sin(t_{j+1}) \right]. $$

For the remaining terms when $ j+2 \not\leq i $, we use the trapezoidal rule: $$ \int_{t_j}^{t_{j+1}} \exp(-(t_i - s)) \sin(s) \, ds \approx \frac{\Delta t}{2} \left[ \exp(-(t_i - t_j)) \sin(t_j) + \exp(-(t_i - t_{j+1})) \sin(t_{j+1}) \right]. $$

#include <stdio.h>
#include <math.h>

#define DELTA_T 0.01
#define MAX_T 10
#define N (int)(MAX_T / DELTA_T) + 1

double analytical_solution(double t);
double k(double t);
double g(double t);

int main() {
    double timeGrid[N];
    double fAnalyticalValues[N];
    double fNumericalValues[N];

    for (int i = 0; i < N; i++) {
        timeGrid[i] = i * DELTA_T;
    }

    for (int i = 0; i < N; i++) {
        fAnalyticalValues[i] = analytical_solution(timeGrid[i]);
    }

    fNumericalValues[0] = 0;
    for (int i = 1; i < N; i++) {
        double sum = 0;
        for (int j = 0; j < i; j += 2) {
            if (j + 2 <= i) {
                double t0 = timeGrid[j];
                double t1 = timeGrid[j + 1];
                double t2 = timeGrid[j + 2];
                double integrand0 = k(timeGrid[i] - t0) * g(t0);
                double integrand1 = k(timeGrid[i] - t1) * g(t1);
                double integrand2 = k(timeGrid[i] - t2) * g(t2);
                sum += (integrand0 + 4 * integrand1 + integrand2) * DELTA_T / 3.0;
            } else if (j + 1 <= i) {
                double t0 = timeGrid[j];
                double t1 = timeGrid[j + 1];
                double integrand0 = k(timeGrid[i] - t0) * g(t0);
                double integrand1 = k(timeGrid[i] - t1) * g(t1);
                sum += 0.5 * (integrand0 + integrand1) * DELTA_T;
            }
        }
        fNumericalValues[i] = sum;
    }

    printf("Time\t\tAnalytical\t\t\tNumerical\n");
    for (int i = 0; i < N; i++) {
        printf("%12.6lf\t%20.12lf\t%20.12lf\n", timeGrid[i], fAnalyticalValues[i], fNumericalValues[i]);
    }

    return 0;
}

double analytical_solution(double t) {
    return 0.5 * (exp(-t) - cos(t) + sin(t));
}

double k(double t) {
    return exp(-t);
}

double g(double t) {
    return sin(t);
}

$$ \left( \begin{array}{ccc} \text{t} & \text{Analytical f(t)} & \text{Numerical f(t)} \\ \hline 0.000000 & 0.000000000000 & 0.000000000000 \\ 0.010000 & 0.000049833333 & 0.000049999167 \\ 0.020000 & 0.000198666667 & 0.000198666667 \\ 0.030000 & 0.000445500001 & 0.000445665785 \\ 0.040000 & 0.000789333339 & 0.000789333339 \\ 0.050000 & 0.001229166688 & 0.001229332356 \\ 0.060000 & 0.001764000064 & 0.001764000064 \\ 0.070000 & 0.002392833495 & 0.002392998980 \\ \cdots & \cdots & \cdots \\ 9.930000 & 0.195556984477 & 0.195556838917 \\ 9.940000 & 0.188751735036 & 0.188751734994 \\ 9.950000 & 0.181927615399 & 0.181927471461 \\ 9.960000 & 0.175085307925 & 0.175085307886 \\ 9.970000 & 0.168225496791 & 0.168225354533 \\ 9.980000 & 0.161348867926 & 0.161348867890 \\ 9.990000 & 0.154456108941 & 0.154455968420 \\ 10.000000 & 0.147547909058 & 0.147547909026 \\ \end{array} \right). $$

Or:

For even $ j $: $$ f_j \approx \frac{\delta_t}{3} \left[ k(t_j - t_0) g_0 + 4 \sum_{\text{odd } i=1}^{j-1} k(t_j - t_i) g_i + 2 \sum_{\text{even } i=2}^{j-2} k(t_j - t_i) g_i + k(0) g_j \right] $$

For odd $ j $: $$ f_j \approx f_{j-1} + \frac{\delta_t}{2} \left[ k(\delta_t) g_{j-1} + k(0) g_j \right] $$

#include <stdio.h>
#include <math.h>

#define DELTA_T 0.01
#define MAX_T 10
#define N ((int)(MAX_T / DELTA_T) + 1)

double k(double t) {
    return exp(-t);
}

double g(double t) {
    return sin(t);
}

double analytical_solution(double t) {
    return 0.5 * (exp(-t) - cos(t) + sin(t));
}

int main() {
    double f[N];
    double t[N];
    double g_values[N];
    double k_values[N];

    for (int i = 0; i < N; i++) {
        t[i] = i * DELTA_T;
        g_values[i] = g(t[i]);
        k_values[i] = k(i * DELTA_T);
    }

    f[0] = 0.0;

    for (int j = 1; j < N; j++) {
        if (j == 1) {
            f[j] = DELTA_T * 0.5 * (g_values[0] * k_values[j] + g_values[j] * k_values[0]);
        } else {
            double sum = g_values[0] * k_values[j] + g_values[j] * k_values[0];
            for (int i = 1; i < j; i += 2) {
                sum += 4 * g_values[i] * k_values[j - i];
            }
            for (int i = 2; i < j; i += 2) {
                sum += 2 * g_values[i] * k_values[j - i];
            }
            f[j] = (DELTA_T / 3.0) * sum;
        }
    }

    printf("t\tApprox. f(t)\tAnalytical f(t)\n");
    for (int j = 0; j < N; j++) {
        double t_current = t[j];
        double f_analytical = analytical_solution(t_current);
        printf("%.3f\t%.12f\t%.12f\n", t_current, f[j], f_analytical);
    }

    return 0;
}

$$\left( \begin{array}{ccc} \text{t} & \text{Analytical f(t)} & \text{Numerical f(t)} \\ \hline 0.000 & 0.000000000000 & 0.000000000000 \\ 0.010 & 0.000049999167 & 0.000049833333 \\ 0.020 & 0.000198666667 & 0.000198666667 \\ 0.030 & 0.000362673823 & 0.000445500001 \\ 0.040 & 0.000789333339 & 0.000789333339 \\ 0.050 & 0.001080048017 & 0.001229166688 \\ 0.060 & 0.001764000064 & 0.001764000064 \\ 0.070 & 0.002177481977 & 0.002392833495 \\ \cdots & \cdots & \cdots \\ 9.940 & 0.188751734994 & 0.188751735036 \\ 9.950 & 0.183576188924 & 0.181927615399 \\ 9.960 & 0.175085307886 & 0.175085307925 \\ 9.970 & 0.169931296344 & 0.168225496791 \\ 9.980 & 0.161348867890 & 0.161348867926 \\ 9.990 & 0.156218452224 & 0.154456108941 \\ 10.000 & 0.147547909026 & 0.147547909058 \\ \end{array} \right).$$

If we use the 10-point Gauss-Legendre quadrature: let $\{x_i\}_{i=1}^{10}$ and $\{w_i\}_{i=1}^{10}$ denote the nodes and weights for the 10-point Gauss-Legendre quadrature on the interval $[-1, 1]$. The quadrature approximates $[a, b]$:

$$ \int_a^b f(x) \, dx \approx \sum_{i=1}^{10} w_i f\left(\frac{b-a}{2} x_i + \frac{a+b}{2}\right) \frac{b-a}{2} $$ $$ \int_{t_m}^{t_{m+1}} k(t_j - s) g(s) \, ds \approx \sum_{i=1}^{10} w_i k\left(t_j - \left(\frac{t_{m+1} - t_m}{2} x_i + \frac{t_{m+1} + t_m}{2}\right)\right) g\left(\frac{t_{m+1} - t_m}{2} x_i + \frac{t_{m+1} + t_m}{2}\right) \frac{t_{m+1} - t_m}{2} $$ $$ \int_{t_m}^{t_{m+1}} k(t_j - s) g(s) \, ds \approx \sum_{i=1}^{10} w_i k\left(t_j - \left(\frac{\delta_t}{2} x_i + t_m + \frac{\delta_t}{2}\right)\right) g\left(\frac{\delta_t}{2} x_i + t_m + \frac{\delta_t}{2}\right) \frac{\delta_t}{2} $$

So, for, $[t_m, t_{m+1}]$, $m = 0, 1, 2, \ldots, j-1$, we obtain:

$$ f_j \approx \sum_{m=0}^{j-1} \sum_{i=1}^{10} w_i k\left(t_j - \left(\frac{\delta_t}{2} x_i + t_m + \frac{\delta_t}{2}\right)\right) g\left(\frac{\delta_t}{2} x_i + t_m + \frac{\delta_t}{2}\right) \frac{\delta_t}{2} $$

#include <stdio.h>
#include <math.h>

#define DELTA_T 0.01
#define MAX_T 10
#define N ((int)(MAX_T / DELTA_T) + 1)

#define GAUSS_POINTS 10

const double gauss_weights[GAUSS_POINTS] = {
    0.2955242247147529, 0.2955242247147529,
    0.2692667193099963, 0.2692667193099963,
    0.2190863625159820, 0.2190863625159820,
    0.1494513491505806, 0.1494513491505806,
    0.0666713443086881, 0.0666713443086881
};

const double gauss_nodes[GAUSS_POINTS] = {
    -0.1488743389816312, 0.1488743389816312,
    -0.4333953941292472, 0.4333953941292472,
    -0.6794095682990244, 0.6794095682990244,
    -0.8650633666889845, 0.8650633666889845,
    -0.9739065285171717, 0.9739065285171717
};

double k(double t) {
    return exp(-t);
}

double g(double t) {
    return sin(t);
}

double analytical_solution(double t) {
    return 0.5 * (exp(-t) - cos(t) + sin(t));
}

int main() {
    double f[N];
    double t[N];
    double g_values[N];

    for (int i = 0; i < N; i++) {
        t[i] = i * DELTA_T;
        g_values[i] = g(t[i]);
    }

    f[0] = 0.0;

    for (int j = 1; j < N; j++) {
        double tj = t[j];
        double sum = 0.0;

        for (int m = 0; m < j; m++) {
            double tm = t[m];
            double tm1 = t[m + 1];

            double mid = 0.5 * (tm + tm1);
            double half_width = 0.5 * (tm1 - tm);

            for (int i = 0; i < GAUSS_POINTS; i++) {
                double xi = mid + half_width * gauss_nodes[i];
                double weight = gauss_weights[i];
                sum += weight * k(tj - xi) * g(xi) * half_width;
            }
        }

        f[j] = sum;
    }

    printf("t\tApprox. f(t)\tAnalytical f(t)\n");
    for (int j = 0; j < N; j++) {
        double t_current = t[j];
        double f_analytical = analytical_solution(t_current);
        printf("%.3f\t%.20f\t%.20f\n", t_current, f[j], f_analytical);
    }

    return 0;
}

$$\left( \begin{array}{ccc} \text{t} & \text{Analytical f(t)} & \text{Numerical f(t)} \\ \hline 0.000 & 0.00000000000000000000 & 0.00000000000000000000 \\ 0.010 & 0.00004983333333472025 & 0.00004983333333475389 \\ 0.020 & 0.00019866666675530159 & 0.00019866666675528825 \\ 0.030 & 0.00044550000100816056 & 0.00044550000100813600 \\ 0.040 & 0.00078933333898971447 & 0.00078933333898971100 \\ 0.050 & 0.00122916668821304565 & 0.00122916668821303199 \\ 0.060 & 0.00176400006424457163 & 0.00176400006424456274 \\ 0.070 & 0.00239283349510070915 & 0.00239283349510072824 \\ \cdots & \cdots & \cdots \\ 9.930 & 0.19555698447728264711 & 0.19555698447728322997 \\ 9.940 & 0.18875173503583422829 & 0.18875173503583472789 \\ 9.950 & 0.18192761539890486833 & 0.18192761539890450750 \\ 9.960 & 0.17508530792480531124 & 0.17508530792480481164 \\ 9.970 & 0.16822549679109002896 & 0.16822549679109022325 \\ 9.980 & 0.16134886792614003270 & 0.16134886792614028250 \\ 9.990 & 0.15445610894056138740 & 0.15445610894056166496 \\ 10.000 & 0.14754790905842288251 & 0.14754790905842257720 \\ \end{array} \right).$$