5
$\begingroup$

Working on a problem in quantum mechanics, I rediscovered a notion of a product integral, which I was not familiar with before. In a simple case, when the integrand is a real-valued scalar function, then by taking a logarithm of it, the product integral can be in a straightforward way reduced to the regular integral that can be evaluated in Mathematica using Integrate and NIntegrate. The problem apparently becomes more complicated if the integrand is a matrix function. In particular, I am interested in a matrix product integral of the form: $$\prod_0^1\exp\!\left(\left[ \begin{array}{cc} 0 & -1 \\ f(x) & 0\\ \end{array} \right]\right)^{\!dx}=\lim_{\Delta x\to 0^+} \left(\,\prod _{x=0,\,\Delta x,\,2\Delta x,\,\dots}^1 \exp\left(\Delta x\left[ \begin{array}{cc} 0 & -1 \\ f(x) & 0\\ \end{array} \right]\right)\right),$$ where $\prod$ denotes the matrix product (Dot), $\exp(\cdot)$ denotes the matrix exponential (MatrixExp), and $f(x)$ is a real-valued, piecewise-smooth function (a particular form of that function depends on the problem we are trying to solve). The notation on the left is borrowed from here; as usual, we write $dx$ to hint at an infinitesimal $\Delta x$.

As far as I can see, this kind of integrals is not directly supported by Mathematica. Symbolic evaluation of such integrals seems a difficult problem even for very simple functions $f(x)$ — I could not find any developed theory for that. For now, I would like to find a way to evaluate matrix product integrals numerically to arbitrary precision. Could you suggest a way to do it?

Note: The matrix exponent from the product can be evaluated as follows: $$\small\exp\left(\left[ \begin{array}{cc} 0 & -1 \\ f(x) & 0\\ \end{array} \right]\right)^{\!h}=\exp\left(h\left[ \begin{array}{cc} 0 & -1 \\ f(x) & 0\\ \end{array} \right]\right)=\left[ \begin{array}{cc} \cos \left(h\sqrt{f(x)}\right) & -h\operatorname{sinc} \left(h\sqrt{f(x)}\right) \\ \sqrt{f(x)} \sin \left(h\sqrt{f(x)}\right) & \cos \left(h \sqrt{f(x)}\right) \\ \end{array} \right]$$ Note: Originally I posted a matrix that had $1$ rather than $-1$ in its top-right corner. Now I’ve fixed this error, and the expression above now correctly contains trigonometric rather than hyperbolic functions.

Update: I found a closed form solution for a quadratic polynomial $f(x)=ax^2+bx+c$; notation $D_\nu(z)$ stands for ParabolicCylinderD. Closed form (see source)

Improve this question
$\endgroup$
10
  • $\begingroup$ My related question at Math.SE: math.stackexchange.com/q/4329847/19661 $\endgroup$ Commented Dec 13, 2021 at 0:54
  • $\begingroup$ For this example, why not use ∏g = exp(∫ log g), so that ∏exp(M) = exp(∫ log exp M) = exp(∫ M)? $\endgroup$ Commented Dec 13, 2021 at 9:33
  • 1
    $\begingroup$ @evanb It is not possible to reduce this matrix product to a sum simply by taking the logarithm because the matrix factors in the product do not commute in general. The correct result can be obtained in form of an infinite series using the Baker–Campbell–Hausdorff formula, which is quite complicated even when there are just two factors in the matrix product: en.wikipedia.org/wiki/… $\endgroup$ Commented Dec 13, 2021 at 21:47
  • 2
    $\begingroup$ Well, this is the $T$-exponent, if I see it right. Indeed, well known in physics. It describes a time-evolution operator. A practical method is to write a differential equation for it. $\endgroup$ Commented Dec 13, 2021 at 22:07
  • 3
    $\begingroup$ You can look in Mahan, Many-particle physics, Sec. 2.1.C (978-0-306-46338-9), or in Stefanucci and van Leeuwen Nonequilibrium Many-Body Theory of Quantum Systems: A Modern Introduction Sec. 3.2 (978-1-139-02397-9). $\endgroup$ Commented Dec 13, 2021 at 23:22

1 Answer 1

