# Why Mathematica is not producing output and taking too much time

I'm trying to solve the given system of ODES but the Mathematica is taking too much time and not producing any output. I was trying to check the error by evaluating one one command but there was no error in any command but the equations EOM2, and EOM3 was taking too much time when I was trying to evaluate the equations.

For simple case aa=0, code works, but when I take non-zero aa, it takes a long time and didn't produce output.

Can anyone please guide me how can I fix this problem? Is there any command in Mathematica that can be used to obtain the fast output?

 R2[r_, θ_] := r^2 + aa^2 Cos[θ]^2;
TR[r_, θ_] := r^2 - 2 M r + aa^2;

gtt[r_, θ_] := -(1 - (2 M r)/R2[r, θ]);
gtϕ[r_, θ_] := -(( 2 r M aa Sin[θ]^2)/
R2[r, θ]);
gϕϕ[
r_, θ_] := (r^2 +
aa^2 + (2  M r (aa^2) )/
R2[r, θ] Sin[θ]^2) Sin[θ]^2;
grr[r_, θ_] := R2[r, θ]/TR[r, θ];
gθθ[r_, θ_] := R2[r, θ];

gUtt[r_, θ_] := -(1/
TR[r, θ]) (r^2 +
aa^2 + (2  M r (aa^2) )/ R2[r, θ] Sin[θ]^2);
gUtϕ[r_, θ_] := -((2 M aa r)/(
TR[r, θ] R2[r, θ]));

gUϕϕ[r_, θ_] := (
TR[r, θ] - aa^2 Sin[θ]^2)/(
TR[r, θ] R2[r, θ] Sin[θ]^2);
gUrr[r_, θ_] := TR[r, θ]/R2[r, θ];
gUθθ[r_, θ_] := 1/R2[r, θ];
M = 1; n = 4;
glo = FullSimplify[{ {gtt[r, θ], 0, 0,
gtϕ[r, θ]}, {0, grr[r, θ], 0, 0}, {0, 0,
gθθ[r, θ], 0}, {gtϕ[r, θ], 0, 0,
gϕϕ[r, θ]}}];
gup = FullSimplify[{ {gUtt[r, θ], 0, 0,
gUtϕ[r, θ]}, {0, gUrr[r, θ], 0, 0}, {0, 0,
gUθθ[r, θ], 0}, {gUtϕ[r, θ], 0,
0, gUϕϕ[r, θ]}}];

dglo = Simplify[Det[glo]];
crd = {t, r, θ, ϕ};

Xup = {t[τ], r[τ], θ[τ], ϕ[τ]};
Vup = {Vt, Vr, Vθ, Vϕ};
Pup = {Pt[τ], Pr[τ], Pθ[τ], Pϕ[τ]};

Sup = {{Stt[τ], Str[τ], Stθ[τ],
Stϕ[τ]},
{Srt[τ], Srr[τ], Srθ[τ], Srϕ[τ]},
{Sθt[τ], Sθr[τ], Sθθ[τ],
Sθϕ[τ]},
{Sϕt[τ], Sϕr[τ], Sϕθ[τ],
Sϕϕ[τ]}};

christoffel =
Table[(1/2)*
Sum[(gup[[i, s]])*(D[glo[[s, k]], crd[[j]] ] +
D[glo[[s, j]], crd[[k]] ] - D[glo[[j, k]], crd[[s]] ]), {s, 1,
n}], {i, 1, n}, {j, 1, n}, {k, 1, n}] ;

riemann =
Table[ D[christoffel[[i, j, l]], crd[[k]] ] -
D[christoffel[[i, j, k]], crd[[l]] ] +
Sum[christoffel[[s, j, l]] christoffel[[i, k, s]] -
christoffel[[s, j, k]] christoffel[[i, l, s]],
{s, 1, n}], {i, 1, n}, {j, 1, n}, {k, 1, n}, {l, 1, n}] ;

