$VersionNumber

Alternate ways to verify the Tomimatsu-Sato solutions for nTS = 1, 2, 3
		(Progr. Theor. Phys. 50, No 1, p 95, 1973)

<< EDCRGTCcode.m

Basic Definitions

nTS is used instead of δ. The scale factor m in the equation relating cylindrical and prolate spheroidal coordinates has been set equal to δ/p.

xCoord = {t, φ, x, y};

coFR = {d[t] - 
            2  nTS q/p (1 - y^2) Cf[x, y]/Af[x, y] d[φ], 
          d[φ], d[x], d[y]};

gTS = 
        DiagonalMatrix[{Af[x, y]/
              Bf[x, y], -Bf[x, y]/
                Af[x, y] (x^2 - 1) (1 - y^2), -Bf[x, y]/
                Kf[x, y]/(x^2 - 1), -Bf[x, y]/
                Kf[x, y]/(1 - y^2)}];

Definitions of Af, Bf, Cf, Kf for nTS = 1, 2, 3

Kf[n_] := p^(2  n) (x^2 - y^2)^(n^2 - 1)

\(
\)

Af[1] = p^2  x^2 - 1 + q^2  y^2;
Bf[1] = (p x + 1)^2 + q^2  y^2;
Cf[1] = -(p x + 1);

\(
\)

$x2 = (x^2 - 1); $y2 = (1 - y^2);

Af[2] = 
        p^4  $x2^4 + q^4  $y2^4 - 
          2  p^2  q^2  $x2 $y2 (3  $x2 $y2 + 2  $x2^2 + 
                2  $y2^2);

Bf[
          2] = (p^2 (x^4 - 1) + q^2 (y^4 - 1) + 2  p x $x2)^2 + 
          4  q^2  y^2 (p x $x2 + (p x + 1) $y2)^2;

Cf[2] = -p^3 x $x2 (2 (x^4 - 1) + (x^2 + 3) $y2) - 
          p^2  $x2 (4  x^2  $x2 + (3  x^2 + 1) $y2) + 
          q^2 (p x + 1) $y2^3;

\(
\)

Af[3] = 
        p^6  $x2^9 - q^6  $y2^9 - 
          3  p^4  q^2  $x2^4  $y2 (3  $x2^4 + 12  $x2^3  $y2 + 
                28  $x2^2  $y2^2 + 30  $x2 $y2^3 + 12  $y2^4) + 
          3  p^2  q^4  $x2 $y2^4 (3  $y2^4 + 12  $y2^3  $x2 + 
                28  $x2^2  $y2^2 + 30  $y2 $x2^3 + 
                12  $x2^4);

Bf[
          3] = (p (x^3 + 3  x) $x2^3 - 
                p q^2 (x^3 + 3  x y^2) ($x2 + $y2)^3 + 
                p^2 (3  x^2 + 1) $x2^3 - 
                q^2 (3  y^2 + 1) $y2^3)^2 + 
          q^2  y^2 (p^2 (3  x^2 + y^2) ($x2 + $y2)^3 - (y^2 + 
                        3) $y2^3 + 
                  12  p x ($x2 + $y2) ($x2*$y2))^2;

Cf[
          3] = -p^5 
            x $x2^4 ($x2^2 (3  x^2 + 1) (x^2 + 3) + 
                4  $x2 $y2 (x^4 + 8  x^2 + 3) + 
                2  $y2^2 (x^4 + 10  x^2 + 5)) + 
          p^3  q^2  
            x $x2 $y2^3 (8  $x2^2 (3  x^4 - 4  x^2 - 3) + 
                4  $x2 $y2 (13  x^4 - 10  x^2 - 15) + 
                4  $y2^2 (9  x^4 - 13) + 6  $y2^3 (x^2 + 3)) - 
          p^4  $x2^4 ($x2^2 (3  x^2 + 1)^2 + 
                2  $x2 $y2 (9  x^4 + 14  x^2 + 1) + 
                2  $y2^2 (5  x^4 + 10  x^2 + 1)) + 
          p^2  q^2  $x2 $y2^3 (54  $x2^3  $y2 + 
                4  $x2^2 (9  x^4 - 18  x^2 + 1) + 
                4  $x2 $y2 (9  x^4 - 8  x^2 - 13) + 
                6  $y2^3 (3  x^2 + 1) + 
                2  $y2^2 (37  x^4 - 42  x^2 - 3)) - 
          q^4 (p x + 1) $y2^8;

