# How to speed up NDSolve when there are ~400k equations?

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!!!

• @LammondW The code contains many typos. I recommend formulating the problem in the form of equations. Dec 19, 2019 at 21:26
• @AlexTrounev Sorry to mislead you, and I will re-edit the description to make it clear right now. Dec 20, 2019 at 2:57
• 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. Dec 20, 2019 at 7:12
• @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! Dec 20, 2019 at 7:36
• 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 Dec 20, 2019 at 10:25

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.

• Additionally, I use all the 24 cores hardly. Does it matter? Of course I'll check it too. Thank you very very much! Dec 20, 2019 at 15:27
• @LammondW Version 12 is much faster on ODE systems than version 11. See mathematica.stackexchange.com/questions/208590/… Dec 20, 2019 at 16:30
• Thanks for your patience and kindness and detailed anwser. I'll figure it out and then reply to you^_^ Dec 21, 2019 at 1:46
• 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)? Dec 24, 2019 at 5:26
• 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 ^_^ Dec 24, 2019 at 7:28