Integration output in expanded form

Question

Apparently Mathematica uses different algorithms to evaluate integrals to those we use. And when it comes to integration by substitution, in most cases, Mathematica gives the outputs in expanded form. For example:

In[1]:= Integrate[12 (y^4 + 4 y^2 + 8)^2 (y^3 + 2 y), y]
Out[1]:= 768 y^2 + 576 y^4 + 256 y^6 + 72 y^8 + 12 y^10 + y^12

When in fact the result can be written as (y^4+4y^2+8)^3. Notice that two expressions are not equal as they differ by the constant

In[2]:= Simplify[768 y^2 + 576 y^4 + 256 y^6 + 72 y^8 + 12 y^10 + 
   y^12 == (y^4 + 4 y^2 + 8)^3]
Out[2]= False

so it's not appropriate to use Factor,FullSimplify etc. to put it in a different form. How can I force Mathematica to give the result in the form as it was evaluated by hand?

Answer 1

Short Answer. We could use unknown coefficients to construct a general polynomial raised to a known power and differentiate that expression. Then we could solve for the coefficients that make the derivative equal to our integrand. For instance, we take $p = a_4 y^4 + a_3 y^3 + ... a_0$. We want to find the $a_4 \ldots a_0$ that make D[p^3,y] == f, where $f$ is our integrand. This works and it can even be generalized to the case where we have a function, $\sin x$, say, instead of $y$. However, if $p$ is a polynomial in $\sin x$, then differentiation will introduce $\cos x$, which may make it difficult or impossible for Solve[] to find the coefficients.

Long Answer. Instead of matching the derivative of a general polynomial to our integrand, we can match a general polynomial to the expression, $g$, returned by Integrate[]. Here is a function that can yield a simpler expression for the anti-derivative than Integrate[]. After the function definition are a few tests, followed by a brief explanation of the function.

ClearAll[myIntegral]
myIntegral[f_, y_] := Block[{g, subexpr, q, q1, sym, ord, expon,
   ndx, \[Mu], m, n, p, coef, a0, soln, lc, isMin, better,
   slen, isShortest, best, msol},
  g = Integrate[f, y];
  subexpr = Union@Level[g, \[Infinity]];

  (* Section 2 *)
  q1 = Select[subexpr, LeafCount[#] == 2 &];
  Which[Length[q1] > 1, Return[{g, {y -> y}}],
   Length[q1] == 1, q = First[q1],
   True, q = y];
  sym = Union[Cases[subexpr, _Symbol], {q}];
  ord = Length[CoefficientList[g, q]] - 1;
  expon = Reverse@Divisors[ord];

  (* Section 3 *)
  For[ndx = 1, ndx < Length[expon], ndx++,
   \[Mu] = expon[[ndx]];
   n = ord/\[Mu];
   coef = ToExpression /@ Table["a" <> ToString[k], {k, 0, n}];
   p = Plus @@ Table[coef[[k + 1]] q^k, {k, 0, n}] // 
     Collect[#, y] &;

   (* Find solution *)
   eqns = Thread[CoefficientList[p^\[Mu] - g - a0^\[Mu], q] == 0];
   soln = SolveAlways[eqns, sym];
   If[Length[soln] == 0,
    soln = Solve[eqns, coef]];
   If[Length[soln] == 0,
    eqns = Thread[CoefficientList[p^\[Mu] - g, q] == 0];
    soln = SolveAlways[eqns, sym]];
   If[Length[soln] == 0,
    soln = Solve[eqns, coef]];

   (* Less complicated is better *)
   If[Length[soln] > 0,
    lc = LeafCount /@ soln;
    isMin = Thread[lc == Min[lc]] // Simplify;
    better = Pick[soln, isMin], better = soln];

   (* Shortest is best *)
   If[Length[better] > 0,
    slen = StringLength[ToString[#]] & /@ better;
    isShortest = Thread[slen == Min[slen]] // Simplify;
    best = Pick[better, isShortest], best = better];

   (* Keep the first solution *)
   If[Length[best] > 0,
    msol = Append[#, m -> \[Mu]] & /@ best;
    Break[]]
   ];

  (* Simplify[ p^m /. First[msol] ] *)
  {p^m, msol}
  ]

Test 1. We take the integrand given in the OP as our first test function. The function myIntegral returns the general polynomial and replacement rules for the coefficients. We use /. on the first set of coefficients to obtain our answer.

f = 3 (y^4 + 4 y^2 + 8)^2 (4 y^3 + 8 y);
{polynom, coefs} = myIntegral[f, y];

answer = polynom /. First@coefs // Simplify

(*  (8 + 4 y^2 + y^4)^3  *)

Test 2. Sometimes Integrate[] returns an expression that already contains a constant of integration, as in this example. Here we see the replacement rules look complicated, but the final answer is same expression that Integrate[] returns.

{polynom, coefs} = myIntegral[4 (4 x + 1)^7, x]

answer = polynom /. First@coefs // Simplify

(* 

{(a0 + a1 x)^m, {{a0 -> 1/2^(3/8), a1 -> 2 2^(5/8), m -> 8}}}

1/8 (1 + 4 x)^8

*)

Test 3. Here we use an integrand that will lead us to a simple expression involving $\sin t$. Try this with various values of the exponent. This test fails if $t$ is replaced with $\pi * t$.

f = 3 b  Cos[t] (a + b Sin[t])^2;
{polynom, coefs} = myIntegral[f, t];

answer = polynom /. First@coefs // Simplify

(* (a + b Sin[t])^3 *)

Test 4. For a final demo / test we use a symbolic function $h(t)$.

f = 3 h'[t] (3 h[t]^2 - 2 h[t]) (6 - h[t]^2 + h[t]^3)^2;
{polynom, coefs} = myIntegral[f, t];

answer = polynom /. First@coefs // Simplify

(*  (6 - h[t]^2 + h[t]^3)^3  *)

For all of these tests it is interesting to compare the results of myIntegral to Integrate. The difference should be constant.

How it works. In the first part of the myIntegrate we integrate $f$ to obtain the integrand $g$. We then use Level to obtain all the subexpresion in $g$. In Section 2 of myIntegrate we first look for functions like $\sin x$ that have a LeafCount of 2. If there is only one, we use it as the base of our polynomial. We find all of the symbols in $g$, determine the polynomial order of $g$ and the exponents we want to try.

In Section 3 of myIntegral we loop over the exponents. We use n-th order polynomials raised to the $\mu$-th power. We try to solve for the coefficients 4 different ways. We try to match our polynomial $p_n^\mu$ to $g$ assuming $g$ has no constant term, $a_0^\mu$ using SolveAlways and Solve. Then we try 2 more times assuming $g$ has the constant term.

When we find solutions for the coefficients, we look for better solution(s) that have the smallest LeafCount. This help to get rid of the complex solutions. Of those solutions we look for best ones, those that produce the shortest string, since we prefer coefficients without negative signs.

Finally, if we have a solution for the coefficients, we return the general polynomial and the coefficients. There could be more than one set of coefficients. We do not look for solutions with lower exponents.

An alternative return value could be replace the coefficients in the polynomial, simplify it, and compare that result's LeafCount to the LeafCount of $g$. Then return the expression with the lower LeafCount.