7

As far as I know spherical gravitational collapse -of central importance to theoretical physics- is thought to be out of the scope of Wolfram Mathematica.

However there are cases, as in this article, where numerical integration with NDSolve seems possible in principle.

If this where to happen Mathematica would prove to have many more applications than is commonly thought.

BOUNTY: Is it possible to reproduce the upper subplot FIG.5 in this article? If not: Can one use eq.5 (or any transformation of it) along with any initial-boundary conditions to produce a plot analogous to FIG.5?

P.S. This question has been modified so as to be more flexible. No substantial change has been done. To be more specific:

  1. In the APPENDIX is posted my first approach towards the problem which has been improved thanks to user21, Alex Trunev and xzczd. It was part of this question. However I feel there may be better approaches to this problem. For example FEM may implement Neumann b.c. automatically. Also one may like to perform an explicit discretization of the PDE system. In any case if anyone feels that my answer should be used as a basis, or contains usefull information there is no problem.

  2. In the initial question I was interested only in the first part of the question: the upper subplot FIG.5 (i.e. the massless-Dirichlet case). However this proved to be quite difficult. So I added the if not part of the question, just in case the first part is impossible: just solve eq.5 with any ic-bc.

P.S.2 I am relatively new to this forum, so any suggestions on how I could imrprove this question are welcome.

Important note: For the massles case $\mu$ in eq.5 is set to $0$!

APPENDIX Here is my own attempt just in case it contains usefull information.

First I costruct initial data according to the Dirichlet-massless case of FIG.5 

    A = 0.04; w = 0.125;
    Pin[r_] := A*Exp[-r^2/w^2] 


    PDE0 = D[u[r], r, r] + 2*D[u[r], r]/r == -Pi*Pin[r]^2*(1 + u[r])^5; 
    (*eqs 23, 24a*)

    rmin = 10^(-40); (* as close to r=0 as possible *)
    BC0 = {u'[rmin] == 0, u'[1] == -u[1]};(*below eq 23*)
    initial = NDSolve[{PDE0, BC0}, {u}, {r, rmin, 1}];


    yin[r_] := First[1 + u[r] /. initial](*since ψ=1+u*)
    ain[r_] := 1
    Fin[r_] := 0
    kin[r_] := 0


    IC = {k[0, r] == kin[r], F[0, r] == Fin[r], a[0, r] == ain[r], 
       P[0, r] == Pin[r], y[0, r] == yin[r]};

Then according to eq.10 I define the following Dirichlet b.c. 

BC1 = {F[t, 1] == 0, P[t, 1] == 0, k[t, 1] == 0, a[t, 1] == 1};

which lead to (see below eq.10) the extra b.c.

