Governing Eq. for Thermodynamic Computing

Following my previous post on thermodynamic computing, I'm going to just write down my summary notes from section II A.1 in ref. [1]. For the most part it is just a rephrasing of the paper, so you could easily just read it, but for me I needed to write it down somewhere. There also is a more through derivation and review of Langevin dynamics in ref. [2] as well as the Fokker-Planck equations.

Langevin Equations

From my perspective the primary focus of thermodynamic computing is the use of continuous variables, that is a signal is given by $x \in \mathbb{R}$. The physics of such signal is given by stochastic dynamics¹ which contains terms related to dissipation and fluctuations. The equations which govern the dynamics are:

$$ \begin{align} \mathrm{d}x &= \frac{p}{M} \mathrm{d}t \\ \mathrm{d}p &= -\nabla U(x) - \frac{\gamma}{M} \, p \,\mathrm{d}t + \mathcal{N}\left(0, \frac{2\gamma}{\beta} \, \mathbb{I} \, \mathrm{d}t\right) \label{eq:langevin} \end{align} $$

The second equation is where the stochastic dynamics comes into play and as a result we call these types of equations stochastic differential equations (SDEs). The first equation is just the kinematic equation which links the change in generalized coordinates.

In the second equation, the first term is the conserved force due to a potential, while the second and third terms correspond to a dissipative force due to dampening-frictional term and a random force due to real-world noise, i.e., white noise that occurs at various time-scales.

The constants $\gamma$ and $\beta$ are positive and real valued, and correspond to the strength of the dissipative force while $\beta$ is usually the temperature-scale which modulates the intensity of the thermal/noise fluctuations. A smaller $\beta$ would mean more intense noise fluctuations. The $M$ corresponds to the mass which is a positive real scalar (or positive definite matrix). Depending on the relative scale of $\gamma$ and $M$ the Langevin equations are referred to as overdamped ($\gamma \gg M$) or underdamped ($\gamma \ll M$).

The eqs. $\ref{eq:langevin}$ correspond to a single trajectory in phase space. However, what we are typically most interested in is the time-evolution of the system's probability density function (PDF). This evolution is described by the Fokker-Planck equation, which is derived from the Langevin equations. The Fokker-Planck equation accounts for both the deterministic dynamics and the stochastic fluctuations captured by the Langevin equations, providing a comprehensive description of the system's statistical behavior over time. For more details on Fokker-Planck equations see ref. [2].

Fokker-Planck solution

The stationary solution of the Fokker-Planck equation, which describes the system's state at thermal equilibrium, is given by the Gibbs distribution:

$$ p(x) = \frac{1}{Z} e^{-\beta U(x)} $$

with $Z$, the partition function, ensuring normalization. This distribution emerges from the interplay between the deterministic forces that drive the system towards lower energy states and the stochastic forces that induce thermal fluctuations, highlighting the fundamental principle of equilibrium statistical mechanics.

The primary interest is now to determine what the exact PDF might be given a specific energy potential $U(x)$. In the case that $U(x)$ is quadratic, i.e., $U(x) = x^T A x + b^T x$ the Gibbs distribution becomes Gaussian²:

$$ \begin{equation} \label{eq:solution} x \sim \mathcal{N}(A^{-1}b, (\beta\,A)^{-1}) \end{equation} $$

It is important to note that $x$, along with corresponding³ momenta $p$, in this context is a vector of generalized coordinates that can represent various physical or analogous quantities, depending on the dimensionality and nature of the system being studied.

If we pay particular attention to eq. $\ref{eq:solution}$ we see that if we can draw enough samples of $x$ then we can determine $A^{-1}$, which is the inverse matrix from $U(x)$. Therefore we have a way using stochastic dynamics to invert a matrix assuming we can prepare and equilibrate the Gibbs distribution given $U(x)$. The excitement is obviously that many problems in ML/AI and science are all about inverting matrices.

Footnotes

---

1. The nuanced point to note about term stochastic dynamics is it means that the dynamics are represented by a probability density function. The distinction is made here because we don't want to just say its a random variable but rather a time-evolving pdf. ↩

