# Parallel computation of a system of differential equation

I need to solve a system of differential equations which contains, at minimum, ten

thousands of differential equations. Actually, I am trying to solve density matrix

equation. It's like something $\frac{\partial}{\partial t}\rho(t) = \frac{-I}{\hbar}(H(t)\rho(t) -\rho(t)H(t))$.

If the dimension of Hilbert space is $N$, then there is $N^2$ coupled differential equations.

I want to know if it is possible to solve this system of differential equations in parallel

in Mathematica? My cpu has six cores, so I guess I can get a six times speed up, if I solve

the equations in parallel.

-
If you have to deal with such a high number of equations, I'd suggest Fortran or C. There was a similar question here: scicomp.stackexchange.com/questions/8855/… –  Gregory Rut Apr 29 at 11:52
From the form of your equation I do not see how the system should be coupled, ist seems like the solution is a simple integration of any matrix element. To the topic: In general mathematica is not the tool to go when it comes to computation speed and memory efficiency. –  Wizard Apr 29 at 13:46
@ Wizrd. I edited the equation. What program do you suggest for such a problem? –  yashar Apr 29 at 14:09

I agree with Wizard's comment that Mathematica isn't the best tool when computation speed and memory efficiency are very important and that you will very likely end up to use specialized code for such a large system.

On the other hand I think NDSolve used as a pure solver is quite impressive and might potentially even be able to solve such a large system on a large enough machine (=having enough memory). NDSolve certainly might be a good choice to learn which methods and features are important for your class of problems if you manage to produce scaled versions of your problem which show the expected behavior of the full system. Such experiments could be a good guide when looking for the right specialized code for the "real" problem if you learn that NDSolve really can't handle it.

The trouble with large systems of ODEs in Mathematica is usually not really the solver (NDSolve) itself but building the system of equations efficiently and handling the large output. Both can be overcome: in your case the generation of the system of ODEs in matrix form seems relatively straightforward, most probably you'll need to look into packed arrays and/or sparse arrays to create an efficient version of H. To avoid NDSolve to generate an interpolating function for each entry of the result vector there are several possibilities and what to do depends a lot on what you really need. You can e.g. only return the final values of these functions or use NDSolveProcessEquations and friends to run the solver without creating the full output and only extract what you need. You can find examples for all these techniques in the exhaustive documentation for NDSolve and also on this site.

As for parallel execution: many methods of NDSolve internally incorporate the solving of linear systems of equations (and potentially other similar algorithms) which is to some extent automatically parallelized, other than that there are AFAIK no possibilities of making use of your 6 processors with NDSolve. The best you probably could do is to run several NDSolves in parallel in case you will need solutions for varying H` or initial/boundary conditions. If you look for another tool you will find that solvers with parallelization at various levels exist, but parallelism for differential equation solvers isn't a feature that can be taken for granted. Whatever you'll find, you should also note that a speedup of 6 with 6 processors is only the theoretical upper limit, it is usually very hard to even come close to that for real applications...

-