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[0]==y00,y1[0]==y01,...,yN[0]==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.
Please advise.
Please feel free to ask for any additional details. Thanks.
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}][0] == 1, Subscript[ρ, {1, 0}][0] == 0,
Subscript[ρ, {-1, 0}][0] == 0, Subscript[ρ, {0, 1}][0] == 0,
Subscript[ρ, {1, 1}][0] == 0, Subscript[ρ, {-1, 1}][0] == 0,
Subscript[ρ, {0, -1}][0] == 0, Subscript[ρ, {1, -1}][0] == 0,
Subscript[ρ, {-1, -1}][0] == 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[#[[2]] == 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.
t->Infinity, andNSolveing that system is very, very slow... What would you suggest? Is there any other way to get the steady state solution? $\endgroup$