loriemann =
Table[Sum[glo[[i, m]]*riemann[[m, j, k, l]], {m, 1, n}], {i, 1,
n}, {j, 1, n}, {k, 1, n}, {l, 1, n}] ;

EOM1 = Table[ D[Xup[[a]], τ] == Vup[[a]] , {a, 1, n}];

EOM2 = Table[
D[Pup[[a]], τ] + \!$$\*UnderoverscriptBox[\(∑$$, $$b = 1$$, $$n$$]$$\*UnderoverscriptBox[\(∑$$, $$c = 1$$, $$n$$]christoffel[$$[$$$$a, b, c$$$$]$$]*
Pup[$$[$$$$b$$$$]$$]*Vup[$$[$$$$c$$$$]$$]\)\) == -(1/2) \!$$\*UnderoverscriptBox[\(∑$$, $$b = 1$$, $$n$$]$$\*UnderoverscriptBox[\(∑$$, $$c = 1$$, $$n$$]$$\*UnderoverscriptBox[\(∑$$, $$d = 1$$, $$n$$]riemann[$$[$$$$a, b, c, d$$$$]$$]*Vup[$$[$$$$b$$$$]$$]*
Sup[$$[$$$$c, d$$$$]$$]\)\)\),
{a, 1, n}];
EOM3 = Table[
D[Sup[[a, b]], τ] + \!$$\*UnderoverscriptBox[\(∑$$, $$c = 1$$, $$n$$]$$\*UnderoverscriptBox[\(∑$$, $$d = 1$$, $$n$$]christoffel[$$[$$$$a, c, d$$$$]$$]*
Sup[$$[$$$$c, b$$$$]$$]*Vup[$$[$$$$d$$$$]$$]\)\) + \!$$\*UnderoverscriptBox[\(∑$$, $$c = 1$$, $$n$$]$$\*UnderoverscriptBox[\(∑$$, $$d = 1$$, $$n$$]christoffel[$$[$$$$b, c, d$$$$]$$]*
Sup[$$[$$$$a, c$$$$]$$]*Vup[$$[$$$$d$$$$]$$]\)\) ==
Pup[[a]]*Vup[[b]] - Pup[[b]]*Vup[[a]],
{a, 1, n}, {b, 1, n}];

Wfactor = 4*μ^2 + \!$$\*UnderoverscriptBox[\(∑$$, $$i = 1$$, $$4$$]$$\*UnderoverscriptBox[\(∑$$, $$j = 1$$, $$4$$]$$\*UnderoverscriptBox[\(∑$$, $$k = 1$$, $$4$$]$$\*UnderoverscriptBox[\(∑$$, $$l = 1$$, $$4$$]$$(loriemann[\([$$$$i, j, k, l$$$$]$$]*$$(Sup[\([$$$$i, j$$$$]$$])\)*\ $$(Sup[\([$$$$k, l$$$$]$$])\))\)\)\)\)\);

Wvec = Table[2/(μ*Wfactor)*(\!$$\*UnderoverscriptBox[\(∑$$, $$i = 1$$, $$4$$]$$\*UnderoverscriptBox[\(∑$$, $$k = 1$$, $$4$$]$$\*UnderoverscriptBox[\(∑$$, $$m = 1$$, $$4$$]$$\*UnderoverscriptBox[\(∑$$, $$l = 1$$, $$4$$]Sup[$$[$$$$j, i$$$$]$$]*
Pup[$$[$$$$k$$$$]$$]*$$(loriemann[\([$$$$i, k, l, m$$$$]$$])\)*$$(Sup[\([$$$$l, m$$$$]$$])\)\)\)\)\)), {j,
1, n}];

NN = 1/Sqrt[1 - \!$$\*UnderoverscriptBox[\(∑$$, $$i = 1$$, $$4$$]$$\*UnderoverscriptBox[\(∑$$, $$k = 1$$, $$4$$]$$(glo[\([$$$$i, k$$$$]$$])\)*Wvec[$$[$$$$i$$$$]$$]*
Wvec[$$[$$$$k$$$$]$$]\)\)];

