I'm trying to integrate the following function:


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?


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.

  • $\begingroup$ You should be aware that your denominator may have some problems ... $\endgroup$ Jan 23 '13 at 23:35
  • $\begingroup$ This post Solving equations involving integrals includes evaluating the same integral from different point of view. $\endgroup$
    – Artes
    Feb 2 at 21:49

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]]

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

enter image description here

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.

  • $\begingroup$ surprisingly, Integrate[ 1/Sqrt[\[CapitalOmega]M (1 + z^3) + \[CapitalOmega]\[CapitalLambda]], {z, 0, x}, GenerateConditions -> False] takes a long time on 9 as well?? $\endgroup$
    – chris
    Jan 24 '13 at 0:21
  • $\begingroup$ I am currently using Mathematica 8. This function is used in cosmology and I've found a few papers dealing with just this issue. By massaging the function a bit and re-expressing it as the difference of two integrals of the form: Integrate[1/(u^4+u),{u,0,s}] - Integrate[1/(u^4+u),{u,0,s/(1+z)}] where s is a function of omega. Produces satisfactory results. $\endgroup$ Jan 24 '13 at 0:26
  • $\begingroup$ It takes a bit but not so long (say less than a minute). $\endgroup$
    – Artes
    Jan 24 '13 at 0:26
  • 1
    $\begingroup$ If OmegaM is negative and x at all large then the argument to the square root eventually becomes negative. So there must be other restrictions. $\endgroup$ Jan 24 '13 at 21:22
  • 1
    $\begingroup$ I apologize, you're right OmegaM and OmegaLambda are both positive and real $\endgroup$ Jan 25 '13 at 10:09

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).

 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.


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]])]

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.