BC2 = {Derivative[0, 1][k][t, 1] == 0, Derivative[0, 1][P][t, 1] == 0,
    Derivative[0, 1][a][t, 1] == ain'[1], y[t, 1] == yin[1], 
   Derivative[0, 1][y][t, 1] == yin'[1]};

Regularity at the center is ensured by implementing Neumann b.c. 

BCreg = {Derivative[0, 1][F][t, rmin] == 0, 
   Derivative[0, 1][P][t, rmin] == 0, 
   Derivative[0, 1][k][t, rmin] == 0, 
   Derivative[0, 1][a][t, rmin] == 0, 
   Derivative[0, 1][y][t, rmin] == 0};

Then according to eq.5 (massless case:$\mu=0$) I define the PDE system that determines the evolution of initial data.

    eqy = D[y[t, r], t] == -a[t, r]*y[t, r]*k[t, r]/6;

    eqk = D[k[t, r], 
    t] == -(1/y[t, r]^4)*(D[a[t, r], r, r] + 2*D[a[t, r], r]/r) - 
    2*D[y[t, r], r]*
     D[a[t, r], r]/y[t, r]^5  + (a[t, r]*k[t, r]^2/3) + (8*Pi*
      P[t, r]^2*a[t, r]);

    eqF = D[F[t, r], t] == -a[t, r]*P[t, r];

    eqP = D[P[t, r], t] == 
    a[t, r]*P[t, r]*
     k[t, r] - (a[t, r]/y[t, r]^4)*(D[F[t, r], r, r] + 
       2*D[F[t, r], r]/r) - D[a[t, r], r]*D[F[t, r], r]/y[t, r]^4 - 
    2*a[t, r]*D[y[t, r], r]*D[F[t, r], r]/y[t, r]^5;

    eqa = D[a[t, r], t] == -2*a[t, r]*k[t, r];

    PDEs = {eqy, eqk, eqF, eqP, eqa};

Finally I integrate with

   tmax = 10^(4)
   evolution = 
   NDSolveValue[{PDEs, IC, BC1, BC2, BCreg}, {y, k, F, P, a}, {t, 0, 
   tmax}, {r, rmin, 1}, "DependentVariables" -> {y, k, F, P, a}, 
   Method ->{"MethodOfLines", TemporalVariable -> t}]
  • 1
    In comparison with the article, you have several errors in the equations. – Alex Trounev Sep 10 '19 at 17:55
  • 2
    I suggest you add a Method->"MethodOfLines" to your NDSolve command and start fixing the errors shown. You can get help on this by clicking on the three ... in front of the NDSolve error message. – user21 Sep 11 '19 at 04:59
  • @AlexTrounev thanks for your response. I checked the PDE system but found no mistakes. Probably I am mistaken. Could it be that I have set $μ=0$ as should be for the massless case? –  Sep 11 '19 at 07:51
  • @user21 adding Method->"MethodOfLines" did help. New error ivone occured which I am trying to resolve. –  Sep 11 '19 at 08:01
  • Are you sure the i.c.s and b.c.s are correctly given? Related: https://math.stackexchange.com/q/450367/58219 “But this is a first-roder system in time and I have given one IC per independent variable.” No, there's a D[a[t, r], t, t] in eqk. – xzczd Sep 11 '19 at 08:15
  • @AlexTrounev,@xzczd you're right there are mistakes in comparison to the article. I wil need a few time to fix it. –  Sep 11 '19 at 08:56
  • @jheidk51 I found 3 typos: 1) D[a[t, r], t, t] + 2*D[a[t, r], r]/r should be D[a[t, r],r,r] + 2*D[a[t, r], r]/r;2)D[F[t, r], t, t] + 2*D[F[t, r], r]/r) should be D[F[t, r],r,r] + 2*D[F[t, r], r]/r;3) 8*Pi*P[t, r]^3*a[t, r] should be 8*Pi*P[t, r]^2*a[t, r]. But I see that you have already corrected two typos. Thank you. – Alex Trounev Sep 11 '19 at 09:01
  • 1
    You have changed your question substantially by removing important information. Please do not make it a habit that you expect people to go and read papers in an effort to answer questions. – user21 Sep 12 '19 at 05:53
  • @user21 I did not have any intention to change my question substantially, but I thought that what I was asking for was too difficult. So tried to be more flexible (i.e. reproduction of the massless-Dirichlet case of FIG.5 is still what I need, but just in case this is impossible, there is an if not alternative). –  Sep 12 '19 at 13:22
  • 2
    Why is this particular problem considered outside the scope of mathematicas ability? – DrMrstheMonarch Sep 14 '19 at 18:43
  • Because gravitational collapse PDE's system in addition to evolution equations usually contain elliptic PDEs with respect to space coordinate r. One may then decide to discretize space by hand, but then there may be a crucial loss of accuracy, given the nonlinear character of the equations. –  Sep 17 '19 at 15:37

1 Answers1

1

I changed the equations in accordance with the article (added terms with mass) and left as many boundary conditions as needed to solve the problem, namely: for first-order equations in r for one boundary condition, for second-order equations for two boundary conditions. The authors of the article write that they have artificial viscosity there. Apparently for this reason they used two boundary conditions in each equation. Without artificial viscosity, up to t = 3 can be calculated.

A = 4/100; w = 125/1000;
Pin[r_] := A*Exp[-r^2/w^2]


PDE0 = D[u[r], r, r] + 2*D[u[r], r]/r == -Pi*Pin[r]^2*(1 + u[r])^5;
(*eqs 23,24a*)

