# Solving a multivariate integral

How can I obtain from mathematica an exact solution to an integral of the form

$\int _0^a x_1^2dx_1\int_0^{\sqrt{a^2-x_1^2}}x_2^2dx_2\int_0^{\sqrt{a^2-x_1^2-x_2^2}}x_3^2dx_3$

And is it possible to obtain a solution for a general case of the same integral

$\int _0^a x_1^2dx_1\int_0^{\sqrt{a^2-x_1^2}}x_2^2dx_2\int_0^{\sqrt{a^2-x_1^2-x_2^2}}x_3^2dx_3...\int_0^{\sqrt{a^2-x_1^2-...-x_n^2}}x_n^2dx_n$

For a given n, one can formulate the integral like this, taking advantage of the symmetry of the integral:

ClearAll[i];
int[n_Integer?Positive] :=
Integrate[Product[x[i]^2, {i, n}],
Array[x, n] ∈ Ball[Table[0, n], a], Assumptions -> a > 0] / 2^n

int /@ Range
(*
{a^3/3, (a^6 π)/96,
(a^9 π)/1890,  <----  n == 3
(a^12 π^2)/184320, (a^15 π^2)/8108100, (a^18 π^3)/1486356480, (a^21 π^3)/109994484600,
(a^24 π^4)/31391848857600, (a^27 π^4)/3415328746830000, (a^30 π^5)/1371195958099968000}
*)


Response to comment:

Bob Hanlon points out below that FindSequenceFunction is a great tool for finding a general formula when you can calculate individual terms for given n. (The OP mentioned in a comment that a solution for a "given n" was desired. I guess that could mean for general n or for any given, particular integer n. The above does the latter.) FindSequenceFunction doesn't always find a formula, but it does in this case.

I spent the last half hour wondering how the integrals above were computed so quickly. I figured a general formula was being applied. Indeed it is. Essentially, the following function is used to compute the integral of a monomial $c\,x_1^{k_1}\cdots x_n^{k_n}$ over a ball of radius r centered at the origin:

monomialBallInt =
Function[{degs, coeff, r, dim},
If[Or @@ OddQ /@ degs,
0,
((2*coeff*r^(Total[degs] + dim)) * Product[Gamma[(1/2)*(1 + degs[[i]])], {i, dim}]) /
((Total[degs] + dim) * Gamma[Sum[(1/2)*(1 + degs[[i]]), {i, dim}]])
]];

monomialBallInt[{2, 2, 2}, 1, a, 3]/2^3
(*  (a^9 π)/1890  *)


We can adapt it to a function that will integrate an arbitrary monomial over the the segment of an $n$-ball in the positive hyperoctant.

ClearAll[monoHyperoctantBallInt];
monoHyperoctantBallInt[degs_, coeff_, r_, dim0_: Automatic] :=
With[{dim = dim0 /. Automatic -> Length[degs]},
((2*coeff*r^(Total[degs] + dim)) * Product[Gamma[(1/2)*(1 + degs[[i]])], {i, dim}]) /
(2^dim*(Total[degs] + dim) * Gamma[Sum[(1/2)*(1 + degs[[i]]), {i, dim}]])
];

monoHyperoctantBallInt[{2, 2, 2}, 1, a]
(*  (a^9 π)/1890  *)


General formula for the OP's problem. We use some coding tricks to make degs appear to be a vector of 2's of length n for the functions in the formula for monoHyperoctantBallInt:

Block[{degs},
degs /: degs[[_]] := 2;
degs /: Total[degs] := 2 n;
degs /: Map[OddQ, degs] := {False};
monoHyperoctantBallInt[degs, 1, a, n]
]
(*  (2^(1 - 2 n) a^(3 n) π^(n/2))/(3 n Gamma[(3 n)/2])  *)


That's all pretty good, I guess, but it would nicer if the problem could be posed to Mathematica as a general $n$-dimensional integral and the formula that is already present be used to return the answer.

• For list = int /@ Range then int2[n_] = FindSequenceFunction[list, n] and check with int == int2 – Bob Hanlon May 19 '18 at 18:52

The Solution of the first Integral is just

Integrate[x1^2 Integrate[x2^2 Integrate[x3^2, {x3, 0, Sqrt[a^2 - x1^2 - x2^2]}], {x2, 0,Sqrt[a^2 - x1^2]}], {x1, 0, a}]


this gives

($a^9 \pi$)/1890

The general solution seems to be

sol[n_] := a^(3 n) Pi^Floor[n/2] / c_n


with $c_n$ an numerical factor. Where the first values are

c={3, 96, 1890, 184320, 8108100, 1486356480, 109994484600,31391848857600}

• Can u generelize it to the nth case? – jarhead May 19 '18 at 15:28
• I am trying right now.. :) – kiara May 19 '18 at 15:44
• FindSequenceFunction[c, n] – Bob Hanlon May 19 '18 at 17:25
• It is not clear how you have obtained the general solution – jarhead May 19 '18 at 18:10
• the integral can be solved analytically , it was done before, but I would like to produce the solution in mathematica for any given n. – jarhead May 19 '18 at 18:11