Deriving an integral reduction formula programmatically

There are a number of integral reduction formulas from basic calculus, including several involving trigonometric or exponential functions. These are of the form where some term in an integrand is raised to the $$n$$th power then gets expressed as constants or other terms but with the same form integral with the power now "reduced" to $$n-1$$.

Here is one such integral reduction formula:

$$\int \frac{1}{(1 + x^2)^n}\ dx = \frac{x}{(2 n - 2)(x^2 + 1)^{n-1}} + \frac{2 n - 3}{2 n - 2} \int \frac{1}{(1 + x^2)^{n - 1}}\ dx$$

I'd like to derive this reduction formula computationally.

The obvious first step is to simply compute the integral:

Assuming[n \[Element] Integers,
Integrate[1/(x^2 + 1)^n, x]]


which yields:

$$x \, _2F_1\left(\frac{1}{2},n;\frac{3}{2};-x^2\right) ,$$

where $$\mbox{}_2F_1 (\cdot)$$ is the Hypergeometric function. This means we can immediately write down the integral which has power $$n-1$$ in the denominator by replacing $$n \to n-1$$ in the solution formula above.

But how to take these results to compute the reduction relation?

One step might be to perform the integral for $$n=1$$, which yields $$\tan^{-1} x$$ and try to form a recursion formula up to higher $$n$$ this way:

f[1, x] = ArcTan[x];
f[n_Integer, x_Real] := x f[n - 1, x]


(I don't understand why this code won't evaluate f[2,x] properly.)

I've tried expanding each result and subtracting to find the difference, or some form of recursion relation, but without success.

Or might there be another approach?

I'm not quite sure I understood what you're looking for.

If you only want to check if your recursion relation is correct, you can use partial memoization of an explicit recursion,

f[n_, x_] = Integrate[1/(x^2 + 1)^n, x]
(*    x Hypergeometric2F1[1/2, n, 3/2, -x^2]    *)

g[1] = ArcTan;
g[n_Integer /; n >= 2] := g[n] = Function[x, Evaluate[
x/((2 n - 2) (x^2 + 1)^(n - 1)) + (2 n - 3)/(2 n - 2) g[n - 1][x]]];
g[n_, x_] := g[n][x]

Table[g[n, x] == f[n, x], {n, 10}] // FullSimplify
(*    {True, True, True, True, True, True, True, True, True, True}    *)


If, on the other hand, you want to discover recursion relations, then the Wolfram functions site could be helpful. Particularly for the Gauss hypergeometric function $$_2F_1$$ we have $${}_2F_1(a,b;c;z)=\frac{c-2b+2+(b-a-1)z}{(b-1)(z-1)}{}_2F_1(a,b-1;c;z)+\frac{b-c-1}{(b-1)(z-1)}{}_2F_1(a,b-2;c;z)$$ which translates to $$f_n(x)=\frac{4n-7+(2n-3)x^2}{(2n-2) (1+x^2)}f_{n-1}(x)+\frac{5-2 n}{(2n-2) (1+x^2)}f_{n-2}(x)$$ Try it out:

f[n, x] == (4 n - 7 + (2 n - 3) x^2)/(2 (n - 1) (1 + x^2)) f[n - 1, x] +
(5 - 2 n)/(2 (n - 1) (1 + x^2)) f[n - 2, x] // FullSimplify
(*    True    *)


This isn't exactly the recurrence relation you're looking for, but it may get you started.

• Your final point is indeed helpful and a step in the right direction. I was hoping, however, to derive such a recursion relation, rather than look it up in the Mathematica documentation. That is, after all, what computer algebra is meant to avoid! Aug 14, 2021 at 19:11

I think this is a solution. The insight is to express the full integral as two terms, one that allows (symbolic) integration by parts to get the "reduced" form. Given the structure of the original integral, the natural variable to define is:

Let $$s(x) = x^2 + 1$$, and thus $$s^\prime (x) = 2 x$$.

s[x]:= x^2 +1;

s'[x]


(* 2 x *)

Suppose we can find polynomials $$a(x)$$ and $$b(x)$$ such that $$a(x)\ s(x) + b(x)\ s^\prime (x) = 1$$. Then our integral could be expressed as:

$$\int \frac{dx}{(x^2+1)^n} = \int \frac{dx}{s^n} = \int \frac{a(x)\ s(x)}{s^n(x)}\ dx + \int \frac{b(x)\ s^\prime (x)}{s^n(x)}\ dx\quad (*)$$

The whole (intermediate) goal here is to get an integrand with $$s^\prime (x)$$ in the numerator so we can perform integration by parts (below).

But what are these polynomials $$a$$ and $$b$$? They must sum to $$1$$:

Simplify@(b /. Solve[a s[x] + b s'[x] == 1 , {a, b}])


$$a(x) = a\quad {\rm and}\quad b(x) = \left\{\frac{1-a \left(x^2+1\right)}{2 x}\right\}$$

We use the general equation for integration by parts:

Assuming[n \[Element] Integers,
Integrate[q'[x]/q[x]^n, x]]


$$\frac{q(x)^{1-n}}{1-n}$$

and thus now the right-most of the integrals on the right-hand-side of (*) becomes:

$$\int \frac{b(x)\ s^\prime (x)}{s^n(x)}\ dx = \frac{-b(x)}{(n-1)s^{n-1}(x)} + \int \frac{b^\prime (x)}{(n-1)s^{n-1}}\ dx$$

We apply that to the term at the right we get:

$$\int \frac{b(x) s^\prime (x)}{s^n(x)}\ dx = \frac{-b(x)}{(n-1) s^{n-1}(x)} + \int \frac{b^\prime (x)}{(n-1) s^{n-1}(x)}\ dx$$

with

b[x_] := (1 - a (1 + x^2))/(2 x);


we find

b'[x]


(* -a - (1 - a (1 + x^2))/(2 x^2) *)

Putting all this together gives the relation to be derived.

If someone can put this into a more-elegant form, or use more symbolic programming, or find a cleaner derivation... GREAT. I'm happy to accept such a better solution. I really thought I could find an induction formula but alas couldn't.