Numerically solve large ODE system

Question

I am new to Mathematica. I am tasked with using NDSolve to numerically solve ~10^6 non-linear, coupled ODEs. However, I am only able to solve up to 100 on my computer. My professor says I should be able to do at least 1000. I'm guessing my code is not efficient, but I'm not sure what else to do. Any ideas? I have attached my code:

latSize = 20;
μ = 3;

Do[Do[
  Subscript[kx, i, j] = i;
  Subscript[ky, i, j] = j;

  Subscript[x0, i, j] = 
   Re[1/2*Sqrt[1 - ((2*μ - (i^2 + j^2))/(2*(i^2 + j^2)))^2] E^(-I*
      ArcTan[j/i])];
  Subscript[y0, i, 
   j] = -1*Im[
     1/2*Sqrt[1 - ((2*μ - (i^2 + j^2))/(2*(i^2 + j^2)))^2] E^(-I*
       ArcTan[j/i])];
  Subscript[z0, i, j] = 1/2*((2*μ - (i^2 + j^2))/(2*(i^2 + j^2)));
  , {j, 1, latSize, 1}]
 , {i, 1, latSize, 1}]

(* Set up delta *)
G = 5;
SuperPlus[Δ1] = -G*
   Sum[Sum[Subscript[kx, i, j]*Subscript[x, i, j][t] + 
      Subscript[ky, i, j]*Subscript[y, i, j][t], {j, 1, latSize, 
      1}], {i, 1, latSize, 1}];
SuperMinus[Δ1] = -G*
   Sum[Sum[Subscript[kx, i, j]*Subscript[x, i, j][t] - 
      Subscript[ky, i, j]*Subscript[y, i, j][t], {j, 1, latSize, 
      1}], {i, 1, latSize, 1}];
SuperPlus[Δ2] = 
  G*Sum[Sum[
     Subscript[ky, i, j]*Subscript[x, i, j][t] - 
      Subscript[kx, i, j]*Subscript[y, i, j][t], {j, 1, latSize, 
      1}], {i, 1, latSize, 1}];
SuperMinus[Δ2] = -G*
   Sum[Sum[Subscript[ky, i, j]*Subscript[x, i, j][t] + 
      Subscript[kx, i, j]*Subscript[y, i, j][t], {j, 1, latSize, 
      1}], {i, 1, latSize, 1}];

(* Set up B *)
Do[Do[
   Subscript[B1, i, j] = 
    2*Subscript[kx, i, 
      j] (SuperPlus[Δ1] + 
        SuperMinus[Δ1]) + 
     2*Subscript[ky, i, 
      j] (SuperMinus[Δ2] - SuperPlus[Δ2]);
   Subscript[B2, i, j] = 
    2*Subscript[kx, i, 
      j] (SuperPlus[Δ2] + 
        SuperMinus[Δ2]) + 
     2*Subscript[ky, i, 
      j] (SuperPlus[Δ1] - SuperMinus[Δ1]);
   Subscript[B3, i, j] = Subscript[kx, i, j]^2 + Subscript[ky, i, j]^2;
   , {j, 1, latSize, 1}], {i, 1, latSize, 1}];

(* Set up DE *)
Do[Do[
   Subscript[x, i, j]'[t] == 
    Subscript[B2, i, j]*Subscript[z, i, j][t] - 
     Subscript[B3, i, j]*Subscript[y, i, j][t];
   Subscript[y, i, j]'[t] == 
    Subscript[B3, i, j]*Subscript[x, i, j][t] - 
     Subscript[B1, i, j]*Subscript[z, i, j][t];
   Subscript[z, i, j]'[t] == 
    Subscript[B1, i, j]*Subscript[y, i, j][t] - 
     Subscript[B2, i, j]*Subscript[x, i, j][t];
   , {j, 1, 2, 1}], {i, 1, latSize, 1}];

