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_, \[Theta]_] := r^2 + aa^2 Cos[\[Theta]]^2; TR[r_, \[Theta]_] := r^2 - 2 M r + aa^2; gtt[r_, \[Theta]_] := -(1 - (2 M r)/R2[r, \[Theta]]); gt\[Phi][r_, \[Theta]_] := -(( 2 r M aa Sin[\[Theta]]^2)/ R2[r, \[Theta]]); g\[Phi]\[Phi][ r_, \[Theta]_] := (r^2 + aa^2 + (2 M r (aa^2) )/ R2[r, \[Theta]] Sin[\[Theta]]^2) Sin[\[Theta]]^2; grr[r_, \[Theta]_] := R2[r, \[Theta]]/TR[r, \[Theta]]; g\[Theta]\[Theta][r_, \[Theta]_] := R2[r, \[Theta]]; gUtt[r_, \[Theta]_] := -(1/ TR[r, \[Theta]]) (r^2 + aa^2 + (2 M r (aa^2) )/ R2[r, \[Theta]] Sin[\[Theta]]^2); gUt\[Phi][r_, \[Theta]_] := -((2 M aa r)/( TR[r, \[Theta]] R2[r, \[Theta]])); gU\[Phi]\[Phi][r_, \[Theta]_] := ( TR[r, \[Theta]] - aa^2 Sin[\[Theta]]^2)/( TR[r, \[Theta]] R2[r, \[Theta]] Sin[\[Theta]]^2); gUrr[r_, \[Theta]_] := TR[r, \[Theta]]/R2[r, \[Theta]]; gU\[Theta]\[Theta][r_, \[Theta]_] := 1/R2[r, \[Theta]]; M = 1; n = 4; glo = FullSimplify[{ {gtt[r, \[Theta]], 0, 0, gt\[Phi][r, \[Theta]]}, {0, grr[r, \[Theta]], 0, 0}, {0, 0, g\[Theta]\[Theta][r, \[Theta]], 0}, {gt\[Phi][r, \[Theta]], 0, 0, g\[Phi]\[Phi][r, \[Theta]]}}]; gup = FullSimplify[{ {gUtt[r, \[Theta]], 0, 0, gUt\[Phi][r, \[Theta]]}, {0, gUrr[r, \[Theta]], 0, 0}, {0, 0, gU\[Theta]\[Theta][r, \[Theta]], 0}, {gUt\[Phi][r, \[Theta]], 0, 0, gU\[Phi]\[Phi][r, \[Theta]]}}]; dglo = Simplify[Det[glo]]; crd = {t, r, \[Theta], \[Phi]}; Xup = {t[\[Tau]], r[\[Tau]], \[Theta][\[Tau]], \[Phi][\[Tau]]}; Vup = {Vt, Vr, V\[Theta], V\[Phi]}; Pup = {Pt[\[Tau]], Pr[\[Tau]], P\[Theta][\[Tau]], P\[Phi][\[Tau]]}; Sup = {{Stt[\[Tau]], Str[\[Tau]], St\[Theta][\[Tau]], St\[Phi][\[Tau]]}, {Srt[\[Tau]], Srr[\[Tau]], Sr\[Theta][\[Tau]], Sr\[Phi][\[Tau]]}, {S\[Theta]t[\[Tau]], S\[Theta]r[\[Tau]], S\[Theta]\[Theta][\[Tau]], S\[Theta]\[Phi][\[Tau]]}, {S\[Phi]t[\[Tau]], S\[Phi]r[\[Tau]], S\[Phi]\[Theta][\[Tau]], S\[Phi]\[Phi][\[Tau]]}}; 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]], \[Tau]] == Vup[[a]] , {a, 1, n}]; EOM2 = Table[ D[Pup[[a]], \[Tau]] + \!\( \*UnderoverscriptBox[\(\[Sum]\), \(b = 1\), \(n\)]\( \*UnderoverscriptBox[\(\[Sum]\), \(c = 1\), \(n\)]christoffel[\([\)\(a, b, c\)\(]\)]* Pup[\([\)\(b\)\(]\)]*Vup[\([\)\(c\)\(]\)]\)\) == -(1/2) \!\( \*UnderoverscriptBox[\(\[Sum]\), \(b = 1\), \(n\)]\( \*UnderoverscriptBox[\(\[Sum]\), \(c = 1\), \(n\)]\( \*UnderoverscriptBox[\(\[Sum]\), \(d = 1\), \(n\)]riemann[\([\)\(a, b, c, d\)\(]\)]*Vup[\([\)\(b\)\(]\)]* Sup[\([\)\(c, d\)\(]\)]\)\)\), {a, 1, n}]; EOM3 = Table[ D[Sup[[a, b]], \[Tau]] + \!\( \*UnderoverscriptBox[\(\[Sum]\), \(c = 1\), \(n\)]\( \*UnderoverscriptBox[\(\[Sum]\), \(d = 1\), \(n\)]christoffel[\([\)\(a, c, d\)\(]\)]* Sup[\([\)\(c, b\)\(]\)]*Vup[\([\)\(d\)\(]\)]\)\) + \!\( \*UnderoverscriptBox[\(\[Sum]\), \(c = 1\), \(n\)]\( \*UnderoverscriptBox[\(\[Sum]\), \(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*\[Mu]^2 + \!\( \*UnderoverscriptBox[\(\[Sum]\), \(i = 1\), \(4\)]\( \*UnderoverscriptBox[\(\[Sum]\), \(j = 1\), \(4\)]\( \*UnderoverscriptBox[\(\[Sum]\), \(k = 1\), \(4\)]\( \*UnderoverscriptBox[\(\[Sum]\), \(l = 1\), \(4\)]\((loriemann[\([\)\(i, j, k, l\)\(]\)]*\((Sup[\([\)\(i, j\)\(]\)])\)*\ \((Sup[\([\)\(k, l\)\(]\)])\))\)\)\)\)\); Wvec = Table[2/(\[Mu]*Wfactor)*(\!\( \*UnderoverscriptBox[\(\[Sum]\), \(i = 1\), \(4\)]\( \*UnderoverscriptBox[\(\[Sum]\), \(k = 1\), \(4\)]\( \*UnderoverscriptBox[\(\[Sum]\), \(m = 1\), \(4\)]\( \*UnderoverscriptBox[\(\[Sum]\), \(l = 1\), \(4\)]Sup[\([\)\(j, i\)\(]\)]* Pup[\([\)\(k\)\(]\)]*\((loriemann[\([\)\(i, k, l, m\)\(]\)])\)*\((Sup[\([\)\(l, m\)\(]\)])\)\)\)\)\)), {j, 1, n}]; NN = 1/Sqrt[1 - \!\( \*UnderoverscriptBox[\(\[Sum]\), \(i = 1\), \(4\)]\( \*UnderoverscriptBox[\(\[Sum]\), \(k = 1\), \(4\)]\((glo[\([\)\(i, k\)\(]\)])\)*Wvec[\([\)\(i\)\(]\)]* Wvec[\([\)\(k\)\(]\)]\)\)]; {Vt, Vr, V\[Theta], V\[Phi]} = NN (Wvec + Pup); EOM = Flatten[ Join[{EOM1, EOM2, EOM3} /. r -> r[\[Tau]] /. \[Theta] -> \[Theta][\[Tau]] /. Derivative[1][r[\[Tau]]][\[Tau]] -> Derivative[1][r][\[Tau]] /. Derivative[1][\[Theta][\[Tau]]][\[Tau]] -> Derivative[1][\[Theta]][\[Tau]]]]; INT1 = {t[0] == 0, r[0] == r0, \[Theta][0] == \[Theta]0, \[Phi][0] == 0}; INT2 = {Pt[0] == 1.32288, Pr[0] == 0, P\[Theta][0] == 0, P\[Phi][0] == 0.07143}; INT3 = {{Stt[0] == 0, Str[0] == 0, St\[Theta][0] == 0, St\[Phi][0] == 0}, {Srt[0] == 0, Srr[0] == 0, Sr\[Theta][0] == 0, Sr\[Phi][0] == 0}, {S\[Theta]t[0] == 0, S\[Theta]r[0] == 0, S\[Theta]\[Theta][0] == 0, S\[Theta]\[Phi][0] == 0}, {S\[Phi]t[0] == 0, S\[Phi]r[0] == 0, S\[Phi]\[Theta][0] == 0, S\[Phi]\[Phi][0] == 0}}; INT = Flatten[Join[{INT1, INT2, INT3}]]; r0 = 7; \[Theta]0 = Pi/2; \[Mu] = 1; aa = 0.5; M = 1; NDSolve[Flatten[Join[{EOM, INT}]], {t, r, \[Theta], \[Phi], Pt, Pr, P\[Theta], P\[Phi], Stt, Str, St\[Theta], St\[Phi], Srt, Srr, Sr\[Theta], Sr\[Phi], S\[Theta]t, S\[Theta]r, S\[Theta]\[Theta], S\[Theta]\[Phi], S\[Phi]t, S\[Phi]r, S\[Phi]\[Theta], S\[Phi]\[Phi]}, {\[Tau], 0, 1000}]
