An explicit numerical algorithm to the solution of Volterra integral equation of the second kind

TL;DR

This paper introduces a numerical algorithm for solving Volterra integral equations of the second kind, accommodating generalized functions like the Dirac delta, by approximating kernels and explicitly computing solutions.

Contribution

It presents a novel algorithm that handles generalized input functions and provides explicit solutions using piecewise constant kernel approximations.

Findings

The algorithm effectively approximates solutions with specified accuracy.

It handles generalized functions such as the Dirac delta.

Explicit solutions are obtained via convolution with known functions.

Abstract

This paper considers a numeric algorithm to solve the equation \begin{align*} y(t)=f(t)+\int^t_0 g(t-\tau)y(\tau)\,d\tau \end{align*} with a kernel $g$ and input $f$ for $y$. In some applications we have a smooth integrable kernel but the input $f$ could be a generalised function, which could involve the Dirac distribution. We call the case when $f=\delta$, the Dirac distribution centred at 0, the fundamental solution $E$, and show that $E=\delta+h$ where $h$ is integrable and solve \begin{align*} h(t)=g(t)+\int^t_0 g(t-\tau)h(\tau)\,d\tau \end{align*} The solution of the general case is then \begin{align*} y(t)=f(t)+(h*f)(t) \end{align*} which involves the convolution of $h$ and $f$. We can approximate $g$ to desired accuracy with piecewise constant kernel for which the solution $h$ is known explicitly. We supply an algorithm for the solution of the integral equation with specified accuracy.

Equations

$$y(t)=f(t)+{\int }_0^{t}g(t-\tau )y(\tau )d\tau$$

$$h(t)=g(t)+{\int }_0^{t}g(t-\tau )h(\tau )d\tau$$

$$y(t)=f(t)+(h\ast f)(t)$$

$$y(t)$$

$$=f(t)+\sum _{i=1}^N{w}_{i}h(t_i)+(h\ast f_1)(t)$$

$$y(t)=f(t)+{\int }_0^{t}y(\tau )g(t-\tau )d\tau ,\text{for}t\ge 0$$

$$g(t)=k\phi (t)$$

$$h(t)=g(t)+{\int }_0^{t}h(\tau )g(t-\tau )d\tau ,t\ge 0$$

$$h(t)$$

$$=\sum _{n=1}^{\infty }{g}^{\otimes n}(t)$$

$$\Vert g^{\otimes n}{\Vert }_{L_{+}^1(\mathbb{R})}=k^n\Vert \phi {\Vert }_{L_{+}^1(\mathbb{R})}^n$$

$$\phi (t)=\frac{{\beta }^{\alpha }}{\Gamma (\alpha )}{t}^{\alpha -1}{e}^{-\beta t}\text{for}t>0,\alpha >0,\beta >0$$

$$\phi (t)$$

$$\Vert h-h_a{\Vert }_1\le \frac{1}{(1-k{)}^2}\Vert g-g_a{\Vert }_1$$

$$\Vert h-h_a{\Vert }_{\infty }\le \frac{1}{(1-k{)}^2}\Vert g-g_a{\Vert }_{\infty }$$

$$g(t)$$

$${\int }_0^{\infty }|g(t)|dt=\delta \sum _{j=0}^{\infty }|{\beta }_j|<1$$

$$g(t)$$

$$\text{rect}(t)$$

$$(Lf)(t)=f(t-1)$$

$$\widehat{\text{rect}}(s)$$

$$=\frac{1-e^{-s}}{s}(s>0)$$

$$\widehat{Lf}(s)$$

$$={\int }_0^{\infty }f(t-1)e^{-st}dt$$

$$={\int }_1^{\infty }f(t-1)e^{-st}dt$$

$$={\int }_0^{\infty }f(u)e^{-s(u+1)}du$$

$$=e^{-s}\hat{f}(s)$$

$$={\hat{M}}_L(s)\hat{f}(s)$$

$$\mathcal{A}=\sum _{j=0}^{\infty }{a}_j{L}^j$$

$$\mathcal{A}f(t)=\sum _{j=0}^{\infty }{a}_{j}f(t-j)$$