s = NDSolve[
   Flatten[Table[{Subscript[x, i, j]'[t] == 
       Subscript[B2, i, j]*Subscript[z, i, j][t] - 
        Subscript[B3, i, j]*Subscript[y, i, j][t], 
      Subscript[y, i, j]'[t] == 
       Subscript[B3, i, j]*Subscript[x, i, j][t] - 
        Subscript[B1, i, j]*Subscript[z, i, j][t], 
      Subscript[z, i, j]'[t] == 
       Subscript[B1, i, j]*Subscript[y, i, j][t] - 
        Subscript[B2, i, j]*Subscript[x, i, j][t], 
      Subscript[x, i, j][0] == Subscript[x0, i, j], 
      Subscript[y, i, j][0] == Subscript[y0, i, j], 
      Subscript[z, i, j][0] == Subscript[z0, i, j]}, {i, 1, 
      latSize}, {j, 1, latSize}]], 
   Flatten[Table[{Subscript[x, i, j], Subscript[y, i, j], Subscript[z,
       i, j]}, {i, 1, latSize}, {j, 1, latSize}]], {t, 0, 1}, 
   Method -> {"EquationSimplification" -> Solve}];

In vector form the equation would be: $\frac{d \vec{r_i}}{dt} = \vec{B_i}\times \vec{r_i}$ Where $\vec{B_i} = \vec{B_i}(\vec{r_1},\vec{r_2},...)$

So if we define $\vec{B} = (B1,B2,B3)$, where the subscript i has been omitted. The B components are of the form: $B_1 = 2*k_x(\Delta1^+ + \Delta1^-) + 2*k_y(\Delta2^- - \Delta2^+)$ Where the deltas have the form: $\Delta1^+ = -G*\sum_{i} k_x*x+k_y*y$

Answer 1

Since OP doesn't give the background information of the equation system, I can just make some limited optimization. Anyway, it's clear that OP's problem is related to this one: large symbolic ODE system is burdensome, not only for the generation and storage of the system, but also for pre-process of NDSolve. (In some cases it even triggers bug, see here for an example. ) So, let's rewrite the system in vectorized form:

b1 = 2 kx (delta[1, plus] + delta[1, minus]) + 2 ky (delta[2, minus] - delta[2, plus]);
b2 = 2 kx (delta[2, minus] + delta[2, plus]) + 2 ky (delta[1, plus] - delta[1, minus]);
b3 = kx^2 + ky^2;

cross = {b1, b2, b3}\[Cross]{x, y, z};

deltarule = {delta[1, plus] :> -G Total@Flatten@(kx x + ky y), 
             delta[1, minus] :> -G Total@Flatten@(kx x - ky y), 
             delta[2, plus] :> G Total@Flatten@(ky x - kx y), 
             delta[2, minus] :> -G Total@Flatten@(ky x + kx y)};

latSize = 20;
μ = 3; G = 5;
xy0 = Table[
   1./2 Sqrt[1 - ((2 μ - (i^2 + j^2))/(2 (i^2 + j^2)))^2] E^(-I ArcTan[j/i]), {i, 
    latSize}, {j, latSize}];
x0 = Re@xy0; y0 = -Im@xy0;
z0 = Table[(2. μ - (i^2 + j^2))/(2 (2 (i^2 + j^2))), {i, latSize}, {j, latSize}];
krule = {kxlst -> Table[i, {i, latSize}, {j, latSize}], 
         kylst -> Table[j, {i, latSize}, {j, latSize}]};

rhsfunc = (Hold@
         Compile[{{r, _Real, 3}}, 
          With[{kx = kxlst, ky = kylst}, Module[{x, y, z}, {x, y, z} = r; #]], 
          CompilationOptions -> {"InlineExternalDefinitions" -> True}, 
          RuntimeOptions -> "EvaluateSymbolically" -> False] /. deltarule /. krule // 
      ReleaseHold) &@cross;
tend = 10^-4;
sol = NDSolveValue[{r'[t] == rhsfunc@r@t, r[0] == {x0, y0, z0}}, 
    r, {t, 0, tend}]; // AbsoluteTiming
(* {0.288241, Null} *)

test = Table[sol[t], {t, 0, tend, tend/100}];
ListLinePlot[test[[All, 1, 1, 1]], DataRange -> {0, tend}]