\(
\)

TSsol[n_] := {nTS -> n, H_[x, y] :> H[n]}

d[p] = 0; d[q] = 0;

p and q satisfy p^2+q^2=1. So use one of the following in simpRules:

pTOq = {p^n_ :> 
            p^Mod[n, 2] (1 - q^2)^Floor[n/2] /; n > 1};

qTOp = {q^n_ :> 
            q^Mod[n, 2] (1 - p^2)^Floor[n/2] /; n > 1};

simpRules = qTOp

If one uses pTOq instead, the calculation takes much longer!

\(
\)

First Method

To verify the Tomimatsu Sato solutions one can first calculate the Ricci tensor for arbitrary functions Af, Bf, Cf, Kf:

RGtensors[FacSimp[gTS], xCoord, FacSimp[coFR], {0, 0}]

ByteCount /@ Rdd

and then use FuncRepRules to replace the functions Af, Bf, Cf, Kf  in terms of the polynomials defined above:

ReplaceACKRules[n_] := 
        Join[FuncRepRules[Af[x, y], Af[n]], FuncRepRules[Cf[x, y], Cf[n]], 
          FuncRepRules[Kf[x, y], Kf[n]], {nTS -> n}];

It is faster if one does not replace all functions at once, but leave one (say Bf[x,y]), and first group terms with respect to Bf and its derivatives:

Rdd1 = 
          Map[Collect[#, {Derivative[__][Bf][x, y], Bf[x, y]}, 
                Factor] &, Rdd, 2]; // Timing

\(
\)

Verification of TS solution n=2:

In the Rdd2 calculation, using FacSimp ensures that the simpRules are applied. Finally, one can replace Bf[x,y]:

Rdd2=Map[Collect[#, {Derivative[__][Bf][x, y], 
                    Bf[x, y]}, FacSimp] &,Rdd1/.ReplaceACKRules[2],2];//Timing

FacSimp[Rdd2 /. FuncRepRules[Bf[x, y], Bf[2]]] // Timing

This way we verify the nTS=2 solution in 1/2 the time of the direct approach. Replacing all functions at once, however, takes much longer (~300 sec):

FacSimp[
        Rdd /. Join[ReplaceACKRules[2], FuncRepRules[Bf[x, y], Bf[2]]]] // 
      Timing

MaxMemoryUsed[]

\(
\)

Verification of TS solution n=3:

Rdd3=Map[Collect[#, {Derivative[__][Bf][x, y], 
                    Bf[x, y]}, FacSimp] &,Rdd1/.ReplaceACKRules[3],2];//Timing

FacSimp[Rdd3 /. FuncRepRules[Bf[x, y], Bf[3]]] // Timing

So, for the nTS=3 solution, this approach is only 10% faster. One could save some more time by examining the distinct Ricci components only (see below).

MaxMemoryUsed[]

\(
\)

Second Method

What if there was a mistake in typing the polynomials A[n], B[n] etc.? One would then wait for hours only to get messy expressions for the Ricci tensor. To avoid such a possibility when dealing with complicated metrics, one should test the expressions using approximations.

Thus, to verify the nTS=3 solution, we will expand each function in a series about some point x0 (or y0) to order n, and show that the Ricci tensor vanishes (to order n-2 because of the factors x^2-1, 1-y^2 entering the metric) for any such expansion.

First we will simplify as much as possible the independent components of the Ricci tensor (drop overall factors):

indRdd = indepTerms[Flatten[Rdd]];

ByteCount /@ indRdd

Head /@ indRdd

Length /@ indRdd

