Ellipsoids in geology

PLEASE DO NOT REDISTRIBUTE WITHOUT PERMISSION
2014 February 14 by Joshua R. Davis, jdavis@carleton.edu

Basic utilities

One can take the logarithm of an ellipsoid matrix using matrixLog, but this specialized ellipsoidLog function is preferred. Once logarithms are defined, we can apply all kinds of standard multivariate statistics to ellipsoid logs, by converting the logs to vectors. Some Sqrt[2] weights are thrown in, so that the Frobenius inner product on the log matrices equals the Euclidean inner product on the log vectors.
	For any ellipsoid tensor, jelinek returns the P-statistic from Jelinek (1981). This statistic measures the anisotropy of an ellipsoid. This P is always at least 1. Spheres have a P of 1. Needles and pancakes have really big P. The orientation of the ellipsoid doesn’t matter --- only the axis lengths do.
	The algorithm of Robin (2002) fits an ellipsoid to a set of elliptical sections of ambiguous relative size. The input secs is a list of sections. A section is a list of four objects: two vectors (3D, Cartesian, unit, perpendicular) and two scalars, indicating the directions and lengths of two semi-axes of an elliptical section. Robin’s tr E = 3 condition is arbitrary and unnecessary, but the subsequent normalization to det E = 1 overrides its arbitrariness. I’ve written a tweaked version of the algorithm, that skips the tr E = 3 condition altogether. My version is faster, but tends to produce worse error than Robin’s original. So the algorithm here is Robin’s original.

Clear[ellipsoidExp, ellipsoidLog, ellipsoidVectorFromLog, ellipsoidLogFromVector, ellipsoidMeanSpread, ellipsoidBHEP, jelinek, robin];

ellipsoidExp[log_] :=
	Module[{vals, vecs},
		{vals, vecs} = Eigensystem[log];
		vecs = Map[unit, vecs];
		Transpose[vecs].DiagonalMatrix[Exp[vals]].vecs];

(* Returns the matrix logarithm, assuming symmetric and positive-definite. *)
ellipsoidLog[ell_] :=
	Module[{vals, vecs},
		{vals, vecs} = Eigensystem[ell];
		vecs = Map[unit, vecs];
		Transpose[vecs].DiagonalMatrix[Log[vals]].vecs];

(* Converts between log-ellipsoid and its 6-vector version, which looks like *)
(* {(1, 1), (2, 2), (3, 3), Sqrt[2](1, 2), Sqrt[2](1, 3), Sqrt[2](2, 3)}. *)
ellipsoidVectorFromLog[ellLog_] :=
	{ellLog[[1, 1]], ellLog[[2, 2]], ellLog[[3, 3]],
		Sqrt[2] ellLog[[1,2]], Sqrt[2] ellLog[[1, 3]], Sqrt[2] ellLog[[2, 3]]};
ellipsoidLogFromVector[ellVec_] :=
	{{ellVec[[1]], ellVec[[4]] / Sqrt[2], ellVec[[5]] / Sqrt[2]},
		{ellVec[[4]] / Sqrt[2], ellVec[[2]], ellVec[[6]] / Sqrt[2]},
		{ellVec[[5]] / Sqrt[2], ellVec[[6]] / Sqrt[2], ellVec[[3]]}};

(* Returns the geometric mean ellipsoid, and a scalar measure of ellipsoid spread. *)
(* The spread is a kind of sample variance, predicated on the sample being normal? *)
ellipsoidMeanSpread[ells_] :=
	Module[{vecs, mean, covar},
		vecs = Table[ellipsoidVectorFromLog[ellipsoidLog[ells[[i]]]], {i, 1, Length[ells]}];
		mean = Mean[vecs];
		covar = Covariance[vecs];
		{MatrixExp[ellipsoidLogFromVector[mean]], Tr[covar]}];

