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}]
EOM
, which has aLeafCount
of a staggeringly large1177079119
. 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. (ApplyingSimplify
toEOM
will take forever itself, however.) $\endgroup$