Definite integral takes a very long time

Question

I'm trying to integrate the following function:

$$r_e=\int_0^z\frac{dz}{\sqrt{\Omega_m(1+z)^3+\Omega_\Lambda}}$$

where $\Omega_\Lambda=1-\Omega_M$ and both $\Omega_\Lambda$ and $\Omega_M$ are real and positive.

Doing the definite integral using Integrate doesn't produce a result even after a few hours of running. I found that just doing the indefinite integral and then putting in the limits later produces a result. Why is this?

Additionally, the result is unexpected and discontinuous for different values of $z$ depending on the value of $\Omega_M$. Any suggestions on the best strategy for integrating this?

FunH[omM_,omL_,z_]:=1/((omM(1+z)^3+omL)^(1./2.))
F[omM_,z_]:=Evaluate[Block[{omL=1-omM},Integrate[FunH[omM,omL,z],z]]]
Table[(1+z)^(-1)*(F[omM,z]-F[omM,0]),{omM,0.000001,1,0.1},{z,0.000001,3,0.001}]

I had to copy the code from another machine so I might have missed a bracket somewhere.

EDIT: If anyone is interested, there is a publicly available paper dealing with a transformation of this integral to make the solution easier.

Essentially one can transform the above integral using a few simple substitutions into something which is proportionate to the difference between two separate integrals of the form: $$\int\frac{du}{\sqrt{u^4+u}}$$ This can be easily integrated by Mathematica.

Answer 1

Start with the transformed integral

$$\begin{align*} \int_0^z\frac{\mathrm du}{\sqrt{a(u+1)^3+b}}&=\int_1^{z+1}\frac{\mathrm du}{\sqrt{au^3+b}}\\ &=\frac{\sqrt[3]{b/a}}{\sqrt{b}}\int_{\frac1{\sqrt[3]{b/a}}}^{\frac{z+1}{\sqrt[3]{b/a}}}\frac{\mathrm du}{\sqrt{u^3+1}} \end{align*}$$

Now, I had already dealt with that elliptic integral in this math.SE answer, so I'm not repeating the work here; suffice it to say that we now have an explicit expression for your integral in terms of EllipticF[]:

With[{a = 4/3, b = 5/2, z = 7},
     {(* integral *)
      NIntegrate[1/Sqrt[a (t + 1)^3 + b], {t, 0, z}, WorkingPrecision -> 20],
      (* closed form solution *)
      N[((b/a)^(1/3)/(3^(1/4) Sqrt[b]))
        (EllipticF[ArcCos[2 Sqrt[3]/(1 + Sqrt[3] + (z + 1)/(b/a)^(1/3)) - 1], (2 + Sqrt[3])/4] -
         EllipticF[ArcCos[2 Sqrt[3]/(1 + Sqrt[3] + 1/(b/a)^(1/3)) - 1], (2 + Sqrt[3])/4]), 20]}]

   {0.97715948427474779292, 0.97715948427474779300}

Alternatively, since EllipticF[ArcCos[u], m] == InverseJacobiCN[u, m], you can also use this function:

With[{a = 4/3, b = 5/2, z = 7},
 N[((b/a)^(1/3)/(3^(1/4) Sqrt[b]))
   (InverseJacobiCN[2 Sqrt[3]/(1 + Sqrt[3] + (z + 1)/(b/a)^(1/3)) - 1, (2 + Sqrt[3])/4] - 
    InverseJacobiCN[2 Sqrt[3]/(1 + Sqrt[3] + 1/(b/a)^(1/3)) - 1, (2 + Sqrt[3])/4]), 20]]

Answer 2

In our case the symbolic definite integral appears to be much more difficult to evaluate. This is because the complex functions underlying the given integral strongly depend on the parameters and may involve branch points and it implies a huge number of possible cases which the system tries to collect in ConditionalExpression. In general one might assume sufficient number of relations between the symbolic parameters, but in our case it seems to fail. As a more convenient way we define this function :

int[ ΩM_, ΩΛ_, x_] := Integrate[ 1/Sqrt[ΩM (1 + z)^3 + ΩΛ], {z, 0, x}]

