# Factoring a separable integral with a product of independent integrands

I would like to input

someFunction[Integrate[p[x] p[y], {x, -1, 1}, {y, -1, 1}]]


and to get the following output

(Integrate[p[x], {x, -1, 1}])^2


Is this possible? (Function p is undefined, I am interested in a symbolic expression only.)

UPD: The question is not that how to define someFunction for this particular case. I am interested in a method which will work for any multiple integral: if there is a product function inside such integral then the answer has to be factorized as well.

Best wishes,

Dmitri.

• For searching purposes: what OP is dealing with is also called a separable integral. May 13, 2015 at 23:41
• Welcome to Mathematica.SE! I suggest that: 1) You take the introductory Tour now! 2) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! 3) As you receive help, try to give it too, by answering questions in your area of expertise. May 14, 2015 at 20:55

This will work with any number of independent integrants. Define:

repl[l_] := # /. Thread[#[[All, 1]] -> Table[x, {Length[#]}]] &[l]

inTfaC[int_] := Times @@ MapThread[Integrate[#1, #2] &,
repl /@ {First[#], Rest[#]} &[
int /. {Integrate -> List, Times -> List}]]


Now verify:

test = Integrate[p[x] p[y] q[z] r[s] r[u],
{x, -1, 1}, {y, -1, 1}, {z, a, 2}, {s, 3, b}, {u, 3, b}]


inTfaC[test]


• Probably, I was not exact. I do not want to teach Mathematica how to factorize the answer. I would like to ask it to do this. The same should be done for a triple integral and so on, May 13, 2015 at 22:18
• @Dmitri updated May 13, 2015 at 22:34
• Great! Thank you very much!!! May 13, 2015 at 22:52
• Be very careful! This will even factor integrals where the variables only "talk to each other" through the integrals' limits! For example, try inTfaC[Integrate[p[x] p[y], {x, -1, y}, {y, -1, 1}]]. May 15, 2015 at 21:17
• @evanb Your integral does not satisfy the question title: "Factoring a separable integral with a product of independent integrands " - and I assume that OP deals with those only. May 15, 2015 at 21:22

Adapting linearExpand from my answer to How to do algebra on unevaluated integrals?, we can come up with some transformations to factor separable multiple integrals. The function someFunction internally deals with and returns Inactive integrals, which can be evaluated with Activate, if appropriate or desired.

Examples

someFunction[Integrate[p[x] p[y], {x, -1, 1}, {y, -1, 1}]]


someFunction[Integrate[p[x] p[y], {x, -1, 1}, {y, -1, 100}]]


someFunction[Integrate[p[x] p[y], {x, -1, 1}, {y, -1, x}]]


someFunction[Integrate[p[x] (p[y] + q[y]), {x, -1, 1}, {y, -1, 1}]^2]


(* Vitaliy Kaurov's example *)
test = Integrate[p[x]*p[y]*q[z]*r[s]*r[u],
{x, -1, 1}, {y, -1, 1}, {z, 0, 2}, {s, 3, 10}, {u, 3, 10}];
someFunction[test]
% /. changeVar[x]


someFunction[Integrate[q[x], {x, 0, 2}] + test]


Outline

We need some auxilliary functions.

• iterated converts a multiple integral into an iterated integral.
• linearExpand applies linearity properties to an integral, distributing the integral over a sum and factoring out constants.
• factorConstants recursively factors constants out of nested integrals using linearExpand.
• changeVar changes the variable of integration.

These can be applied by Simplify via the TransformationFunctions option. The trick, and it can be tricky, is devising a ComplexityFunction that will prefer a result in a factored out form. Note that we use changeVar only to combine two integrals that are equivalent. At the end we can change the variable in the independent integrals to the same one.

Code dump

ClearAll[linearExpand, iterated, factorConstants, changeVar, someFunction];

e //. {op : head[arg_Plus, __] :> Distribute[op],
head[arg1_Times, var_, rest___] /; ! FreeQ[var, x] :>
With[{dependencies = InternalDependsOnQ[#, x] & /@ List @@ arg1},
Pick[arg1, dependencies, False] *

iterated[Integrate[f_, dom : {_Symbol, _, _} ..]] :=
Fold[Inactive[Integrate], f, Reverse@{dom}];
iterated[Inactive[Integrate][f_, dom : {_Symbol, _, _} ..]] :=
Fold[Inactive[Integrate], f, Reverse@{dom}];

factorConstants[Inactive[Integrate][f_, {v_Symbol, a_, b_}]] :=
Inactive[Integrate][
f /. j : Inactive[Integrate][_, _] :> factorConstants[j], {v, a, b}] /.
i : Inactive[Integrate][_, {v, _, _}] :>
linearExpand[i, v, Inactive[Integrate]];
factorConstants[x_] := x;

changeVar[Inactive[Integrate][f_, {v_, a_, b_}]] :=
i : Inactive[Integrate][g_, {w_, a, b}] /;
Simplify[f == g /. v -> w] :> (i /. w -> v);
changeVar[v_] :=
i : Inactive[Integrate][g_, {w_Symbol, _, _}] :> (i /. w /; FreeQ[g, Integrate] -> v);

someFunction[integral_] :=
With[{v =
FirstCase[Hold[integral],
HoldPattern[(Integrate | Inactive[Integrate])[_, {x_, _, _}, ___]] :> x,
Missing[],
Infinity]},
(Simplify[
integral /.
{i : (Integrate | Inactive[Integrate])[_, {_Symbol, _, _}, {_Symbol, _, _} ..] :>
iterated@Inactivate[i, Integrate],
i : (Integrate | Inactive[Integrate])[_, {_Symbol, _, _}] :>
Inactivate[i, Integrate]},
TransformationFunctions -> {Automatic,
# /. i : Inactive[Integrate][f_, {_, _, _}] :> factorConstants[i] &,
# /. changeVar /@
DeleteDuplicates@
Cases[#, Inactive[Integrate][f_, {v_, a_, b_}], Infinity] &},
ComplexityFunction -> (LeafCount[#] +
10 Count[{#},
Inactive[Integrate][Inactive[Integrate][__], __], Infinity] +
5 Count[{#},
Inactive[Integrate][f_, __] /; ! FreeQ[f, Integrate], Infinity] -
Count[{#}, Power[Inactive[Integrate][f_, __], _], Infinity] +
Length@DeleteDuplicates@
Cases[{#}, Inactive[Integrate][f_, {x_, _, _}] :> x, Infinity] &)
] /. changeVar[v]) /; FreeQ[v, Missing]
];
someFunction[expr_] := expr;


The first two terms after LeafCount of the ComplexityFunction penalize nested integrals, the first one penalizing nested integrals in which the constants have not been factored out. The next term, being subtracted, rewards gathering powers of integrals. The last term penalizes extra variables of integration, which drives Simplify to prefer changing variables to the same thing.

Maybe for this case use:

someFunction[int_] := (int /. _[y] :> 1/2)^2
`