# Is there any way to save time using of “Cross” or “Table”?

I need to solve a equation by NDSlove, so I write the equation first, but it takes me a lot of time. The follow is my code

energy = Subdivide[0, 50, 5];
theta0 = Reverse[N[ArcSin[Subdivide[0, 1, 5]]]]
theta = Cos[
ArcSin[(R Sin[theta0])/Sqrt[
t^2 + 2 R t Cos[theta0] + R^2 Cos[theta0]^2 + R^2 Sin[theta0]^2]]];
ni = Length[theta];nj = Length[energy];
Ebare = 10;Ebaree = 15;Ebarx = 24;betae = 0.315;betaee = 0.21;betax = 0.131;R = 10;
Fj = Compile[{{x, _Real}}, (
betae (betae x)^2)/((E^(betae x) + 1) Ebare) + (
betax (betax x)^2)/((E^(betax x) + 1) Ebarx)];
Fjbar = Compile[{{x, _Real}}, ((
betaee (betaee x)^2)/((E^(betaee x) + 1) Ebaree) + (
betax (betax x)^2)/((E^(betax x) + 1) Ebarx)) ];
For[j = 2, j <= nj, j++,
For[i = 1, i <= ni, i++,
P[theta0[[i]], energy[[j]], t] = {Px[theta0[[i]], energy[[j]], t],
Py[theta0[[i]], energy[[j]], t],
Pz[theta0[[i]], energy[[j]], t]}]];
For[j = 2, j <= nj, j++,
For[i = 1, i <= ni, i++,
Pbar[theta0[[i]], energy[[j]],
t] = {Pbarx[theta0[[i]], energy[[j]], t],
Pbary[theta0[[i]], energy[[j]], t],
Pbarz[theta0[[i]], energy[[j]], t]}]];
Clear["i", "j"];
timing = AbsoluteTiming[
MPP = Flatten[
ParallelTable[D[P[theta0[[i]], energy[[j]], t], t] ==
Cross[Sum[(1 - theta[[ii]] theta[[i]]) (theta[[ii]] -
theta[[ii -
1]]) Sum[(Fj[energy[[jj]]] P[theta0[[ii]],
energy[[jj]], t] -
Fjbar[energy[[jj]]] Pbar[theta0[[ii]], energy[[jj]],
t]) (energy[[jj]] - energy[[jj - 1]]), {jj, 2,
nj}], {ii, 2, ni}], P[theta0[[i]], energy[[j]], t]], {i,
1, ni}, {j, 2, nj}]]];
timing[[1]]


I run the former code in a HPC which have 16 kernels, when ni=nj=6, the output is 301..., when ni=6, nj=8 the output is 1229...s, as you can find the time used arise quiet quickly. And what final I want is ni=80, nj=32, I can't accept the time used of this.

So, is there anyway to save time to do Table or Cross. Thanks for any suggestions.

• Honestly, this is not the way to write down a system of differential equations of this size. You try to press everything into symbolic computations. But these are slow. Better submit the system to NDSolve in the form D[X[t], t] == F[X[t]] with a function F[x_?NumericQ]:= that eats a numeric vector and throws out the forcing term. Then you can make use of fast operations like Dot for matrices and vectors. – Henrik Schumacher Dec 1 '18 at 10:09
• Otherwise it does actually not matter if you work on an "HPC" or not; you have to reformulate your problem in a way such that an HPC can play out it strengths -- and that is mere number crunching power. – Henrik Schumacher Dec 1 '18 at 10:10
• @HenrikSchumacher I'm not really understand what you mean, do you mean that I should not do D[P[t], t] but do D[Px[t]], D[Py[t]], D[Pz[t]] separately? – 袁子奕 Dec 1 '18 at 11:12
• I would go more into detail if I were able to understand what the right hand sides of the equations is supposed to do. For instance, I am puzzled what Cross is supposed to do there. Are the P supposed to be vector-valued? In total, there is really not enough information for providing a detailed answer. Certainly, your system is very expensive as each degree of freedom seems to couple nontrivially with all the other ones. – Henrik Schumacher Dec 1 '18 at 11:24
• If possible, describe your problem in natural language. – Αλέξανδρος Ζεγγ Dec 1 '18 at 11:26

Here's a 51 x 31 computation taking around a minute without parallelization. All I really did was precompute the cross product and substitute the components. Perhaps Cross[] does some symbolic operations that are avoided this way. The 6 x 6 calculation took around 0.03 seconds.

energy = Subdivide[0., 50., 50];
theta0 = Reverse[N[ArcSin[Subdivide[0., 1., 30]]]];
(* ... *)

timing = AbsoluteTiming[
MPP = Flatten[
Table[D[P[theta0[[i]], energy[[j]], t],
t] == (Cross[{a, b, c}, {d, e, f}] /.
Thread[{a, b, c, d, e, f} ->
Join[Sum[(1 - theta[[ii]] theta[[i]]) (theta[[ii]] -
theta[[ii -
1]]) Sum[(Fj[energy[[jj]]] P[theta0[[ii]],
energy[[jj]], t] -
Fjbar[energy[[jj]]] Pbar[theta0[[ii]], energy[[jj]],
t]) (energy[[jj]] - energy[[jj - 1]]), {jj, 2,
nj}], {ii, 2, ni}],
P[theta0[[i]], energy[[j]], t]]]), {i, 1, ni}, {j, 2,
nj}]]];
timing[[1]]

(*  67.5372  *)


Compile isn't always faster, especially when the code depends heavily on external variables. The following changes reduces the computation to 42.3379 seconds:

Fj = Compile[{{x, _Real}}, (betae (betae x)^2)/((E^(betae x) +
1) Ebare) + (betax (betax x)^2)/((E^(betax x) + 1) Ebarx) //
Evaluate];
Fjbar = Compile[{{x, _Real}}, ((betaee (betaee x)^2)/((E^(betaee x) +
1) Ebaree) + (betax (betax x)^2)/((E^(betax x) +
1) Ebarx)) // Evaluate];


If you don't compile the expressions at all, it takes 45.0974 sec., which shows how little Compile helps in this application.