(* yukterez.net *)

Quit[]

"strahlungsdominiertes Universum mit ΩR=1"
Quiet[FullSimplify[DSolve[{((a'[t])/(a[t]))==H0 (Sqrt[((ΩR)/((a[t])^4))]), a[0]==0},a[t],t],H0>0&&t>0]]

"Skalenfaktor"
a[t_]:=(Sqrt[2]) (Sqrt[H0 t]) ((ΩR)^(1/4));
T=(Solve[a[t]==α, t]/.α->a)[[1,1,2]];
"a[t]"->FullSimplify[a[t],Reals]->Block[{t=T},FullSimplify[a[t],Reals]]
"a'[t]"->FullSimplify[a'[t],Reals]->Block[{t=T},FullSimplify[a'[t],Reals]]
"a''[t]"->FullSimplify[a''[t],Reals]->Block[{t=T},FullSimplify[a''[t],Reals]]
"a'''[t]"->FullSimplify[a'''[t],Reals]->Block[{t=T},FullSimplify[a'''[t],Reals]]

"Hubbleparameter"
h[t_]=((a'[t])/(a[t]));
"H[t]"->Simplify[h[t]]->FullSimplify[Block[{t=T},h[t]],Reals]
"H'[t]"->Simplify[h'[t]]->Block[{t=T},FullSimplify[h'[t],Reals]]
"H''[t]"->Simplify[h''[t]]->Block[{t=T},FullSimplify[h''[t],Reals]]
"H'''[t]"->Simplify[h'''[t]]->Block[{t=T},FullSimplify[h'''[t],Reals]]

"Zeit"
"t[a]"->T

Quit[]

"materiedominiertes Universum, ΩM=1"
Quiet[FullSimplify[DSolve[{((a'[t])/(a[t]))==H0 (Sqrt[((ΩM)/((a[t])^3))]), a[0]==0},a[t],t],H0>0&&t>0]]

"Skalenfaktor"
a[t_]:=((((3/2)))^(2/3)) (((H0 t))^(2/3)) ((ΩM)^(1/3)); 
T=(Solve[a[t]==α, t]/.α->a)[[1,1,2]];
"a[t]"->Simplify[a[t]]->Block[{t=T},FullSimplify[a[t],Reals]]
"a'[t]"->Simplify[a'[t]]->Block[{t=T},FullSimplify[a'[t],Reals]]
"a''[t]"->Simplify[a''[t]]->Block[{t=T},FullSimplify[a''[t],Reals]]
"a'''[t]"->Simplify[a'''[t]]->Block[{t=T},FullSimplify[a'''[t],Reals]]

"Hubbleparameter"
h[t_]=((a'[t])/(a[t]));
"H[t]"->Simplify[h[t]]->FullSimplify[Block[{t=T},h[t]],Reals]
"H'[t]"->Simplify[h'[t]]->FullSimplify[Block[{t=T},h'[t]],Reals]
"H''[t]"->Simplify[h''[t]]->FullSimplify[Block[{t=T},h''[t]],Reals]
"H'''[t]"->Simplify[h'''[t]]->FullSimplify[Block[{t=T},h'''[t]],Reals]

"Zeit"
"t[a]"->T

Quit[]

"DE dominiertes Universum mit ΩΛ=1"
sol=Quiet[FullSimplify[DSolve[{((a'[t])/(a[t]))==H0 (Sqrt[ΩΛ]) , a[0]==1},a[t],t],H0>0&&t>0]]

"Skalenfaktor"
a[t_]:=(E^(H0 t (Sqrt[ΩΛ])));
T=(Solve[a[t]==α, t]/.α->a)[[1,1,2]];
"a[t]"->Simplify[a[t]]->Block[{t=T},Normal@FullSimplify[a[t],Reals]]
"a'[t]"->Simplify[a'[t]]->Block[{t=T},Normal@FullSimplify[a'[t],Reals]]
"a''[t]"->Simplify[a''[t]]->Block[{t=T},Normal@FullSimplify[a''[t],Reals]]
"a'''[t]"->Simplify[a'''[t]]->Block[{t=T},Normal@FullSimplify[a'''[t],Reals]]

"Hubbleparameter"
h[t_]=((a'[t])/(a[t]));
"H[t]"->Simplify[h[t]]
"H'[t]"->Simplify[h'[t]]
"H''[t]"->Simplify[h''[t]]
"H'''[t]"->Simplify[h'''[t]]

"Zeit"
"t[a]"->Normal@(T/.(C[1])->0)

Quit[]

"flaches Universum mit dunkler Energie, ΩΛ=1-ΩM"
Quiet[FullSimplify[DSolve[{((a'[t])/(a[t]))==H0 (Sqrt[((1-ΩΛ)/((a[t])^3))+ΩΛ]), a[0]==0},a[t],t],H0>0&&t>0&&0<ΩΛ<1]]

"Skalenfaktor"
a[t_]:=(((-1+(1/(ΩΛ))))^(1/3)) ((Sinh[(3/2) H0 t (Sqrt[ΩΛ])])^(2/3));
T=(Solve[a[t]==α, t]/.α->a)[[1,1,2]];
"a[t]"->Simplify[a[t]]->Block[{t=T},FullSimplify[a[t],H0>0&&a>0&&1>ΩΛ>0]]
"a'[t]"->Simplify[a'[t]]->Block[{t=T},FullSimplify[a'[t],H0>0&&a>0&&1>ΩΛ>0]]
"a''[t]"->Simplify[a''[t]]->Block[{t=T},FullSimplify[a''[t],H0>0&&a>0&&1>ΩΛ>0]]
"a'''[t]"->Simplify[a'''[t]]->Block[{t=T},FullSimplify[a'''[t],H0>0&&a>0&&1>ΩΛ>0]]

"Hubbleparameter"
h[t_]=((a'[t])/(a[t]));
"H[t]"->Simplify[h[t]]->FullSimplify[Block[{t=T},h[t]],H0>0&&a>0&&1>ΩΛ>0]
"H'[t]"->Simplify[h'[t]]->FullSimplify[Block[{t=T},h'[t]],H0>0&&a>0&&1>ΩΛ>0]
"H''[t]"->Simplify[h''[t]]->FullSimplify[Block[{t=T},h''[t]],H0>0&&a>0&&1>ΩΛ>0]
"H'''[t]"->Simplify[h'''[t]]->FullSimplify[Block[{t=T},h'''[t]],H0>0&&a>0&&1>ΩΛ>0]

"Zeit"
"t[a]"->FullSimplify[T,H0>0&&a>0&&1>ΩΛ>0]

Quit[]

"Geschlossenes Universum mit ΩR>1"
sol=Quiet[FullSimplify[DSolve[{((a'[t])/(a[t]))==H0 (Sqrt[((ΩR)/((a[t])^4))+((1-ΩR)/((a[t])^2))]) , a[0]==a0},a[t],t],ΩR>1&&t>0]]

"Skalenfaktor"
a[t_]:=(Sqrt[(a0^2)+H0 t (2 (Sqrt[(a0^2)+ΩR-(a0^2) ΩR])+H0 (t-t ΩR))])

"Zeit"
T=(Solve[a[t]==α, t]/.α->a)[[1,1]];

"Beispiel"
ΩR=11/10;H0=1;a0=1/1000; 
tMax=FindMaximum[a[t],{t,1}, WorkingPrecision->32][[2,1,2]]; T:=If[t<tMax,t,2tMax-t]
Η[t_]:=a'[t]/a[t]; H[t_]:=If[t<tMax,Η[t],-Η[2tMax-t]]
Plot[a[T],{t,0,2tMax}, Frame->True,GridLines->{{tMax},{}}, AspectRatio->1/2]
Plot[H[t],{t,0,2tMax}, Frame->True,GridLines->{{tMax},{}}, AspectRatio->1/2]

Quit[]

"Geschlossenes Universum mit ΩM>1"
sol=Quiet[FullSimplify[DSolve[{((a'[t])/(a[t]))==H0 (Sqrt[((ΩM)/((a[t])^3))+((1-ΩM)/((a[t])^2))]) , a[0]==a0},a[t],t],ΩM>1&&t>0]]

"Funktion"
f[x_]:=((ΩM Log[-(Sqrt[1-ΩM]) (Sqrt[x])+(Sqrt[ΩM+x-ΩM x])] (Sqrt[ΩM+x-ΩM x])+(Sqrt[1-ΩM]) (Sqrt[x]) (ΩM+x-ΩM x))/((((1-ΩM))^(3/2)) (x^(3/2)) (Sqrt[((ΩM+x-ΩM x)/(x^3))])))

"Argument"
x[t_]:=((-(Sqrt[a0 (1-ΩM)]) (-ΩM+a0 (-1+ΩM) (1+H0 t (Sqrt[((a0+ΩM-a0 ΩM)/(a0^3))])))+ΩM (Sqrt[a0+ΩM-a0 ΩM]) Log[-(Sqrt[a0 (1-ΩM)])+(Sqrt[a0+ΩM-a0 ΩM])])/((a0^(3/2)) (((1-ΩM))^(3/2)) (Sqrt[((a0+ΩM-a0 ΩM)/(a0^3))])))

"Skalenfaktor"
a[t_]:=InverseFunction[f[#1]&][x[t]]

"Beispiel"
ΩM=11/10;H0=1;a0=1/1000; 
tMax=FindMaximum[a[t],{t,1}, WorkingPrecision->32][[2,1,2]]; T:=If[t<tMax,t,2tMax-t]
Η[t_]:=a'[t]/a[t]; H[t_]:=If[t<tMax,Η[t],-Η[2tMax-t]]
Plot[a[T],{t,0,2tMax}, Frame->True,GridLines->{{tMax},{}}, AspectRatio->1/2]
Plot[H[t],{t,0,2tMax}, Frame->True,GridLines->{{tMax},{}}, AspectRatio->1/2]

Quit[]

"Geschlossenes Universum mit ΩΛ>1"
sol=Quiet[FullSimplify[DSolve[{((a'[t])/(a[t]))==H0 (Sqrt[ΩΛ+((1-ΩΛ)/((a[t])^2))]) , a[0]==a0},a[t],t],ΩΛ>1&&a0>(Sqrt[1-(1/(ΩΛ))])]]

"Skalenfaktor"
a[t_]:=((Tanh[H0 t (Sqrt[ΩΛ])+ArcTanh[((a0 ΩΛ)/(Sqrt[ΩΛ (1+(-1+(a0^2)) ΩΛ)]))]])/(Sqrt[-((ΩΛ ((Sech[H0 t (Sqrt[ΩΛ])+ArcTanh[((a0 ΩΛ)/(Sqrt[ΩΛ (1+(-1+(a0^2)) ΩΛ)]))]])^2))/(-1+ΩΛ))]))

"Beispiel"
ΩΛ=11/10;H0=1;a0=1;tMax=2;
H[t_]:=a'[t]/a[t];
Plot[a[t],{t,0,2tMax}, Frame->True,GridLines->{{tMax},{}}, AspectRatio->1/2]
Plot[H[t],{t,0,2tMax}, Frame->True,GridLines->{{tMax},{}}, AspectRatio->1/2]

"Minimal a"
H->0 @ a->(Sqrt[1-1/ΩΛ])

Quit[]

"Geschlossenes Universum mit ΩΛ<1"
sol=Quiet[FullSimplify[DSolve[{((a'[t])/(a[t]))==H0 (Sqrt[ΩΛ+((1-ΩΛ)/((a[t])^2))]) , a[0]==a0},a[t],t],ΩΛ>1&&a0>(Sqrt[1-(1/(ΩΛ))])]]

"Skalenfaktor"
a[t_]:=((Tanh[H0 t (Sqrt[ΩΛ])+ArcTanh[((a0 ΩΛ)/(Sqrt[ΩΛ (1+(-1+(a0^2)) ΩΛ)]))]])/(Sqrt[-((ΩΛ ((Sech[H0 t (Sqrt[ΩΛ])+ArcTanh[((a0 ΩΛ)/(Sqrt[ΩΛ (1+(-1+(a0^2)) ΩΛ)]))]])^2))/(-1+ΩΛ))]))

"Beispiel"
ΩΛ=1/10;H0=1;a0=0;tMax=2;
H[t_]:=a'[t]/a[t];
Plot[a[t],{t,0,2tMax}, Frame->True,GridLines->{{tMax},{}}, AspectRatio->1/2]
Plot[H[t],{t,0,2tMax}, Frame->True,GridLines->{{tMax},{}}, AspectRatio->1/2]

Quit[]

"vollständige Lösung, numerisch"

ΩR=92136*^-9;
ΩM=315/1000;
ΩΛ=1-ΩM-ΩR;
ΩK=1-ΩM-ΩR-ΩΛ; 
Gyr=(10^7)*36525*24*3600;
H0=217618*^-23;

plot[f_,text_]:=
Grid[{{text},
{Plot[f,{t,0,20},
Frame->True,GridLines->{{t0},{0,1}}, 
ImageSize->280, PlotPoints->280-56-5,
PlotStyle->Thick,
ImagePadding->{{56,5},{18,5}}]}}]

sol=NDSolve[{((a'[t])/(a[t]))==H0 (Sqrt[((ΩR)/((a[t])^4))+((ΩM)/((a[t])^3))+((ΩK)/((a[t])^2))+ΩΛ]), a[0]==1*^-15},
a[t],{t,0,20Gyr}, WorkingPrecision->48,
Method->{"ImplicitRungeKutta","DifferenceOrder"->20}];
t0:=t/.FindRoot[а[t Gyr]-1,{t,13}];

"Skalenfaktor"
а[t_]:=Evaluate[a[t]/.sol]
Grid[{{
plot[а[t Gyr],"a"],
plot[а'[t Gyr],"a'"],
plot[а''[t Gyr],"a''"],
plot[а'''[t Gyr],"a'''"]}}]

"Hubbleparameter"
H[t_]:=((а'[t])/(а[t]))
Grid[{{
plot[H[t Gyr],"H"],
plot[H'[t Gyr],"H'"],
plot[H''[t Gyr],"H''"],
plot[H'''[t Gyr],"H'''"]}}]

"t0"->t0 "Gyr"