# 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...

-
How does one obtain only the final value? @Kuba asked a question regarding this topic (mathematica.stackexchange.com/questions/31094/…) and the provided answer wasn't exactly a solution. –  Gregory Rut Sep 12 at 10:18
@GregoryRut: Why do you think the answer given there was not a solution? I think it basically does what was asked for and also what you are intersted in: run NDSolve and only store the solutions at some points, you could of course also only store at one point, the end point. Neverthless you might find the additional answer that I just posted for that question also of interest, it gives a solution which is somewhat simpler if you really only are interested in the final value. –  Albert Retey Sep 12 at 22:46
Your solution is simply great. I didn't expect that it be done in such a simple way. My point was that NDSolveProcessSolutions actually interpolates the solutions so actually this is exactly the same as if one created a loop. From my understanding, that method might be also faulty in such a sense that such division may lead to increasing error (since you lose the information on the previous points). –  Gregory Rut Sep 13 at 10:17
@GregoryRut: I see, that is of course a difference that I didn't have in mind: if you reinitialize you'll loose previous points, if you don't ProcessSolution will do the interpolation. I can't remember exactly but I think in older versions there was no way to get avoid the interpolation but in newer version it might be possible to extract results from StateData without building the full interpolation (or maybe it was just planned?) and continue without reinitialization. Anyway, fortunately there are also these other possibilities that I did mention... –  Albert Retey Sep 13 at 10:32