$$\widehat{\mathcal{A}f}(s)$$

$$={\hat{M}}_{\mathcal{A}}(s)\hat{f}(s)$$

$$\Delta =I-L$$

$$(\Delta f)(t)=f(t)-f(t-1)$$

$$Jf(t)={\int }_0^{t}f(\tau )d\tau$$

Full text

[TABLE]

with a kernel $g$ and input $f$ for $y$. In some applications we have a smooth integrable kernel but the input $f$ could be a generalised function, which could involve the Dirac distribution. We call the case when $f=\delta$, the Dirac distribution centred at 0, the fundamental solution $E$, and show that $E=\delta+h$ where $h$ is integrable and solve

[TABLE]

The solution of the general case is then

[TABLE]

which involves the convolution of $h$ and $f$. We can approximate $g$ to desired accuracy with piecewise constant kernel for which the solution $h$ is known explicitly. We supply an algorithm for the solution of the integral equation with specified accuracy.

1. Volterra Integral Equation of the Second Kind

Applications of Hawkes process in various grounds, such as in quantitative finance and machine learning (See for instance [2, 4, 3]) requires one to study Volterra Equation of the second kind when $f=\delta$, where $\delta$ is the Dirac distribution centred at [math]. We call the solution $y(t)$ the fundamental solution of the second order Volterra equation. Using the theory of distributions by Schwartz. One can show that the fundamental solution have the form $\delta+h$ where $h$ is a $L^{1}$ function solving the equation $h=g+h*g$. It also follow from this study that the general solution of volterra equation is $y=f+h*f$ and this convolution is well defined for many examples used in the studies of Hawkes.

[TABLE]

We seek solution of

[TABLE]

In the Hawkes’s setup,

[TABLE]

where $\varphi$ is $L^{1}_{+}(\mathbb{R})$ with norm $\|\cdot\|_{1}$, that is $\varphi(t)\geq 0$ if $t\geq 0$ and $\varphi(t)=0$ if $t<0$ and $\int^{\infty}_{0}\varphi(t)\,dt=1$, $0<k<1$. We seek $h$

[TABLE]

In fact

[TABLE]

where $g^{\otimes n}(t)=g\otimes g\otimes g\cdots\otimes g$ is $n$-fold convolution, $g^{\otimes 1}=g.$

Remark*.*

• •

The $n$-fold convolution always exists in $L^{1}_{+}(\mathbb{R})$ and

[TABLE]

• •

In very few cases analytic expressions for $h$ are available. If $\varphi(t)=\theta e^{-\theta t}$ for $t\geq 0$, then $h(t)=k\theta e^{-(1-k)\theta t}$ for $t\geq 0$. Analytic expressions are also available for

[TABLE]

[TABLE]

For other choices, numerical procedure are needed.

• •

If $g_{a}$ is an approximation of $g$ and $\|g\|_{1}$, $\|g_{a}\|_{1}\leq k<1$, then

[TABLE]

[TABLE]

where $\|h\|_{1}=\int^{\infty}_{0}|h(t)|\,dt$ and $\|h\|_{\infty}=\sup_{t\geq 0}|h(t)|$. The proofs of (1.6) and (1.7) are in appendix.

The key result here would be to find $h$ explicitly where

[TABLE]

and

[TABLE]

In our explicit calculation here, the $\beta_{j}\geq 0$. We provide an algorithm for $h$, which involves no approximation.

2. Part I : Special Case

We study the case

[TABLE]

and then derive the case in (1.9) from it.

Define

[TABLE]

[TABLE]

when $f\in L^{1}_{+}(\mathbb{R})$. In (2.3) we note that $(Lf)(t)=0$ if $t<1$. We use the Laplace transform

[TABLE]

[TABLE]

where $\widehat{M}_{L}(s)$ is called a multiplier for obvious reasons. Let us define

[TABLE]

where $L^{j}=L\circ L\circ\cdots\circ L$ ($j$ times) and $L^{j}f(t)=f(t-j)$ for $t\geq 0$. So

[TABLE]

[TABLE]

where $\widehat{M}_{\mathcal{A}}$ is the multiple of $\mathcal{A}$.

Define

