I want to solve the time-dependent Schrödinger equation:

$$ i\partial_t \psi(t) = H(t)\psi(t)$$

for matrix, time-dependent $H(t)$ and vector $\psi$.

What is an efficient way of doing this so that it efficiently scales to high-dimensional spaces?

  • $\begingroup$ what are the valuse of b and omega for the second plot? I want to check my code to see is working or not! Thanks Jiyan. $\endgroup$
    – user21640
    Commented Oct 23, 2014 at 17:26

3 Answers 3


Time-dependent case

in the time-dependent case, $[H(t),H(t')]\neq0$ in general and we need to time-order, ie, the operator taking a state from $t=0$ to $t=\tau$ is $U(0,\tau)=\mathcal{T}\exp(-i\int_0^\tau dt\, H(t))$ with $\mathcal{T}$ the time-ordering operator. In practice we just split the time interval into lots of small pieces (basically using the Baker-Campbell-Hausdorff thing).

So, consider the time-dependent Hamiltonian for a two-level system:

$$ H = \left( \begin{array}{cc} \epsilon_1 && b \cos(\omega t) \\ b\cos(\omega t) && \epsilon_2 \\ \end{array} \right) $$

i.e. two level coupled by a time-periodic driving (see here). Even this simplest possible periodically-driven system can't be solved analytically in general.

Anyway, here's a function to construct the hamiltonian:

ham[e1_, e2_, b_, omega_, 
  t_] := {{e1, b*Cos[omega*t]}, {b*Cos[omega*t], e2}}

and here's one to construct the propagator from some initial time to some final time, given a function to construct the Hamiltonian matrix at each point in time (and splitting the interval into $n$ slices--you should try with increasing $n$ until your results stop changing):

constructU::usage = "constructU[h,tinit,tfinal,n]";
constructU[h_, tinit_, tfinal_, n_] := 
 Module[{dt = N[(tfinal - tinit)/n], 
   curVal = IdentityMatrix[Length@h[0]]}, 
  Do[curVal = MatrixExp[-I*h[t]*dt].curVal, {t, tinit, tfinal - dt, 

This constructs the operator $U(0,\tau)=\mathcal{T}\exp(-i\int_0^\tau dt\,H(t))$ as $$ U(0,\tau)\approx\prod_{n=0}^{N}\exp\left( -iH(ndt)dt \right) $$ with $N=\tau/dt-1$ (or its ceiling anyway). This is an approximation to the correct $U$.

And now here is how to look at the time-dependent expectation of $\sigma_z$ for different coupling strengths $b$:

ClearAll[cU, psi0];
psi0 = {1., 0};
   Chop[#\[Conjugate].PauliMatrix[3].#] &@(constructU[
       ham[-1., 1., b, 1., #] &, 0, upt, 100].psi0),
   {upt, .01, 20, .1}
  Joined -> True,
  PlotRange -> {-1, 1}
 {b, 0, 2}

Mathematica graphics

Alternatively, you could calculate the wavefunction at some time tfinal given the wavefunction at time tinit with this:

propPsi[h_, psi0_, tinit_, tfinal_, n_] := 
 Module[{dt = N[(tfinal - tinit)/n],
   psi = psi0},
   psi = MatrixExp[-I*h[t]*dt, psi], {t, tinit, tfinal - dt, dt}

which uses the form MatrixExp[-I*h*t,v]. For large sparse matrices (eg, for h a many-body Hamiltonian), this can be much faster at the cost of losing access to $U$.

  • $\begingroup$ Thanks a lot for all this. However as you mentioned in your previous comment, my problem is a time-dependent schrodinger equation. In this case, the Hamiltonian doesn't commute in different times, and it can't just be a simple exponential; it should be a path ordered exponential. This is the reason that I can't do this. Such this problems don't have an analytical solution, they ave to be solved numerically!! $\endgroup$ Commented Mar 21, 2013 at 13:49
  • $\begingroup$ @Zahra For a time-dep H, you can simply construct the propagator from $t=0$ to some time $t=\tau$, say. I can explain how if you want. but let me know if you actually want so I don't waste my time if you insist on doing it with NDSolve--but ask yourself which is the most practical way if you have a hilbert space of dimension 20000, for instance (so you'd need to solve 20000 coupled ODEs with your approach) $\endgroup$
    – acl
    Commented Mar 21, 2013 at 15:35
  • $\begingroup$ Thanks a lot for the time that you give to ask my question. I really appreciate that. I'm not insisting to solve my problem with NDSolve. It would be great if I can solve it the way you are explaining. I just thought it is not possible to solve it non-numerically. Can you please tell me how to do that? I appreciate it. $\endgroup$ Commented Mar 21, 2013 at 15:55
  • 1
    $\begingroup$ @Zahra oh I see, no, what I suggest is fully numerical. You're absolutely right that such problems cannot be solved analytically in general. OK let me write it up quickly and you can see if it's useful (I routinely use this on systems with much bigger Hilbert spaces than yours, up to 20-30000) $\endgroup$
    – acl
    Commented Mar 21, 2013 at 15:58
  • 1
    $\begingroup$ here you go (I had actually done this the first time, just posted only the time-independent limit because I did not realize you had a time-dependent hamiltonian). Note that the way I construct the Manipulate is not efficient because I recalculate $U$ from scratch all the time, but it's fast enough... $\endgroup$
    – acl
    Commented Mar 21, 2013 at 16:04

Since there hasn't been any discussion of NDSOlve yet, let me point out that for a finite-dimensional Hilbert space where the Schrödinger equation is merely a first-order equation in time, it's easiest to just do this (using the two-dimensional Hamiltonian ham from acl's answer):

ham[e1_, e2_, b_, omega_, 
   t_] := {{e1, b*Cos[omega*t]}, {b*Cos[omega*t], e2}};
Manipulate[Module[{ψ, sol, tMax = 20},
  sol = First@NDSolve[{I D[ψ[t], t] == 
       ham[-1, 1, b, 1, t] .ψ[t], ψ[0] == {1,0}}, ψ, {t, 0, tMax}];
  Plot[Chop[#\[Conjugate].PauliMatrix[3].#] &@(ψ /. sol)[t], 
     {t, 0, tMax}, PlotRange -> {-1, 1}]
 {{b, 1}, 0, 2}


I copied the parameters from acl's answer too, to show the direct comparison in the Manipulate. Here the vector $\psi$ is recognized by NDSolve as two-dimensional, so the formulation of the problem is quite concise, and we can leave the time step choice up to Mathematica instead of choosing a discretization ourselves.

  • $\begingroup$ In fact the original question explicitly mentioned NDSolve (I edited it to make it less localized). There's nothing wrong with NDSolve for up to a few thousand states, but the approach I gave scales much better (I use it for systems with dimensions in the tens of thousands; NDSolve seems to tank much earlier); of course the way I wrote the code it's inefficient. $\endgroup$
    – acl
    Commented Mar 22, 2013 at 9:33
  • $\begingroup$ as an additional comment, this approach (using NDSolve directly) works also for cases where the "Hamiltonian" depends on the wavefunction, so that we have a set of nonlinear coupled ODEs. This kind of problem appears in various mean-field approaches to many-body systems (eg, Gutzwiller-ansatz approach to many-body dynamics of bosons, see eg eq 3 here ). I've used NDSolve for precisely this problem with up to a couple of thousand coupled ODEs; it's really not practical at those sizes, but there's no alternative (in mma) for nonlinear ODEs. $\endgroup$
    – acl
    Commented Mar 22, 2013 at 14:10
  • $\begingroup$ @acl Thanks for pointing that out (I already upvoted your answer). I think it would have been better to edit the question in such a way as to retain some information on what the OP originally tried already. $\endgroup$
    – Jens
    Commented Mar 22, 2013 at 14:19
  • $\begingroup$ feel free to change it, I would not object (and I imagine neither would Zahra). I thought the question as it was was way too localized (eg, it was asked about a specific Hamiltonian, defined only in hard-to read code--take a look at the original form if you haven't) and wanted to make it as general as possible so it's useful. I think that phrasing it the current way it admits as many approaches as possible. $\endgroup$
    – acl
    Commented Mar 22, 2013 at 14:47
  • $\begingroup$ @acl I see your point. No worries, it's fine the way it is. $\endgroup$
    – Jens
    Commented Mar 22, 2013 at 15:16

Frame it as a set of Linear ODEs and solve it somehow. I usually use Implicit Runge Kutta in the interaction picture.

solver[H_, a_] := 
soln = Module[{d, init, eq, vars, solargs, t, t0, tf}, 
d = Dimensions[H][[1]];
t0 = a[[2]];
tf = a[[3]];
t = a[[1]];
u[t_] := Table[Subscript[u, i, j][t], {i, 1, d}, {j, 1, d}];
init = (u[t0] == IdentityMatrix[d]);
eq = (I u'[t] == H.u[t]);
vars = Flatten[Table[Subscript[u, i, j], {i, 1, d}, {j, 1, d}]];
solargs = LogicalExpand[eq && init];
 NDSolve[solargs, vars, a, 
  Method -> {"FixedStep", 
    Method -> {"ImplicitRungeKutta", "DifferenceOrder" -> 10}}, 
  StartingStepSize -> tf/100, MaxSteps -> Infinity]]];

U[t_] := u[t] /. soln[[1]]

Alternatively, you could solve for the $\psi(t)$ and obtain U as $|\psi(t)\rangle\langle \psi(0)| $ and use an appropriate normalisation to preserve probability.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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