Exercises lecture 5

ClearAll["Global`*"];SetDirectory[NotebookDirectory[]]

Exercise 1: Consider this list of points on a plane: points= {{1,6},{4,1.5},{0.5,2.5},{3,3},{1,5}}. Find by trying out systematically all permutations of the points the minimal path that starts from the origin, touches all the points once and comes back to the origin.
Check your result by using the built-in function  FindShortestTour. Plot the resulting path through the points.

points= {{1,6},{4,1.5},{0.5,2.5},{3,3},{1,5}};
distance[p1_,p2_]:=Sqrt[(p1[[1]]-p2[[1]])^2+(p1[[2]]-p2[[2]])^2];

newmod[a_,b_]=If[Mod[a,b]==0,b,Mod[a,b]];
np=Length[points];
list=Permutations[points];
len=Length[list];tour={};
start={0,0};Do[AppendTo[tour,distance[start,list[[i,1]]]+Sum[distance[list[[i,j]],list[[i,j+1]]],{j,np-1}]+distance[list[[i,np]],start]],{i,len}]
len
minpath=Min[tour]
pos=Flatten[Position[tour,minpath]]
Table[list[[pos[[i]]]],{i,Length[pos]}]

tour1=FindShortestTour[Join[{start},points]]
tour1[[2]]

full=Join[{start},points,{start}];fullbest=Table[full[[tour1[[2,i]]]],{i,Length[full]}]

ListPlot[fullbest,Joined->True]

Exercise 2: draw a contour plot of a 2 D Gaussian function Exp[-1/2*(x^2+y^2)], Then normalize numerically the Gaussian by integrating it in a region such that x^2 + y^2 < a, where a is big enough that the results is close to the normalization 2Pi of the Gaussian, then find by trial and error the value of a such that the integral of the normalized Gaussian contain 68% of the probability

ContourPlot[Exp[-1/2*(x^2+y^2)],{x, -3, 3},{y,-3,3}]

norm=NIntegrate[Exp[-1/2*(x^2+y^2)],{x,y}∈ImplicitRegion[x^2+y^2+x*y/2<16,{x,y}]]

NIntegrate[Exp[-1/2*(x^2+y^2)]/norm,{x,y}∈ImplicitRegion[x^2+y^2<1.51^2,{x,y}]]

