# NIntegrate with a complicated function over the 4-sphere

I would like to maximize the following expression: $$\sum_{k = 1}^4\int_{\mathbb{S^4}}|u - \mathbf{e}_k.u_k||\nabla u_k|.\Big(1 - \frac{1}{8}\Big(\langle{v, \mathbf{e}_k}\rangle^2 + \langle{ e_k,\frac{\nabla u_k}{|\nabla u_k|}}\rangle^2 + \langle{\frac{\nabla u_k}{|\nabla u_k|}, \frac{u - u_k.\mathbf{e}_k}{u - u_k.\mathbf{e}_k}}\rangle^2 + \langle{\frac{u - u_k.\mathbf{e}_k}{u - u_k.\mathbf{e}_k}, v}\rangle^2\Big)\Big).dS$$ where $$dS$$ denotes the normalised measure of $$\mathbb{S}^4$$, $$u = u_1.\mathbf{e}_1 + u_2.\mathbf{e}_2 + u_3.\mathbf{e}_3 + u_4.\mathbf{e}_4$$, where $$\mathbf{e}_i$$ are the standard vector and $$\nabla$$ is the gradient on $$\mathbb{S}^4$$. For the maximization we ask that $$|u|^2 = u_1^2 + u_2^2 + u_3^2 + u_4^2$$ integrates to $$1$$ over $$\mathbb{S}^4$$ AND that $$u_i$$ are spherical harmonics of degree 3 over $$\mathbb{S}^4$$.

To do this with Mathematica, I am using a package by Sheldon Axler (https://www.axler.net/HFT_Math.html) that generates an orthonormal basis of spherical harmonics of degree 3 for me. The first step is to integrate the term for $$k = 1$$.

m = 4;
n = 3;
setDimension[x, n];
d = dimHarmonic[m, n];
basis = basisH[m, x,
Sphere] /. \
{\[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar]^2 ->
x^2 + y^2 + z^2 + t^2, x[1] -> x, x[2] -> y, x[3] -> z, x[4] -> t};
weight[a_?(VectorQ[#, NumericQ] &), b_?(VectorQ[#, NumericQ] &),
c_?(VectorQ[#, NumericQ] &), d_?(VectorQ[#, NumericQ] &)] :=
1 - 1/8*(x^2
+ ({x, y, z, t} . unit[{0, b . basis, c . basis, d . basis}])^2
+ (unit[spheregrad[a . basis]] . {1, 0, 0, 0})^2
unit[{0, b . basis, c . basis, d . basis}])^2);
term[a_?(VectorQ[#, NumericQ] &), b_?(VectorQ[#, NumericQ] &),
c_?(VectorQ[#, NumericQ] &),
d_?(VectorQ[#, NumericQ] &)] := (1/surfaceArea[n])*
NIntegrate[
Sqrt[(b . basis)^2 + (c . basis)^2 + (d . basis)^2]*
weighta[a, b, c, d], {x, y, z, t} \[Element] Sphere[n],
WorkingPrecision -> 15, MaxPoints -> 10000];


Here unit denotes the re-scaling of a vector to norm $$1$$, spheregrad is the function that return the spherical gradient. Evaluating the function term at a random point on $$\mathbb{S}^{64}$$ ($$d = 16$$ is the dimension of spherical harmonics of degree $$3$$ on $$\mathbb{S}^4$$, so writing $$(a,b,c,d) \in \mathbb{S}^{64}$$ with $$u_1 = a.basis$$, $$u_2 = b.basis$$ and so on) simply gives Warnings  and the function returns no answer. I suspect the reason is that the function weight is complicated. Any help in making term return some value is much appreciated!!

• To integrate a function $f(x,y,z,w)$ over the unit 4-sphere's surface you can integrate it in explicit hyperspherical coordinates: NIntegrate[f[Cos[θ], Cos[χ]*Sin[θ], Cos[φ]*Sin[θ]*Sin[χ], Sin[θ]*Sin[φ]*Sin[χ]] * Sin[θ]^2*Sin[χ], {θ, 0, π}, {χ, 0, π}, {φ, 0, 2 π}]. Aug 21, 2021 at 16:11
• Thanks! That might help...
– Ceka
Aug 21, 2021 at 16:53
• Thank you, that indeed works and I can accept that as an aswer. Do you perhaps know how to run NMaximize now with this? It seems that NMaximize[improvement[a[[1 ;; d]], a[[d + 1 ;; 2*d]], a[[2*d + 1 ;; 3*d]], a[[3*d + 1 ;; 4*d]]], a [Element] Sphere[4*d]] goes through only one iteration and returns NMaximize::ivar: True is not a valid variable., and then NMaximize[4.90928724127370, True]. Improvement is the sum of term functions as above for $k = 1, 2, 3, 4$. Thank you again!
– Ceka
Aug 22, 2021 at 10:42
• Please ask a separate question for this follow-up! Aug 22, 2021 at 15:18

As mentioned in the comments, hyperspherical coordinates give a direct parametrization of the surface of the unit 4-sphere. To integrate a function $$f(x,y,z,w)$$ on this hyperspherical surface $$x^2+y^2+z^2+w^2=1$$, we use the coordinates

$$x = \cos (\theta )\\ y = \sin (\theta ) \cos (\chi )\\ z = \sin (\theta ) \sin (\chi ) \cos (\phi)\\ w = \sin (\theta ) \sin (\chi ) \sin (\phi )$$

to give the 4-surface integral

NIntegrate[f[Cos[θ], Sin[θ]*Cos[χ], Sin[θ]*Sin[χ]*Cos[φ], Sin[θ]*Sin[χ]*Sin[φ]] *
Sin[θ]^2*Sin[χ], {θ, 0, π}, {χ, 0, π}, {φ, 0, 2π}]


To get more uniform sampling, use $$(x,y,\phi)$$ as parameters and thus eliminate the Jacobian:

NIntegrate[f[x, y, Sqrt[1-x^2-y^2]*Cos[φ], Sqrt[1-x^2-y^2]*Sin[φ]],
{x, -1, 1}, {y, -Sqrt[1-x^2], Sqrt[1-x^2]}, {φ, 0, 2π}]


There is a lot of symmetry in this integral, and you are allowed to permute all arguments if it makes your integral simpler.