# Solve numerical differential equations at t->Infinity

I have a first order non-homogeneous system of differential equation (100+ equations, so no hope to solve them analytically, due to the Abel–Ruffini theorem). If I solve them using NDSolve, I have to use the form:

NDSove[Join[MyEquations,{y0==y00,y1==y01,...,yN==yN0}],{y0[t],y1[t],...,yN[t]},
{t,timeStart,timeEnd}]


Now this will give the solution as an interpolation function between timeStart and timeEnd. This is, however, not really what I need. I need the solution of this system at t->Infinity only (in physics it's called the steady state solution).

The question is: How can I get the solution of y0[Infinity] numerically?

The problem with doing it the "easy" way, i.e., by choosing a high value of t, is that it consumes so much memory that the kernel crashes.

Update:

Since only help can be provided with an example, I created a simplified system showing the problem. I hate going into details because it's the wrong way to ask a question here, but there seems to be no way around it. Now this system is an atomic system with 3 levels, leading to 9 density matrix equations (3 populations and 6 coherences). This system simulates the famous problem called EIT (Electromagnetically Induced Transparency, and the steady state solution shows the red curve from the Wikipedia page). I replaced all the parameters, and all that's left is the parameter $\Delta$, which represents the detuning (in GHz). The task is: Get the steady state solution of this system for about a 500 values of $\Delta$ to see that red curve. This is equivalent to a scan of light frequency in an experiment. Now this is doable for 3 levels with a simple Table[NDSolve[...],{Δ,...,...}].

Here's the system the full NDSolve call (please just copy/paste to your Mathematica notebook):

MyEquations = {(I Subscript[ρ, {0, 0}]'[t])/(2 π) == -0.5 E^(-2 I π t Δ) Subscript[ρ, {-1, 0}][t] + 0.5 E^(2 I π t Δ) Subscript[ρ, {0, -1}][t] - (0. + 1.273255460229472 I) Subscript[ρ, {0, 0}][t] + 0.5 E^(2 I π t Δ) Subscript[ρ, {0, 1}][t] - 0.5 E^(-2 I π t Δ) Subscript[ρ, {1, 0}][t], (I Subscript[ρ, {1, 0}]'[t])/(2 π) == -0.5 E^(2 I π t Δ) Subscript[ρ, {0, 0}][t] + 0.5 E^(2 I π t Δ) Subscript[ρ, {1, -1}][t] - (0. + 0.6366356878618905 I) Subscript[ρ, {1, 0}][t] + Δ Subscript[ρ, {1, 0}][t] + 0.5 E^(2 I π t Δ) Subscript[ρ, {1, 1}][t], (I Subscript[ρ, {-1, 0}]'[t])/(2 π) == 0.5 E^(2 I π t Δ) Subscript[ρ, {-1, -1}][t] - (0. + 0.6366356878618905 I) Subscript[ρ, {-1, 0}][t] + Δ Subscript[ρ, {-1, 0}][t] + 0.5 E^(2 I π t Δ) Subscript[ρ, {-1, 1}][t] - 0.5 E^(2 I π t Δ) Subscript[ρ, {0, 0}][t], (I Subscript[ρ, {0, 1}]'[t])/(2 π) == -0.5 E^(-2 I π t Δ) Subscript[ρ, {-1, 1}][t] + 0.5 E^(-2 I π t Δ) Subscript[ρ, {0, 0}][t] - (0. + 0.6366356878618905 I) Subscript[ρ, {0, 1}][t] - Δ Subscript[ρ, {0, 1}][t] - 0.5 E^(-2 I π t Δ) Subscript[ρ, {1, 1}][t], (I Subscript[ρ, {1, 1}]'[t])/(2 π) == (I (0.00005 + 4 Subscript[ρ, {0, 0}][t]))/(2 π) - 0.5 E^(2 I π t Δ) Subscript[ρ, {0, 1}][t] + 0.5 E^(-2 I π t Δ) Subscript[ρ, {1, 0}][t] - (0. + 0.000015915494309189534 I) Subscript[ρ, {1, 1}][t], (I Subscript[ρ, {-1, 1}]'[t])/(2 π) == 0.5 E^(-2 I π t Δ) Subscript[ρ, {-1, 0}][t] - (0. + 0.000015915494309189534 I) Subscript[ρ, {-1, 1}][t] - 0.5 E^(2 I π t Δ) Subscript[ρ, {0, 1}][t], (I Subscript[ρ, {0, -1}]'[t])/(2 π) == -0.5 E^(-2 I π t Δ) Subscript[ρ, {-1, -1}][t] - (0. + 0.6366356878618905 I) Subscript[ρ, {0, -1}][t] - Δ Subscript[ρ, {0, -1}][t] + 0.5 E^(-2 I π t Δ) Subscript[ρ, {0, 0}][t] - 0.5 E^(-2 I π t Δ) Subscript[ρ, {1, -1}][t], (I Subscript[ρ, {1, -1}]'[t])/(2 π) == -0.5 E^(2 I π t Δ) Subscript[ρ, {0, -1}][t] - (0. + 0.000015915494309189534 I) Subscript[ρ, {1, -1}][t] + 0.5 E^(-2 I π t Δ) Subscript[ρ, {1, 0}][t], (I Subscript[ρ, {-1, -1}]'[t])/(2 π) == (0. - 0.000015915494309189534 I) Subscript[ρ, {-1, -1}][t] + 0.5 E^(-2 I π t Δ) Subscript[ρ, {-1, 0}][t] - 0.5 E^(2 I π t Δ) Subscript[ρ, {0, -1}][t] + (I (0.00005 + 4 Subscript[ρ, {0, 0}][t]))/(2 π)} /. Δ -> 3;
MyBC = {Subscript[ρ, {0, 0}] == 1, Subscript[ρ, {1, 0}] == 0,
Subscript[ρ, {-1, 0}] == 0, Subscript[ρ, {0, 1}] == 0,
Subscript[ρ, {1, 1}] == 0, Subscript[ρ, {-1, 1}] == 0,
Subscript[ρ, {0, -1}] == 0, Subscript[ρ, {1, -1}] == 0,
Subscript[ρ, {-1, -1}] == 0};
MyVariables = {Subscript[ρ, {0, 0}][t], Subscript[ρ, {1, 0}][t], Subscript[ρ, {-1, 0}][t],
Subscript[ρ, {0, 1}][t], Subscript[ρ, {1, 1}][t],
Subscript[ρ, {-1, 1}][t], Subscript[ρ, {0, -1}][t],
Subscript[ρ, {1, -1}][t], Subscript[ρ, {-1, -1}][t]};
NDSolve[{MyEquations, MyBC}, MyVariables, {t, 0, 500}];
Plot[Subscript[ρ, {0, 0}][t] /. s, {t, 0, 500}]


Now as you see, going to time t=500 gives a steady state solution, but:

1- Takes quite a while to solve

2- For my large system, I need t=10^6 in order to reach the steady state solution.

If I just replace 500 with 10^6, not only that solving will take forever, but also the kernel of Mathematica will crash.

What do I need? I probably need some old fashioned Runge-Kutta solver, where I solve this set of differential equations progressively until rho[t]-(rho[t-dt] is comparable to machine precision (or to some predefined precision). I don't need interpolation!

Now if I try to solve this with NSolve:

NSolve[#[] == 0 & /@ MyEquations /. Δ -> 3, MyVariables]


Then:

1- t will still appear on the other side of the equation.

2- This will take forever with a huge system. There's no way to take a limit of t->Infinity before solving the system.

Just for completeness, I would like to point out that these equations for this simplified system can be further simplified using the famous RWA (Rotating Wave Approximation), but this is not possible in my larger system because there's multiple generators of rotations (multiple angular momenta) involved there.

• Often, setting derivatives to zero gives a steady-state solution without actually solving the differential equations. Readers probably can provide good answers, if you include a few sample equations. Dec 31 '15 at 10:13
• @bbgodfrey Thank you for the tip, I'm trying that now. It's difficult to provide examples because it's 100+ equations. It's actually the Liouville equations of a density matrix system I'm trying to solve. Hope that helps. Dec 31 '15 at 10:14
• @bbgodfrey By setting the derivatives to zero, I still have complex exponentials on the other side of the equations that depend on time. I don't seem to find a way to put a value for them. Mathematica is unable to take the limit of that when t->Infinity, and NSolveing that system is very, very slow... What would you suggest? Is there any other way to get the steady state solution? Dec 31 '15 at 10:46
• @bbgodfrey Thank you. I updated my question with a simplified set of equations... hope that provides an exmaple of the problem. Dec 31 '15 at 17:09
• Thanks to adaptive step sizes, running to t=10^6 won't take 2000 times as long as running to t=500 if the system does settle down to a steady state. Jan 2 '16 at 3:02

It is convenient to gather the components of this ODE system into MyEquations, MyBC, and MyVariables, as I did while editing the question to correct a few transcription errors, so that the NDSolve call can be written as

s = NDSolve[{MyEquations, MyBC}, MyVariables, {t, 0, 500}];
Grid[Table[Plot[Evaluate[ReIm[Chop[Subscript[ρ, {i1, i2}][t] /. s]]],
{t, 0, 100}], {i1, -1, 1}, {i2, -1, 1}]] Although four of the solutions remain oscillatory, the oscillations are becoming very small, and Subscript[ρ, {0, 0}] has reached a steady state.

LogPlot[Chop[Subscript[ρ, {0, 0}][t] /. s], {t, 0, 500}] with a value of

Subscript[ρ, {0, 0}][t] /. s /. t -> 500
(* {0.0000124782 - 1.20307*10^-19 I} *)


It seems plausible to seek a steady state by setting derivatives equal to zero in MyEquations, as I suggested in a comment above.

steady = MyEquations /. Equal[z1_, z2_] -> Equal[0, z2];
ss = Solve[steady, MyVariables] // Simplify // Flatten;
Subscript[ρ, {0, 0}][t] /. ss
(* 0.0000124989 + 2.38692*10^-9 I *)


which agrees well with the numerical result at large t. For this ODE system, at least, the approach just outlined indeed yields the steady state.

Alternative Solution

Only four of the nine dependent variables exhibit rapid oscillation, suggesting that this oscillation can be eliminated by a change of variables. Before doing so, let us eliminate the function Subscript, which often introduces needless difficulties.

eqs = MyEquations /. Subscript[ρ, {i1_, i2_}] -> ρ[i1, i2]
(* {((I/2)*Derivative[ρ[0, 0]][t])/Pi == (-0.5*ρ[-1, 0][t])/E^((2*I)*Pi*t*Δ) + 0.5*E^((2*I)*Pi*t*Δ)*ρ[0, -1][t] -
(0. + 1.273255460229472*I)*ρ[0, 0][t] + 0.5*E^((2*I)*Pi*t*Δ)*ρ[0, 1][t] - (0.5*ρ[1, 0][t])/E^((2*I)*Pi*t*Δ),
((I/2)*Derivative[ρ[1, 0]][t])/Pi == -0.5*E^((2*I)*Pi*t*Δ)*ρ[0, 0][t] + 0.5*E^((2*I)*Pi*t*Δ)*ρ[1, -1][t] -
(0. + 0.6366356878618905*I)*ρ[1, 0][t] + Δ*ρ[1, 0][t] + 0.5*E^((2*I)*Pi*t*Δ)*ρ[1, 1][t],
((I/2)*Derivative[ρ[-1, 0]][t])/Pi == 0.5*E^((2*I)*Pi*t*Δ)*ρ[-1, -1][t] - (0. + 0.6366356878618905*I)*ρ[-1, 0][t] + Δ*ρ[-1, 0][t] +
0.5*E^((2*I)*Pi*t*Δ)*ρ[-1, 1][t] - 0.5*E^((2*I)*Pi*t*Δ)*ρ[0, 0][t],
((I/2)*Derivative[ρ[0, 1]][t])/Pi ==
(-0.5*ρ[-1, 1][t])/E^((2*I)*Pi*t*Δ) + (0.5*ρ[0, 0][t])/E^((2*I)*Pi*t*Δ) - (0. + 0.6366356878618905*I)*ρ[0, 1][t] - Δ*ρ[0, 1][t] -
(0.5*ρ[1, 1][t])/E^((2*I)*Pi*t*Δ),
((I/2)*Derivative[ρ[1, 1]][t])/Pi == ((I/2)*(0.00005 + 4*ρ[0, 0][t]))/Pi -
0.5*E^((2*I)*Pi*t*Δ)*ρ[0, 1][t] + (0.5*ρ[1, 0][t])/E^((2*I)*Pi*t*Δ) - (0. + 0.000015915494309189534*I)*ρ[1, 1][t],
((I/2)*Derivative[ρ[-1, 1]][t])/Pi == (0.5*ρ[-1, 0][t])/E^((2*I)*Pi*t*Δ) - (0. + 0.000015915494309189534*I)*ρ[-1, 1][t] -
0.5*E^((2*I)*Pi*t*Δ)*ρ[0, 1][t],
((I/2)*Derivative[ρ[0, -1]][t])/Pi == (-0.5*ρ[-1, -1][t])/E^((2*I)*Pi*t*Δ) -
(0. + 0.6366356878618905*I)*ρ[0, -1][t] - Δ*ρ[0, -1][t] + (0.5*ρ[0, 0][t])/E^((2*I)*Pi*t*Δ) - (0.5*ρ[1, -1][t])/E^((2*I)*Pi*t*Δ),
((I/2)*Derivative[ρ[1, -1]][t])/Pi == -0.5*E^((2*I)*Pi*t*Δ)*ρ[0, -1][t] - (0. + 0.000015915494309189534*I)*ρ[1, -1][t] +
(0.5*ρ[1, 0][t])/E^((2*I)*Pi*t*Δ),
((I/2)*Derivative[ρ[-1, -1]][t])/Pi == (0. - 0.000015915494309189534*I)*ρ[-1, -1][t] +
(0.5*ρ[-1, 0][t])/E^((2*I)*Pi*t*Δ) - 0.5*E^((2*I)*Pi*t*Δ)*ρ[0, -1][t] + ((I/2)*(0.00005 + 4*ρ[0, 0][t]))/Pi} *)

bc = MyBC /. Subscript[ρ, {i1_, i2_}] -> ρ[i1, i2]
(* {ρ[0, 0] == 1, ρ[1, 0] == 0, ρ[-1, 0] == 0, ρ[0, 1] == 0, ρ[1, 1] == 0,
ρ[-1, 1] == 0, ρ[0, -1] == 0, ρ[1, -1] == 0, ρ[-1, -1] == 0} *)

var = MyVariables /. Subscript[ρ, {i1_, i2_}] -> ρ[i1, i2]
(* {ρ[0, 0][t], ρ[1, 0][t], ρ[-1, 0][t], ρ[0, 1][t], ρ[1, 1][t], ρ[-1, 1][t],
ρ[0, -1][t], ρ[1, -1][t], ρ[-1, -1][t]} *)


Note that the value of Δ has not yet been specified in eqs. Next, make the substitution

{ρ[-1, 0][t] -> E^(2 I π t Δ) σ[-1, 0][t], ρ[1, 0][t] -> E^(2 I π t Δ) σ[1, 0][t],
ρ[0, -1][t] -> E^(-2 I π t Δ) σ[0, -1][t], ρ[0, 1][t] -> E^(-2 I π t Δ) σ[0, 1][t]}


which after additional manipulations eliminates E^(2 I π t Δ) from the equations.

eqstt = First[#] == Simplify[Last[#] /.
{ρ[-1, 0][t] -> E^(2 I π t Δ) σ[-1, 0][t], ρ[1, 0][t] -> E^(2 I π t Δ) σ[1, 0][t],
ρ[0, -1][t] -> E^(-2 I π t Δ) σ[0, -1][t], ρ[0, 1][t] -> E^(-2 I π t Δ) σ[0, 1][t]}] &
/@ eqs;
eqstt[] = Thread[E^(-2 I π t Δ) eqstt[] /. Derivative[ρ[1, 0]][t] ->
(Unevaluated[D[ρ[1, 0][t], t]] /. ρ[1, 0][t] -> E^(2 I π t Δ) σ[1, 0] [t]), Equal]
// Expand;
eqstt[] = Thread[E^(-2 I π t Δ) eqstt[] /. Derivative[ρ[-1, 0]][t] ->
(Unevaluated[D[ρ[-1, 0][t], t]] /. ρ[-1, 0][t] -> E^(2 I π t Δ) σ[-1, 0] [t]), Equal]
// Expand;
eqstt[] = Thread[E^(2 I π t Δ) eqstt[] /. Derivative[ρ[0, 1]][t] ->
(Unevaluated[D[ρ[0, 1][t], t]] /. ρ[0, 1][t] -> E^(-2 I π t Δ) σ[0, 1] [t]), Equal]
// Expand;
eqstt[] = Thread[E^(2 I π t Δ) eqstt[] /. Derivative[ρ[0, -1]][t] ->
(Unevaluated[D[ρ[0, -1][t], t]] /. ρ[0, -1][t] -> E^(-2 I π t Δ) σ[0, -1] [t]), Equal]
// Expand;

sstt = NDSolve[{eqstt /. Δ -> 3, bctt}, vartt, {t, 0, 500}];
Grid[Table[Plot[Evaluate[ReIm[ρ[i1, i2][t] /.
{ρ[-1, 0][t] -> σ[-1, 0][t], ρ[1, 0][t] -> σ[1, 0][t], ρ[0, -1][t] -> σ[0, -1][t],
ρ[0, 1][t] -> σ[0, 1][t]} /. sstt]], {t, 0, 500}], {i1, -1, 1}, {i2, -1, 1}]] The AbsoluteTiming of this computation is almost two orders of magnitude less than that of the solution provided a few days earlier.

The steady state solution is obtained by setting derivatives to zero.

steady = eqstt /. Δ -> 3 /. {ρ[0, 0]'[t] -> 0, σ[1, 0]'[t] -> 0, σ[-1, 0]'[t] -> 0,
σ[0, 1]'[t] -> 0, ρ[1, 1]'[t] -> 0, ρ[-1, 1]'[t] -> 0, σ[0, -1]'[t] -> 0,
ρ[1, -1]'[t] -> 0, ρ[-1, -1]'[t] -> 0};
ss = Solve[steady, vartt] // Simplify // Flatten // Chop
(* {ρ[0, 0][t] -> 0.0000124767, σ[1, 0][t] -> -0.0000748591 - 7.94299*10^-6 I,
σ[-1, 0][t] -> -0.0000748591 - 7.94299*10^-6 I, σ[0, 1][t] -> .0000748591 +
7.94299*10^-6 I,
ρ[1, 1][t] -> 0.499994, ρ[-1, 1][t] -> -0.499073,
σ[0, -1][t] -> -0.0000748591 + 7.94299*10^-6 I, ρ[1, -1][t] -> -0.499073,
ρ[-1, -1][t] -> 0.499994} *)


which agree with the corresponding NDSolve solutions at t -> 500 to at least six significant figures.

A cursory examination of the OP's full 81 x 81 system of equations suggests that the transformation used here will work there as well, although the computations are likely to be one to two orders of magnitude slower. However, there appear to be additional symmetries in the equations that, if taken advantage of, may further decrease computational time.

• Thank you for the answer. What about going to t=10^6. That's actually the real problem? My huge set of equations oscillate until that time... and using NSolve there is very slow... Jan 1 '16 at 0:49
• @TheQuantumPhysicist The point of my answer was to demonstrate that setting derivatives to zero and using Solve on the resulting equations gives the steady-state answer. Presumably, you can do the same for your complete set of equations too. To be sure, using Solve or NSolve is likely to be slow, but that seems inevitable for such a large set of equations. I suggest you try it on some intermediate sized sets of equations to verify that the method used above for 9 equations works for, say, 20 equations, then 50, finally 100. Jan 1 '16 at 4:09
• Thanks for the help. To be honest it's quite disappointing that Mathematica deals with differential equations in this weird way... While it's helpful for some applications, but who says all people require a range to solve their differential equations instead of just points?! This is deliberate wasting of memory and resources for no good reason... Jan 1 '16 at 10:07
• @TheQuantumPhysicist I take it you would like Mathematica to integrate from zero to some very large number, discarding earlier results as it goes in order to conserve memory. Actually, it is possible to do this by using the final result of computing to t1 as initialization for a computation to t2, etc. Also, it is possible to run NDSolve using Runge-Kutta and a fixed time step, allow the check you described in your question. Do you need a steady-state in all variables or only in some? I shall take another look at your equations later today. Best wishes for the New Year. Jan 1 '16 at 10:29
• Thank you very much for caring! Yes exactly. I need the computation up to some very large number. The problem with evaluating just the last number is that while doing it Mathematica assumes that I need all the values in between (thus creating an interpolation function), which gets my kernel to crash because it requires so much memory. I need the steady-state solutions only for diagonal components of the density matrix; i.e., components where the two subscript indices are equal. For the simple example I provided it's 3 components. Jan 1 '16 at 12:29

To help with the memory problem, you could save only the last value by using

s = NDSolve[{MyEquations, MyBC}, MyVariables, {t, 500, 500}];


and extract the final values using

Table[Subscript[ρ, {i1, i2}][t], {i1, -1, 1}, {i2, -1, 1}] /.s/.t->tmax


As an alternative to NSolve for finding the equilibrium, how about FindRoot?

• Would this work even though I have the initial conditions at t=0? Jan 2 '16 at 2:41
• Yes -- I was surprised by this trick as well. See mathematica.stackexchange.com/a/59577/6358 Jan 2 '16 at 2:46
• Wow... I'm really surprised! All I can think about right now is: How the hell didn't I think of this before!! I'm +1ing this but the problem still exists that the process is very slow. Thanks for the info :) Jan 2 '16 at 3:41