I would some help on showing that the following integral is positive: Fix $1\leq p < \infty$ and let $(a,b) \in (0.5,1)^2$
$$ \int_0^1 \int_0^1 a(p+2)|(a,b)|^p \Big(\sqrt{(x-a)^2+(y-b)^2}\Big)^{p} - p|(a,b)|^{p+2}(a-x)\Big(\sqrt{(x-a)^2+(y-b)^2}\Big)^{p-2}dydx $$
Context: I am trying to show that the quotient $$\frac{|(a,b)|^{p+2}}{\int_0^1\int_0^1 |(x,y) - (a,b)|^pdydx}$$ achieves its maximum at $(a,b) = (1,1)$. So, I took its derivative with respect to $a$, which resulted in the integral above (numerator only), and now I'm trying to show it is positive, i.e. the quotient is increasing as $a$ goes to 1. Something analogous for $b$ will allow me to conclude.
I know that the maximum happens at $(1,1)$ for $(a,b) \in [0.5,1]^2$ by using a graphing calculator. See Here
You can't prove that the derivative w.r.t. $a$ is positive, because it has zeroes and regions where it is negative. It's still possible that your quotient has its maximum in the desired corner $(a,b)=(1,1)$. For an extremum somewhere inside the unit square we would need the derivatives w.r.t. $a$ and $b$ simultaneously zero, and I don't immediately see such points... But since the question was about the $d/da$ integral, we'll show here that it can become negative. This happens for $p$ larger than approximately $4.49$. We first split it into two parts: $$ I = I_1+I_2, \quad \mbox{with:} $$ $$ I_1 = a\ (p+2)\ |(a,b)|^p \ \int_0^1 \int_0^1\Big(\sqrt{(x-a)^2+(y-b)^2}\Big)^{p} dy\ dx $$ $$ I_2 = -p\ |(a,b)|^{p+2} \int_0^1 \int_0^1(a-x)\Big(\sqrt{(x-a)^2+(y-b)^2}\Big)^{p-2}dy\ dx $$
For $I_2$ we can write the intgrand as a derivative w.r.t. $x$ (because it obviously comes from a derivative w.r.t. $a$) which makes the $x$ integration trivial: $$ I2 = -p\ |(a,b)|^{p+2}\ \int_0^1 \int_0^1\frac{-d}{dx}\Big(\sqrt{(x-a)^2+(y-b)^2}\Big)^{p}dy\ dx $$ $$ = p\ |(a,b)|^{p+2}\ \left. \int_0^1\Big(\sqrt{(x-a)^2+(y-b)^2}\Big)^{p}dy \ \right|_{x=0}^1 $$ Now we can shift the variables and write: $$ I_1= a\ (p+2)\ |(a,b)|^p \int_{-a}^{1-a} dx \int_{-b}^{1-b} dy\ \left(x^2+y^2\right)^{p/2} $$ $$ I_2 = p\ |(a,b)|^{p+2} \left. \int_{-b}^{1-b} dy\ \left(x^2+y^2\right)^{p/2} \ \right|_{x=-a}^{1-a} $$ We first do the $y$ integral, which is needed for both $I_1$ and $I_2$: $$ I_y = \int_{-b}^{1-b} dy\ \left(x^2+y^2\right)^{p/2} = $$ $$ = | x|^{\ p} \left\{b\ \ _2F_1\Big(\frac{1}{2},-\frac{p}{2};\frac{3}{2};-\frac{b^2}{x^2}\Big) - (b-1) \ \ _2F_1\Big(\frac{1}{2},-\frac{p}{2};\frac{3}{2};-\frac{(b-1)^2}{x^2}\Big) \right\} $$ The hypergeometric functions are only needed for the general case with non-integer $p$. To get a simpler result we'll now look at $p=6$, which is in the problem region, and gives: $$ I_y^{(6)} =\left((b-1) b+x^2+1\right) \left(2 (b-1) b x^2+(b-1) b ((b-1) b+1)+x^4\right)+\frac{3 x^2}{5}+\frac{1}{7} $$ By evaluating this in $x=1-a$ and $x=-a$ we find $I_2$ as: $$ I_2^{(6)} = -6\ |(a,b)|^{8} \ \frac{2 a-1}{5} \Big[15 \Big((2 (a-1) a+3) b^2-2 (a-1) a b+ $$ $$ \quad\quad\quad\quad + (a-1) a ((a-1) a+2)+b^4-2 b^3\Big)-30 b+13\Big] $$ To get $I_1$ we do the $x$-integral, which for $p=6$ gives us: $$ I_1^{(6)}= 8 a\ |(a,b)|^6 \ \Big\{ a^6-3 a^5+3 a^4 b^2-3 a^4 b+6 a^4-6 a^3 b^2+6 a^3 b-7 a^3+ $$ $$ +3 a^2 b^4-6 a^2 b^3+12 a^2 b^2-9 a^2 b+\frac{28 a^2}{5}-3 a b^4+6 a b^3-9 a b^2+ $$ $$ + 6 a b-\frac{13 a}{5}+b^6-3 b^5+6 b^4-7 b^3+\frac{28 b^2}{5}-\frac{13 b}{5}+\frac{24}{35}\Big\} $$ Now we can plot the original $I=I_1+I_2$ for $p=6$:
In the plot, the blue plane is added to indicate zero, so it is clear that a whole region has a negative value. To be more precise we can look at the function at the $a=1$ edge, which gives just a 12th-order polynomial in $b$: $$ I^{(p=6)}_{a\rightarrow 1} = \frac{(b^2+1)^3}{35} \left(175 b^6-630 b^5+1260 b^4-1540 b^3+1162 b^2-518 b+101\right) $$ The negative region is created by the 6th order factor of this polynomial, it is the region where $0.544892 < b < 0.880778$. As mentioned, this failure to have positive $d/da$ everywhere does not rule out that we'll still find the maximum at $(a,b)=(1,1)$ for the original quotient of the question: $$ Q =\frac{|(a,b)|^{p+2}}{\int_0^1\int_0^1 |(x,y) - (a,b)|^pdydx}.$$ Here's a plot of $Q$ for $p=8$, where the derivative has an even larger negative region, but still the maximum appears to be in the $(1,1)$ corner:
PS: Derivation with Wolfram Engine of the $y$-integral goes as follows:
with the result:
And it can then be integrated w.r.t. $x$, using
Which allows plotting for any value of $p$, but gives a very long expression for the integral $I_{xy}$: $$ \frac{1}{(p+1) (p+2)}\ \left\{ b (p+1) a^{p+1} \ _2F_1\Big(\frac{1}{2},-\frac{p}{2};\frac{3}{2};-\frac{b^2}{a^2}\Big)+b a^{p+1} \ _2F_1\Big(\frac{1}{2} (-p-1),-\frac{p}{2};\frac{1-p}{2};-\frac{b^2}{a^2}\Big)-(b-1) (p+1) a^{p+1} \ _2F_1\Big(\frac{1}{2},-\frac{p}{2};\frac{3}{2};-\frac{(b-1)^2}{a^2}\Big)-(b-1) a^{p+1} \ _2F_1\Big(\frac{1}{2} (-p-1),-\frac{p}{2};\frac{1-p}{2};-\frac{(b-1)^2}{a^2}\Big)+b (p+1) (1-a)^{p+1} \ _2F_1\Big(\frac{1}{2},-\frac{p}{2};\frac{3}{2};-\frac{b^2}{(a-1)^2}\Big)+b (1-a)^{p+1} \ _2F_1\Big(\frac{1}{2} (-p-1),-\frac{p}{2};\frac{1-p}{2};-\frac{b^2}{(a-1)^2}\Big)+(a-1) (b-1) (p+1) (1-a)^p \ _2F_1\Big(\frac{1}{2},-\frac{p}{2};\frac{3}{2};-\frac{(b-1)^2}{(a-1)^2}\Big)+(a-1) (b-1) (1-a)^p \ _2F_1\Big(\frac{1}{2} (-p-1),-\frac{p}{2};\frac{1-p}{2};-\frac{(b-1)^2}{(a-1)^2}\Big)- 2 \sqrt{\pi } \left(b^{p+2}+(1-b)^{p+2}\right) \Gamma \Big(\frac{1}{2}-\frac{p}{2}\Big) \left/ \ \Gamma \Big(-\frac{p}{2}\Big) \right. \right\} $$ while the plot for $p=9.12345$ in the example looks more or less identical to the $p=8$ case from the earlier plot.