The goal of this notebook is to plot the function \((\lambda,\gamma)\mapsto A_k(\lambda,\gamma)\) defined in Lemma 4 of https://arxiv.org/abs/1906.08530
functionA <- function(k,lambda,gamma)
{
require(pracma)
library(pracma)
A <- sqrt(lambda - 1)/lambda
A <- A * (2*sqrt(lambda)/log(lambda-1))^k * k
gamma_arg <- 0.5*sqrt(gamma)*log(lambda-1)
A <- A * gammainc(gamma_arg,k)[[2]]
A <- A + (k * (gamma*lambda)^(0.5* k - 1) - 2)/(k-2)
}
We will now compute the values of \(A_4(\cdot,\cdot)\) on a sufficiently fine grid of the rectangle \((\lambda,\gamma)\in [12,18]\times [4,6]\).
m <- 30
gamma_vec <- 4 + (6-4)*((1:m)/m)
lambda_vec <- 12 + (18 - 12) * ((1:m)/m)
for (ncol in 1:m)
{
gamma = gamma_vec[ncol]
for (nrow in 1:m)
{
lambda = lambda_vec[nrow]
Matrix_A4[nrow,ncol] <- functionA(4, lambda, gamma)
}
}
We will now plot the corresponding surface.
lambda <- lambda_vec
gamma <- gamma_vec
A4 <- Matrix_A4
fig <- plot_ly(x = ~lambda, y = ~gamma, z = ~A4) %>% add_surface(
contours = list(
z = list(
show=TRUE,
start = 441,
end = 443, size = 0.3,
usecolormap=TRUE,
highlightcolor="#ff0000",
project=list(z=TRUE)
)
)
)
fig <- fig %>% layout(
scene = list(
camera=list(
eye = list(x=2, y=0.98, z=0.24)
)
)
)
fig
#print(functionA(4,15,5))
We do now the same for \(k=3\).
m <- 30
Matrix_A3 = matrix(0,m,m)
gamma_vec <- 4 + (6-4)*((1:m)/m)
lambda_vec <- 12 + (18 - 12) * ((1:m)/m)
for (ncol in 1:m)
{
gam = gamma_vec[ncol]
for (nrow in 1:m)
{
lamb = lambda_vec[nrow]
Matrix_A3[nrow,ncol] <- functionA(3, lamb, gam)
}
}
A3 <- Matrix_A3
fig <- plot_ly(x = ~lambda, y = ~gamma, z = ~A3) %>% add_surface(
contours = list(
z = list(
show=TRUE,
start = 40,
end = 41, size = 0.1,
usecolormap=TRUE,
highlightcolor="#ff0000",
project=list(z=TRUE)
)
)
)
fig <- fig %>% layout(
scene = list(
camera=list(
eye = list(x=2, y=0.98, z=0.24)
)
)
)
fig
The graphically determined nearly optimal values are computed here.
print(functionA(3,13,5))
[1] 40.73317
print(functionA(4,15,5))
[1] 441.4962
LS0tCnRpdGxlOiAiUGxvdHMgZm9yIGNvbXB1dGluZyB0aGUgY29uc3RhbnQgJEFfayQiCm91dHB1dDogaHRtbF9ub3RlYm9vawotLS0KClRoZSBnb2FsIG9mIHRoaXMgbm90ZWJvb2sgaXMgdG8gcGxvdCB0aGUgZnVuY3Rpb24gJChcbGFtYmRhLFxnYW1tYSlcbWFwc3RvIApBX2soXGxhbWJkYSxcZ2FtbWEpJCBkZWZpbmVkIGluIExlbW1hIDQgb2YgaHR0cHM6Ly9hcnhpdi5vcmcvYWJzLzE5MDYuMDg1MzAKCmBgYHtyfQpmdW5jdGlvbkEgPC0gZnVuY3Rpb24oayxsYW1iZGEsZ2FtbWEpCnsKICBpZiAoIXJlcXVpcmUoInByYWNtYSIpKSBpbnN0YWxsLnBhY2thZ2VzKCJwcmFjbWEiKQogIGxpYnJhcnkocHJhY21hKQogIEEgPC0gc3FydChsYW1iZGEgLSAxKS9sYW1iZGEKICBBIDwtIEEgKiAoMipzcXJ0KGxhbWJkYSkvbG9nKGxhbWJkYS0xKSleayAqIGsKICBnYW1tYV9hcmcgPC0gMC41KnNxcnQoZ2FtbWEpKmxvZyhsYW1iZGEtMSkKICBBIDwtIEEgKiBnYW1tYWluYyhnYW1tYV9hcmcsaylbWzJdXQogIEEgPC0gQSArIChrICogKGdhbW1hKmxhbWJkYSleKDAuNSogayAtIDEpIC0gMikvKGstMikKfQpgYGAKCldlIHdpbGwgbm93IGNvbXB1dGUgdGhlIHZhbHVlcyBvZiAkQV80KFxjZG90LFxjZG90KSQgb24gYSBzdWZmaWNpZW50bHkgCmZpbmUgZ3JpZCBvZiB0aGUgcmVjdGFuZ2xlICQoXGxhbWJkYSxcZ2FtbWEpXGluIFsxMiwxOF1cdGltZXMgWzQsNl0kLgpgYGB7cn0KCm0gPC0gMzAKTWF0cml4X0E0ID0gbWF0cml4KDAsbSxtKQpnYW1tYV92ZWMgPC0gNCArICg2LTQpKigoMTptKS9tKQpsYW1iZGFfdmVjIDwtIDEyICsgKDE4IC0gMTIpICogKCgxOm0pL20pCmZvciAobmNvbCBpbiAxOm0pCnsKICBnYW1tYSA9IGdhbW1hX3ZlY1tuY29sXSAKICBmb3IgKG5yb3cgaW4gMTptKQogIHsKICAgIGxhbWJkYSA9IGxhbWJkYV92ZWNbbnJvd10gCiAgICBNYXRyaXhfQTRbbnJvdyxuY29sXSA8LSBmdW5jdGlvbkEoNCwgbGFtYmRhLCBnYW1tYSkKICB9Cn0KYGBgCgpXZSB3aWxsIG5vdyBwbG90IHRoZSBjb3JyZXNwb25kaW5nIHN1cmZhY2UuCmBgYHtyfQppZiAoIXJlcXVpcmUoInBsb3RseSIpKSBpbnN0YWxsLnBhY2thZ2VzKCJwbG90bHkiKQpsaWJyYXJ5KHBsb3RseSkKbGFtYmRhIDwtIGxhbWJkYV92ZWMKZ2FtbWEgPC0gZ2FtbWFfdmVjCkE0IDwtIE1hdHJpeF9BNApmaWcgPC0gcGxvdF9seSh4ID0gfmxhbWJkYSwgeSA9IH5nYW1tYSwgeiA9IH5BNCkgJT4lIGFkZF9zdXJmYWNlKAogIGNvbnRvdXJzID0gbGlzdCgKICAgIHogPSBsaXN0KAogICAgICBzaG93PVRSVUUsCiAgICAgIHN0YXJ0ID0gNDQxLCAKICAgICAgZW5kID0gNDQzLCBzaXplID0gMC4zLAogICAgICB1c2Vjb2xvcm1hcD1UUlVFLAogICAgICBoaWdobGlnaHRjb2xvcj0iI2ZmMDAwMCIsCiAgICAgIHByb2plY3Q9bGlzdCh6PVRSVUUpCiAgICApCiAgKQopCmZpZyA8LSBmaWcgJT4lIGxheW91dCgKICBzY2VuZSA9IGxpc3QoCiAgICBjYW1lcmE9bGlzdCgKICAgICAgZXllID0gbGlzdCh4PTIsIHk9MC45OCwgej0wLjI0KQogICAgKQogICkKKQoKZmlnCiNwcmludChmdW5jdGlvbkEoNCwxNSw1KSkKYGBgCgpXZSBkbyBub3cgdGhlIHNhbWUgZm9yICRrPTMkLiAKYGBge3J9Cm0gPC0gMzAKTWF0cml4X0EzID0gbWF0cml4KDAsbSxtKQpnYW1tYV92ZWMgPC0gNCArICg2LTQpKigoMTptKS9tKQpsYW1iZGFfdmVjIDwtIDEyICsgKDE4IC0gMTIpICogKCgxOm0pL20pCmZvciAobmNvbCBpbiAxOm0pCnsKICBnYW0gPSBnYW1tYV92ZWNbbmNvbF0gCiAgZm9yIChucm93IGluIDE6bSkKICB7CiAgICBsYW1iID0gbGFtYmRhX3ZlY1tucm93XSAKICAgIE1hdHJpeF9BM1tucm93LG5jb2xdIDwtIGZ1bmN0aW9uQSgzLCBsYW1iLCBnYW0pCiAgfQp9CkEzIDwtIE1hdHJpeF9BMwpmaWcgPC0gcGxvdF9seSh4ID0gfmxhbWJkYSwgeSA9IH5nYW1tYSwgeiA9IH5BMykgJT4lIGFkZF9zdXJmYWNlKAogIGNvbnRvdXJzID0gbGlzdCgKICAgIHogPSBsaXN0KAogICAgICBzaG93PVRSVUUsCiAgICAgIHN0YXJ0ID0gNDAsIAogICAgICBlbmQgPSA0MSwgc2l6ZSA9IDAuMSwKICAgICAgdXNlY29sb3JtYXA9VFJVRSwKICAgICAgaGlnaGxpZ2h0Y29sb3I9IiNmZjAwMDAiLAogICAgICBwcm9qZWN0PWxpc3Qoej1UUlVFKQogICAgKQogICkKKQpmaWcgPC0gZmlnICU+JSBsYXlvdXQoCiAgc2NlbmUgPSBsaXN0KAogICAgY2FtZXJhPWxpc3QoCiAgICAgIGV5ZSA9IGxpc3QoeD0yLCB5PTAuOTgsIHo9MC4yNCkKICAgICkKICApCikKZmlnCmBgYAoKVGhlIGdyYXBoaWNhbGx5IGRldGVybWluZWQgbmVhcmx5IG9wdGltYWwgdmFsdWVzIGFyZSBjb21wdXRlZCBoZXJlLgpgYGB7cn0KcHJpbnQoZnVuY3Rpb25BKDMsMTMsNSkpCnByaW50KGZ1bmN0aW9uQSg0LDE1LDUpKQpgYGA=