0
$\begingroup$

I need to find the value of $z$ for a particular value of $D_c$ (eg. $500$), but $z$ is inside an integral, and I'm not able to use Solve since the integral is giving Hypergeometric2F1 function as the output.

OmegaM = 0.3111;
OmegaLambda = 0.6889;
Dc = 500;

eqn = Integrate[(OmegaM (1 + z1)^3 + OmegaLambda)^(-1/2), {z1, 0, z},
                  Assumptions -> z > 0]
-1.1473+(1.20482+1.20482z)Hypergeometric2F1[0.333333,0.5,1.33333,-0.451589(1.+z)^3]
zvalue = Solve[eqn == Dc, z]
Solve was unable to solve the system with inexact coefficients or the 
system obtained by direct rationalization of inexact numbers present
in the system. Since many of the methods used by Solve require exact
input, providing Solve with an exact version of the system may help.

Is there any other way I can solve this equation?

Also, Integrate is taking some time and I'd like it to be fast since I need to put it in a loop with lots of $z$ values to be computed for corresponding $D_c$ values.

$\endgroup$
  • $\begingroup$ There is a factor of c/H0 multiplied with the integral where c=300000, H0=68. I think that will make the equation consistent. Thanks mathematica.stackexchange.com/users/35368/superciocia for suggesting FindRoot, it solves the issue. Now I just want to know if I can speed up this integral somehow. $\endgroup$ – Divyajyoti May 4 at 5:45
4
$\begingroup$

From your integral I get the following equation (same as yours when you plot it):

