6
$\begingroup$

I'm trying to use NDSolve to solve a set of ~400k equations. The problem I met is: when I ran the evaluation, Mathematica (version 11.0) raised an error:

NDSolve::ntdv: Cannot solve to find an explicit formula for the derivatives. Consider using the option Method->{"EquationSimplification"->"Residual"}.

However, after I added the option Method->{"EquationSimplification"->"Residual"} into the evaluation, it took almost 24h to get the anwsers.

So, what does the option Method->{"EquationSimplification"->"Residual"} really means? Are there other ways I can get the equations solved correctly and quickly, maybe in 1h?

Below are my simlified codes:

(*Equations*)
EqD[nC_,nN_] := c[nC,nN]'[t];
Eqk[nC_,nN_] := Sum[b[nC-i,nN-(1-i),i,1-i]*c[nC-i,nN-(1-i)][t]-(b[nC,nN,i,1-i]+a[nC,nN,i,1-i])*c[nC,nN][t]+a[nC+i,nN+(1-i),i,1-i]*c[nC+i,nN+(1-i)][t],{i,{0,1}}];
Eq[nC_,nN_] := EqD[nC,nN] == Eqk[nC,nN];
(*Coefficients*)
b[nC_,nN_,i_,j_] := 4\[Pi]*(R[nC,nN]+R[i,j])/\[CapitalOmega]*D\[Alpha][i,j]*c[i,j][t];
a[nC_,nN_,i_,j_] := 4\[Pi]*(R[nC-i,nN-j]+R[i,j])/\[CapitalOmega]*D\[Alpha][i,j]*Exp[-((F[nC-i,nN-j]+F[i,j]-F[nC,nN])/(k*T))]
(*Other definitions needed*)
R[nC_,nN_] := (3(nC+nN)*\[CapitalOmega]/(4\[Pi]))^(1/3)
D\[Alpha][i_,j_] := 0;D\[Alpha][1,0] := 8.0*10^-4;D\[Alpha][0,1] := 1.0*10^-4;
F[nC_,nN_] := -k*T*Log[(nC+nN)!/(nC!*nN!)]+(nC*Fc[nC+nN]+nN*Fn[nC+nN])/(nC+nN);
(*\[Pi],\[CapitalOmega],k and T are constants;Fc[nC_] and Fn[nN_] are two functions related to R[nC,0] and R[0,nN]*)

Are there any useful suggestions to speed up my evaluation? Thanks a lot!!!

$\endgroup$
  • $\begingroup$ @LammondW The code contains many typos. I recommend formulating the problem in the form of equations. $\endgroup$ – Alex Trounev Dec 19 '19 at 21:26
  • $\begingroup$ @AlexTrounev Sorry to mislead you, and I will re-edit the description to make it clear right now. $\endgroup$ – Lammond W Dec 20 '19 at 2:57
  • $\begingroup$ You code involves things like Cinit = 0.439673 % inside a Module, which is obviously incorrect. Please notice % doesn't represent percentage in Mathematica, it's the short form of Out. $\endgroup$ – xzczd Dec 20 '19 at 7:12
  • $\begingroup$ @xzczd Thanks for your reminding me. I've simplified the code to make it easy to understand for those who aren't familiar with material computing, so there may be some mistakes. Anyway, I've corrected the error. Thanks again! $\endgroup$ – Lammond W Dec 20 '19 at 7:36
  • 1
    $\begingroup$ How much time does {Sys,Bc} = PrepareEquations[tMin,tMax,NoC,NoN]; take? "So, what does the option Method->{"EquationSimplification"->"Residual"} really means?" We already have this post: mathematica.stackexchange.com/a/158519/1871 $\endgroup$ – xzczd Dec 20 '19 at 10:25
5
$\begingroup$

I tested the code as it is without optimization. I give a report and general recommendations: Test 1

In[11]:= Equations[tMin, tMax, 200, 20]; // AbsoluteTiming

Out[11]= {7.76804, Null}

Test 2

In[12]:= Equations[tMin, tMax, 200, 200]; // AbsoluteTiming

Out[12]= {86.8025, Null}

Test 3