(* A normality test statistic for ellipsoids. *)
(* Applies BHEP test statistic to the ellipsoid vectors. *)
(* beta > 0 is a smoothing parameter. Try beta = 0.5? *)
(* You can pass a list of betas to get a list of BHEP test statistics. *)
ellipsoidBHEP[ells_, beta_] :=
	Module[{n, vecs, s},
		n = Length[ells];
		vecs = Table[ellipsoidVectorFromLog[ellipsoidLog[ells[[i]]]], {i, 1, n}];
		bhep[vecs, beta]];

(* Returns the Jelinek (1981) P-statistic, measuring the anisotropy of the *)
(* given ellipsoid tensor. Spheres have P = 1. Needles and pancakes have large P. *)
jelinek[ell_] :=
	Module[{as, etas, ave},
		(* as = Eigenvalues[ell]^(-1 / 2); *)
		(* etas = 2 Log[as]; *)
		etas = -Log[Eigenvalues[ell]];
		ave = (etas[[1]] + etas[[2]] + etas[[3]]) / 3;
		E^Sqrt[2 ((etas[[1]] - ave)^2 + (etas[[2]] - ave)^2 + (etas[[3]] - ave)^2)]];

(* Returns the best-fit ellipsoid for the given ellipsoid sections, *)
(* according to the algorithm of Robin (2002). Each section is a list of *)
(* four objects: two vectors (unit, 3D, Cartesian, orthogonal) and two *)
(* scalars, indicating the directions and lengths of the two semi-axes. *)
robin[secs_] :=
	Module[{e11, e12, e13, e22, e23, c, form, f, formVars, vars, gradF, matF, vecB, sol, e},
		form = {{e11, e12, e13}, {e12, e22, e23}, {e13, e23, 3 - e11 - e22}};
		f = Total[Table[
			(secs[[i, 1]].form.secs[[i, 1]] - c[i] secs[[i, 3]]^-2)^2 + (secs[[i, 1]].form.secs[[i, 2]])^2 + (secs[[i, 2]].form.secs[[i, 2]] - c[i] secs[[i, 4]]^-2)^2,
			{i, 1, Length[secs]}]];
		formVars = Variables[form];
		vars = Join[formVars, Array[c, {Length[secs]}]];
		gradF = D[f, {vars}];
		matF = Table[D[gradF[[i]], {vars}], {i, 1, Length[vars]}];
		vecB = gradF /. Table[vars[[i]] -> 0, {i, 1, Length[vars]}];
		sol = -Inverse[matF].vecB;
		e = form /. Table[formVars[[i]] -> sol[[i]], {i, 1, Length[formVars]}];
		determinantOne[e]];

Here’s some stuff taken from elsewhere.

(* Returns the octahedral shear strain from the Finger tensor F Ft. *)
octahedralShearStrainOLDOLD[ffT_] :=
	Module[{l},
		l = Eigenvalues[ffT];
		Sqrt[3] / 8 * (Log[l[[1]] / l[[2]]]^2 + Log[l[[2]] / l[[3]]]^2 + Log[l[[1]] / l[[3]]]^2)];

octahedralShearStrainOLD[semiaxes_] :=
	Module[{logs},
		logs = Log[semiaxes];
		Sqrt[3] (logs.logs - logs[[1]] logs[[2]] - logs[[2]] logs[[3]] - logs[[3]] logs[[1]])];

octahedralShearStrainOLDNEW[semiaxes_] :=
	Module[{logs},
		logs = Log[semiaxes];
		Sqrt[3] Sqrt[logs.logs - logs[[1]] logs[[2]] - logs[[2]] logs[[3]] - logs[[3]] logs[[1]]]];

Clear[ellipsoidVectorFromEllipsoid, ellipsoidFromVector, nectorFromNellipsoid, nellipsoidFromNector];
ellipsoidVectorFromEllipsoid[ell_] :=
	ellipsoidVectorFromLog[ellipsoidLog[ell]];
ellipsoidFromVector[vec_] :=
	ellipsoidExp[ellipsoidLogFromVector[vec]];
nectorFromNellipsoid[nell_] :=
	Module[{vec},
		vec = ellipsoidVectorFromEllipsoid[nell];
		{vec[[1]], vec[[2]], vec[[4]], vec[[5]], vec[[6]]}];