{Vt, Vr, Vθ, Vϕ} = NN (Wvec + Pup);

EOM = Flatten[
Join[{EOM1, EOM2, EOM3} /.
r -> r[τ] /. θ -> θ[τ] /.
Derivative[1][r[τ]][τ] -> Derivative[1][r][τ] /.
Derivative[1][θ[τ]][τ] ->
Derivative[1][θ][τ]]];

INT1 = {t[0] == 0,
r[0] == r0, θ[0] == θ0, ϕ[0] == 0};
INT2 = {Pt[0] == 1.32288, Pr[0] == 0, Pθ[0] == 0,
Pϕ[0] == 0.07143};
INT3 = {{Stt[0] == 0, Str[0] == 0, Stθ[0] == 0,
Stϕ[0] == 0},
{Srt[0] == 0, Srr[0] == 0, Srθ[0] == 0, Srϕ[0] == 0},
{Sθt[0] == 0, Sθr[0] == 0, Sθθ[0] == 0,
Sθϕ[0] == 0},
{Sϕt[0] == 0, Sϕr[0] == 0, Sϕθ[0] == 0,
Sϕϕ[0] == 0}};
INT = Flatten[Join[{INT1, INT2, INT3}]];
r0 = 7; θ0 = Pi/2; μ = 1; aa = 0.5; M = 1;

NDSolve[Flatten[Join[{EOM, INT}]], {t, r, θ, ϕ, Pt, Pr,
Pθ, Pϕ, Stt, Str, Stθ, Stϕ, Srt, Srr,
Srθ, Srϕ,
Sθt, Sθr, Sθθ, Sθϕ,
Sϕt, Sϕr, Sϕθ, Sϕϕ}, {τ, 0,
1000}]

• Not answer, but the suit of packages xAct might be useful to you: xact.es Oct 31, 2020 at 14:15
• @mattiav27 How can I use this package?
– MMS
Oct 31, 2020 at 20:11
• It took over 30 minutes just to obtain EOM, which has a LeafCount of a staggeringly large 1177079119. The problem quickly saturated my 16 GB computer memory. On this basis, NDSolve will take an enormously long time to evaluate. If you know that any of the dependent variables are identically zero, set them to zero before the computation. Also, look for ways to simplify intermediate results. (Applying Simplify to EOM will take forever itself, however.) Nov 1, 2020 at 0:09

As I noted in a comment above, EOM as computed in the question takes over 30 minutes and has a LeafCount of 1177079119. These values can be reduced by FullSimplifying christoffel. riemann. and loriemann. For instance,

christoffel = FullSimplify@Table[...]


I also replaced Simplify by FullSimplify in the definition of dglo, although doing so did not have a significant impact. With these changes EOM can be computed in a few minutes, and

LeafCount[EOM]
(* 23063610 *)


which, will still large, is much smaller than before. Next, for convenience, define

var = Join[Xup, Pup, Flatten[Sup]];


and solve the 24 ODEs.

SetSystemOptions["NDSolveOptions" -> "DefaultSolveTimeConstraint" -> 100.];
NDSolveValue[Flatten[Join[{EOM, INT}]], var, {τ, 0, 1000}];


SetSystemOptions is needed to keep NDSolve from timing out with the initialization error message, "NDSolve::ntdv", as explained here. A plot of the solutions is given by

Plot[Evaluate[%[[;; 8]]], {τ, 0, 250}, PlotRange -> {Automatic, 10}, ImageSize -> Large,
PlotLegends -> Placed[ToString /@ var, {.9, .5}], LabelStyle -> {15, Black, Bold}]


The remaining variables are identically zero.

%% /. τ -> 1000
(* {1269.97, 8.32087, 1.5708, 56.0725, 1.2494, -0.0268073, 7.30725*10^-17,
0.0500981, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.} *)
`