0
$\begingroup$

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]
```
Improve this question
$\endgroup$
2
  • 1
    $\begingroup$ 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. $\endgroup$
    – tueda
    Apr 30, 2022 at 8:47
  • $\begingroup$ 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. $\endgroup$
    – josh
    Apr 30, 2022 at 13:10

2 Answers 2

Reset to default
3
$\begingroup$

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, π}]

enter image description here

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}}]

enter image description here

Improve this answer
$\endgroup$
3
  • $\begingroup$ 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}$ $\endgroup$
    – LionCereals
    May 7, 2022 at 8:36
  • $\begingroup$ Also, I made a sign error: It's r32'[t] = -1/2*a* (....), not +. $\endgroup$
    – LionCereals
    May 7, 2022 at 8:57
  • $\begingroup$ Nevermind, an application of FullSimplify[others] did the job. $\endgroup$
    – LionCereals
    May 7, 2022 at 9:11
2
$\begingroup$

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];
solc = FullSimplify[ToRadicals@Flatten@Solve[%%, %], 
    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.

Improve this answer
$\endgroup$
2
  • $\begingroup$ 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. $\endgroup$
    – Akku14
    May 1, 2022 at 6:25
  • $\begingroup$ @Akku14 Thank you. $\endgroup$
    – bbgodfrey
    May 1, 2022 at 20:26

Your Answer

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

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