Mathematica notebook implementing the Ricci flow and Newton method algorithm discussed in:

Chapter 10 of the book `Black Holes in Higher Dimensions' to be published by Cambridge University Press (editor: G. Horowitz)  [ arXiv:1107.5513 [gr-qc] ]

Author: Toby Wiseman (Imperial College)

   [Copyright Toby Wiseman, Theoretical Physics, Imperial College, Jan 2012]
   

Please feel free to email me if you have found this notebook useful, if you spot any errors in it or if you have any comments on it;   t.wiseman _at_ imperial.ac.uk

This notebook is divided into 4 sections;

	Sec 1)  Computes the continuum Ricci-DeTurck flow equations for the metric components for the ansatz in eqn 10.33 of chapter 10. This notebook uses simple tensor manipulation routines that I have written, and many more powerful packages for Mathematica are available to do these calculations [eg. GRTensorM, xTensor ...].
	
	Sec 2)  Following the discretization scheme outlined in eqns 10.34 & 10.35 the discrete vector of equations in eqn 10.37 is constructed.
	
	Sec 3)  Ricci-DeTurck flow is performed to various initial data.
	
	Sec 4)  Newton's method is applied to find the Schwarzschild solution.

Note:

The discretization resolution used in this notebook may be adjusted. By default it is very low, only 20 points, so that the notebook runs in a reasonable time. The user may change this to higher values, but Ricci flow in particular will run rather slowly. 

I emphasize that this notebook is written to give a very simple implementation of the explicit example discussed in chapter 10  for the two algorithms, and is very far from optimized for speed. In particular it is written in an extremely inefficient way for the Ricci flow algorithm. 

For more details on implementation of these algorithms in more complicated settings and with greater efficiency please refer to the research articles:  hep-th/0606086,  arXiv:0905.1822,  arXiv:1104.4489

Section 1

Code to generate Ricci-DeTurck flow equations

Required routines:

Please note that I have written these routines myself, but one could quite simply use other mathematica packages to perform the same manipulations.

BuildMetric:=(Gu=FullSimplify[Inverse[Gd]];G={{u,d},IdentityMatrix[Nind]};)

getind[i1_]:=(If[IntegerQ[i1],i1,Position[ind,i1][[1,1]]])

lstind[i1_]:=(lsttmp2={};Do[AppendTo[lsttmp2,getind[i1[[lsttmp1,1]]]],{lsttmp1,1,Length[i1]}];lsttmp2)

up[i1_,i2_]:=(lsttmp2={};Do[AppendTo[lsttmp2,If[i1[[lsttmp1,2]]=!=i2[[lsttmp1]],If[i1[[lsttmp1,2]]===d,d,u],0]],{lsttmp1,1,Length[i1]}];lsttmp2)

MT[t1_,dim_,indpos_]:=(MTtmp1=0;If[Length[indpos]=!=dim,Print["Error in index position"]];Do[MTtmp1=Table[MTtmp1,{MTtmp3,1,Nind}],{MTtmp2,1,dim}];t1={indpos,MTtmp1})

ST[t1_,tind_,val_]:=(STtmp1=t1;If[up[tind,t1[[1]]]=!=Table[0,{Length[tind]}],Print["Error in Index Assignment!!"];];STtmp1[[2]]=ReplacePart[t1[[2]],val,lstind[tind]];STtmp1)

GT[t1_,tind_]:=(GTtmp2=up[tind,t1[[1]]];GTtmp3="";GTtmp4="";GTtmp5="";GTtmp6=t1;
		Do[
			Switch[GTtmp2[[GTtmp1]],
				0,
				(GTtmp3=GTtmp3<>","<>ToString[getind[tind[[GTtmp1,1]]]];),
				u,
				(GTtmp3=GTtmp3<>","<>"tmpind"<>ToString[GTtmp1];
					GTtmp4=GTtmp4<>" Gu[["<>ToString[getind[tind[[GTtmp1,1]]]]<>", tmpind"<>ToString[GTtmp1]<>"]] ";
					GTtmp5=GTtmp5<>",{tmpind"<>ToString[GTtmp1]<>",1,Nind}"),
				d,
				(GTtmp3=GTtmp3<>","<>"tmpind"<>ToString[GTtmp1];
					GTtmp4=GTtmp4<>" Gd[["<>ToString[getind[tind[[GTtmp1,1]]]]<>", tmpind"<>ToString[GTtmp1]<>"]] ";
					GTtmp5=GTtmp5<>",{tmpind"<>ToString[GTtmp1]<>",1,Nind}";)];
			,{GTtmp1,1,Length[t1[[1]]]}];If[GTtmp4===""&&GTtmp5==="",ToExpression["GTtmp6[[2"<>GTtmp3<>"]]"],ToExpression["Sum[GTtmp6[[2"<>GTtmp3<>"]]"<>GTtmp4<>GTtmp5<>"]"]])

Load Metric ansatz (Eqn 10.33):  (ds^2) = (r^2) A ((dτ)^2) + 4 (f^2) B (dr^2) + f C ((dΩ)^2)   where f = 1/(((1-(r^2)))^2)

Metric functions A, B, C are functions of the radial coordinate r and Ricci flow time λ.

f[r_]:=(1/((1-r^2)^2))

Gd=(({{r^2A[λ,r],   0, 0, 0}, {0, 4f[r]^2B[λ,r], 0, 0}, {0, 0, f[r]C[λ,r], 0}, {0, 0, 0, f[r]C[λ,r]Sin[θ]^2}}))

ind={τ,r,θ,ϕ}

Nind=Length[ind]

BuildMetric;

Calculate Christoffel symbol (Γ) and Ricci Tensor (R2)

Clear[Γ,R4,R2]
MT[Γ,3,{d,d,u}];
MT[R4,4,{d,d,d,u}];
MT[R2,2,{d,d}];

(* Compute the Christoffel symbol components *)

Do[Γ=ST[Γ,{{ind[[i1]],d},{ind[[i2]],d},{ind[[i3]],u}},Sum[(1/2)GT[G,{{ind[[i3]],u},{ind[[ss1]],u}}](D[GT[G,{{ind[[i1]],d},{ind[[ss1]],d}}],ind[[i2]]]+D[GT[G,{{ind[[ss1]],d},{ind[[i2]],d}}],ind[[i1]]]-D[GT[G,{{ind[[i1]],d},{ind[[i2]],d}}],ind[[ss1]]]),{ss1,1,Nind}]],{i1,1,Nind},{i2,1,Nind},{i3,1,Nind}]

Γ=FullSimplify[Γ];

(* Compute the Riemann tensor components *)

Do[
	R4=ST[R4,
			{{ind[[i1]],d},{ind[[i2]],d},{ind[[i3]],d},{ind[[i4]],u}},
	
			D[GT[Γ,{{ind[[i1]],d},{ind[[i3]],d},{ind[[i4]],u}}],ind[[i2]]]-D[GT[Γ,{{ind[[i2]],d},{ind[[i3]],d},{ind[[i4]],u}}],ind[[i1]]]+
			Sum[GT[Γ,{{ind[[i1]],d},{ind[[i3]],d},{ind[[ss1]],u}}]GT[Γ,{{ind[[i2]],d},{ind[[ss1]],d},{ind[[i4]],u}}]-GT[Γ,{{ind[[i2]],d},{ind[[i3]],d},{ind[[ss1]],u}}]GT[Γ,{{ind[[i1]],d},{ind[[ss1]],d},{ind[[i4]],u}}],{ss1,1,Nind}]
		
		]
		,
	{i1,1,Nind},{i2,1,Nind},{i3,1,Nind},{i4,1,Nind}]

(* Compute the Ricci tensor components *)

Do[R2=ST[R2,
			{{ind[[i1]],d},{ind[[i2]],d}},
			Sum[GT[R4,{{ind[[i1]],d},{ind[[ss1]],u},{ind[[i2]],d},{ind[[ss1]],d}}],{ss1,1,Nind}]
		]
	,{i1,1,Nind},{i2,1,Nind}]

R2=Simplify[R2];

Compute the Ricci-DeTurck flow equations, with flow time λ.

Choose the background metric to be:   ((OverBar[ds])^2) = (r^2) A0 ((dτ)^2) + 4 (f^2) B0 (dr^2) + f C0 ((dΩ)^2)

(* Compute the reference connection *)

Γ0=Γ/.{A->A0,B->B0,C->C0};

(* Compute the ξ vector *)

Clear[ξ,dΓ]

MT[dΓ,3,{d,d,u}];
MT[ξ,1,{d}];

Do[dΓ=ST[dΓ,{{ind[[i1]],d},{ind[[i2]],d},{ind[[i3]],u}},
GT[Γ,{{ind[[i1]],d},{ind[[i2]],d},{ind[[i3]],u}}]-GT[Γ0,{{ind[[i1]],d},{ind[[i2]],d},{ind[[i3]],u}}]
],{i1,1,Nind},{i2,1,Nind},{i3,1,Nind}];

dΓ=Simplify[dΓ];

Do[
ξ=ST[ξ,{{i1,d}},
 Sum[GT[dΓ,{{s1,d},{s1,u},{i1,d}}],{s1,1,Nind}]
];,
{i1,1,Nind}];

ξ=Simplify[ξ];

(* Compute ∇ξ *)

Clear[dξ]

MT[dξ,2,{d,d}];

Do[
dξ=ST[dξ,{{i1,d},{i2,d}},
(1/2)(D[GT[ξ,{{i1,d}}],ind[[i2]]]-Sum[GT[Γ,{{i1,d},{i2,d},{s1,u}}]GT[ξ,{{s1,d}}],{s1,1,Nind}]
+D[GT[ξ,{{i2,d}}],ind[[i1]]]-Sum[GT[Γ,{{i2,d},{i1,d},{s1,u}}]GT[ξ,{{s1,d}}],{s1,1,Nind}])
];,
{i1,1,Nind},{i2,1,Nind}];

dξ=Simplify[dξ];

Clear[RFeq]

MT[RFeq,2,{d,d}];

Do[
RFeq=ST[RFeq,{{i1,d},{i2,d}},
D[GT[G,{{i1,d},{i2,d}}],λ]+2( GT[R2,{{i1,d},{i2,d}}]- GT[dξ,{{i1,d},{i2,d}}] )
];,
{i1,1,Nind},{i2,1,Nind}];

RFsol=Simplify[Solve[RFeq[[2]]==0,{Derivative[1,0][A][λ,r],Derivative[1,0][B][λ,r],Derivative[1,0][C][λ,r]}]][[1]];

Collect the Ricci-DeTurck flow equations into the list RF

RF={Derivative[1,0][A][λ,r],Derivative[1,0][B][λ,r],Derivative[1,0][C][λ,r]}/.RFsol/.{A[λ,r]->A,B[λ,r]->B,C[λ,r]->C,Derivative[0,1][A][λ,r]->Ar,Derivative[0,1][B][λ,r]->Br,Derivative[0,1][C][λ,r]->Cr,Derivative[0,2][A][λ,r]->Arr,Derivative[0,2][B][λ,r]->Brr,Derivative[0,2][C][λ,r]->Crr}

Section 2

Code to discretize the flow equation

NR=20;

Note that by default 20 points are used to discretize the equation. This may be modified by appropriately adjusting the value above. Be warned that in this very simple Mathematica notebook taking too many points will cause very slow running, particularly for the Ricci flow which is implemented here in an extremely inefficient manner to ensure this notebook is  straightforward to understand.

DR=1/(NR-1);

rarr=Table[r,{r,0,1,DR}];

size=Length[rarr]

The lines below initialize the background metric to the form in eqn 10.38

initialA[r_]=1-α(1-r^2)^2;

initialB[r_]=1-α(1-r^2)^2;

initialC[r_]=(1/2)(1+r^2);

A0[λ_,r_]=initialA[r];

B0[λ_,r_]=initialB[r];

C0[λ_,r_]=initialC[r];

The table 'data' holds the metric component values at the grid points.

metriccmpts=Table[{a[ii],b[ii],c[ii]},{ii,1,size}];

metriccmpts

In the following lines the boundary conditions are imposed on this table at the horizon and asymptotic boundaries.

metriccmpts[[1,1]]=(1/3)(4metriccmpts[[2,1]]-metriccmpts[[3,1]]);

metriccmpts[[1,2]]=(1/3)(4metriccmpts[[2,1]]-metriccmpts[[3,1]]);

metriccmpts[[1,3]]=(1/3)(4metriccmpts[[2,3]]-metriccmpts[[3,3]]);

metriccmpts[[2,2]]=metriccmpts[[2,1]]+(1/4)(metriccmpts[[3,2]]-metriccmpts[[3,1]]);

metriccmpts[[size,1]]=1;metriccmpts[[size,2]]=1;metriccmpts[[size,3]]=1;

metriccmpts=Simplify[metriccmpts];

metriccmpts

initialmetricvalues and metricvariables give a list with the initial and formal values of the metric component functions arranged as in the eqn just below eqn 10.34

initialmetricvalues=Flatten[Join[Table[initialA[rarr[[ii]]],{ii,2,size-1}],Table[initialB[rarr[[ii]]],{ii,3,size-1}],Table[initialC[rarr[[ii]]],{ii,2,size-1}]]];

metricvariables=Flatten[Join[Table[a[ii],{ii,2,size-1}],Table[b[ii],{ii,3,size-1}],Table[c[ii],{ii,2,size-1}]]];

derivs gives replacement rules to discretize the continuum equations

derivs=Join[{"End"},Table[Flatten[Table[{{A,B,C}[[jj]]->metriccmpts[[ii,jj]],{Ar,Br,Cr}[[jj]]->((metriccmpts[[ii+1,jj]]-metriccmpts[[ii-1,jj]])/(2DR)),{Arr,Brr,Crr}[[jj]]->((metriccmpts[[ii+1,jj]]+metriccmpts[[ii-1,jj]]-2metriccmpts[[ii,jj]])/(DR*DR))},{jj,1,3}]],{ii,2,size-1}],{"End"}];

equationvector will hold the discretized equations as in eqn 10.37.

equationvector=Simplify[Flatten[Join[Table[RF[[1]]/.derivs[[ii]]/.{r->rarr[[ii]]},{ii,2,size-1}],Table[RF[[2]]/.derivs[[ii]]/.{r->rarr[[ii]]},{ii,3,size-1}],Table[RF[[3]]/.derivs[[ii]]/.{r->rarr[[ii]]},{ii,2,size-1}]]]];

vecsize=Length[equationvector]

The linearized operator (linearOp) required by Newton's method is computed by explicit differentiation of the equations (equationvector) with respect to the discretized variables (metricvariables).

linearOp=Simplify[Table[D[equationvector[[ii]],metricvariables[[jj]]],{ii,1,vecsize},{jj,1,vecsize}]];

Section 3

Code to perform Ricci-DeTurck flow

Perform a Ricci-DeTurck flow using forward Euler differencing ( = Jacobi relaxation, see discussion around eqn 10.27) for an initial metric & reference metric with parameter α = 0.6 ;  this should collapse to a singularity after finite flow time  l ~ 2

α=0.6;

metricvalues=initialmetricvalues;

time=0;

Dλ=0.5*DR*DR;

radlist={};

Do[

Do[

dmetricvaluesBydλ=N[equationvector/.Thread[metricvariables->metricvalues]];

metricvalues+=Dλ*dmetricvaluesBydλ;
time+=Dλ;

,{iter2,1,50}];

sphereradius=Sqrt[(metriccmpts[[1,3]]/.Thread[metricvariables->metricvalues])];

AppendTo[radlist,{time,sphereradius}];

Print[" flow time λ = ",time,"  radius of horizon 2-sphere ( = (Sqrt[C (0)]) ) = ",sphereradius];

,{iter1,1,27}];

plt1=ListPlot[radlist,AxesLabel->{"λ","horizon radius"}]

Perform a Ricci-DeTurck flow for an initial metric & reference metric with parameter α = +0.8 ;  this should expand for ever.

α=0.8;

metricvalues=initialmetricvalues;

time=0;

Dλ=0.5*DR*DR;

radlist={};

Do[

Do[

dmetricvaluesBydλ=N[equationvector/.Thread[metricvariables->metricvalues]];

metricvalues+=Dλ*dmetricvaluesBydλ;
time+=Dλ;

,{iter2,1,50}];

sphereradius=Sqrt[(metriccmpts[[1,3]]/.Thread[metricvariables->metricvalues])];

AppendTo[radlist,{time,sphereradius}];

Print[" flow time λ = ",time,"  radius of horizon 2-sphere ( = (Sqrt[C (0)]) ) = ",sphereradius];

,{iter1,1,35}];

plt2=ListPlot[radlist,AxesLabel->{"λ","horizon radius"}]

Perform a Ricci-DeTurck flow for an initial metric & reference metric with parameter α = +0.732 ;  this is close to the value α* for 20 points.

α=0.732;

metricvalues=initialmetricvalues;

time=0;

Dλ=0.5*DR*DR;

radlist={};

Do[

Do[

dmetricvaluesBydλ=N[equationvector/.Thread[metricvariables->metricvalues]];

metricvalues+=Dλ*dmetricvaluesBydλ;
time+=Dλ;

,{iter2,1,50}];

sphereradius=Sqrt[(metriccmpts[[1,3]]/.Thread[metricvariables->metricvalues])];

AppendTo[radlist,{time,sphereradius}];

Print[" flow time λ = ",time,"  radius of horizon 2-sphere ( = (Sqrt[C (0)]) ) = ",sphereradius];

,{iter1,1,35}];

plt3=ListPlot[radlist,AxesLabel->{"λ","horizon radius"}]

Show horizon radius against flow time for the 3 flows above;

Show[plt1,plt2,plt3,PlotRange->{{0,All},{0,1.5}}]

Section 4

Code to perform Newton's method

Solve using Newton's method for an initial guess (& background metric) with parameter α = 0.6

α=0.6;

metricvalues=initialmetricvalues;

Do[

equationvectorValue=N[equationvector/.Thread[metricvariables->metricvalues]];

linearOpValue=N[linearOp/.Thread[metricvariables->metricvalues]];

Δmetricvalues=-LinearSolve[linearOpValue,equationvectorValue];

Print[iter," -> norm = ",Norm[Δmetricvalues]/vecsize,"  horizon 2-sphere radius ( = (Sqrt[C (0)]) ) = ",N[(Sqrt[metriccmpts[[1,3]]]/.Thread[metricvariables->metricvalues])]];

metricvalues+=Δmetricvalues;

,{iter,1,10}];

Note that the correct continuum value for the horizon radius is exactly one.

This simple routine allows the metric functions to be visualized.

disp:=(

dataval=metriccmpts/.Thread[metricvariables->metricvalues];

pltA=ListPlot[Transpose[{rarr,dataval[[All,1]]}],PlotStyle->{Hue[0.]}];

pltB=ListPlot[Transpose[{rarr,dataval[[All,2]]}]];

pltC=ListPlot[Transpose[{rarr,dataval[[All,3]]}],PlotStyle->{Hue[0.6]}];

)

disp;

Print[Show[pltA,pltC,PlotRange->{{0,1.01},{0,1.5}},AxesLabel->{"r","A,C"}]]