1
$\begingroup$

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$

$\endgroup$
  • $\begingroup$ Please do not post images of your work, especially when the images display at a size that make them difficult to read. Please post your actual Mathematica code in the form of text that can be copied and pasted into a Mathematica notebook. Without such, it will be difficult to reproduce your problem and to experiment with possible solutions. $\endgroup$ – m_goldberg Jun 1 '17 at 0:11
  • 1
    $\begingroup$ Ah, thanks for the tip. I've replaced the image with the code in text. $\endgroup$ – olabaz Jun 1 '17 at 0:19
  • $\begingroup$ The last Do loop is redundant. Then, related: mathematica.stackexchange.com/q/131411/1871 . I think the best solution is to rewrite the system in vector form, but it's tedious to figure out the vector form of the system from your code. Do you already have the vector form of the system at hand? Or at least the system expressed with traditional math notation? $\endgroup$ – xzczd Jun 1 '17 at 2:29
  • $\begingroup$ Hi, thanks for the help. The system is basically: dr_i/dt = B_i x r_i $\endgroup$ – olabaz Jun 1 '17 at 2:52
  • $\begingroup$ You need to add @xzczd in your comment, or I won't get the reminder. Can you add the specific formula in $\LaTeX$ form to your question? $\endgroup$ – xzczd Jun 1 '17 at 3:35
3
$\begingroup$

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

Mathematica graphics

| improve this answer | |
$\endgroup$
  • $\begingroup$ Thank you for the quick reply! I have been going through the code and I was wondering if you could clarify a few things: Could you explain what rhsfunc is doing? I have no experience with compiled functions and how they interact with "with", "module", and "&@cross". Also, NDSolveValue gives me the value of r(t) from 0 to tend and then it is plotted? How can I choose which kx,ky I want to plot? $\endgroup$ – olabaz Jun 2 '17 at 16:37
  • $\begingroup$ rhsfunc[r[t]] (in my code, I wrote rhsfunc@r@t, they're equivalent) is $\vec{B_i}\times \vec{r_i}$. Well, I admit compiling is advanced technique, to fully understand the rhsfunc, consider starting from the following uncompiled version: Clear@rhsfunc; core = (Hold@ Function[r, With[{kx = kxlst, ky = kylst}, Module[{x, y, z}, {x, y, z} = r; #]]] /. deltarule /. krule // ReleaseHold) &@cross; rhsfunc[r_?ArrayQ] := core[r]. … $\endgroup$ – xzczd Jun 3 '17 at 1:29
  • $\begingroup$ …You may also want to read this and this. sol[(*time*)] will evaluate to a 3 dimensional list whose size is 3×latSize×latSize. (Try sol[(*time*)]//Dimensions. ) To choose the specific $k_x$, $k_y$, you just need sol[(*time*)][[All, (*index of kx*), (*index of ky*)]]. $\endgroup$ – xzczd Jun 3 '17 at 1:33
  • $\begingroup$ Oooh, I see this is very cool. How would I plot delta[1, plus] as a function of time? Is there any easy way? $\endgroup$ – olabaz Jun 3 '17 at 2:06
  • $\begingroup$ Ah, the delta should account for all of them for example, in a 2x2 lattice: $delta[1, plus] :> -G [(1) x_1 + (1)*y_1 +(2)x_2 + (1)*y_2 + (1)*x_3 + (2)*y_3 +(2)x_4 + (2)*y_4]$ $\endgroup$ – olabaz Jun 3 '17 at 2:15

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.