How to use Compile with NDSolve function for large no of differential equations?

Question

Matrix M[t] and and Hamiltonian of same order L^n x L^n.

L=3; n=3;tmax = 2;
Rho[t_] = Table[M[i,j][t], {i,L^n}, {j,L^n}]; //Timing

Rho0=ParallelTable[M[i,j][0],{i,L^n},{j,L^n}];
Hop[t_, i_] =
Transpose[
  Table[If[i == j, 1, 1+Cos[t]], {i, {1}}, {j, L}]].Table[If[i == j, 1, 1+Sin[t]^2], {i, {2}}, {j, L}];

I1[i_Integer] := I1[i]=IdentityMatrix[L^(i-1)];
I2[i_Integer] := I2[i]=IdentityMatrix[L^(n-i)];

H[t_] = Sum[FixedPoint[ArrayFlatten, I1[i] \[TensorProduct] Hop[t,i] 
\[TensorProduct]I2[i]],{i,n}];

I have made a code for large no. of time dependent differential equations with initial values. I want to compile this for using NDSolve. Here is my code...

Rho1[t_] =
Table[
    D[Rho[t],t][[i,j]] == H[t][[i,j]],
    {i,L^n}, {j,L^n}]; // Timing

Fun[t_]= -I (#1.#2 - #2.#1)&[H[t],Rho[t]];
Eq0[i_,j_]:=If[i == 1 && j == 1,#==1.,#==0. ]&[Rho0[[i,j]]];
Eq0a=ParallelTable[Eq0[i,j],{i,L^n},{j,L^n}];
Eqsopwo[t_]=Join[Thread /@ Thread[D[Rho[t], t] == Fun[t]],Eq0a]//Flatten;
\[Rho]Termsop[t_]=Rho[t]//Flatten;

sol= NDSolve[Eqsopwo[t],Diagonal[Rho[t]],{t,0,tmax},DependentVariables->\[Rho]Termsop[t],MaxSteps->10^5,MaxStepFraction-> 1/50, Method-> {"EquationSimplification"->"Solve"}]//Chop;
Plot[Evaluate[{Re[Diagonal[Rho[t]]]}/.sol],{t,0,tmax},PlotRange->All]

Related http://mathematica.stackexchange.com/questions/1096/list-of-compilable-functions — Dr. belisarius, Jan 28 '14 at 19:47
I remembered once I read some material (probably in somewhere inside the document) that said NDSolve will use Compile automatically… — xzczd, Jan 29 '14 at 08:47
Really? that means higher than double precision of workingprecision will be impossible — LCFactorization, Jan 25 '16 at 08:41

score 12 · Accepted Answer · edited Apr 13 '17 at 12:56

I won't use your specific set of equations, but will answer more generally. Suppose we have the following set of first order equations with forcing:

$$ \frac{d\vec{x}}{dt}=A(t,\vec{x}(t))\vec{x}(t)+f(t) $$

Suppose further that we have $A$ and $f$ already as Mathematica functions, for example, here is something random:

d = 2;
A[t_, x_] := {{-1, -1 Exp[-t]}, {1 + Cos[x[[2]]], -1}};
y[t_] := {Sin[t], 0};

Now we wish to use compiled versions of these functions in NDSolve, and we won't anywhere make an assumption that d=2. For each row of our equation vector we will need a separate compilation. We have to jump through some hoops because Compile has the HoldAll attribute, and we can also not let the expression x[[n]] evaluate since x is just a symbol with no parts.

expr = A[t, Table[x[n], {n, d}]].Table[x[n], {n, d}] + y[t];
Table[
    comp[n] = ReleaseHold[
         Hold[Compile[{{t, _Real}, {x, _Real, 1}}, inside]] //. 
         {inside -> expr[[n]], x[m_] :> x[[m]]}
    ];,
    {n, d}
];

At this point is is a good idea to make sure that the compilation actually worked, or if it is just calling the main evaluator in some trivial way. See an answer here to learn about this.

Now Compile will be angry with us if we try to use comp directly as the RHS in NDSolve because we would be giving it symbolic inputs. As a work around, we define rhs which only evaluates when time is numeric:

rhs[n_, t_?NumericQ, x_]:=comp[n][t, x];

Now we can define our set of equations and initial conditions to be put into NDSolve:

eqns=Join[
    Table[
        x[n]'[t]==rhs[n,t,Table[x[m][t],{m,d}]],
        {n,d}
    ],
    Table[x[n][-1]==0,{n,d}]
];

And finally solve them:

fcns[t_] = Table[x[n][t], {n, d}] /. NDSolve[eqns, Table[x[n], {n, d}], {t, 0, 10}];
Plot[fcns[t], {t, 0, 10}]

Final Result Plot

Very under-rated answer (from 6.5 years ago)! – Anton Antonov Dec 12 '20 at 02:19 — Anton Antonov, Dec 12 '20 at 02:19

How to use Compile with NDSolve function for large no of differential equations?

1 Answers1

Linked