Numerical Integration Bounded by Two Singularities

Question

I would like to solve the following definite integral numerically using Simpson's Rule, however it has singularities at both ends. I was told it's possible to perform a simple change of variable in order to transform it into something usable, but I've been unable to find the correct substitution. How would one replace the variables in this integral so that there would no longer be a singularity in the domain of the integral?

$$\int_0^1 x^{-2/3}(1-x)^{-1/3} \,$$

Answer 1

$\newcommand{\bbx}[1]{\,\bbox[15px,border:1px groove navy]{\displaystyle{#1}}\,} \newcommand{\braces}[1]{\left\lbrace\,{#1}\,\right\rbrace} \newcommand{\bracks}[1]{\left\lbrack\,{#1}\,\right\rbrack} \newcommand{\dd}{\mathrm{d}} \newcommand{\ds}[1]{\displaystyle{#1}} \newcommand{\expo}[1]{\,\mathrm{e}^{#1}\,} \newcommand{\ic}{\mathrm{i}} \newcommand{\mc}[1]{\mathcal{#1}} \newcommand{\mrm}[1]{\mathrm{#1}} \newcommand{\pars}[1]{\left(\,{#1}\,\right)} \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}} \newcommand{\root}[2][]{\,\sqrt[#1]{\,{#2}\,}\,} \newcommand{\totald}[3][]{\frac{\mathrm{d}^{#1} #2}{\mathrm{d} #3^{#1}}} \newcommand{\verts}[1]{\left\vert\,{#1}\,\right\vert}$ \begin{align} &\bbox[10px,#ffd]{\int_{0}^{1}x^{-2/3}\pars{1 - x}^{-1/3}\,\dd x} = \int_{0}^{1/2}x^{-2/3}\pars{1 - x}^{-1/3}\,\dd x + \int_{1/2}^{1}x^{-2/3}\pars{1 - x}^{-1/3}\,\dd x \end{align}

Set $\ds{x = t^{3}}$ in the First Integral and $\ds{x = 1 - t^{3}}$ in the Second Integral:

\begin{align} &\bbox[10px,#ffd]{\int_{0}^{1}x^{-2/3}\pars{1 - x}^{-1/3}\,\dd x} = 3\int_{0}^{2^{\large -1/3}}\pars{1 - t^{3}}^{-1/3}\,\dd t + 3\int_{0}^{2^{\large -1/3}}t\pars{1 - t^{3}}^{-2/3}\,\dd t \end{align}

Answer 2

Simpson's rule and most other numerical methods won't work in general for a function with an infinite singularity. In this case, it's guaranteed to fail - two of the points we evaluate at get us a literal $\infty$.

So, then, we need to transform the integral somehow into something we can use the rule on. How to we transform an integral to eliminate a weak singularity? With a singularity that looks like a multiple of $x^{-r}$, its antiderivative looks like a multiple of $x^{1-r}$, and that antiderivative will transform with the substitution. Substituting $t=x^{1-r}$ will make that antiderivative into a multiple of $t$, which matches to a function that tends to something nonzero at zero. This, of course, will only work for $r<1$.

So, let's apply that here:

\begin{align*}I &= \int_0^1 x^{-\frac23}(1-x)^{-\frac13}\,dx\\ &\phantom{|}^{y=x^{1/3}}_{dt=\frac13x^{-2/3}\,dx}\\ I &= \int_0^1 3(1-y^3)^{-\frac13}\,dy\\ &\phantom{|}^{z=(1-y)^{\frac23},\, y=1-z^{3/2}}_{dy=-\frac32 z^{1/2}}\\ I &= \int_1^0 -\frac92\left(1-(1-z^{\frac32})^3\right)^{-\frac13}z^{\frac12}\,dz\\ &= \frac92\int_0^1 \left(3z^{\frac32}-3z^3+z^{\frac92}\right)^{-\frac13}z^{\frac12}\,dz = \frac92\int_0^1 \left(3-3z^{\frac32}+z^3\right)^{-\frac13}\,dz\end{align*} Note that for the singularity at $1$, we use a power of $1-y$ rather than the $1-y^3$ in the parenthesis. Using $1-y^3$ would just reintroduce the singularity at the other end. As with any chain of substitutions, we could combine them into a single substitution - but it's just a lot clearer this way.

The new integral isn't pretty, but it's at least a proper Riemann integral. Numerical integration methods will work. Since it's not smooth at zero, we can't trust the numerical methods to converge as fast as they would for a nicer function, but they will at least converge.