# Why does DSolve not yield a solution to this coupled differential equation?

I am having 16 coupled equations that I try to solve using DSolve. However, when running my notebook, the Kernel keeps running for over 30 mins, so there might be an error somewhere. I know that there is an analytic solution to the set of equations, and I also know that I can't solve them by hand. I checked the equations for typos for multiple hours and I really think that there is no typo in there anymore.

What I tried:

1. Specified that my input parameters are real
2. Solved the coupled equations without the imaginary part, which gives the solution I am looking for.
3. If I add the imaginary part, something goes wrong.

Maybe someone here immediately sees an error!

$Assumptions = \[Omega] \[Element] Reals && a \[Element] Reals eqns = {r11'[t] == -1/2*a*(2*r11[t] - r22[t] - r33[t]), r22'[t] == -1/2*a*(2*r22[t] - r11[t] - r44[t]), r33'[t] == -1/2*a*(2*r33[t] - r44[t] - r11[t]), r44'[t] == -1/2*a*(2*r44[t] - r33[t] - r22[t]), r12'[t] == I*\[Omega]*r12[t] - 1/2*a*(2*r12[t] - r21[t] - r34[t]), r13'[t] == I*\[Omega]*r13[t] - 1/2*a*(2*r13[t] - r24[t] - r31[t]), r14'[t] == 2*I*\[Omega]*r14[t] - 1/2*a*(2*r14[t] - r23[t] - r32[t]), r21'[t] == -I*\[Omega]*r21[t] - 1/2*a*(2*r21[t] - r12[t] - r43[t]), r23'[t] == -1/2*a*(2*r23[t] - r41[t] - r14[t]), r24'[t] == I*\[Omega]*r24[t] - 1/2*a*(2*r24[t] - r13[t] - r42[t]), r31'[t] == -I*\[Omega]*r31[t] - 1/2*a*(2*r31[t] - r42[t] - r13[t]), r32'[t] == 1/2*a*(2*r32[t] - r41[t] - r14[t]), r34'[t] == I*\[Omega]*r34[t] - 1/2*a*(2*r34[t] - r43[t] - r12[t]), r41'[t] == -2*I*\[Omega]*r41[t] - 1/2*a*(2*r41[t] - r32[t] - r23[t]), r42'[t] == -I*\[Omega]*r42[t] - 1/2*a*(2*r42[t] - r31[t] - r24[t]), r43'[t] == -I*\[Omega]*r43[t] - 1/2*a*(2*r43[t] - r34[t] - r21[t]), r11[0] == 0, r22[0] == 1/2, r33[0] == 1/2, r44[0] == 0, r12[0] == 0, r13[0] == 0, r14[0] == 0, r21[0] == 0, r23[0] == 1/2, r24[0] == 0, r31[0] == 0, r32[0] == 1/2, r34[0] == 0, r41[0] == 0, r42[0] == 0, r43[0] == 0}; sol = DSolve[ eqns, {r11, r12, r13, r14, r22, r33, r44, r21, r23, r24, r31, r32, r34, r41, r42, r43}, t] $$$$  • Though this may not be an answer you want to get, simply the system may be too complicated for Mathematica to automatically handle. If you set \[Omega] and a to certain numbers (say, both are 1), then DSolve and NDSolve work well. Commented Apr 30, 2022 at 8:47 • I don't think 30 minutes is long enough if there are no typos and there is an analytic solution that can be solved manually. I'd wait at least 4 hours. – josh Commented Apr 30, 2022 at 13:10 ## 2 Answers This system can be solved using LaplaceTransform. The first step is to separate the system into a list of ODEs and a list of initial conditions: ClearAll["Global*"] eqnt = {r11'[t] == -1/2*a*(2*r11[t] - r22[t] - r33[t]), r22'[t] == -1/2*a*(2*r22[t] - r11[t] - r44[t]), r33'[t] == -1/2*a*(2*r33[t] - r44[t] - r11[t]), r44'[t] == -1/2*a*(2*r44[t] - r33[t] - r22[t]), r12'[t] == I*\[Omega]*r12[t] - 1/2*a*(2*r12[t] - r21[t] - r34[t]), r13'[t] == I*\[Omega]*r13[t] - 1/2*a*(2*r13[t] - r24[t] - r31[t]), r14'[t] == 2*I*\[Omega]*r14[t] - 1/2*a*(2*r14[t] - r23[t] - r32[t]), r21'[t] == -I*\[Omega]*r21[t] - 1/2*a*(2*r21[t] - r12[t] - r43[t]), r23'[t] == -1/2*a*(2*r23[t] - r41[t] - r14[t]), r24'[t] == I*\[Omega]*r24[t] - 1/2*a*(2*r24[t] - r13[t] - r42[t]), r31'[t] == -I*\[Omega]*r31[t] - 1/2*a*(2*r31[t] - r42[t] - r13[t]), r32'[t] == 1/2*a*(2*r32[t] - r41[t] - r14[t]), r34'[t] == I*\[Omega]*r34[t] - 1/2*a*(2*r34[t] - r43[t] - r12[t]), r41'[t] == -2*I*\[Omega]*r41[t] - 1/2*a*(2*r41[t] - r32[t] - r23[t]), r42'[t] == -I*\[Omega]*r42[t] - 1/2*a*(2*r42[t] - r31[t] - r24[t]), r43'[t] == -I*\[Omega]*r43[t] - 1/2*a*(2*r43[t] - r34[t] - r21[t])}; ics = {r11[0] == 0, r22[0] == 1/2, r33[0] == 1/2, r44[0] == 0, r12[0] == 0, r13[0] == 0, r14[0] == 0, r21[0] == 0, r23[0] == 1/2, r24[0] == 0, r31[0] == 0, r32[0] == 1/2, r34[0] == 0, r41[0] == 0, r42[0] == 0, r43[0] == 0};  The next steps are to transform the equations, to apply the initial conditions, to solve for the transforms (vars) and to transform back. eqns = LaplaceTransform[eqnt, t, s] /. Rule @@@ ics; vars = Cases[eqns, _LaplaceTransform, ∞] // DeleteDuplicates; soln = Solve[eqns, vars] // Flatten; rsoln = InverseLaplaceTransform[soln, s, t];  The 16 solutions for the functions $$r_{ij}(t)$$ given in the list rsoln fall into three categories. Eight of the solutions are trivial. Four of the solutions are independent of $$\omega$$. The remaining four solutions may have imaginary parts. trivial = Select[rsoln, Last[#] == 0 &]; freeω = Select[Complement[rsoln, trivial], FreeQ[#, ω] &]; others = Complement[rsoln, trivial, freeω];  For real, positive $$a$$, the solutions that are independent of $$\omega$$ decay exponentially. f = freeω /. a -> 1 // Values; Plot[f, {t, 0, π}]  The solutions with imaginary parts can also be plotted for specific values of $$\omega$$ and $$a$$ g = others /. ω -> 1 /. a -> -1 // Values; ReImPlot[g, {t, 0, 2 π}, PlotRange -> {All, {-20, 20}}]  • Thanks. One follow-up: Is there a way to simplifiy the 4 equations in others? For example,$r_{14}$has the solution$r_{14} = \frac{a}{2}e^{-at}\left(\frac{2 i \omega[\cosh(\Omega t) - 1]}{\Omega^2} + \frac{\sinh(\Omega t)}{\Omega}\right)$,$\Omega = \sqrt{a^2 - 4\omega^2}$. In particular, I care about$r_{32}(t) = \frac{1}{2}e^{-a t} \frac{a^2 \cosh(\Omega t) - 4\omega^2}{\Omega^2}\$ Commented May 7, 2022 at 8:36
• Also, I made a sign error: It's r32'[t] = -1/2*a* (....), not +. Commented May 7, 2022 at 8:57
• Nevermind, an application of FullSimplify[others] did the job. Commented May 7, 2022 at 9:11

