Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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

Hop[t_, i_] =
  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] 

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_] =
    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]]];
Eqsopwo[t_]=Join[Thread /@ Thread[D[Rho[t], t] == Fun[t]],Eq0a]//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;
share|improve this question
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 8:47
Really? that means higher than double precision of workingprecision will be impossible – LCFactorization Jan 25 at 8:41
up vote 7 down vote accepted

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];
    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:


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

share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.