2. You can show that this potential indeed results in a Gaussian distribution if you use the technique to complete the square, which is just basic algebra manipulation. You'll end up with something like $p(x) = \frac{1}{(2\pi)^{n/2} (\beta^{-1} |A^{-1}|)^{1/2}} e^{-\frac{1}{2} \beta (x - \mu)^T A (x - \mu)}$, with $\mu = A^{-1}b$ and $\Sigma = (\beta\,A)^{-1}$. The assumption is $A$ is symmetric and invertible. ↩

3. This is usually referred to as the canonical conjugate momenta (see Hamiltonian mechanics). ↩

References

[1] D. Melanson, M.A. Khater, M. Aifer, K. Donatella, M.H. Gordon, T. Ahle, G. Crooks, A.J. Martinez, F. Sbahi, P.J. Coles, Thermodynamic Computing System for AI Applications, arXiv Preprint arXiv:2312.04836 (2023).

[2] L. Peliti, S. Pigolotti, Stochastic Thermodynamics: An Introduction, Princeton University Press, 2021. URL.

Reuse and Attribution

Refresher: Bessel Functions

 Let's talk about Bessel functions! First, who was Bessel? Friedrich Bessel was a German astronomer and mathematician whose namesake was given to the solutions used to solve a specific type of differential equation. The differential equation in question is the following:

\begin{equation}x^2\frac{d^2y}{dx^2} + x\frac{dy}{dx} + (x^2-n^2)y = 0\label{eq:ode}\end{equation}

