2
$\begingroup$

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}]
$\endgroup$
4
  • $\begingroup$ Not answer, but the suit of packages xAct might be useful to you: xact.es $\endgroup$
    – mattiav27
    Oct 31, 2020 at 14:15
  • $\begingroup$ @mattiav27 How can I use this package? $\endgroup$
    – MMS
    Oct 31, 2020 at 20:11
  • $\begingroup$ Download the package from xact.es/download.html and follow the intruction xact.es/download/install $\endgroup$
    – mattiav27
    Oct 31, 2020 at 20:21
  • 3
    $\begingroup$ 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.) $\endgroup$
    – bbgodfrey
    Nov 1, 2020 at 0:09

1 Answer 1

3
$\begingroup$

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}]

enter image description here

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.} *)
$\endgroup$
1
  • $\begingroup$ Thanks @bbgodfrey for your answer, it was really helpful. $\endgroup$
    – MMS
    Nov 1, 2020 at 8:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.