[TABLE]

So

[TABLE]

and $\widehat{\Delta}f(s)=(1-e^{-s})\hat{f}(s)=\widehat{M}_{\Delta}(s)\hat{f}(s)$.

Define

[TABLE]

Hence

[TABLE]

Lemma 2.1**.**

[TABLE]

Proof.

[TABLE]

and

[TABLE]

∎

We now compute $g\otimes\cdots\otimes g=g^{\otimes n}$: we note that

[TABLE]

becomes

[TABLE]

and

[TABLE]

Lemma 2.2**.**

Let

[TABLE]

then

[TABLE]

Proof.

We have

[TABLE]

∎

where $\Gamma$ is the usual Gamma function, namely $\Gamma(x)=\int^{\infty}_{0}t^{x-1}e^{-t}dt$ for $x>0$. Using binomial expansion $(a+b)^{n}=a^{n}+\binom{n}{1}a^{n-1}b+\binom{n}{2}a^{n-2}b^{2}+\cdots+b^{n}$, we found that

[TABLE]

This leads to:

Lemma 2.3**.**

[TABLE]

for $n\geq 1$

Proof.

This is an immediate result from our previous calculation. ∎

Here are some graphics for $n=1$, $2$ and $3$.

In general it can be shown that $\text{rect}^{\otimes n}$ are unimodal functions and have maxima at $n/2$ and they are $(n-2)$ times differentiable on $(0,n)$ for $n\geq 3$. We can write

[TABLE]

and is given by the expression in (2.3). These are universal functions. They may have other good applications in numerical analysis.

We can give an alternative expression for $\gamma_{n}$. Note that

[TABLE]

and by Lemma 2.2 with $a=0$ and $\gamma=n-1$.

[TABLE]

If we select $f\equiv 1$, then

[TABLE]

and so

[TABLE]

and

[TABLE]

Lemma 2.4**.**

We have

[TABLE]

Proof.

This follows from

[TABLE]

As

[TABLE]

and so (2.20) is equivalent to (2.3). The last step to $g^{\otimes n}$ is to compute $\mathcal{A}^{n}$. ∎

Lemma 2.5**.**

[TABLE]

where $\{\beta^{n}_{r}\}$ are calculated as follows.

[TABLE]

[TABLE]

Proof.

We start with

[TABLE]

and equate coefficient of $L^{r}$ on both sides. ∎

Corollary** (Corollary to 2.5).**

We note that $\beta^{n}_{r}$ is the coefficient of $x^{r}$ in

[TABLE]

and if we put $x=1$,

[TABLE]

Theorem 2.6**.**

[TABLE]

Remark*.*

The intuition of Theorem 2.6 is that, if we want $n$-fold convolution, the way to do it is to apply $A$ $n$-times to the $\gamma$. Then $\gamma_{n}(t)$ is the $n$-fold convolution of the rectanglular function and itself. Then the operator $\mathcal{A}$ applies to $\gamma_{n}$, $n$ times.

Theorem 2.7**.**

[TABLE]

Instead of taking unit interval, we approximate in step of $\delta$. If you solve for the case $\delta=1$. Then Theorem 2.7 gives us formula for $h$.

3. Part II : Actual case for Part I

We now let

[TABLE]

Let us now define some useful operators. First, define $L_{\delta}$ ($\delta>0$)

[TABLE]

[TABLE]

[TABLE]

Lemma 3.1**.**

[TABLE]

Proof.

For the left hand side,

[TABLE]

For the right hand side,

[TABLE]

∎

Let $g$ be as in Part I. Then

[TABLE]

So

[TABLE]

when $a_{j}=\beta_{j}\delta$ for $j=0,1,\cdots$. Suppose we want to find the solution of

[TABLE]

We show:

Lemma 3.2**.**

[TABLE]

Proof.

Note that $h\delta$ is the solution of (3.7) but $g_{\delta}(t)=\frac{1}{\delta}g\left(\frac{t}{\delta}\right)$, then

[TABLE]

Therefore,

[TABLE]

or

[TABLE]

and by uniqueness of solution of the Volterra equation (1.3),

[TABLE]

which implies

[TABLE]

