# Numerically solve large ODE system

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$

• 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. Commented Jun 1, 2017 at 0:11
• Ah, thanks for the tip. I've replaced the image with the code in text. Commented Jun 1, 2017 at 0:19
• 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? Commented Jun 1, 2017 at 2:29
• Hi, thanks for the help. The system is basically: dr_i/dt = B_i x r_i Commented Jun 1, 2017 at 2:52
• 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? Commented Jun 1, 2017 at 3:35

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


• 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? Commented Jun 2, 2017 at 16:37
• 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]. … Commented Jun 3, 2017 at 1:29
• …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*)]]. Commented Jun 3, 2017 at 1:33
• Oooh, I see this is very cool. How would I plot delta[1, plus] as a function of time? Is there any easy way? Commented Jun 3, 2017 at 2:06
• 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]$ Commented Jun 3, 2017 at 2:15