# Solving a partial integro-differential equation numerically

I am trying to find the solutions for a probability density $$p(x,t)$$, governed by,

$$\frac{\partial p(x,t)}{\partial t} =\int_{-\infty}^\infty dx' \; \Lambda(x-x')\frac{\partial^2 p(x',t)}{\partial x'^2}$$ for any function $$\Lambda(x)$$, which can even be given numerically.

• It's a linear problem, and any linear problem can be converted to linear algebra. In this case it will be of the form $\partial_t \vec{p} = \hat{M} \vec{p}$, where $\vec{p}$ is the state vector, and matrix $\hat{M}$ represents your linear operator in the RHS. Dec 11, 2022 at 17:40
• How do you intend to represent a function that is defined all the way to plus and minus infinity? Dec 12, 2022 at 5:47
• In this case $p(x)$ is known to be probability density, so it is positive definite and $\int_{-\infty}^{\infty} p dx$=1. Since this improper integral converges one can always find a finite domain $[-L,L]$ such that it contains the integral with desired accuracy, i.e., $\int_{-L}^{L} p dx \ge 1-\epsilon$.So the infinite domain can be truncated to a finite size. Another possibility is changing the $x$ coordinate by mapping the infinite domain to a finite domain. Dec 12, 2022 at 6:29
• @MaximUmansky Can you please elaborate or provide a link about how to write this as a linear algebra problem? Is the kernel in the integral like the matrix element? If so, I also have a second derivative with respect to position, how do I handle that? Dec 12, 2022 at 6:35
• @WolfgangBangerth, you can also use Hermite polynomials as a basis to deal with the infiniteness of the domain. Dec 12, 2022 at 20:44

Let's make a simplifying assumption that the function $$p(x)$$ vanishes outside of some finite domain $$[-L,L]$$. This can be justified since $$p(x)$$ is known to be a probability density function, so it is positive definite and $$\int_{-\infty}^{\infty} pdx=1$$, and since this improper integral converges one can always find a finite domain $$[-L,L]$$ such that it contains the integral with desired accuracy, i.e., $$\int_{-L}^{L} pdx \ge 1-\epsilon$$. So the infinite domain can be truncated to a finite size.

Next, let's put a uniform grid on the domain $$[-L,L]$$, and then the function $$p(x)$$ is represented by a vector $$\vec{p}$$.

The right-hand side of the integro-differential equation is a combination of two linear operators, one for the second derivative, and the other one for the integral. For the discretized solution, each of those operators becomes a matrix operator, and the product of those matrices represents the right-hand side of the discretized equation.

For the second derivative, the simplest finite difference representation is

$$\partial^2_{xx} p \rightarrow (p_{i+1}+p_{i-1}-2 p_i)/(\delta x)^2$$

so the second derivative operator is represented by a three-diagonal matrix, let's call it $$\hat{D}$$.

The integral operator with the kernel $$\Lambda$$ acting on a function $$f(x)$$, is approximated by a sum, e.g., by the method of rectangles,

$$\int_{-L}^{L} \Lambda(x-x') f(x') dx' \rightarrow \sum_j \Lambda(x_i - x_j) f_j \delta x$$

So if the function $$f$$ on the grid is represented by a vector $$\vec{f}$$, the integral operator becomes a matrix, let's call it $$\hat{L}$$, such that

$$\hat{L}_{ij} = \Lambda(x_i - x_j) \delta x$$

In the end, the integro-differential equation discretized on the grid is cast to the form,

$$\partial_t \vec{p} = \hat{M} \vec{p},$$

where $$\hat{M}$$ is the product of matrices $$\hat{L}$$ and $$\hat{D}$$.

One can formally solve the time-evolution equation for $$\vec{p}$$ as a matrix exponent. Or one can write the solution in terms of eigenvalues and eigenvectors of matrix $$\hat{M}$$. Or one can solve it numerically as a system of ODEs with a standard package. Or one can do time-stepping by discretizing in time, e.g., by the explicit first-order Euler method, which leads to a matrix-vector product,

$$\vec{p}(t+\tau) = (\tau \hat{M} + \hat{I}) \vec{p}(t),$$

or by the implicit first-order Euler method, which leads to a linear system,

$$(\hat{I}-\tau \hat{M}) \vec{p}(t+\tau) = \vec{p}(t),$$

One should note that the solution method outlined here is not the only way to approach it. Instead, one could, for example, expand the functions $$p(x)$$ and $$\Lambda(x)$$ in Fourier series and in the end obtain a linear system for the coefficients of the $$p(x)$$ expansion.

• In order to find $\vec{p}$, do I have to solve the last equation iteratively? Dec 13, 2022 at 5:16
• It is a linear system, so if the size of it is not large you can use a direct (non-iterative) solver just fine. Dec 13, 2022 at 5:29