Head /@ Map[Take[#, -1] &, indRdd]

Head /@ Map[Take[#, {-2}] &, indRdd]

indRdd = 
        Sort[indepTerms[Map[Take[#, -1] &, indRdd]], 
          ByteCount[#1] < ByteCount[#2] &];

Head /@ indRdd

ByteCount /@ indRdd

Next, remove powers of q higher than 1 from Af, Bf, Cf, and define a general procedure for generating series replacement rules

Af[3] = FacSimp[Af[3] //. qTOp]; Bf[3] = FacSimp[Bf[3] //. qTOp]; 
    Cf[3] = FacSimp[Cf[3] //. qTOp];

\(\(ReplaceAllSeriesRules[x1_, x0_, 
          n_] := \(ReplaceAllSeriesRules[x1, x0, n] = 
          FacSimp[
		Join[
              FuncRepRules[Af[x, y], Series[Af[3], {x1, x0, n}]], 
              FuncRepRules[Bf[x, y], Series[Bf[3], {x1, x0, n}]], 
              FuncRepRules[Cf[x, y], Series[Cf[3], {x1, x0, n}]], 
              FuncRepRules[Kf[x, y], 
                Series[Kf[3], {x1, x0, n}]], {nTS -> 
                  3}]]\);\)\)

Warning: There are significant differences between different versions of Mathematica in doing the following calculations

FacSimp[indRdd /. ReplaceAllSeriesRules[y, 0, 3]] // Timing

The above calculation takes a few seconds with Mathematica 4! The following were evaluated using Mathematica 4.0. (Mathematica 3 and, remarkably, Mathematica 5 can take more than one hour to complete some of these substitutions!!!)

If[$VersionNumber > 3.5, 
      FacSimp[indRdd /. ReplaceAllSeriesRules[x, 0, 3]] // Timing]

If[$VersionNumber > 3.5 && $VersionNumber < 5, 
      FacSimp[indRdd /. ReplaceAllSeriesRules[x, 1/2, 3]] // 
        Timing]

If[$VersionNumber > 3.5 && $VersionNumber < 5, 
      FacSimp[indRdd /. ReplaceAllSeriesRules[x, 1, 3]] // Timing]

If[$VersionNumber > 3.5, 
      FacSimp[indRdd /. ReplaceAllSeriesRules[y, 0, 4]] // Timing]

If[$VersionNumber > 3.5 && $VersionNumber < 5, 
      FacSimp[indRdd /. ReplaceAllSeriesRules[y, 1/2, 4]] // 
        Timing]

If[$VersionNumber > 3.5 && $VersionNumber < 5, 
      FacSimp[indRdd /. ReplaceAllSeriesRules[y, -1, 5]] // 
        Timing]

One can also expand with respect to the parameter p:

If[$VersionNumber > 3.5, 
      FacSimp[indRdd /. ReplaceAllSeriesRules[p, 0, 4]] // Timing]

If[$VersionNumber > 3.5 && $VersionNumber < 5, 
      FacSimp[indRdd /. ReplaceAllSeriesRules[p, 1, 3]] // Timing]

If[$VersionNumber > 3.5 && $VersionNumber < 5, 
      FacSimp[indRdd /. ReplaceAllSeriesRules[p, 1, 4]] // Timing]

If[$VersionNumber > 3.5 && $VersionNumber < 5, 
      FacSimp[indRdd /. ReplaceAllSeriesRules[p, 1/3, 3]] // 
        Timing]

\(
\)

It is faster, instead of making substitutions, to evaluate RGtensors with the approximate frame and metric:

Clear$dx

gTSser = FacSimp[gTS /. ReplaceAllSeriesRules[y, 0, 3]]

coFRser = FacSimp[coFR /. ReplaceAllSeriesRules[y, 0, 3]]

RGtensors[gTSser, xCoord, coFRser, {0, 0}]

Rdd

Same calculation with one order higher:

Clear$dx

gTSser = FacSimp[gTS /. ReplaceAllSeriesRules[y, 0, 4]]

coFRser = FacSimp[coFR /. ReplaceAllSeriesRules[y, 0, 4]]

RGtensors[gTSser, xCoord, coFRser, {0, 0}]

Rdd