1

I'm trying to speed up this code that seems to take forever to evaluate. I'm solving a system of 9 coupled differential equations using a Laplace transform (a version of this code that works perfectly for a system of 4 equations was created in an answer to this question previously posted by me). 

h2 = {{0, 0, -Ω1}, {0, -Δ1 + Δ2 + k1 v - k2 v, -Ω2}, {-Ω1, -Ω2, -Δ1 + k1 v}};
ρ = {{ρ11[t], ρ12[t], ρ13[t]}, {ρ21[t], ρ22[t], ρ23[t]}, {ρ31[t], ρ32[t], ρ33[t]}};
ρprime = -I (h2.ρ - ρ.h2) + {{γ31 ρ33[t] + γ21 ρ22[t], -(1/2) γ21 ρ12[t], -(1/2) (γ31 + γ32) ρ13[t]}, {-(1/2) γ21 ρ21[t], -γ21 ρ22[t] + γ32 ρ33[t], -(1/2) (γ21 + γ31 + γ32) ρ23[t]}, {-(1/2) (γ31 + γ32) ρ31[t], -(1/2) (γ21 + γ31 + γ32) ρ32[t], -ρ33[t] (γ31 + γ32)}};
replace3 = {Δ1 -> ( 2 π)/(500*10^-9)*10^3, Δ2 -> ( 2 π)/(500*10^-9)*10^3, γ21 -> 1/(16*10^-9), γ31 -> 1/(16*10^-9), γ32 -> 1/(16*10^-9), Ω1 -> 10^9,Ω2 -> 10^9, k1 -> ( 2 π)/(500*10^-9), k2 -> ( 2 π)/(500*10^-9)};

var = Flatten@ρ;

{eq, ic} = {D[var, t] == Flatten@ρprime // Thread, 
var == {1, 0, 0, 0, 0, 0, 0, 0, 0} /. t -> 0 // Thread} /. replace3;

tvar = LaplaceTransform[var, t, s];
tsol = tvar /. First@Solve[LaplaceTransform[eq, t, s] /. Rule @@@ ic, tvar] // Simplify;
(sol = InverseLaplaceTransform[tsol, s, t]) // AbsoluteTiming;

Thanks for the help.

Caleb Horwitz
  • 53
  • 4

1 Answers1

1

Since the system is in the form $v'(t)=a.v(t)$ where $a$ ia a constant matrix and $v(t)$ is a vector, it can be solved with MatrixExp:

{barray, marray} = CoefficientArrays[Flatten@ρprime /. replace3, var]

sol = MatrixExp[marray t, {1, 0, 0, 0, 0, 0, 0, 0, 0}]; // AbsoluteTiming
(* {23.7716, Null} *)

Be careful the symbolic solution is rather large:

LeafCount@sol
(* 74603173 *)

So it's not easy to check its correctness with Simplify, but we can compare it to the numeric solution:

sol /. {v -> 3`16, t -> 10^-7} // AbsoluteTiming
(* {265.675, 
    {0.52868393 + 0.*10^-9 I, -0.01450599 + 0.18208956 I, 
    -0.0403699 - 0.0147429 I, -0.01450599 - 0.18208956 I, 
     0.46386734 + 0.*10^-9 I, -0.03540213 + 0.01425282 I,
    -0.0403699 + 0.0147429 I, -0.0354021 - 0.0142528 I, 
     0.00744873 + 0.*10^-9 I}} *)

NDSolveValue[{eq, ic} /. replace3 /. v -> 3, var /. t -> 10^-7, {t, 0, 10^-7}, 
  WorkingPrecision -> 16, MaxSteps -> Infinity] // AbsoluteTiming
(*{3.10284, 
  {0.5286838344358492, -0.01450596346156398 + 0.1820894384493129 I, 
  -0.04037047804336591 - 0.01474265129379949 I, -0.01450596346156398 - 0.1820894384493129 I, 
   0.4638674822129191, -0.03540122024132100 + 0.01425467017942954I, 
  -0.04037047804336591 + 0.01474265129379949 I, -0.03540122024132100 - 0.01425467017942954I,  
   0.00744868335123165}} *)
xzczd
  • 65,995
  • 9
  • 163
  • 468
  • There seems to be a discrepancy between this method and the Laplace Transform one. If you apply the Matrix Exp method to this and plot Re[sol[[4]] for an arbitrary t, from {v,0,1000}, and then you do the same but with the Laplace Transform method, you get different plots. The Matrix Exp one is negative, which for physical reasons should not be the case. – Caleb Horwitz Oct 06 '19 at 00:20
  • @CalebHorwitz With Plot[sol[[4]] /. t -> 1 // Evaluate, {v, 0, 1000}, PlotRange -> All] I obtain the same result. What parameter do you choose? – xzczd Oct 06 '19 at 05:24
  • Upon fixing a typo, I get the same result for both. – Caleb Horwitz Oct 07 '19 at 00:31
  • this method takes infinitely long to complete for matrices 4x4 and greater. Is there any way to do this with a 4x4 matrix that doesn't take forever to evaluate? – Caleb Horwitz Feb 24 '20 at 12:42
  • @caleb Hard to give advice without seeing the specific system. BTW, I'm now with a mobile phone so can't check, but have you LeafCounted the solution to 2×2 system? If it's much smaller compared to that of a 3×3 system, then it'll probably be difficult to achieve your goal. – xzczd Feb 24 '20 at 15:01
  • I've posted a new question here – Caleb Horwitz Feb 25 '20 at 08:25