nellipsoidFromNector[nec_] :=
	ellipsoidExp[ellipsoidLogFromVector[
		{nec[[1]], nec[[2]], -nec[[1]] - nec[[2]], nec[[3]], nec[[4]], nec[[5]]}]]

Clear[octahedralShearStrain, lodesParameter];

(* Octahedral shear strain, which is one measure of ellipsoid shape. *)
(* The literature may disagree about the scaling 1 / Sqrt[3]. *)
(* This scaling agrees with Sarah and Vasilis. *)
(* According to Ramsay and Huber, p. 200: *)
(* The principal strains are e1, e2, e3. The semi-axis lengths are 1 + e1, etc. *)
(* The natural or logarithmic strains are epsilon1 = ln(1 + e1), *)
(* and it is these natural strains that enter into the calculations below. *)
(* According to Nadai (1963) this measures the work of a coaxial deformation. *)
(* See also Mookerjee and Peek (2014). *)
octahedralShearStrain[semiaxes_] :=
	Module[{logs},
		logs = Log[semiaxes];
		Sqrt[((logs[[1]] - logs[[2]])^2 + (logs[[2]] - logs[[3]])^2 + (logs[[3]] - logs[[1]])^2) / 3]];

(* Lode's parameter nu, which is one measure of the shape of an ellipsoid. *)
lodesParameter[semiaxes_] :=
	Module[{logs, l1, l2, l3},
		logs = Log[semiaxes];
		l1 = Max[logs];
		l3 = Min[logs];
		l2 = Total[logs] - l1 - l3;
		(2 l2 - l1 - l3) / (l1 - l3)];

Here are some functions to generate random ellipsoids according to various distributions of questionable utility.

Clear[ellipsoidNormalSemiaxes, ellipsoidUniformSemiaxes, ellipsoidUniformNorm, ellipsoidNormal];

(* Returns a randomly generated ellipsoid tensor. *)
(* Orientation is uniform on SO(3). *)
(* Semiaxis logs are normal with prescribed mean, stddev. *)
ellipsoidNormalSemiaxes[meanLog_, stddevLog_] :=
	Module[{rot, as},
		rot = rotationUniform[];
		as = E^Table[RandomVariate[NormalDistribution[meanLog, stddevLog]], {3}];
		Transpose[rot].DiagonalMatrix[as^-2].rot];

(* Returns a randomly generated ellipsoid tensor. *)
(* Orientation is uniform on SO(3). *)
(* Semiaxis logs are uniform with prescribed bounds. *)
ellipsoidUniformSemiaxes[logA_, logB_] :=
	Module[{rot, as},
		rot = rotationUniform[];
		as = E^Table[RandomReal[{logA, logB}], {3}];
		Transpose[rot].DiagonalMatrix[as^-2].rot];

(* Returns a randomly generated ellipsoid tensor. *)
(* Ellipsoid vector is uniform in ball of radius b in R^6. *)
ellipsoidUniformNorm[b_] :=
	Module[{vec},
		vec = Table[RandomReal[{-b, b}], {6}];
		While[vec.vec > b^2,
			vec = Table[RandomReal[{-b, b}], {6}]];
		MatrixExp[ellipsoidLogFromVector[vec]]];

(* Returns a randomly generated ellipsoid tensor. *)
(* Ellipsoid vector is multivariate normal in R^6 *)
(* with prescribed mean and covariance matrix. *)
ellipsoidNormal[mean_, covariance_] :=
	MatrixExp[ellipsoidLogFromVector[RandomVariate[MultinormalDistribution[mean, covariance]]]];

Deformable and rigid ellipsoids

The functions in this section compute the rigid (Jeffery, 1922) and deformable (Eshelby, 1957; Bilby et al., 1975) ellipsoid dynamics. For example, you give defEByFInverse an ellipsoid tensor E, an ambient velocity gradient tensor L, a viscosity ratio r, and a step count n. The function deforms the ellipsoid for one unit of time, and outputs the resulting ellipsoid tensor. We begin with common functions and rigid functions.

Clear[ellipsoidCorrectedQ, ellipsoidQA, ellipsoidRigidWTilde, ellipsoidRigid];

(* Helper functions for ellipsoidRigid. See Davis et al. (2013). *)
ellipsoidCorrectedQ[q_, a_, d_] :=
	Module[{newQ, b, t, r},
		If[Not[Element[a[[1]], Reals] && a[[1]] > 0 && Element[a[[2]], Reals]
				&& a[[2]] > 0 && Element[a[[3]], Reals] && a[[3]] > 0],
			Throw[{"ellCorrectedQ: q, a:", q, a}],
			newQ = q;
			(* Perform the a1 = a2 correction. *)
			If[a[[1]] == a[[2]],
				b = newQ.d.Transpose[newQ];
				If[b[[2, 2]] == b[[1, 1]],
					t = Pi / 4,
					t = 0.5 ArcTan[2 b[[1, 2]] / (b[[2, 2]] - b[[1, 1]])]];
				r = {{Cos[t], -Sin[t], 0}, {Sin[t], Cos[t], 0}, {0, 0, 1}};
				newQ = r.newQ];
			(* Perform the a1 = a3 correction. *)
			If[a[[1]] == a[[3]],
				b = newQ.d.Transpose[newQ];
				If[b[[1, 1]] == b[[3, 3]],
					t = Pi / 4,
					t = 0.5 ArcTan[2 b[[1, 3]] / (b[[1, 1]] - b[[3, 3]])]];
				r = {{Cos[t], 0, Sin[t]}, {0, 1, 0}, {-Sin[t], 0, Cos[t]}};
				newQ = r.newQ];
			(* Perform the a2 = a3 correction. *)
			If[a[[2]] == a[[3]],
				b = newQ.d.Transpose[newQ];
				If[b[[3, 3]] == b[[2, 2]],
					t = Pi / 4,
					t = 0.5 ArcTan[2 b[[2, 3]] / (b[[3, 3]] - b[[2, 2]])]];
				r = {{1, 0, 0}, {0, Cos[t], -Sin[t]}, {0, Sin[t], Cos[t]}};
				newQ = r.newQ];
			newQ]];
ellipsoidQA[e_, d_] :=
	Module[{aNeg2, q, a},
		(* Diagonalize Q using a rotation matrix. The unit Map is important when E is exact! *)
		{aNeg2, q} = Eigensystem[e];
		a = aNeg2^(-1 / 2);
		q = Map[unit, q];
		If[Det[q] < 0.0, q[[3]] = -q[[3]]];
		{ellipsoidCorrectedQ[q, a, d], a}];
ellipsoidRigidWTilde[q_, a_, d_] :=
	Module[{dTilde, wTilde},
		dTilde = q.d.Transpose[q];
		wTilde = {{0, (a[[1]]^2 - a[[2]]^2) / (a[[1]]^2 + a[[2]]^2) dTilde[[1, 2]],
			(a[[1]]^2 - a[[3]]^2) / (a[[1]]^2 + a[[3]]^2) dTilde[[1, 3]]}};
		wTilde = Append[wTilde, {-wTilde[[1, 2]], 0,
			(a[[2]]^2 - a[[3]]^2) / (a[[2]]^2 + a[[3]]^2) dTilde[[2, 3]]}];
		wTilde = Append[wTilde, {-wTilde[[1, 3]], -wTilde[[2, 3]], 0}];
		wTilde];

(* Computes the Jeffery (1922) rigid ellipsoid dynamics. *)
(* e0 is the ellipsoid tensor at time t = 0. *)
(* l is the ambient velocity gradient tensor, constant. *)
(* n is the number of steps in the underlying left Lie group RK4. *)
(* Returns the ellipsoid tensor at time t = 1. *)
ellipsoidRigid[e0_, l_, n_] :=
	Module[{d, w, q, a, x},
		(* Compute constants and the starting Q. *)
		d = 0.5 (l + Transpose[l]);
		w = 0.5 (l - Transpose[l]);
		{q, a} = ellipsoidQA[e0, d];
		x = Function[qq, Module[{corr, negK},
			corr = ellipsoidCorrectedQ[qq, a, d];
			negK = Transpose[corr].ellipsoidRigidWTilde[corr, a, d].corr - w;
			negK]];
		q = passError[leftLieRungeKutta, {q, x, n}, {"ellipsoidRigid: e0, l:", e0, l}];
		Transpose[q].DiagonalMatrix[a^-2].q];

The defJ function computes Eshelby’s (1957) magic coefficients J_ij and J_i. Computing these naively requires nine numerical integrations. My old implementation used six integrations. In the typical case, defJ uses two integrations, and is 2 or 3 times as fast as my old implementation. In special cases, defJ uses no integrations, and is hundreds of times as fast as my old implementation. It may be possible to optimize the typical case further, by reducing the two integrations down to the evaluation of two elliptic integrals, for which Mathematica has optimized routines. Eshelby (1957) mentions this optimization, and his text agrees with Kellogg (1922 aka 1953), but I cannot make it produce the same results as my other methods, so I do not use it.

Clear[ellipsoidDeformableJSphere, ellipsoidDeformableJOblate, ellipsoidDeformableJProlate, ellipsoidDeformableJTypical, ellipsoidDeformableJ, ellipsoidDeformableCTilde, ellipsoidDeformableK, ellipsoidDeformable];

(* Computes the Eshelby J-tensor. See Davis et al. (2013). *)
(* Assumes a == b == c. *)
ellipsoidDeformableJSphere[a_] :=
	Module[{jTwo, jOne},
		jTwo = 4 Pi / (15 a^2) * {{3, 1, 1}, {1, 3, 1}, {1, 1, 3}};
		jOne = {4 Pi / 3, 4 Pi / 3, 4 Pi / 3};
		{jTwo, jOne}];
(* Assumes a == b > c. *)
ellipsoidDeformableJOblate[a_, c_] :=
	Module[{jTwo, jOne, j1, j2, j3, j13, j12, j11, j23, j22, j33},
		j1 = 2 Pi a^2 c / (a^2 - c^2)^(3 / 2) (ArcCos[c / a] - c / a (1 - c^2 / a^2)^(1 / 2));
		j2 = j1;
		j3 = 4 Pi - j1 - j2;
		j13 = (j3 - j1) / (3 (a^2 - c^2));
		j12 = Pi / (3 a^2) - j13 / 4;
		j11 = 3 j12;
		j23 = (j3 - j2) / (3 (a^2 - c^2));
		j22 = 4 Pi / (3 a^2) - j12 - j23;
		j33 = 4 Pi / (3 c^2) - j13 - j23;
		jTwo = {{j11, j12, j13}, {j12, j22, j23}, {j13, j23, j33}};
		jOne = {j1, j2, j3};
		{jTwo, jOne}];
(* Assumes a > b == c. *)
ellipsoidDeformableJProlate[a_, b_] :=
	Module[{jTwo, jOne, j1, j2, j3, j13, j12, j11, j23, j22, j33},
		j3 = 2 Pi a b^2 / (a^2 - b^2)^(3 / 2) (a / b (a^2 / b^2 - 1)^(1 / 2) - ArcCosh[a / b]);
		j2 = j3;
		j1 = 4 Pi - j2 - j3;
		j13 = (j3 - j1) / (3 (a^2 - b^2));
		j12 = (j2 - j1) / (3 (a^2 - b^2));
		j11 = 4 Pi / (3 a^2) - j12 - j13;
		j23 = Pi / (3 b^2) - j12 / 4;
		j22 = 3 j23;
		j33 = 4 Pi / (3 b^2) - j13 - j23;
		jTwo = {{j11, j12, j13}, {j12, j22, j23}, {j13, j23, j33}};
		jOne = {j1, j2, j3};
		{jTwo, jOne}];
(* Assumes a ≠ b ≠ c ≠ a. Doesn't use the elliptic integral optimization, so doesn't assume a > b > c. *)
ellipsoidDeformableJTypical[a_, b_, c_] :=
	Module[{jTwo, jOne, j1, j2, j3, j13, j12, j11, j23, j22, j33, theta, k, e, f, u, delta},
		(* This optimization is from Eshelby (1957). I checked it with Kellogg (1953).
		But I couldn't get it to agree with the other methods. N.B. It assumes a > b > c.
		theta = ArcSin[(1 - c^2 / a^2)^(1 / 2)];
		k = ((a^2 - b^2) / (a^2 - c^2))^(1 / 2);
		f = EllipticF[theta, k];
		e = EllipticE[theta, k];
		j1 = 4 Pi a b c / ((a^2 - b^2) (a^2 - c^2)^(1 / 2)) * (f - e);
		j3 = 4 Pi a b c / ((b^2 - c^2) (a^2 - c^2)^(1 / 2)) * (b (a^2 - c^2)^(1 / 2) / (a c) - e);*)
		delta = ((a^2 + u) (b^2 + u) (c^2 + u))^(1 / 2);
		j1 = 2 Pi a b c NIntegrate[1 / ((a^2 + u) delta), {u, 0, Infinity}];
		j3 = 2 Pi a b c NIntegrate[1 / ((c^2 + u) delta), {u, 0, Infinity}];
		j2 = 4 Pi - j1 - j3;
		j13 = (j3 - j1) / (3 (a^2 - c^2));
		j12 = (j2 - j1) / (3 (a^2 - b^2));
		j11 = 4 Pi / (3 a^2) - j12 - j13;
		j23 = (j3 - j2) / (3 (b^2 - c^2));
		j22 = 4 Pi / (3 b^2) - j12 - j23;
		j33 = 4 Pi / (3 c^2) - j13 - j23;
		jTwo = {{j11, j12, j13}, {j12, j22, j23}, {j13, j23, j33}};
		jOne = {j1, j2, j3};
		{jTwo, jOne}];
ellipsoidDeformableJ[a_] :=
	Module[{j, p},
		Which[
			a[[1]] ≠ a[[2]] && a[[2]] ≠ a[[3]] && a[[3]] ≠ a[[1]],
				ellipsoidDeformableJTypical[a[[1]], a[[2]], a[[3]]],
			a[[1]] == a[[2]] == a[[3]],
				ellipsoidDeformableJSphere[a[[1]]],
			a[[1]] == a[[2]] > a[[3]],
				ellipsoidDeformableJOblate[a[[1]], a[[3]]],
			a[[1]] > a[[2]] == a[[3]],
				ellipsoidDeformableJProlate[a[[1]], a[[3]]],
			a[[2]] == a[[3]] > a[[1]],
				j = ellipsoidDeformableJOblate[a[[2]], a[[1]]];
				p = {{0, 0, 1}, {0, 1, 0}, {1, 0, 0}};
				{p.j[[1]].p, p.j[[2]]},
			a[[2]] > a[[3]] == a[[1]],
				j = ellipsoidDeformableJProlate[a[[2]], a[[1]]];
				p = {{0, 1, 0}, {1, 0, 0}, {0, 0, 1}};
				{p.j[[1]].p, p.j[[2]]},
			a[[3]] == a[[1]] > a[[2]],
				j = ellipsoidDeformableJOblate[a[[3]], a[[2]]];
				p = {{1, 0, 0}, {0, 0, 1}, {0, 1, 0}};
				{p.j[[1]].p, p.j[[2]]},
			a[[3]] > a[[1]] == a[[2]],
				j = ellipsoidDeformableJProlate[a[[3]], a[[1]]];
				p = {{0, 0, 1}, {0, 1, 0}, {1, 0, 0}};
				{p.j[[1]].p, p.j[[2]]},
			True,
				{"error: defJ:", a}]];

(* Computes the Eshelby VGT K. See Davis et al. (2013). *)
ellipsoidDeformableCTilde[q_, a_, d_, r_, jTwo_] :=
	Module[{dTilde, b, cTilde11, cTilde22, cTilde33, cTilde12, cTilde13, cTilde23},
		dTilde = q.d.Transpose[q];
		b = IdentityMatrix[3] + (3 / (4 Pi)) (r - 1) jTwo.DiagonalMatrix[a^2];
		{cTilde11, cTilde22, cTilde33} = Inverse[b].{dTilde[[1, 1]], dTilde[[2, 2]], dTilde[[3, 3]]};
		cTilde12 = dTilde[[1, 2]] / (1 + (r - 1) (3 / (4 Pi)) (a[[1]]^2 + a[[2]]^2) jTwo[[1, 2]]);
		cTilde13 = dTilde[[1, 3]] / (1 + (r - 1) (3 / (4 Pi)) (a[[1]]^2 + a[[3]]^2) jTwo[[1, 3]]);
		cTilde23 = dTilde[[2, 3]] / (1 + (r - 1) (3 / (4 Pi)) (a[[2]]^2 + a[[3]]^2) jTwo[[2, 3]]);
		{{cTilde11, cTilde12, cTilde13}, {cTilde12, cTilde22, cTilde23}, {cTilde13, cTilde23, cTilde33}}];
ellipsoidDeformableK[q_, a_, d_, w_, r_] :=
	Module[{jTwo, jOne, cTilde, kTildeMinusWTilde},
		{jTwo, jOne} = ellipsoidDeformableJ[a];
		cTilde = ellipsoidDeformableCTilde[q, a, d, r, jTwo];
		(* Compute K~ - W~ without computing K~ or W~. *)
		kTildeMinusWTilde = Table[Table[
				(1.0 + (r - 1.0) / (4.0 Pi) (jOne[[i]] - jOne[[j]])) cTilde[[i, j]],
			{j, 1, 3}], {i, 1, 3}];
		(* Return K = Q^T (K~ - W~) Q + W. *)
		N[Transpose[q].kTildeMinusWTilde.q + w]];

(* Computes the Eshelby (1957) deformable ellipsoid dynamics. *)
(* e0 is the ellipsoid tensor at time t = 0. *)
(* l is the ambient velocity gradient tensor, constant. *)
(* r is the viscosity ratio between inclusion and surroundings. *)
(* n is the number of steps in the underlying left Lie group RK4. *)
(* Returns the ellipsoid tensor at time t = 1. *)
ellipsoidDeformable[e0_, l_, r_, n_] :=
	Module[{d, w, q, aNeg2, f0I, x, fnI},
		d = 0.5 (l + Transpose[l]);
		w = 0.5 (l - Transpose[l]);
		x = Function[fI, Module[{qa},
			qa = ellipsoidQA[Transpose[fI].e0.fI, d];
			-ellipsoidDeformableK[qa[[1]], qa[[2]], d, w, r]]];
		fnI = passError[leftLieRungeKutta, {IdentityMatrix[3], x, n}, {"ellipsoidDeformable: e0, l, r:", e0, l, r}];
		Transpose[fnI].e0.fnI];

Experimental

This is the E-dot method of Davis et al. (2013). I coded it for comparison to the same method in R. They work identically. So do the two versions of the left Lie group method.

eshelbyTest[e0_, l_, r_, n_] :=
	Module[{d, w, x, e1},
		d = 0.5 (l + Transpose[l]);
		w = 0.5 (l - Transpose[l]);
		x = Function[e,  Module[{qa, k},
			qa = ellipsoidQA[e, d];
			k = ellipsoidDeformableK[qa[[1]], qa[[2]], d, w, r];
			-Transpose[k].e - e.k]];
		e1 = passError[euclideanRungeKutta, {e0, x, n}, {"ellipsoidDeformable: e0, l, r:", e0, l, r}];
		e1];

eshelbyTest[IdentityMatrix[3], {{1, 5, -2}, {1, 1, 0}, {0, -3, -1}}, 10, 4] // MatrixForm

ellipsoidDeformable[IdentityMatrix[3], {{1, 5, -2}, {1, 1, 0}, {0, -3, -1}}, 10, 4] // MatrixForm