In[13]:= Equations[tMin, tMax, 2000, 20]; // AbsoluteTiming

Out[13]= {109.039, Null}

Test 4

In[14]:= Equations[tMin, tMax, 2000, 200]; // AbsoluteTiming

Out[14]= {2174.46, Null}

And so it took 36 minutes, 30-50% of the CPU and 12GB of memory (maximum). Most of the time 1 kernel worked. I think that we can reduce the time to 3-5 minutes. Version report:

 $Version

Out[]= "12.0.0 for Microsoft Windows (64-bit) (April 6, 2019)"

Let's start the optimization. Here we are dealing with the Riccati matrix equation $$c'=c.A.c+B1.c+c.B2...(1) $$ where c is a 2000x200 matrix. Therefore we used Method -> {"EquationSimplification" -> "Residual"to solve it. Then it takes about 36 min on my ASUS laptop. We can reduce time by simple replacement t->tmax t. The corresponding code is

PrepareEquations[tMin_, tMax_, NoC_, NoN_] := 
  Module[{i, 
    j,(*C and N species*)\[Alpha]C = 1,(*C monomer*)\[Alpha]N = 
     1,(*N monomer*)nCmin = 0, nCmax = NoC,(*Number of C atoms*)
    nNmin = 0, nNmax = NoN,(*Number of N atoms*)
    c0 = 10^-30(*For numerical stability,avoid null value*)}, 
   Print["Function: PrepareEquations"];
   Clear[c, a, b, F, Fc, Fn];
   Cinit = 0.00439673;(*Initial C concentration*)
   Ninit = 0.00761477;(*Initial N concentration*)
   k = 1.380649*10^-23;(*Boltzmann constant,unit J/K*)
   T = 723;(*Temperature,unit K*)
   R[nC_, nN_] := ((3*(nC + nN)*0.0118)/(4 \[Pi]))^(1/3);(*Radius*)
   D\[Alpha][i_, j_] := 0;(*Diffusion coefficients*)
   D\[Alpha][\[Alpha]C, 0] := 8.0*10^-4;
   D\[Alpha][0, \[Alpha]N] := 1.0*10^-4;
   k\[Alpha][nC_, nN_, i_, j_] := 
    4 \[Pi]*(R[nC, nN] + R[i, j])/0.0118*D\[Alpha][i, j];
   (*Calculation of energy*)RsC = 0.02780299;
   \[Sigma]C = 4.09324*10^-19;
   Fc[nC_] := 4 \[Pi]*(R[nC, 0])^2*\[Sigma]C*(1 - RsC/R[nC, 0]);
   \[Sigma]N = 5.54492*10^-20;
   RsN = -0.14758424;
   Fn[nN_] := 4 \[Pi]*(R[0, nN])^2*\[Sigma]N*(1 - RsN/R[0, nN]);
   F[nC_, 
     nN_] := -k*T*
      Log[(nC + nN)!/(nC!*nN!)] + (nC*Fc[nC + nN] + 
        nN*Fn[nC + nN])/(nC + nN);
   b[nC_, nN_, i_, j_] := 
    k\[Alpha][nC, nN, i, j]*c[i, j][t];(*Condensation rate*)
   a[nC_, nN_, i_, j_] := 
    k\[Alpha][nC - i, nN - j, i, j]*
     Exp[-((F[nC - i, nN - j] + F[i, j] - F[nC, nN])/(k*
           T))];(*Emission rate*)a[i_, j_, i_, j_] := 0;(*Monomer*)
   a[0, nN_, 1, j_] := 0;
   a[nC_, 0, i_, 1] := 0;
   Clear[CL, Eqk, EqD, Eq];
   With[{},(*Initial conditions*)CL[nC_, nN_] := c[nC, nN][tMin] == 0;
    (*Definition of the equilibrium equations*)
    EqD[nC_, nN_] := c[nC, nN]'[t];
    Eqk[nC_, nN_] := 
     Sum[b[nC - i, nN - (1 - i), i, 1 - i]*
        c[nC - i, nN - (1 - i)][
         t] - (b[nC, nN, i, 1 - i] + a[nC, nN, i, 1 - i])*
        c[nC, nN][t] + 
       a[nC + i, nN + (1 - i), i, 1 - i]*
        c[nC + i, nN + (1 - i)][t], {i, {0, 1}}];
    Eq[nC_, nN_] := EqD[nC, nN] == tMax Eqk[nC, nN];
    (*Definition of equilibrium concentrations for each special size*)
    Module[{nC, nN}, nC = nCmax;
     For[nN = 1, nN < nNmax, nN++, 
      Eqk[nC, nN] = (b[nC, nN - 1, 0, 1]*
           c[nC, nN - 1][t] - (b[nC, nN, 0, 1] + a[nC, nN, 0, 1])*
           c[nC, nN][t] + 
          a[nC, nN + 1, 0, 1]*
           c[nC, nN + 1][t]) + (b[nC - 1, nN, 1, 0]*c[nC - 1, nN][t] -
           a[nC, nN, 1, 0]*c[nC, nN][t])];
     nN = nNmax;
     For[nC = 1, nC < nCmax, nC++, 
      Eqk[nC, nN] = (b[nC - 1, nN, 1, 0]*
           c[nC - 1, nN][t] - (b[nC, nN, 1, 0] + a[nC, nN, 1, 0])*
           c[nC, nN][t] + 
          a[nC + 1, nN, 1, 0]*
           c[nC + 1, nN][t]) + (b[nC, nN - 1, 0, 1]*c[nC, nN - 1][t] -
           a[nC, nN, 0, 1]*c[nC, nN][t])];
     nC = 0;
     For[nN = 2, nN < nNmax, nN++, 
      Eqk[nC, nN] = (b[nC, nN - 1, 0, 1]*
           c[nC, nN - 1][t] - (b[nC, nN, 0, 1] + a[nC, nN, 0, 1])*
           c[nC, nN][t] + 
          a[nC, nN + 1, 0, 1]*
           c[nC, nN + 1][t]) + (-b[nC, nN, 1, 0] c[nC, nN][t] + 
          a[nC + 1, nN, 1, 0]*c[nC + 1, nN][t])];
     nN = 0;
     For[nC = 2, nC < nCmax, nC++, 
      Eqk[nC, nN] = (b[nC - 1, nN, 1, 0]*
           c[nC - 1, nN][t] - (b[nC, nN, 1, 0] + a[nC, nN, 1, 0])*
           c[nC, nN][t] + 
          a[nC + 1, nN, 1, 0]*
           c[nC + 1, nN][t]) + (-b[nC, nN, 0, 1] c[nC, nN][t] + 
          a[nC, nN + 1, 0, 1]*c[nC, nN + 1][t])];
     nC = nCmax; nN = 0;
     Eqk[nC, 
       nN] = (-b[nC, nN, 0, 1]*c[nC, nN][t] + 
         a[nC, nN + 1, 0, 1]*c[nC, nN + 1][t]) + (b[nC - 1, nN, 1, 0]*
          c[nC - 1, nN][t] - a[nC, nN, 1, 0]*c[nC, nN][t]);
     nC = 0; nN = nNmax;
     Eqk[nC, 
       nN] = (-b[nC, nN, 1, 0]*c[nC, nN][t] + 
         a[nC + 1, nN, 1, 0]*c[nC + 1, nN][t]) + (b[nC, nN - 1, 0, 1]*
          c[nC, nN - 1][t] - a[nC, nN, 0, 1]*c[nC, nN][t]);
     nC = nCmax; nN = nNmax;
     Eqk[nC, 
       nN] = (b[nC, nN - 1, 0, 1]*c[nC, nN - 1][t] - 
         a[nC, nN, 0, 1]*c[nC, nN][t]) + (b[nC - 1, nN, 1, 0]*
          c[nC - 1, nN][t] - a[nC, nN, 1, 0]*c[nC, nN][t]);
     nC = 0; nN = 0;
     CL[nC, nN] = c[nC, nN][tMin] == 0;
     Eqk[nC, nN] = 0;
     nC = \[Alpha]C; nN = 0;
     CL[nC, nN] = c[nC, nN][tMin] == Cinit + c0;
     Eqk[nC, 
       nN] = (-b[nC, nN, 1, 0]*c[nC, nN][t] + 
         a[nC + 1, nN, 1, 0]*c[nC + 1, nN][t]) + (-b[nC, nN, 0, 1]*
          c[nC, nN][t] + 
         a[nC, nN + 1, 0, 1]*
          c[nC, nN + 1][
           t]) - (Sum[(b[x, y, nC, nN] - a[x, y, nC, nN])*
           c[x, y][t], {x, 0, nCmax - 1}, {y, 0, nNmax}] + 
         Sum[(-a[nCmax, i, nC, nN]*c[nCmax, i][t]), {i, 0, nNmax}]);
     nC = 0; nN = \[Alpha]N;
     CL[nC, nN] = c[nC, nN][tMin] == Ninit + c0;
     Eqk[nC, 
       nN] = (-b[nC, nN, 1, 0]*c[nC, nN][t] + 
         a[nC + 1, nN, 1, 0]*c[nC + 1, nN][t]) + (-b[nC, nN, 0, 1]*
          c[nC, nN][t] + 
         a[nC, nN + 1, 0, 1]*
          c[nC, nN + 1][
           t]) - (Sum[(b[x, y, nC, nN] - a[x, y, nC, nN])*
           c[x, y][t], {x, 0, nCmax}, {y, 0, nNmax - 1}] + 
         Sum[(-a[i, nNmax, nC, nN]*c[i, nNmax][t]), {i, 0, nCmax}]);];
    Bc := 
     Flatten[ParallelTable[
        c[nC, nN], {nC, 0, nCmax}, {nN, 0, nNmax}]] // Evaluate];
   DistributeDefinitions[Eq, CL];
   Sys = Join[
       Flatten[ParallelTable[
         Eq[nC, nN], {nC, 0, nCmax}, {nN, 0, nNmax}]], 
       Flatten[ParallelTable[
         CL[nC, nN], {nC, 0, nCmax}, {nN, 0, nNmax}]]] // Flatten // 
     Evaluate;
   {Sys, Bc}];
tMin = 0; tMax = 7.75*10^5;

Here we have replaced only one equation Eq[nC_, nN_] := EqD[nC, nN] == tMax Eqk[nC, nN];.Now we call

PrepareEquations[tMin, tMax, 2000, 200]; // AbsoluteTiming

And it takes {536.24, Null}. Now we solve the system on {t, 0, 1}:

RSol = NDSolve[Sys, Bc, {t, tMin, 1}, 
    Method -> {"EquationSimplification" -> 
       "Residual"}]; // AbsoluteTiming

It takes only {1148.03, Null}. Thus, we saved 9 minutes with a simple replacement t->tmax t. However, if we can write the equation in matrix form, then we can use exponential Rosenbrock-type integrators. Then we can reduce the time to 1 minute.

| improve this answer | |
$\endgroup$
  • $\begingroup$ Additionally, I use all the 24 cores hardly. Does it matter? Of course I'll check it too. Thank you very very much! $\endgroup$ – Lammond W Dec 20 '19 at 15:27
  • $\begingroup$ @LammondW Version 12 is much faster on ODE systems than version 11. See mathematica.stackexchange.com/questions/208590/… $\endgroup$ – Alex Trounev Dec 20 '19 at 16:30
  • $\begingroup$ Thanks for your patience and kindness and detailed anwser. I'll figure it out and then reply to you^_^ $\endgroup$ – Lammond W Dec 21 '19 at 1:46
  • $\begingroup$ Thanks for your anwser again!! Someone else also suggested I should use a vector form of equations. Honestly I'm not very well with Mathematica, so I'm wondering whether the method applies to the situation where I will change the number of equations (namely nC and nN)? $\endgroup$ – Lammond W Dec 24 '19 at 5:26
  • $\begingroup$ I've checked your optimization and here's another problem. The solution gained using your method is a series of InterpolatingFunction with a domain of {0.,1.}, while I need them in the domain of {0.,tMax}. Any suggestions? Thanks a lot ^_^ $\endgroup$ – Lammond W Dec 24 '19 at 7:28

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.