Exercise 3: Open outside of Mathematica the external file Union2.1.dat (located at https://www.thphys.uni-heidelberg.de/~amendola/teaching/Union2.1.dat, which contains supernovae Ia data and examine its structure.
Read the file Union2.1.dat and make a plot of distance modulus dm, defined as m-M=25+5Log[dL], versus redshift, including errorbars from the fourth column (sigma). Then use the expressions below to find the luminosity distance in a LCDM model, and evaluate the chisquare (dm-dT)^2/sigma^2 where dT=25+5Log[dL] summed over every SNIa. Then find the best value of Omega_matter (assuming Omega_Lambda=1-Omega_matter) by minimizing the chisquare over a grid of Omega_matter values. Optional: find the confidence region for Omega_matter.

h=0.7;h0=1/3000*h;
hub[z_,Ωm0_,ΩΛ0_]:=h0 (Sqrt[ΩΛ0+(1-Ωm0-ΩΛ0)(((1+z))^2)+Ωm0(((1+z))^3)]);
(* luminosity distance in Mpc: *)
dLum[z_,Ωm0_,ΩΛ0_]:=(1+z)NIntegrate[(1/(hub[zz,Ωm0,ΩΛ0])),{zz,0,z},PrecisionGoal->5,AccuracyGoal->5]; 
(* theoretical distance modulus as a function of z, omega_matter, omega_lambda *)
dModTeo[z_,Ωm0_,ΩΛ0_]:=25+5*Log10[dLum[z,Ωm0,ΩΛ0]];

SetDirectory[NotebookDirectory[]];

cat=Import["Union2.1.dat"];
sncat=cat[[2;;-1]]; (* remove first line *)
Dimensions[sncat]
ztab=sncat[[All,2]];
distvec=sncat[[All,3]];
σDistvec=sncat[[All,4]];
len=Length[ztab]

h=0.7;h0=1/3000*h;
hub[z_,Ωm0_,ΩΛ0_]:=h0 (Sqrt[ΩΛ0+(1-Ωm0-ΩΛ0)(((1+z))^2)+Ωm0(((1+z))^3)]);
dLum[z_,Ωm0_,ΩΛ0_]:=(1+z)NIntegrate[(1/(hub[zz,Ωm0,ΩΛ0])),{zz,0,z},PrecisionGoal->5,AccuracyGoal->5]; (* in Mpc *)
dModTeo[z_,Ωm0_,ΩΛ0_]:=25+5*Log10[dLum[z,Ωm0,ΩΛ0]];

lp=ListPlot[Table[{ztab[[i]],Around[distvec[[i]],σDistvec[[i]]]},{i,len}],Frame->True];
tp=Plot[dModTeo[z,.3,.7],{z,.01,1.4}];
Show[lp,tp]

chisquared[Ωm0_,ΩΛ0_]:=Module[{},
dModTeovec=Table[dModTeo[ztab[[i]],Ωm0,ΩΛ0],{i,Length[ztab]}];
Total[(((((dModTeovec-distvec)/(σDistvec))))^2)]];

AbsoluteTiming[chisquared[.3,.7]]

chitab=Table[{om,chisquared[om,1-om]},{om,.1,.6,.01}];
chimin=Min[chitab[[All,2]]]
pos=Position[chitab[[All,2]],chimin][[1,1]]
chitab[[pos]]

Exercise 4: Define a function f[lx,ly,lab1,lab2,s] that produces a plot of two functions (eg two polynomials) with a legend placed at position lx,ly carrying labels lab1, lab2 and with frame fonts of size s

f[lx_,ly_,lab1_,lab2_,s_]:=Plot[{BesselY[0,x],BesselY[1,x]},{x,0,12},PlotLegends->Placed[{"a","b"},{lx,ly}],Frame->True,FrameStyle->s];f[.5,.15,"a","b",14]

Exercise 5: Find the Ricci scalar for 3D spatial non-flat metric. You can use either the metric in Cartesian form

Clear[metric];
metric:=DiagonalMatrix[
a^2*{1,1,1}/(1+k/4*(x^2+y^2+z^2))^2]

coord={x,y,z};

or in spherical form

Clear[metric];
metric:=DiagonalMatrix[
a^2*{1/(1-k*r^2),r^2,r^2*Sin[theta]^2}]

coord={r,theta,phi};

Clear[metric];
metric:=DiagonalMatrix[
a^2*{1,1,1}/(1+k/4*(x^2+y^2+z^2))^2]

coord={x,y,z};

raise[metric_, s_] := Inverse[metric] . s;
lower[metric_, s_] := metric . s;
Clear[grad]; grad[s_List, var_] := Outer[D, s, var];

Clear[gradScal]; gradScal[s_, var_] := Map[D[s, #] &, var]


div[s_, var_] := Inner[D, s, var, Plus];

Clear[christoffel]; christoffel[metric_, var_] := Inverse[metric] . (grad[metric, var] +
     			Transpose[grad[metric, var], {1, 3, 2}] - Transpose[grad[metric, var], {2, 3, 1}])/2
     			
chrContracted[metric_, var_] := gradScal[Log[(-Det[metric])^(1/2)], var]
     			

Clear[ricci]; ricci[metric_, var_] :=
 div[Transpose[christoffel[metric, var], {3, 2, 1}], var] - grad[chrContracted[metric, var], var] +
  		chrContracted[metric, var] . christoffel[metric, var] -
  		Table[Sum[(christoffel[metric, var] . christoffel[metric, var])[[i, j, kk, i]],
    				{i, Length[var]}], {j, Length[var]}, {kk, Length[var]}]
		
    				
Simplify[Tr[raise[metric, ricci[metric, coord]]]]

Clear[metric];
metric:=DiagonalMatrix[
a^2*{1/(1-k*r^2),r^2,r^2*Sin[theta]^2}]

coord={r,theta,phi};

raise[metric_, s_] := Inverse[metric] . s;
lower[metric_, s_] := metric . s;
Clear[grad]; grad[s_List, var_] := Outer[D, s, var];

Clear[gradScal]; gradScal[s_, var_] := Map[D[s, #] &, var]


div[s_, var_] := Inner[D, s, var, Plus];

Clear[christoffel]; christoffel[metric_, var_] := Inverse[metric] . (grad[metric, var] +
     			Transpose[grad[metric, var], {1, 3, 2}] - Transpose[grad[metric, var], {2, 3, 1}])/2
     			
chrContracted[metric_, var_] := gradScal[Log[(-Det[metric])^(1/2)], var]
     			

Clear[ricci]; ricci[metric_, var_] :=
 div[Transpose[christoffel[metric, var], {3, 2, 1}], var] - grad[chrContracted[metric, var], var] +
  		chrContracted[metric, var] . christoffel[metric, var] -
  		Table[Sum[(christoffel[metric, var] . christoffel[metric, var])[[i, j, kk, i]],
    				{i, Length[var]}], {j, Length[var]}, {kk, Length[var]}]
		
    				
Simplify[Tr[raise[metric, ricci[metric, coord]]]]