This is a second-order linear ordinary differential equation and is known as a canonical ODE because its solutions are special types of functions. For those unfamiliar with second-order linear ODEs let me explain ( I often forget so it's helpful for me as well). The "second-order" just means the highest derivative in the equation, in this case, its two. The ordinary is because there is only a single independent variable $x$. Linear indicates that we don't have terms like $y^2\frac{dy}{dx}$, in other words,  the dependent variable and its derivatives only appear to the first power and the coefficients are not in terms of the dependent variable. There is also the case that this equation is homogenous because $y$ is involved with every term.

The solutions to the Bessel equation are called Bessel functions. The functions form a basis because they are linearly independent and there are two flavors of functions called the "first-kind" and the "second-kind". The first-kind, usually denoted as $J_n(x)$, has the form:

\begin{equation}J_n(x) = \frac{1}{\pi} \int_0^{\pi} \cos(n*t - x \sin(t)) dt\end{equation}

where $n$ is an integer. This is actually a special case known as Bessel integrals. The solutions of the first kind can be multiplied by a scalar value and remain solutions to the ODE. Given Bessel functions of the first kind are written in terms of sine and cosine they are oscillatory functions. They are linearly dependent and form an orthogonal basis, given that:

$$\int_0^x J_n(x) \cdot J_m(x) dx = 0 \;  m \neq n $$

They also have recurrence relations whereby negative values for $n$ are the same as positive values for $n$. The Bessel functions of the first kind are used to obtain solutions in cylindrical coordinate representations because the solutions are the radial part. In spherical coordinates, the radial part is given by the Bessel functions of the second kind. So what do the Bessel functions of the first kind look like? Let me show it with some Julia code that would compute the Bessel function integral of the first kind for integer values of $n$.

"""
	Jₙ(x;n=1,ν=1000)

First kind Bessel numerical (i.e., trapizodial) integrals for
integer values of n.
""" 
function Jₙ(x,n=1;ν=1000)
    ∫f = 0.0
    for i=0:ν
    	t = π/ν * i
        f = cos(n * t - x * sin(t))
        # trapezoidal weights
        w = (i == 0 || i == ν) ? 0.5 : 1.0 
        ∫f += w * f
    end
    return 1/ν * ∫f
end

The code simply calculates the integral form for the first kind of Bessel function, which is possible when $n$ is an integer. I believe this is also a valid function for non-integer values of $n$ but only when $x$ is real and $\gt 0$. When we plot this function for different values of $n$ we get the following solutions:

Bessel functions (integrals) of the first kind where $n$ is an integer.

So what about Bessel functions of the second kind? When $n$ is not an integer these take the form:

\begin{equation}Y_m(x) = J_m(x) * \cos(m\pi) - \frac{J_{-m}(x)}{\sin(m\pi)}\label{eq:bessel1}\end{equation}

with $m$ being a non-integer value. As you can see, $Y_m(x)$ can be defined in terms of Bessel functions of the first kind, $J_m(x)$ and $J_{-m}(x)$. The functions in eq. (\ref{eq:bessel1}) are clear valid because they are linear combinations of first-kind solutions. One thing of importance is that Bessel functions of the second kind have a singularity at $x=0$ for all values of $m$.  Also, keep in mind that $J_m(x)$ is not the function $J_n(x)$ since $m$ is a non-integer. If $m\to +n$, i.e., is an integer and non-negative, then the solution [1] is given by:

\begin{multline}Y_n(x) =-\frac{\left(\frac{x}{2}\right)^{-n}}{\pi}\sum_{k=0}^{n-1} \frac{(n-k-1)!}{k!}\left(\frac{x^2}{4}\right)^k +\frac{2}{\pi} J_n(x) \ln \frac{x}{2} \\ -\frac{\left(\frac{x}{2}\right)^n}{\pi}\sum_{k=0}^\infty (\psi(k+1)+\psi(n+k+1)) \frac{\left(-\frac{x^2}{4}\right)^k}{k!(n+k)!}\label{eq:bessel2}\end{multline}

where $\psi(x)$ is the logarithmic derivative of the gamma function. You can review the Wikipedia entry for Bessel functions which provides more details on the derivation of this solution. This solutions $Y_n(x)$ are linearly independent from $J_n(x)$ hence why we have "second-kind" solutions Bessels equations. There is also an integral form that exists when $x$ is positive and real-valued.  So what does the plot look like for eq. (\ref{eq:bessel2})? Below shows the plot for the limiting case discussed.

Bessel functions of the second kind where $x>0$ and $n$ is an integer.

If you want to look at the Julia implementation for the Bessel functions of the second kind you can take a look at the Pluto.jl notebook here. You'll need to set up Julia+Pluto.jl to run the notebook locally. You can also use the "run in cloud" but its not worth it. For Windows users who just want to run the notebook without any setup, check out the PlutoDesktop exe. For python users you can probably just copy the raw notebook file content into an LLM and ask for it to rewrite the code in Python.

In general, the use of Bessel functions of the first kind vs the second kind often depends on the specific problem at hand and the boundary conditions involved. Both types of Bessel functions are solutions to eq. (\ref{eq:ode}), but they behave differently and hence have different applications. Here's a quick rundown:

1. First Kind ($J_n$): These functions are finite at the origin ($x=0$), and this makes them suitable for problems with boundary conditions that require the solution to be finite at the origin. For example, they often appear in problems of wave propagation and heat conduction in cylindrical or spherical symmetry where the solution must be finite at $r=0$.
2. Second Kind ($Y_n$): These functions, on the other hand, have a singularity at the origin. They are often useful in problems where the physical system has a singularity at the origin or where the boundary conditions do not require the solution to be finite at the origin. They are also used when the full set of linearly independent solutions to eq. \ref{eq:ode} are needed.

The complete solution to eq. $\ref{eq:ode}$ is a linear combination of Bessel functions, where the coefficients (i.e., amplitudes) are determined by the boundary conditions. In the context of partial differential equations such as Laplace's equation in cylindrical or spherical coordinates, the Bessel function provides solutions to the radial part whereas the spherical harmonics represent the angular parts.

References

[1] Y. L. Luke, Integrals of Bessel Functions (Courier Corporation, 2014).

[2] DLMF: Chapter 10 Bessel Functions, https://dlmf.nist.gov/10.

Reuse and Attribution

Linear Differential Equation With Constant Prefactors

I decided I wanted to revisit some techniques for solving ordinary differential equations. It's always good to through in some refreshers as I continue to write more for my blog. Here I'll focus on  homogenous differential equations that have the following form:

$$ A_1 \frac{d^2 y}{dx^2} + A_2\frac{dy}{dx} - A_3y =0 .$$

The coefficients $A_1$, $A_2$, and $A_3$ are all constant. This equation can be solved using the auxiliary equation, which is derived from the fact that such differential equations have a solution of the form $y=e^{mx}$. The auxiliary equation is:

$$ f(m) = A_0 m^{n} + A_{1} m^{n-1} + \dots+ A_{n-1}m + A_{n} = 0$$

The goal is to find the roots of the auxiliary equation, and if they are real, we can use the general solution to our homogenous differential equation which is:

$$ y(x) = C_{1}e^{m_1 x}+C_{2}e^{m_2 x} + \dots + C_n e^{m_n x} $$

Let us now solve our initial differential equation at the beginning of this post, where $A_1=10$, $A_2=1$, and $A_3=-2$. We start by rewriting the equation in terms of operator power, i.e., the degree of differentiation, this would be:

$$y(10L^2+L-2)=0,$$

$L$ is the differential operator. Now we want to plug the coefficients and operator powers into the auxiliary equation above:

$$10m^{2}+m-2$$

We can find the roots for this polynomial, which are $m=-\frac{1}{2}$ and $m=\frac{2}{5}$. Substitution of these roots into our general solution as given above yields:

$$y(x)=C_1e^{\frac{-1}{2}x}+C_2e^{\frac{2}{5}x}$$

the undetermined coefficients to our general condition, and therefore the particular solution, can be found from initial or boundary conditions to our differential equation. It is also worth mentioning that it is common to have complex roots and therefore solutions. 

References  & Additional Reading

[1] Advanced Engineering Mathematics, Wiley, 8th ed., E. Kreyszig.

[2] The Mathematics Companion, IoP, 1st ed., A.C. Fischer-Cripps.

Reuse and Attribution

Exact Differential Equations

The general form of a first-order differential equation is given by the following:

$$ M(x,y)dx + N(x,y)dy = 0$$

our differential equation is said to be exact if it satisfies the following exactness test:

$$\frac{\partial M\left(x,y\right)}{\partial y} = \frac{\partial N\left(x,y\right)}{\partial x}$$

The goal is to determine a function $f(x,y)$ that satisfies the following:

$$df = M(x,y)dx + N(x,y)dy$$

$$ \frac{\partial f\left(x,y\right)}{\partial x} = M\left(x,y\right)$$

$$ \frac{\partial f\left(x,y\right)}{\partial y} = N\left(x,y\right)$$

Let us look at the following example differential equation:

$$\left(y^{2}-2x\right)dx + \left(2xy+1\right)dy=0$$

Taking the partial derivatives of the functions corresponding to $M\left(x,y\right)$ and $N\left(x,y\right)$, we get:

$$ \frac{\partial M}{\partial y} = 2y $$

$$ \frac{\partial N}{\partial x} = 2y $$

So our differential equation is indeed exact and we now can find the total function, $f(x,y)$, whose derivative is equal to our differential equation. This is done by integrating  the functions $M\left(x,y\right)$ and $N\left(x,y\right)$,

$$M\left(x,y\right) = \frac{\partial f\left(x,y\right)}{\partial x}$$

$$f = \int{\left(y^{2}-2x\right) dx} = xy^{2}-x^{2} $$

similarly for $N\left(x,y\right)$,

$$N\left(x,y\right) = \frac{\partial f\left(x,y\right)}{\partial y}$$
$$f = \int{\left(2xy+1\right)dy} = xy^{2}+y $$

In both cases, we ignore the constant of integration. We now can identify unique terms and construct the function, $f(x,y)$, by summing these terms:

$$f\left(x,y\right) = xy^{2}-x^{2}+y=\text{constant}$$

So we have identified a function, $f\left(x,y\right)$, that is a solution to our exact differential equation.

Now for our quote:

I became an atheist because, as a graduate student studying quantum physics, life seemed to be reducible to second-order differential equations. It thus became apparent to me that mathematics, physics, and chemistry had it all and I didn't see any need to go beyond that.
-Attributed to Francis Collins but unconfirmed.

Reuse and Attribution