Help speeding up solution of system of differential equations using Laplace transform

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.

• Just leave away the Simplify. The rest terminates after half a minute. Sep 26, 2019 at 4:52
• You could try Refine instead: reference.wolfram.com/language/ref/Refine.html Sep 26, 2019 at 5:21
• @HenrikSchumacher I tried this and it still takes forever. I tried breaking it down line by line and the last line is what appears to be causing the problem. Oct 1, 2019 at 23:32
• There're 2 rules for Δ1 and no rule for Δ2 inside replace3, is this intended or a typo? Oct 2, 2019 at 13:44
• @xzczd Definitely a typo, although I fixed it and it appears to make no difference. Oct 2, 2019 at 15:35

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 -> 316, 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}} *)

• 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. Oct 6, 2019 at 0:20
• @CalebHorwitz With Plot[sol[[4]] /. t -> 1 // Evaluate, {v, 0, 1000}, PlotRange -> All] I obtain the same result. What parameter do you choose? Oct 6, 2019 at 5:24
• Upon fixing a typo, I get the same result for both. Oct 7, 2019 at 0: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? Feb 24, 2020 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 LeafCount`ed 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. Feb 24, 2020 at 15:01