|
%I #33 Nov 20 2022 18:31:10
%S 1,5,16,41,102,242,558,1263,2817,6214,13583,29471,63548,136305,291019,
%T 618849,1311314,2769847,5834119,12257072,25691785,53738815,112188059,
%U 233796875,486435094,1010552580,2096469429,4343666482
%N Dimensions where the volume of an L^p unit ball is maximized.
%C The volume of an n-dimensional L^p-ball of radius R is vol(R, p, n) = ((2^n)*(gamma((1/p)+1))^n /gamma((n/p)+1)*R^n for all real p > 0. For the Hilbert space L^2 we recover the familiar Euclidean volume of a unit n-ball (Pi^(n/2))/gamma((n/2)+1) which attains a maximum value at dimension n=5. The present sequence describes the same phenomenon for positive integer values of p.
%C From _Robert L. Diersing_, May 26 2020: (Start)
%C We can find very tight bounds for a(n):
%C define c := 2*gamma(1 + 1/p), then
%C floor(p*(c^p - 1/2) - 1/2) <= a(n) <= ceiling(p*(c^p - 1/2)).
%C First we treat V(n, p, 1) as a real-valued smooth function w.r.t. n, then we can maximize it. Since log(x) is an increasing function, the critical points of log(V(n, p)) will be the same as the critical points for V(n, p).
%C (d/dn)log(V(n, p)) = log(c) - 1/(p(digamma((n+p)/p))).
%C (d^2/dn^2)log(V(n, p)) = -polygamma(1, (n+p)/p))/p^2 < 0 for all n > 0, p > 0.
%C This is because polygamma(1, x) = Sum_{k>=1} 1/(x + k)^2 > 0 for all x > 0.
%C Therefore (d/dx)V(n, p) is a decreasing function w.r.t. n, which means there is only one local extremum. We can bound that local extremum similar to the way that we find bounds for the case where p=2 and the radius varies: see Wikipedia link.
%C Let's call this extremum n_.
%C dV/dn = log(c) - 1/(p*(digamma((n_+p)/p))) = 0;
%C p*log(c) = digamma((n_ + p) / p);
%C p*(inverse_digamma(plog(c)) - 1) = n_.
%C We can see from [2] that
%C log(y - 1/2) < digammma(y) < log(y + e^-m - 1),
%C where m is the Euler-Mascheroni constant.
%C Thus
%C log((n_+ p)/p - 1/2) < p*log(c) < log((n_ + p)/p + e^-m - 1);
%C c^p - e^-m + 1 < (n_ + p)/p < c^p + 1/2;
%C p*(c^p - e^-m) < n_ < p*(c^p - 1/2).
%C We can now improve on these bounds; define:
%C n` = p*(c^p - 1/2).
%C We can prove
%C n` - n_ < 1/2;
%C p*(c^p - 1/2 - inverse_digamma(p*log(c))) < p/c^p < 1/2;
%C c^p + 1/2 - 1/c^p < inverse_digamma(log(c^p)), say x := c^p;
%C x + 1/2 - 1/x < inverse_digamma(log(x));
%C digamma(x + 1/2 - 1/x) < log(x)
%C (inverse_digamma is an increasing function).
%C We see from [1] that digamma(1/log(1 + 1/x)) < log(x), therefore we only need to show
%C digamma(x + 1/2 - 1/x) < digamma(1/log(1 + 1/x)),
%C and since inverse_digamma is an increasing function, it suffices to show
%C x + 1/2 - 1/x < 1/log(1 + 1/x),
%C which is true for x > 0, and can be seen from an inequality solver such as Wolfram Alpha: https://www.wolframalpha.com/input/?i=x+%2B+%C2%BD+-+1%2Fx+%3C+1%2Flog%281+%2B+1%2Fx%29.
%C thus
%C n` - 1/2 < n_ < n`.
%C Now we have to confirm that n_ is the global maximum by checking the boundaries.
%C V(p, 0) = c^p < n_;
%C Limit_{n->oo} V(p, n) = lim_{n->oo} (2*sqrt(2*Pi/p))^n / (sqrt(2*Pi*n/p) n^n) = 0 < n_.
%C Therefore n_ is the global maximum, and thus we can solve for a(n).
%C floor(n` - 1/2) <= a(n) <= ceiling(n`).
%C So to find a(n) we only need to iterate over at most three values. (End)
%H Necdet Batir, <a href="https://arxiv.org/abs/1705.06547">Inequalities for the Inverses of the Polygamma Functions</a>, arXiv:1705.06547 [math.CA], 2017.
%H N. Elezovic, C. Giordano and J. Pecaric, <a href="https://dx.doi.org/10.7153/mia-03-26">The best bounds in Gautschi's inequality</a>, Math. Inequal. Appl. 3 (2000), 239-252.
%H Wikipedia, <a href="http://en.wikipedia.org/wiki/Volume_of_an_n-ball#Balls_in_Lp_norms">Balls in L^p norms</a>
%F From _Robert L. Diersing_, May 26 2020: (Start)
%F Define c = 2*gamma(1 + 1/p), and n` = p*(c^p - 1/2)
%F for x in [floor(n` - 1/2), floor(n` + 1/2), floor(n` + 3/2)]
%F a(n) = the value x that gives the largest volume. (End)
%e a(2) = 5 because the sequence vol(R=1, p=2, n) increases monotonically up to n=5 and then decreases monotonically.
%e The approximate values for vol(1,p,a(p)) are 2, 5.263789015, 50.05958637, 5970.510613, p = 1..4. - _Wolfdieter Lang_, Mar 20 2015
%t f[p_] = Round[FindArgMax[((2^n)*(Gamma[1+(1/p)])^n)/Gamma[1+(n/p)], n]]
%K nonn
%O 1,2
%A _Dimitris Cardaris_, Feb 14 2015
%E Edited by _Wolfdieter Lang_, Mar 20 2015
%E More terms from _Robert L. Diersing_, May 26 2020
| 