3

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} \,$$

GrainOfSalt
  • 33

2 Answers2

3

$\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}

Felix Marin
  • 94,079
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.

jmerry
  • 19,943
  • This is wrong: tanh-sinh, Gaussian, and Gauss-Kronrod quadrature all handle infinite singularities just fine. This is why "strictly interior nodes" is a desirable property of quadrature methods. – user14717 Jan 27 '19 at 23:34
  • Unless you happen to put it in exactly the wrong place. Gaussian quadrature, for example, evaluates at the exact center of the interval half the time (I'm not familiar with the other methods you mentioned). Anyway, that's not the point - the point is that having an infinite singularity we might evaluate near completely wrecks any error estimates we might have, and evaluating at the infinite singularity leads to complete failure. – jmerry Jan 27 '19 at 23:40
  • Not true. You should try it. There is an entire class of endpoint singularities that Gaussian, tanh-sinh, Gauss-Kronrod evaluates correctly. And the problem at hand is not interior singularities, the problem is endpoint singularities. – user14717 Jan 27 '19 at 23:46
  • So there's a specific class of methods that handle (absolutely integrable) endpoint singularities. They don't handle singularities anywhere else. I stand by my original statement. – jmerry Jan 27 '19 at 23:49
  • It's not a specific class of quadrature methods that handle endpoint singularities, it's basically every method that you will find in provided in a modern numerical software package. As to interior singularities, you have to define it's meaning via (say) Cauchy principle value. And modern numerical packages handle that as well! – user14717 Jan 27 '19 at 23:57
  • Oh, the solution to an interior singularity is obvious in principle - you find out exactly where it is, and split there before applying something that handles endpoint singularities. If you're talking about principal values - well, I'm not. That way lies bogus results. And all of this is way off-topic to a question that was asking about Simpson's rule. – jmerry Jan 28 '19 at 00:16
  • For future reference, since it seems both this and the method above produce a proper integral, is one of them typically more efficient for computation? Is it better to keep your integral together or split it when you have the choice? – GrainOfSalt Jan 28 '19 at 01:22