0
$\begingroup$

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.

$\endgroup$
5
  • $\begingroup$ 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. $\endgroup$ Commented Dec 11, 2022 at 17:40
  • $\begingroup$ How do you intend to represent a function that is defined all the way to plus and minus infinity? $\endgroup$ Commented Dec 12, 2022 at 5:47
  • 1
    $\begingroup$ 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. $\endgroup$ Commented Dec 12, 2022 at 6:29
  • $\begingroup$ @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? $\endgroup$
    – user35952
    Commented Dec 12, 2022 at 6:35
  • 1
    $\begingroup$ @WolfgangBangerth, you can also use Hermite polynomials as a basis to deal with the infiniteness of the domain. $\endgroup$
    – nicoguaro
    Commented Dec 12, 2022 at 20:44

1 Answer 1

3
$\begingroup$

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.

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.