Reset to default
6
$\begingroup$

This is a little on the exploratory/speculative side, so caveat emptor.

Letting $f(x)=\log(1+x)$, if we try three different Riemann-like discretizations:

(* left *) lSeq = With[{prec = 35}, Table[Apply[Dot, MapThread[Function[{h, x}, MatrixExp[h {{0, 1}, {Log[1 + x], 0}}]], Through[{Differences, Most}[N[Subdivide[n], prec]]]]], {n, LinearRecurrence[{0, 2}, {8, 12}, 14]}]]; (* right *) rSeq = With[{prec = 35}, Table[Apply[Dot, MapThread[Function[{h, x}, MatrixExp[h {{0, 1}, {Log[1 + x], 0}}]], Through[{Differences, Rest}[N[Subdivide[n], prec]]]]], {n, LinearRecurrence[{0, 2}, {8, 12}, 14]}]]; (* midpoint *) mSeq = With[{prec = 35}, Table[Apply[Dot, MapThread[Function[{h, x}, MatrixExp[h {{0, 1}, {Log[1 + x], 0}}]], Through[{Differences, MovingAverage[#, 2] &}[N[Subdivide[n], prec]]]]], {n, LinearRecurrence[{0, 2}, {8, 12}, 14]}]]

(where I have used the so-called "Bulirsch sequence" $2,3,4,6,8,12,\dots$ with high-precision evaluation), and then feed the results to NumericalMath`NSequenceLimit[], I get the following results:

NumericalMath`NSequenceLimit[lSeq, Method -> {"WynnEpsilon", "Degree" -> 2}] {{1.258267041020026417604, 1.06692480900025473528}, {0.40712785463903170298, 1.139960566770598057424}} NumericalMath`NSequenceLimit[rSeq, Method -> {"WynnEpsilon", "Degree" -> 2}] {{1.258267041014204419775, 1.066924811507401634993}, {0.40712785465310075194, 1.139960566387343039431}} NumericalMath`NSequenceLimit[mSeq, Method -> {"WynnEpsilon", "Degree" -> 2}] {{1.25826704102519941889787971, 1.06692481134707652125411183}, {0.40712785448828023753731087, 1.13996056685658832633988270}} NumericalMath`NSequenceLimit[(lSeq + rSeq)/2, Method -> {"WynnEpsilon", "Degree" -> 2}] {{1.25826704102519941643151917, 1.06692481134707650233505942}, {0.407127854488280241496255464, 1.13996056685658832978970988}}

which seem to point to an answer that is approximately {{1.2582670410, 1.0669248}, {0.407127854, 1.139960566}}. Note that NumericalMath`NSequenceLimit[] uses an extrapolation method that is intended for scalar sequences, and one might get better results if one uses methods specially adapted for matrix sequences (e.g. this one). One could also play around with the discretization used; I chose the Bulirsch sequence since the exponential growth is not as terrible as the sequence $2,4,8,16,\dots$ (sometimes referred to as the "Romberg sequence" in these contexts).

Improve this answer
$\endgroup$
6
  • $\begingroup$ Thanks. Yes, I tried a similar approach with NSequenceLimit, but the convergence is painfully slow (but the result seems to agree with the one obtained by solving the whole problem numerically rather than symbolically, using NDSolve). $\endgroup$ Commented Dec 16, 2021 at 2:54
  • $\begingroup$ @Vladimir, yes, I expected this kind of discretization to have relatively slow convergence even after acceleration. You mention in a comment that you are considering linear and quadratic $f(x)$; are these the particular cases you are interested in, or are there others? $\endgroup$ Commented Dec 16, 2021 at 6:40
  • $\begingroup$ I am mostly interested in polynomials because they can readily be used to approximate other smooth functions. $\endgroup$ Commented Dec 16, 2021 at 20:09
  • 1
    $\begingroup$ I got much better convergence by using powers of $2$ to subdivide the interval. $\endgroup$ Commented Dec 18, 2021 at 4:46
  • 1
    $\begingroup$ I've added a closed form for quadratic polynomials to my question above. It is pretty heavy. $\endgroup$ Commented Dec 19, 2021 at 3:54

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.