or

[TABLE]

and so

[TABLE]

using the definition of $S_{\delta}$ in (3.3).

∎

4. Applications

The solution to equation (1) is

[TABLE]

Example 1*.*

Suppose

[TABLE]

where $\{w_{1},\cdots,w_{N}\}$ are some weights and $t_{1}<t_{2}<\cdots<t_{N}$, $\delta$ is the diract delta function and $f_{1}\in L^{1}_{\text{loc}}(\mathbb{R}_{+})$ (This means that $f_{1}(t)=-$ for $t<0$ and for any $T>0$, $\int^{T}_{0}|f_{1}(t)|\,dt<\infty$). If $f_{1}\in L^{1}(\mathbb{R}_{+})$, that is $\int^{\infty}_{0}|f(t)|\,dt<\infty$, then $f_{1}\in L^{1}_{\text{loc}}(\mathbb{R}_{+})$. If $f(t)=2+\sin(t)$ for $t\geq 0$ and $=0$ for $t<0$, then $w_{i}=0$, $i=1,\cdots,N$ and $f_{1}(t)=2+\sin{t}\in L^{1}_{\text{loc}}(\mathbb{R}_{+})$, $f_{1}\notin L^{1}(\mathbb{R}_{+})$. Thus

[TABLE]

as $f_{1}\in L^{1}_{\text{loc}}(\mathbb{R}_{+})$, then $(h*f_{1})(t)$ is well defined (as $h\in L^{1}(\mathbb{R}_{+})$) and $(h*f_{1})\in L^{1}_{\text{loc}}(\mathbb{R}_{+})$ and

[TABLE]

This means that we need a formula for $(h*f_{1})(t)$ when $f_{1}\in L^{1}{\text{loc}}(R_{+}).$

Lemma 4.1**.**

Let $f\in L^{1}_{\text{loc}}(\mathbb{R}^{+})$, then

[TABLE]

Proof.

Let

[TABLE]

So

[TABLE]

and so

[TABLE]

[TABLE]

and so

[TABLE]

∎

This means that $h_{\delta}*f$ can be calculated via $h$.

**Comment

**

We might want to vary $\delta>0$, we could vary $h_{\delta}$ or we could keep $h$ constant (once calculated) and vary $S_{\delta}f$. If

[TABLE]

then

[TABLE]

Theorem 4.2**.**

If $f\in L^{1}_{\text{loc}}(\mathbb{R}^{+})$

[TABLE]

Proof.

We note that

[TABLE]

and the theorem is proved. ∎

Corollary**.**

If one requires $(h*f)(t)$ for $t\leq n$ (for some integer $n\geq 1$) then one may use $L^{r}$ with $r\leq n-1$.

Corollary**.**

We have

[TABLE]

and if we only need values for $t\leq n\delta$, then may use $L^{r}$ with $r\leq n-1$.

Lemma 4.3**.**

If $f_{1}$, $f_{2}$ are in $L^{1}_{\text{loc}}(\mathbb{R}^{+})$, then

[TABLE]

Proof.

We have

[TABLE]

as $f_{1}*f_{2}=f_{2}*f_{1}$, the theorem follows. ∎

Remark*.*

We have

[TABLE]

We could have $L^{r}\gamma_{n}$ functions tabulated (universally calculated) for $r$ and $n$, we only need $r\leq n-1$ in (2.24) and (4.5) as $\gamma_{n}(t)=0$ for $t\geq n$

Lemma 4.4**.**

Let $f_{1}$, $f_{2}\in L^{1}_{\text{loc}}(\mathbb{R}^{+})$.

[TABLE]

Proof.

We have

[TABLE]

and as before. ∎

We now try to simplify (4.5) stated earlier by utilising the operators defined in Lemma 3.2

[TABLE]

Lemma 4.5**.**

Let $f\in L^{1}(\mathbb{R}^{+})$,

[TABLE]

Proof.

We have for $s>0$,

[TABLE]

∎

5. Error Analysis

We note that the results in (1.6) and (1.7). If $f\in L^{1}(\mathbb{R}^{+})$, then111See proof in appendix.

[TABLE]

[TABLE]

