# 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. – Henrik Schumacher Sep 26 '19 at 4:52
• You could try Refine instead: reference.wolfram.com/language/ref/Refine.html – infinitezero Sep 26 '19 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. – Caleb Horwitz Oct 1 '19 at 23:32
• There're 2 rules for Δ1 and no rule for Δ2 inside replace3, is this intended or a typo? – xzczd Oct 2 '19 at 13:44
• @xzczd Definitely a typo, although I fixed it and it appears to make no difference. – Caleb Horwitz Oct 2 '19 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. – Caleb Horwitz Oct 6 '19 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? – xzczd Oct 6 '19 at 5:24
• Upon fixing a typo, I get the same result for both. – Caleb Horwitz Oct 7 '19 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? – 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 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. – xzczd Feb 24 '20 at 15:01