Because the ODEs are linear with constant coefficients, DSolve must be able to solve them, although perhaps slowly. Here we provide a quick approach.

With the array of 16 ODEs (but not the initial conditions) set to eqns, DSolve finds a solution in seconds.

sol = First@DSolve[eqns, {r11[t], r12[t], r13[t], r14[t], r22[t], r33[t],
r44[t], r21[t], r23[t], r24[t], r31[t], r32[t], r34[t], r41[t],
r42[t], r43[t]}, t, Assumptions -> ω ∈ Reals && a ∈ Reals];


although sol is far too long to reproduce here. Next evaluate the initial conditions,

sol /. t -> 0;
{r11[0] == 0, r22[0] == 1/2, r33[0] == 1/2, r44[0] == 0,
r12[0] == 0, r13[0] == 0, r14[0] == 0, r21[0] == 0, r23[0] == 1/2,
r24[0] == 0, r31[0] == 0, r32[0] == 1/2, r34[0] == 0, r41[0] == 0,
r42[0] == 0, r43[0] == 0} /. %;
% // Simplify;
Array[C, 16];
Assumptions -> ω ∈ Reals && a ∈ Reals]
(* {C[1] -> 0, C[2] -> 1/2, C[3] -> 1/2, C[4] -> 0, C[5] -> 0, C[6] -> 0,
C[7] -> 0, C[8] -> 0, C[9] -> 0, C[10] -> 0, C[11] -> 0, C[12] -> 0,
C[13] -> 0, C[14] -> 1/2, C[15] -> 1/2, C[16] -> 0} *)


(My thanks to Akku14 for pointing out the simplification above in solc.) Finally,

sode = Simplify[sol /. solc]


gives the desired solution. The simpler solutions are given by

{r11[t], r22[t], r33[t], r44[t]} /. sode
(* {1/4 - 1/4 E^(-2 a t), 1/4 (1 + E^(-2 a t)), 1/4 (1 + E^(-2 a t)),
1/4 - 1/4 E^(-2 a t)} *)
{r12[t], r21[t], r34[t], r43[t], r13[t], r24[t], r31[t], r42[t]} /. sode
(* {0, 0, 0, 0, 0, 0, 0, 0} *)


and {r13, r14, r15, r16} involve lengthy Root functions.

• to @bbgodfrey, with  // ToRadicals // FullSimplify[#, Assumptions -> [Omega] [Element] Reals && a [Element] Reals] &  you can show C[14] and C[15] to be == 1/2 both. Commented May 1, 2022 at 6:25
• @Akku14 Thank you. Commented May 1, 2022 at 20:26