5
$\begingroup$

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}]

enter image description here

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}]

enter image description here

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$.

$\endgroup$
  • $\begingroup$ To clarify: do you need the sum of the positive eigenvalues (as your text says: positive real part and zero imaginary part) or the sum of the eigenvalues with positive real part (as your code says: positive real part and arbitrary imaginary part)? $\endgroup$ – Roman May 15 at 9:09
  • $\begingroup$ @Roman, sorry, positive real part yes. I'll clarify that, thanks. $\endgroup$ – KraZug May 15 at 10:17
  • 2
    $\begingroup$ There is no need to Sow/Reap the points at which totalPosEigs is called, as these are stored as DownValues anyway. You can extract them with pts = Cases[DownValues[totalPosEigs], RuleDelayed[_[totalPosEigs[x_]], y_] /; NumericQ[x] && NumericQ[y] :> {x, y}] after the code finishes. $\endgroup$ – Roman May 15 at 12:17
5
$\begingroup$

Hmmm. I would have set up this system as follows and this way, it is about 5 times as fast:

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}};

n = Length[Q];
cA = With[{code = N@A}, Compile[{{x, _Real}}, code, CompilationTarget -> "C"]];
cQ = With[{code = N@Q}, Compile[{{x, _Real}}, code, CompilationTarget -> "C"]];

sysmat[x_?NumericQ] := With[{
    λ = (1. - UnitStep[-Re[#]]).# &@Eigenvalues[cA[x]]
    },
   cQ[x] - DiagonalMatrix[ConstantArray[λ, n]]
   ];

xa = -4; xb = 4;    
Y0 = {0, 0, 0, 0, -1, 0};
Ysol = NDSolveValue[{
      Y'[x] == sysmat[x].Y[x],
      Y[xa] == Y0
      },
     Y, {x, xa, xb}
     ]; // AbsoluteTiming // First

0.043133

Most essential steps:

  • Compiling definitions of A and Q.
  • Teplacing Select by UnitStep and Total by Dot.
  • Performing all relevant computation in sysmat with scoped variables so that memoization is actually not required.

Addendum

OP asked whether this would work also with ParametricNDSolveValue if A and Q depend on a further parameter, say z. My first idea was to add z as a further argument sysmat to sysmat, but upon calling ParametricNDSolveValue, this resulted in an error message

"Dependent variables {Y,sysmat[x,z]} cannot depend on parameters {z}."

I hve no idea for what reason this happens. I consider it a bug.

One can circumvent this by submitting also Y[x] as variable to the function F on the right hand side of the ODE:

F[x_?NumericQ, z_?NumericQ, Y_?VectorQ] := 
  With[{λ = (1. - UnitStep[-Re[#]]).# &@Eigenvalues[cA[x]]},
   (cQ[x] - DiagonalMatrix[ConstantArray[λ, n]]).Y
   ];

Y0 = {0, 0, 0, 0, -1, 0};
Ysol = ParametricNDSolveValue[
     {
      Y'[x] == F[x, z, Y[x]],
      Y[xa] == Y0
      },
     Y,
     {x, xa, xb},
     z
     ];

(I am well aware that the parameter z does not effect the result.)

$\endgroup$
  • $\begingroup$ Thank you Henrik. I don't think the compilation is making any real difference here, but the UnitStep/Dot magic instead of Select makes such a difference. $\endgroup$ – KraZug May 15 at 12:54
  • $\begingroup$ Actually, compiling the functions cA and cQ here takes ~0.1 seconds each, the compilation is helping the speed (when there are no interpolation functions in the matrix, else compiling needs MainEvaluate anyway) but at the trade-off of the time to compile it. I can probably rewrite my code to only need to compile once for multiple parameters. $\endgroup$ – KraZug May 15 at 15:05
  • $\begingroup$ Yes, of course: You may add just further arguments to the compiled functions to supply further parameters. $\endgroup$ – Henrik Schumacher May 15 at 15:16
  • $\begingroup$ Thanks Henrik, i rewrote my code to only compile once for all x and another parameter and to use the vectorized form. $\endgroup$ – KraZug May 16 at 18:59
  • $\begingroup$ You're welcome! $\endgroup$ – Henrik Schumacher May 16 at 19:00

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.