Is there any way to save time using of "Cross" or "Table"?

Question

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 01 '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 01 '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 01 '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 01 '18 at 11:24
@HenrikSchumacher Really sorry about that, I miss the definition of P and Pbar. Now you could just copy my code and it can run to get the time used. I follow your advice do D[P[t]] separately, it almost finished right away, I'm checking my code of the separately definition, and try to make it more clearly. As you can see if it write separately, it will be too long. — 袁子奕, Dec 01 '18 at 11:53
Marginally related: Why should I avoid the For loop in Mathematica? and Avoiding procedural loops — Michael E2, Dec 01 '18 at 13:56
How is this related to https://mathematica.stackexchange.com/questions/185067/nlnum-problem-of-ndsolve?noredirect=1#comment482667_185067 ? — Alex Trounev, Dec 01 '18 at 14:05
@AlexTrounev It's not entirely unrelated, the Cross in that question won't take such a long time. In other words, I didn't meet this problem in that question. But both question are proposed to solve the same physics problem. I ask this new question is because I did't find any related question mentioned the time used by Cross in this forum. — 袁子奕, Dec 01 '18 at 14:33

Michael E2 · Accepted Answer · 2018-12-01T16:23:58.610

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.

Is there any way to save time using of "Cross" or "Table"?

1 Answers1