if $f\in L^{1}_{\text{loc}}(\mathbb{R}^{+})$. Then (5.2) becomes

[TABLE]

which was proved in earlier notes.

6. Examples: Power Law and Rayleigh Kernel

Example 2*.*

Let

[TABLE]

We set

[TABLE]

then for $j\delta\leq t<(j+1)\delta$

[TABLE]

[TABLE]

and so $|g^{\prime}(t)|\leq k\theta(1+\theta)$ for all $t\geq 0$.

So

[TABLE]

So $g_{\delta}(t)$ with this choice of $\beta_{j}$ satisfies

[TABLE]

So

[TABLE]

and for $t\leq T$.

[TABLE]

Example 3*.*

[TABLE]

[TABLE]

So

[TABLE]

Again choose

[TABLE]

and then

[TABLE]

and the error analysis follows the same argument before.

Example 4*.*

Let

[TABLE]

Thus this is the case where $a_{0}=k$, $a_{j}=0$ for $j\geq 1$, so $\beta^{n}_{r}=0$ if $r\geq 1$ and $\beta^{n}_{0}=k^{n}$

[TABLE]

as $\gamma_{n}(t)=0$ for $t\geq n$. We only need a finite number of terms in (6.8) to calculate $h(t)$ for $t\leq T$. In fact for $t\leq T$

[TABLE]

7. Numerical implementation

Compute $\beta^{n}_{r}$ coefficients for $r\leq n-1$, all rest equal 0 for $n=1,2,\cdots$. Solve the function $\gamma_{n}(t)$ for each $n$ (MATLAB or C++) so that their values can be called. See equation (2.3). If $k$ is small, we can drop out terms in (2.24) with a small error that can be estimated. may lead to dropping further terms in calculation. No approximation is required, because for $0\leq t\leq T$ only a finite number of terms of the series are non zero.

8. Future outlook

We will implement our explicit algorithm into software and experiment with its behaviour. We will also experiment with real-world data by using our generative integral equation model to predict for arbitrary time point. (See the ODE version in [1]. This will contribute to not only the current scientific computing literature but also benefit the Machine Learning community.

9. Acknowledgments

This material was motivated from a problem in computational social science. We thank Behavioral Data Science group, especially Dr Marian-Andrei Rizoiu in facilitating discussions and supporting us with research environment.

10. Appendix

We now prove (1.6) and (1.7) If we write (1.1) as:

[TABLE]

[TABLE]

So

[TABLE]

So

[TABLE]

where all norms are $L^{1}$ norms for which

[TABLE]

If $g_{1}=g_{2}$ then $h_{1}=h_{2}$. This implies (1.1) has an unique solution.

Also if $g_{n}\to g$ in $L^{1}$ then $h_{n}\to h$ in $L^{1}$. If we assume $\|g_{n}\|\leq k<1$ for all $n\geq 1$, which leads to

[TABLE]

(We do not want $\|g_{n}\|\to 1$ as $n\to\infty$)

We can also show convergence in sup norms. Let us define

[TABLE]

Let $t>0$ and $s\in[0,t]$, then

[TABLE]

as

[TABLE]

So

[TABLE]

or

[TABLE]

and also

[TABLE]

Now we would like to derive (6) under $\text{sup}_{t}$, or $\|\cdot\|_{t}$ now.

For $0\leq s\leq t$, we have

[TABLE]

Thus as before

[TABLE]

and if $\|g_{1}\|$, $\|g_{2}\|\leq k<1$, then

[TABLE]

and so

[TABLE]

and so

[TABLE]

and so

[TABLE]

and letting $t\to\infty$

[TABLE]

where $\|f\|_{\mathcal{D}}=\text{sup}_{t\geq 0}|f(t)|$ This implies that if $g_{n}\to g$ in sup norm then so does $h_{n}\to h$ in sup norm, provided $\|g_{n}\|\leq k<1$ for all $n\geq 1$.

10.1. Proofs of equations (5.1) and (5.2)

[TABLE]

where $\|\cdot\|$ is either $L^{1}$ and $L^{\infty}$ norm. Then following the same arguments of proving (1.6), one has

[TABLE]