2
$\begingroup$

I was wondering if one could benefit from Mathematica's rich linear algebra methods for diagonalizing 2nd rank tensors. Namely, in the context of systems (fluids) comprised of capsule-like particles, in order to quantify the amount of orientational order of the particles, the Nematic order parameter $p$ is chosen, which is $\approx 1$ when the particles are aligned towards a common vector (but positions uncorrelated) and $\approx 0$ when they are orientated randomly. The definition of $p$ is usually given as the largest eigenvalue of the 2nd rank tensor $Q$ given by:

$$ Q_{ab} = \frac{1}{2n}\sum_{i=1}^n (3u_a^i u_b^i - \delta_{ab}) \tag{1} $$ where $u_a^i$ and $u_b^i$ are the a-th and b-th components of the orientation vector (normalized) of particle $i$ and $\delta$ is here the Kronecker delta.

An example to generate such a tensor with random oriented vectors (normalized):

randomVec = #*Normalize@RandomReal[{-1, 1}, 3] & (*[taken from here][1]*)
n = 2; (*particles count*)
u = Table[randomVec@1, {n}] 
matQ = Table[
   1/(2.*n)*
    Sum[3*Part[u, i, a]*Part[u, i, b] - KroneckerDelta[a, b], {i, 1, 
      n}], {a, 3}, {b, 3}];

Edited after comments:

As pointed out by Henrik Schumacher in the comments, $Q$ is expected to be always diagonalizable as it is in fact a 3-by-3 symmetric matrix, therefore the eigendecomposition is not computationally challenging. Instead, the computation of the matrix itself may pose a performance problem, and can therefore be optimized (my approach given in the question is rather the most basic of going about it.). So the question is, given we know our target, the largest eigenvalue and its eigenvector, how could we optimize the computation and diagonalization of Q?

[1]: https://mathematica.stackexchange.com/a/13040/52181

$\endgroup$
  • 1
    $\begingroup$ Could you add an example Q in Mathematica form for people to work with? $\endgroup$ – Chris K Apr 25 at 8:41
  • $\begingroup$ @ChrisK Sure, I've added one way of generating such matrix Q. $\endgroup$ – user929304 Apr 25 at 9:24
  • $\begingroup$ Q should always be diagonalizable because it is symmetric. $\endgroup$ – Henrik Schumacher Apr 25 at 9:35
  • $\begingroup$ The performance problem here is not to diagonalize the matrix Q but to compute in the first place. $\endgroup$ – Henrik Schumacher Apr 25 at 9:41
  • $\begingroup$ @HenrikSchumacher Thanks, I've re-iterated the question to incorporate your points. $\endgroup$ – user929304 Apr 25 at 9:48
3
$\begingroup$

Creating random orientations:

SeedRandom[123];
n = 1000000;
u = RandomPoint[Sphere[{0, 0, 0}], n];

A helper function to compute many KroneckerProducts in a threaded way (utilizing the undocumented function NDSolve`FEM`MapThreadDot):

f = u \[Function] NDSolve`FEM`MapThreadDot[
 ArrayReshape[u, Append[Dimensions[u], 1]],
 ArrayReshape[u, Insert[Dimensions[u], 1, 2]]
 ];

That's how I interpret your definitions of Q and p:

Q = 0.5 (ConstantArray[3./n, n].f[u] - N@IdentityMatrix[3]);
{λ, U} = Eigensystem[Q];
order = Ordering[λ];
λ = λ[[order]];
U = U[[order]];
p = λ[[-1]]

0.00141718

The corresponding eigenvector is

U[[-1]]

{-0.134387, -0.74292, 0.655751}

$\endgroup$
  • $\begingroup$ Can we still perform the eigendecomposition in an approximate manner even when $Q$ is not diagonalizable? (for particles in 3D, $Q$ is always a 3-by-3 matrix). $\endgroup$ – user929304 Apr 25 at 9:35
  • $\begingroup$ Q should always be diagonalizable because it is symmetric. $\endgroup$ – Henrik Schumacher Apr 25 at 9:35
  • $\begingroup$ I completely missed that, nice observation! I checked your definition of $Q$ matches exactly mine given in the post now. Could you kindly add how you obtain the eigenvector corresponding to that eigenvalue? $\endgroup$ – user929304 Apr 25 at 9:41
  • $\begingroup$ Thanks again, very instructive as always! $\endgroup$ – user929304 Apr 25 at 10:05
  • $\begingroup$ You're welcome! $\endgroup$ – Henrik Schumacher Apr 25 at 10:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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