eqn[z_] := 
 3.2566440560469836` - (
  3.5857498598223954` Hypergeometric2F1[1/6, 1/2, 7/
    6, -(2.2144005143040824`/(1 + z)^3)])/Sqrt[1 + z]

Both Solveand NSolve fail.

So I tried FindRoot:

Dc = 3.1;
FindRoot[ eqn[z] - Dc, {z, 0}]
{z -> 523.001}

which agrees with a graphical solution: enter image description here

I don't think it has a solution for Dc=500 as the eqn flattens out to 3.25664 as $z\rightarrow \infty$:

Limit[ eqn[z], {z -> ∞}]
3.25664
| improve this answer | |
$\endgroup$
  • $\begingroup$ How's your eqn is so different than OP's? $\endgroup$ – m0nhawk May 3 at 20:23
  • 1
    $\begingroup$ It's the same function. I had to change the integrand from $z1+1$ to $z2$ in order for Mathematica to give me an analytical result. Maybe it's because of different versions. $\endgroup$ – SuperCiocia May 3 at 20:25
  • $\begingroup$ Strange, I can reproduce this only for z1+1 as new variable, and Assumptions -> z>0. Otherwise it's OP's function. $\endgroup$ – m0nhawk May 3 at 20:33
  • $\begingroup$ What version do you have? $\endgroup$ – SuperCiocia May 3 at 20:34
  • $\begingroup$ 12.1. Yours? And did you specify Assumptions for eqn[z_]? $\endgroup$ – m0nhawk May 3 at 20:35
5
$\begingroup$

We will demonstrate that the exact formula for $z$ reads: $$z=\wp\bigg(\frac{\sqrt{\Omega_M}}{2}D_c+\wp^{-1}\big(1;0,-\frac{4\Omega_\Lambda}{\Omega_M}\big);0,-\frac{4\Omega_\Lambda}{\Omega_M}\bigg)-1$$ where $\wp(x;g_2,g_3)$ is the Weierstrass elliptic function, which yields a value $w$ in the elliptic integral $$x=\int^{w}_{\infty}\frac{d t}{\sqrt{4t^3-g_2\;t-g_3}}$$ and so generalizing the answer to the original question for any integrand of the form $\frac{1}{\sqrt{R(t)}}$, where $R(t)$ is a fourth or a third order polynomial in $t$. This formula can be implemented in the following way:

z[ Dc_, OM_, OL_]:= WeierstrassP[ Sqrt[OM/4] Dc+ InverseWeierstrassP[ 1, { 0,-4OL/OM}],
                                  { 0, -4OL/OM}]-1

We rationalize numeric constants to play with the system seamlessly (although this step is not neccesary):

{ OM, OL} = Rationalize[{ OmegaM = 0.3111, OmegaLambda = 0.6889}];

Let's derive $z$: $$D_c=\int^{z}_{0}\frac{d s}{\sqrt{\Omega_M (s+1)^3+\Omega_{\Lambda}}}=\frac{2}{\sqrt{\Omega_M}}\int^{z+1}_{1}\frac{d s}{\sqrt{4 s^3+\frac{4\Omega_{\Lambda}}{\Omega_M}}}=\\=\frac{2}{\sqrt{\Omega_M}}\Bigg(\int^{\infty}_{1}\frac{d s}{\sqrt{4 s^3+\frac{4\Omega_{\Lambda}}{\Omega_M}}}-\int^{\infty}_{z+1}\frac{d s}{\sqrt{4 s^3+\frac{4\Omega_{\Lambda}}{\Omega_M}}}\Bigg)=\\=\frac{2}{\sqrt{\Omega_M}}\Bigg(-\wp^{-1}\big(1;0,-\frac{4\Omega_\Lambda}{\Omega_M}\big)+\wp^{-1}\big(z+1;0,-\frac{4\Omega_\Lambda}{\Omega_M}\big)\Bigg)$$ and this implies our formula for $z$.

The formula for $z$ is valid in the range $0<D_c<D_{m}=3.25664$ and we can also derive an exact formula for $D_m$: $$D_m=\frac{2}{\sqrt{\Omega_M}} \Re\Big( 2\;\omega_{1}(0,g_3)-\wp^{-1}\big(1;0,g_3\big)\Big)$$ where $\Re$ is the real part, $\omega_{1}(0,g_3)$ is the Weierstrass half period and $g_3$ is the Weierstrass invariant, in our case $g_3=-\frac{4\Omega_{\Lambda}}{\Omega_M}$, implementing it:

g3=-4OL/OM;
Dm = 2/Sqrt[OM]( 2WeierstrassHalfPeriodW1[{0, g3}]-InverseWeierstrassP[1,{0, g3}])//Re//N
3.25664

Dm been calculated in version 12.1, however in earlier versions one should evaluate simply Dm = -2/Sqrt[OM] InverseWeierstrassP[1,{0, g3}]. This is because InverseWeierstrassP[1,{0, g3}] is computed in a adjacent parallelogram (see e.g this discussion). One should also note better handling of symbolic input in WeierstrssHalfPeriodW1 etc. For the presentation of the structure of $z$ being an elliptic function (shifted and rescaled $\wp$) we define:

wHP = Through @ { WeierstrassHalfPeriodW1,WeierstrassHalfPeriodW2,
                  WeierstrassHalfPeriodW3} @{ 0,-4OL/OM}//ReIm
                 // FullSimplify;
GraphicsRow @ Table[ ContourPlot[ Evaluate @ Table[p[z[x+I y,OM,OL]] ==k,       
                            {k, wHP[[#1,#2]]& @@@ {{2,1},{2,2},{3,1},{3,2}}}],
                        {x, -8, 8}, {y, -8, 8}, ContourStyle ->Thread[
                          {Thick,{Red,Darker@Cyan,Darker@Green,Orange}}]],
                     {p, {Re, Im}}]

enter image description here

There was an assumption that $z>0$, however $D_c=500$ can be reached for negative $z$, e.g.

z[ 500,OM, OL]//N//Chop
-1.73134

and for $0< z<D_m$ e.g.

z[ 2, OM, OL]//N//Chop
7.13731
| improve this answer | |
$\endgroup$
0
$\begingroup$
OmegaM = 0.3111 // Rationalize;
OmegaLambda = 0.6889 // Rationalize;
Dc = 500;

eqn = Integrate[(OmegaM (1 + z1)^3 + OmegaLambda)^(-1/2), {z1, 0, z}, 
  Assumptions -> z > 0]

enter image description here

For z > 0, eqn is monotonically increasing

Assuming[z > 0, D[eqn, z] > 0 // Simplify]

(* True *)

The maximum value of eqn is

(lim = Limit[eqn, z -> Infinity]) // N

(* 3.25664 *)

LogLinearPlot[{lim, eqn}, {z, 10^-2, 10^4},
 PlotLegends -> Placed["Expressions", {.3, .7}]]

enter image description here

Consequently, eqn can never equal the specified value of Dc

Using instead

Dc = 2;

Use NSolve

zvalue = NSolve[{eqn == Dc, z > 0}, z]

(* {{z -> 7.13731}} *)

Or FindRoot

zvalue = FindRoot[eqn == Dc, {z, 1}]

{z -> 7.13731}

Or Reduce (provides the exact value as a Root expression)

zvalue = Reduce[{eqn == Dc, z > 0}, z]

enter image description here

zvalue // N

(* z == 7.13731 *)

Similarly with Solve

zvalue = Solve[{eqn == Dc, z > 0}, z][[1]]

enter image description here

| improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.