1
$\begingroup$

The below code runs fine for anything up to a 64^2 matrix. However, moving to 128^2 throws a segmentation fault, which, despite my best efforts to debug, I cannot find a cause for. Or, to be honest, even any output whatsoever, the kernel just crashes, debugger and all. I am aware such a fault occurs when a program tries to write to memory it shouldn't be trying to access. I have also tried to run this on a larger external cluster, and have gotten the same result, so I think I can rule out a resources issue.

id= IdentityMatrix[2]; Z=PauliMatrix[3];
n = 7;
(*Create state with entries dependent on time*)
\[Rho]iDSecond = 
  Table[Table[Subscript[\[Rho], i, j][t], {j, 1, 2^n}], {i, 1, 2^n}];
(*Hamiltonian for interaction between states*)

HamilDSecond = 
  J (-1/2)*(Flatten[ 
      Map[KroneckerProduct @@ # &, {{id, Z, Z, id, id, id, id}}], 1] +
      Flatten[
      Map[KroneckerProduct @@ # &, {{id, Z, id, Z, id, id, id}}], 1] +
      Flatten[
      Map[KroneckerProduct @@ # &, {{id, Z, id, id, Z, id, id}}], 1] +
      Flatten[
      Map[KroneckerProduct @@ # &, {{id, Z, id, id, id, Z, id}}], 1] +
      Flatten[
      Map[KroneckerProduct @@ # &, {{id, Z, id, id, id, id, Z}}], 
      1]);
(*diss[\[Rho]i_]:=-\[Gamma]*excitations*(-Flatten[Map[\
KroneckerProduct@@#&,{{Z,id,id,id}}],1].\[Rho]i.Flatten[Map[\
KroneckerProduct@@#&,{{Z,id,id,id}}],1]+\[Rho]i);*)

LDephasingSecond[\[Rho]i_] := -I (HamilDSecond.\[Rho]i - \
\[Rho]i.HamilDSecond);
matrixDSecond = LDephasingSecond[\[Rho]iDSecond];
Clear[J, \[Gamma], excitations];

eqsDSecond = 
  Flatten[Table[
    D[\[Rho]iDSecond[[i, j]], t] == matrixDSecond[[i, j]], {i, 1, 
     2^n}, {j, 1, 2^n}], 1];
varsDSecond = Flatten[\[Rho]iDSecond];
initialDSecond = 
  Table[varsDSecond[[i]] == Flatten[totalSysInitial][[i]], {i, 1, 
     Length[varsDSecond]}] /. {t -> 0};
equationsDSecond = Join[eqsDSecond, initialDSecond];
solutionDSecond = 
  DSolve[equationsDSecond, varsDSecond, t] // FullSimplify // 
   Flatten;
stateAtDSecond[t_, J_, \[Gamma]_, excitations_] = 
  N[\[Rho]iDSecond] /. solutionDSecond // Chop;


It runs fine up until DSolve, which I assume is causing the issues.

eqsDSecond is the series of differential equations we wish to solve, set equal to the corresponding entry resulting from the output matrix from the function LDephasingSecond.

initialDSecond is the series of initial values for the matrix.

Both are passed into DSolve, where they are solved in terms of t.

id is an Identity matrix, Z is the pauli-Z operator, but I don't think any of that is relevant to this error.

I can't give you an input matrix example, due to the massive size, obviously. However, all input matrices are density matrices. Different input states produce the same issue. I am hoping that this is enough, and I am just missing something obvious. Either that, or it's just too large a number of equations, or a Mathematica bug, which has happened to me more than once.

$\endgroup$
3
  • $\begingroup$ If I understand correctly, id=IdentityMatrix[2]; Z=PauliMatrix[3]; (you should edit that into your code to make it easy for people to run your code). But then the matrix HamilDSecond is diagonal and the differential equations eqsDSecond are a list of 2^(2*n) uncoupled differential equations. I would not feed all of them at once into DSolve, that is asking for trouble. $\endgroup$
    – user293787
    Commented Oct 14, 2022 at 14:37
  • $\begingroup$ So this is an issue of DSolve having too much to process at once then? How would one work around this? $\endgroup$ Commented Oct 14, 2022 at 14:50
  • $\begingroup$ I do not know exactly why you get a segmentation error, but I would not call DSolve for a system of this size. Maybe my answer below does what you want. $\endgroup$
    – user293787
    Commented Oct 14, 2022 at 15:19

1 Answer 1

1
$\begingroup$
DiagonalMatrixQ[HamilDSecond]
(* True *)

(* OP did not provide initial data, let me pick it randomly *)
SeedRandom[1];
totalSysInitial=RandomInteger[{-10,10},{2^n,2^n}];

(* the equations are uncoupled and linear,
   the solution can be written down without DSolve *)
sol1=Exp[t*matrixDSecond/ρiDSecond]*totalSysInitial;
sol=Thread[Flatten[ρiDSecond]->Flatten[sol1]]

(* {Subscript[ρ, 1,1][t]->-5,
    Subscript[ρ, 1,2][t]->-10 E^(I J t),
    Subscript[ρ, 1,3][t]->-3 E^(I J t),
    Subscript[ρ, 1,4][t]->-10 E^(2 I J t),
    ...} *)

Note: Exp is not the matrix exponential.

$\endgroup$
1
  • 1
    $\begingroup$ Thank you for this. I've previously avoided doing this due to most hamiltonians not being diagonal $\endgroup$ Commented Oct 17, 2022 at 13:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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