I have a smallish (e.g. 2x2 or 4x4, but ideally up to 10x10) non-symmetric square matrix $\mathbf{A}(x)$. I need to define a function $f(x)$ which is the sum of the eigenvalues with positive real part of $\mathbf{A}(x)$, for $x \in [x_a,x_b]$.
I then wish to use this function as a normalisation factor within a system of differential equations, $\mathbf{y}= (\mathbf{Q}(x) - f(x) \cdot \mathbf{I}) \,\mathbf{y}$.
For a small initial example, for the 2x2 case $\mathbf{Q}=\mathbf{A}$.
A = {{0, 1}, {-1 + x^2, 0}};
Q = A;
The eigenvalues here become imaginary for $|x|<1$, so I need to be able to detect those transitions cleanly.
Plot[Evaluate@Eigenvalues[A], {x, -4, 4}]
xa = -4; xb = 4;
yVector = {y[1][x], y[2][x]};
ICs = {y[1][xa] == 0, y[2][xa] == 1};
Clear[eee, totalPosEigs];
eee[z_] = Eigenvalues[N[A/.x->z]];
totalPosEigs[z_?NumericQ] :=
totalPosEigs[z] = (Sow[z]; Total@Select[eee[z], Re[#] > 0 &]);
yeqn = Thread[D[yVector, x] == (Q - totalPosEigs[x]*IdentityMatrix[Length[Q]]).yVector];
NDSolve[{yeqn, ICs},
Array[y, {Length[yVector]}], {x, xa, xb}]; // AbsoluteTiming
(* {0.01045, Null} *)
This works, and is fast enough for the 2x2 case. However, as the matrix $\mathbf{A}$ gets larger (here 4x4), this approach gets a bit too clunky (and Q is now 6x6):
A = {{0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1},
{-(1/ϵ^4), -((4 Cos[x])/ϵ^2), -((4 Sin[x])/ϵ^2), 0}} /. ϵ -> 1/10;
Q = {{0, 1, 0, 0, 0, 0}, {0, 0, 1, 1, 0, 0}, {-220 Cos[x], -220 Sin[x], 0, 0, 1, 0},
{0, 0, 0, 0, 1, 0}, {10000, 0, 0, -220 Sin[x], 0, 1}, {0, 10000, 0, 220 Cos[x], 0, 0}};
The eigenvalue structure is also more complicated in this case:
Plot[Evaluate@Re@Eigenvalues[A], {x, xa, xb}]
xa = -4; xb = 4;
yVector = Through[Array[y, {Length[Q]}][x]];
ICs = {y[1][xa] == 0, y[2][xa] == 0, y[3][xa] == 0, y[4][xa] == 0,
y[5][xa] == -1, y[6][xa] == 0};
Clear[eee, totalPosEigs];
eee[z_] = Eigenvalues[N[A /. x -> z]];
totalPosEigs[z_?NumericQ] :=
totalPosEigs[z] = (Sow[z]; Chop@Total@Select[eee[z], Re[#] > 0 &]);
yeqn = Thread[D[yVector, x] == (Q - totalPosEigs[x]*IdentityMatrix[Length[Q]]).yVector];
AbsoluteTiming[{sol, {pts}} =
Reap[NDSolve[{yeqn, ICs}, Array[y, {Length[yVector]}], {x, xa, xb}]];]
{0.229799, Null}
Here totalPosEigs
has been evaluated 567 times (even with caching), each evaluation is fast but they are adding up to contribute most of the time of the calculation. If I was to just use the sum of the absolute value of the eigenvalues for instance it takes a quarter of the time, so the Select
is a large part of the slowdown.
Its worth noting that the change points may have an imaginary part when they become positive, and the matrix may contain interpolation functions of $x$.
Sow
/Reap
the points at whichtotalPosEigs
is called, as these are stored asDownValues
anyway. You can extract them withpts = Cases[DownValues[totalPosEigs], RuleDelayed[_[totalPosEigs[x_]], y_] /; NumericQ[x] && NumericQ[y] :> {x, y}]
after the code finishes. $\endgroup$