now we can evaluate definite integrals in various cases :

int[2, 1, 1]

  (1/(3^(1/4)))(-1)^(1/12) 2^(2/3) ( 
     EllipticF[ ArcSin[((-1)^(11/12) Sqrt[1 + 2^(1/3)])/3^(1/4)], (-1)^(1/3)]
     - EllipticF[ ArcSin[((-1)^(11/12) Sqrt[1 + 2 2^(1/3)])/3^(1/4)], (-1)^(1/3)])

Nevertheless built-in rewrite rules are not always capable to give interesting results, e.g. in ver.9 the first and the second integrals are unevaluated unlike in ver.8

{int[1, 2, 3], int[2, Pi/2, 1], int[Sqrt[3], Pi, -1]} // TraditionalForm

Many assumptions underlying the symbolic definite integral appear to be a no-go issue unlike in a case by case basis of the definite integral.

Answer 3

For reasons noted in a comment, I don't think oM should be allowed to go negative. Not sure if it can be arbitrarily large. Let's say we are operating under the assumption that 0 < oM.

This is what I would like to evaluate:

ii = Integrate[1/Sqrt[oM (1 + z)^3 + (1 - oM)], {z, 0, x}, 
  Assumptions -> {x > 0, 0 < oM < 1}]

Unfortunately, this does not succeed. The reason is that the integration code becomes very sulky when it sees elliptics, especially with parameters, due to the difficulty of discerning path singularities.

So here is an old trick. Give it an explicit value for oM, one that you can use to later recover a (possibly correct) parametrized result. I use $1/e$ for this purpose (since I don't see $e$ in the indefinite integral).

Timing[
 ii = Integrate[1/Sqrt[1/E (1 + z)^3 + (1 - 1/E)], {z, 0, x}, 
    Assumptions -> 0 < x];]

(* Out[1]= {21.960000, Null} *)

Numeric check:

Chop[Table[ii, {x, 1., 11., 1}]]

(* Out[2]= {0.742608949781, 1.14591119954, 1.39565792161, 1.56802998702, 
1.69583835424, 1.79538009669, 1.87570826106, 1.94228852684, 
1.9986361674, 2.04712646089, 2.08942907726} *)

Table[NIntegrate[1/Sqrt[1/E (1 + z)^3 + (1 - 1/E)], {z, 0, x}], {x, 1., 11., 1}]

(* Out[3]= {0.742608949781, 1.14591119954, 1.39565792161, 1.56802998702, 
1.69583835424, 1.79538009669, 1.87570826106, 1.94228852684, 
1.9986361674, 2.04712646089, 2.08942907726} *)

Okay, so it's fine at least for one value of the parameter oM, and a handful of values for x. (That proves that it never fails. NEVER.) We can now do

ii2 = ii /. E -> 1/oM

Use at your own risk, though. It could become garbage depending on settings for x and/or oM.

Answer 4

With "12.2.0 for Mac OS X x86 (64-bit) (December 12, 2020)"

Integrate[1/Sqrt[Subscript[\[CapitalOmega], \[CapitalLambda]] + 
Subscript[\[CapitalOmega], m]*(1 + x)^3], {x, 0, z}]

gives a solution in terms of Hypergeometric Functions:

(Sqrt[Subscript[\[CapitalOmega], \[CapitalLambda]] + (1 + z)^3*
   Subscript[\[CapitalOmega], m]]/
Subscript[\[CapitalOmega], \[CapitalLambda]])*(1 + z)*
Hypergeometric2F1[5/6, 1, 
4/3, (-(Subscript[\[CapitalOmega], m]/
    Subscript[\[CapitalOmega], \[CapitalLambda]]))*(1 + z)^3] - 
(Sqrt[Subscript[\[CapitalOmega], m] + 
  Subscript[\[CapitalOmega], \[CapitalLambda]]]/
Subscript[\[CapitalOmega], \[CapitalLambda]])*
Hypergeometric2F1[5/6, 1, 
4/3, -(Subscript[\[CapitalOmega], m]/
  Subscript[\[CapitalOmega], \[CapitalLambda]])]