1
$\begingroup$

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.

$\endgroup$
  • 3
    $\begingroup$ 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. $\endgroup$ – Henrik Schumacher Dec 1 '18 at 10:09
  • 1
    $\begingroup$ 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. $\endgroup$ – Henrik Schumacher Dec 1 '18 at 10:10
  • $\begingroup$ @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? $\endgroup$ – 袁子奕 Dec 1 '18 at 11:12
  • $\begingroup$ 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. $\endgroup$ – Henrik Schumacher Dec 1 '18 at 11:24
  • $\begingroup$ If possible, describe your problem in natural language. $\endgroup$ – Αλέξανδρος Ζεγγ Dec 1 '18 at 11:26
3
$\begingroup$

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.

$\endgroup$

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.