Maximizing the fraction of two integrals using matlab

Question

Good day,

I've run into an issue in solving a certain problem with matlab, and I was hoping if anyone could help me out. I am relatively new to matlab, so I don't even know if this is possible, but I figured that in principle it should be.

What I want to do is maximize with respect to $a$

$$\frac{\displaystyle\int_a^\infty\frac{y^2}{e^y-1}\,\mathrm dy}{\displaystyle\int_a^\infty\frac{y^3}{e^y-1}\,\mathrm dy}$$

for positive values of $a$. With the calculus I know, I can't evaluate these integrals, so I have no applicable formula's so I am pretty stuck at this point. I tried using mathematica, and although it can evaluate the integrals, it cannot maximize the fraction. So instead I figured I should use matlab (which I also have to my disposal) to do so instead, but I simply don't know how.

I don't get much further than defining the integrals, as fun1 = @(x) x.^2./(exp(x)-1); fun2 = @(x) x.^3./(exp(x)-1);

Could anyone help me out?

Answer 1

Let the derivative of your ratio be $0$, and you find $a\int y^2/(e^y-1)dy=\int y^3/(e^y-1)dy$.
But since every piece of the second integral is $y$ times that piece of the first integral, and $y$ is always more than $a$, I don't think it is possible. So I think the ratio will be decreasing for all $y>0$.

Here is one way to find an expression for the integrals. For large values of $a$, $$\frac{y^2}{e^y-1}=e^{-y}y^2(1+e^{-y}+e^{-2y}+...)$$ and you should be able to integrate that by parts.

Answer 2

Let $f(a)$ be the expression you want to maximize. Then $$\begin{align} f'(a)&=\frac{-\dfrac{a^2}{e^a-1}\displaystyle\int_0^\infty\dfrac{y^3}{e^y-1}dy+\dfrac{a^3}{e^a-1}\displaystyle\int_0^\infty\dfrac{y^2}{e^y-1}dy}{\Bigl(\displaystyle\int_0^\infty\dfrac{y^3}{e^y-1}dy\Bigr)^2}\\ &=\frac{a^2}{e^a-1}\Bigl(\int_0^\infty\frac{y^3}{e^y-1}dy\Bigr)^{-2}\int_a^\infty\frac{y^2(a-y)}{e^y-1}\,dy\\ &<0. \end{align}$$ It follows that $f$ is decreasing and that the maximum is achieved at $a=0$. According to Mathematica, $$ f(0)=\frac{30}{4}\zeta(3)=0.370209\dots $$

Answer 3

If you must use MATLAB throw the biggest hammer at it: fmincon. Set your cost function up and as a constraint make $a\geq0$. Beware, this will be slow.

As a side note, have you looked at Leibniz's integral rule?