rmin = 10^(-30);(*as close to r=0 as possible*)BC0 = {u'[rmin] == 0, 
  u'[1] == -u[1]};(*below eq 23*){initial, initial1} = 
 NDSolveValue[{PDE0, BC0}, {u, u'}, {r, rmin, 1}, 
  WorkingPrecision -> 30];


yin[r_] := 1 + initial[r](*since \[Psi]=1+u*)
ain[r_] := 1
Fin[r_] := 0
kin[r_] := 0

{Plot[initial[r], {r, rmin, 1}], Plot[initial1[r], {r, rmin, 1}]}


mu = 0;

rmin = 10^-3; IC = {k[0, r] == kin[r], F[0, r] == Fin[r], 
  a[0, r] == ain[r], P[0, r] == Pin[r], y[0, r] == yin[r]};
BC1 = {F[t, 1] == 0, P[t, 1] == 0, k[t, 1] == 0, a[t, 1] == 1};
BC2 = {Derivative[0, 1][a][t, 1] == 0, 
  y[t, 1] == yin[1]}; BCreg = {Derivative[0, 1][F][t, rmin] == 0, 
  Derivative[0, 1][a][t, rmin] == 0};
eqy = D[y[t, r], t] == (-a[t, r]*y[t, r]*k[t, r]/6);

eqk = D[k[t, r], 
    t] == (-(1/y[t, r]^4)*(D[a[t, r], r, r] + 2*D[a[t, r], r]/r) - 
     2*D[y[t, r], r]*D[a[t, r], r]/y[t, r]^5 + (a[t, r]*k[t, r]^2)/
      3 + 4*Pi*a[t, r] (2 P[t, r]^2 - mu^2 F[t, r]^2));

eqF = D[F[t, r], t] == (-a[t, r]*P[t, r]);

eqP = D[P[t, r], 
    t] == (a[t, r]*P[t, r]*
      k[t, r] - (a[t, r]/y[t, r]^4)*(D[F[t, r], r, r] + 
        2*D[F[t, r], r]/r) - D[a[t, r], r]*D[F[t, r], r]/y[t, r]^4 - 
     2*a[t, r]*D[y[t, r], r]*D[F[t, r], r]/y[t, r]^5 + 
     mu^2 a[t, r] F[t, r]);

eqa = D[a[t, r], t] == (-2*a[t, r]*k[t, r]);

PDEs = {eqy, eqk, eqF, eqP, eqa};
tmax = 3;
evolution = 
  NDSolveValue[{PDEs, IC, BC1, BC2, BCreg}, {y, k, F, P, a}, {t, 0, 
    tmax}, {r, rmin, 1}, 
   Method -> {"MethodOfLines", 
     "SpatialDiscretization" -> {"TensorProductGrid", 
       "MinPoints" -> 40, "MaxPoints" -> 100, 
       "DifferenceOrder" -> "Pseudospectral"}}, MaxSteps -> 10^6];

It can be seen from Fig. 1 that even at the very beginning of evolution, characteristic oscillations appeared. In this example, artificial viscosity is not yet used, and the mass $\mu = 0$ Figure 1

I managed to build a stable code for calculating evolution over a period of $tmax=10^3-10^4$. For this, I added two more equations to the system that describe the conservation of the Hamiltonian h[t,r] and the momentum m[t,r] (eq (6)-(7) in the paper). In addition, I added artificial viscosity (not to all equations) and the equation for calculating the scalar, which is shown in Fig. 5 (Kretschmann scalar). The result is such a code

mu = 4; {av1, av2, av3, av4, av5, av6, 
  av7} = {0, 1, 1, 0, 1, 1, 1} 10^-3; nn = 999;
rmin = 1/nn; IC = {k[0, r] == kin[r], F[0, r] == Fin[r], 
  a[0, r] == ain[r], P[0, r] == Pin[r], y[0, r] == yin[r], 
  h[0, r] == 0, m[0, r] == 0};
BC1 = {F[t, 1] == 0, P[t, 1] == 0, k[t, 1] == 0, a[t, 1] == 1, 
   y[t, 1] == yin[1], h[t, 1] == 0, m[t, 1] == 0};
BCreg = {Derivative[0, 1][F][t, rmin] == 0, 
  Derivative[0, 1][a][t, rmin] == 0, 
  Derivative[0, 1][k][t, rmin] == 0, h[t, rmin] == 0, m[t, rmin] == 0};
eqy = D[y[t, r], t] == (-a[t, r]*y[t, r]*k[t, r]/6) + 
    av1 D[y[t, r], r, r];

eqk = D[k[t, r], 
    t] == (-(1/y[t, r]^4)*(D[a[t, r], r, r] + 2*D[a[t, r], r]/r) - 
      2*D[y[t, r], r]*D[a[t, r], r]/y[t, r]^5 + (a[t, r]*k[t, r]^2)/
       3 + 4*Pi*a[t, r] (2 P[t, r]^2 - mu^2 F[t, r]^2)) + 
    av2 D[k[t, r], r, r];

eqF = D[F[t, r], t] == (-a[t, r]*P[t, r]) + av3 D[F[t, r], r, r];

eqP = D[P[t, r], 
    t] == (a[t, r]*P[t, r]*
       k[t, r] - (a[t, r]/y[t, r]^4)*(D[F[t, r], r, r] + 
         2*D[F[t, r], r]/r) - D[a[t, r], r]*D[F[t, r], r]/y[t, r]^4 - 
      2*a[t, r]*D[y[t, r], r]*D[F[t, r], r]/y[t, r]^5 + 
      mu^2 a[t, r] F[t, r]) + av4 D[P[t, r], r, r];

eqa = D[a[t, r], t] == (-2*a[t, r]*k[t, r]) + av5 D[a[t, r], r, r];
eqh = D[h[t, r], 
    t] == ((D[y[t, r], r, r] + 2/r D[y[t, r], r])/y[t, r]^5 - 
      k[t, r]^2/12 + 
      Pi (P[t, r]^2 + D[F[t, r], r]^2/y[t, r]^4 + mu^2 F[t, r]^2)) + 
    av6 D[h[t, r], r, r];
eqm = D[m[t, r], 
    t] == (2/3 D[k[t, r], r] + 8 Pi P[t, r] D[F[t, r], r]) + 
    av7 D[m[t, r], r, r];

PDEs = {eqy, eqk, eqF, eqP, eqa, eqh, eqm};
tmax = 1000;
evolution = 
  NDSolveValue[{PDEs, IC, BC1, BCreg}, {y, k, F, P, a, h, m}, {t, 0, 
    tmax}, {r, rmin, 1}, 
   Method -> {"MethodOfLines", 
     "SpatialDiscretization" -> {"TensorProductGrid", 
       "MinPoints" -> nn, "MaxPoints" -> nn, "DifferenceOrder" -> 2}},
    MaxSteps -> 10^6];

lb = {y, k, F, P, a, h, m};

Table[Plot3D[evolution[[i]][t, r], {t, 0, tmax}, {r, rmin, 1}, 
  Mesh -> None, ColorFunction -> "Rainbow", AxesLabel -> Automatic, 
  PlotLabel -> lb[[i]], PlotRange -> All], {i, 1, 7}]
(*Kretschmann scalar*)
ks = (2/27 (k[t, r]^4 - 
        24 Pi k[t, r]^2 (P[t, r]^2 + mu^2 F[t, r]^2)) + 
     8 D[a[t, r], r, r]^2/(3 a[t, r]^2 y[t, r]^8) + 
     8/3 (4 Pi^2 (11 P[t, r]^4 - 2 mu^2 P[t, r]^2 F[t, r]^2 + 
          5 mu^4 F[t, r]^4))) /. 
   Flatten[Table[lb[[i]] -> evolution[[i]], {i, 1, 5}]];

Plot[ks /. r -> rmin, {t, 0, tmax}, PlotRange -> All]

Figure 2 shows the results for $\mu = 4$. It can be seen that nonlinear oscillations are observed only at the very beginning of evolution. Moreover, on small grids with nn = 3200, these oscillations disappear altogether. Figure 2

There is another solution method at large time intervals. Here I did not include the Hamiltonian and momentum in the system of equations and put $\mu =0$. In this case, self-oscillations also occur at t< 100, even at nn = 3200 (this number was used in the construction of figure 5).

A = 4/100; w = 125/1000;
Pin[r_] := A*Exp[-r^2/w^2]


PDE0 = D[u[r], r, r] + 2*D[u[r], r]/r == -Pi*Pin[r]^2*(1 + u[r])^5;
(*eqs 23,24a*)

rmin = 10^(-30);(*as close to r=0 as possible*)BC0 = {u'[rmin] == 0, 
  u'[1] == -u[1]};(*below eq 23*){initial, initial1} = 
 NDSolveValue[{PDE0, BC0}, {u, u'}, {r, rmin, 1}, 
  WorkingPrecision -> 30];


yin[r_] := 1 + initial[r](*since \[Psi]=1+u*)
ain[r_] := 1
Fin[r_] := 0
kin[r_] := 0

{Plot[initial[r], {r, rmin, 1}], Plot[initial1[r], {r, rmin, 1}]}


mu = 0; {av1, av2, av3, av4, av5, av6, 
  av7} = {0, 1, 1, 0, 1, 1, 1} 10^-3; nn = 3200;
rmin = 1/nn; IC = {k[0, r] == kin[r], F[0, r] == Fin[r], 
  a[0, r] == ain[r], P[0, r] == Pin[r], y[0, r] == yin[r]};
BC1 = {F[t, 1] == 0, P[t, 1] == 0, k[t, 1] == 0, a[t, 1] == 1, 
   y[t, 1] == yin[1]};
BC2 = {Derivative[0, 1][a][t, 1] == 0, Derivative[0, 1][y][t, 1] == 0,
   Derivative[0, 1][k][t, 1] == 0, Derivative[0, 1][P][t, 1] == 0, 
  Derivative[0, 1][F][t, 1] == 
   0}; BCreg = {Derivative[0, 1][F][t, rmin] == 0, 
  Derivative[0, 1][a][t, rmin] == 0, 
  Derivative[0, 1][k][t, rmin] == 0};
eqy = D[y[t, r], t] == (-a[t, r]*y[t, r]*k[t, r]/6) + 
    av1 D[y[t, r], r, r];

eqk = D[k[t, r], 
    t] == (-(1/y[t, r]^4)*(D[a[t, r], r, r] + 2*D[a[t, r], r]/r) - 
      2*D[y[t, r], r]*D[a[t, r], r]/y[t, r]^5 + (a[t, r]*k[t, r]^2)/
       3 + 4*Pi*a[t, r] (2 P[t, r]^2 - mu^2 F[t, r]^2)) + 
    av2 D[k[t, r], r, r];

eqF = D[F[t, r], t] == (-a[t, r]*P[t, r]) + av3 D[F[t, r], r, r];

eqP = D[P[t, r], 
    t] == (a[t, r]*P[t, r]*
       k[t, r] - (a[t, r]/y[t, r]^4)*(D[F[t, r], r, r] + 
         2*D[F[t, r], r]/r) - D[a[t, r], r]*D[F[t, r], r]/y[t, r]^4 - 
      2*a[t, r]*D[y[t, r], r]*D[F[t, r], r]/y[t, r]^5 + 
      mu^2 a[t, r] F[t, r]) + av4 D[P[t, r], r, r];

eqa = D[a[t, r], t] == (-2*a[t, r]*k[t, r]) + av5 D[a[t, r], r, r];
PDEs = {eqy, eqk, eqF, eqP, eqa};
tmax = 10000;
evolution = 
  NDSolveValue[{PDEs, IC, BC1, BCreg}, {y, k, F, P, a}, {t, 0, 
    tmax}, {r, rmin, 1}, 
   Method -> {"MethodOfLines", 
     "SpatialDiscretization" -> {"TensorProductGrid", 
       "MinPoints" -> nn, "MaxPoints" -> nn, "DifferenceOrder" -> 4}},
    MaxSteps -> 10^6];
lb = {y, k, F, P, a};

Table[Plot3D[evolution[[i]][t, r], {t, 0, tmax}, {r, rmin, 1}, 
  Mesh -> None, ColorFunction -> "Rainbow", AxesLabel -> Automatic, 
  PlotLabel -> lb[[i]], PlotRange -> All], {i, 1, 5}]

(*Kretschmann scalar*)
ks = (2/27 (k[t, r]^4 - 
        24 Pi k[t, r]^2 (P[t, r]^2 + mu^2 F[t, r]^2)) + 
     8 D[a[t, r], r, r]^2/(3 a[t, r]^2 y[t, r]^8) + 
     8/3 (4 Pi^2 (11 P[t, r]^4 - 2 mu^2 P[t, r]^2 F[t, r]^2 + 
          5 mu^4 F[t, r]^4))) /. 
   Flatten[Table[lb[[i]] -> evolution[[i]], {i, 1, 5}]];

LogLogPlot[ks /. r -> rmin, {t, 0, tmax}, PlotRange -> All, 
 PlotLabel -> "Kretschmann scalar", AxesLabel -> Automatic]

Figure 3

Alex Trounev
  • 44,369
  • 3
  • 48
  • 106
  • Alex Trounev how do rmin = 10^(-30);(*as clos and rmin = 10^-3; IC = { work together? –  Sep 11 '19 at 10:53
  • There is no point in using rmin=10^-30 in PDE, since the instability develops at the boundary r=1, and not at the boundary r=0. In any case, this is working code and you can check all the options including rmin=10^-30.Please note that I have removed the extra boundary conditions. – Alex Trounev Sep 11 '19 at 10:54
  • Alex Trounev you're right. But then how do the authors actualy implement Neumann boundary conditions for k,y,P? Will the result be the same if I omit these conditions? –  Sep 11 '19 at 11:16
  • 1
    I don’t think that they understand there to the end what they write. They, for example, do not describe the method of numerical solution. It is possible that they add artificial viscosity. Then the number of boundary conditions must correspond to five parabolic equations. – Alex Trounev Sep 11 '19 at 11:27
  • Does this mean that results stated in the article cannot be addressed with mathematica? –  Sep 11 '19 at 12:32
  • 3
    @jheidk51 To be precise, the results cannot be addressed with any software or program without knowing the exact method they use. As to the boundary part, you can still force the b.c. at $r=0$ to be the given ones using the method in e.g. https://mathematica.stackexchange.com/a/133140/1871 but I doubt if it's the right way to go, because this'll probably lead to something like shock wave in my limited experience. – xzczd Sep 11 '19 at 13:26