Mismatch between theory and experiments

Question

Consider the Lagrangians (for the problem definition given below):

$$L_0=\int_{\mathbb{R}} mg^u{f}^{1-u}\mathrm{d}\mu+\int_{\mathbb{R}}\lambda_0(f-f_L)\mathrm{d}\mu+\mu_0\left(\int_{\mathbb{R}} f\mathrm{d}\mu-1\right)$$

$$L_1=\int_{\mathbb{R}} g^u{f}^{1-u}\mathrm{d}\mu+\int_{\mathbb{R}}\lambda_1(g-g_L)\mathrm{d}\mu+\mu_1\left(\int_{\mathbb{R}} g\mathrm{d}\mu-1\right)$$

Taking the Gateux derivatives of the Lagrangians, at the direction of funtions $\psi_0$ and $\psi_1$, respectively, leads to

$$\frac{\partial L_0}{\partial f}=\int\left((1-u)\left(\frac{g}{f}\right)^u+\lambda_0+\mu_0\right)\psi_0\mathrm{d}\mu=0\quad\quad\quad (1)$$

$$\frac{\partial L_1}{\partial g}=\int\left(u\left(\frac{g}{f}\right)^{u-1}+\lambda_1+\mu_1\right)\psi_1\mathrm{d}\mu=0\quad\quad\quad\quad\,\, (2)$$

Here according to Gateux derivative, $\psi_0$ and $\psi_1$ are arbitrary functions. I take them as integrable functions with $\int \psi_0 \mathrm{d}\mu=1$ and $\int \psi_1 \mathrm{d}\mu=1$. The conditions above $\partial L_0/\partial f=0$ and $\partial L_1/\partial g=0$ make sense only for the conditions $f>f_L$ and $g>g_L$.

Question: As a result of many experiments I know that the result of Gateux derivatives ($(1)$ and $(2)$ above) must give $g/f=k$, where $k$ is a constant. This also implies that $\lambda_*$ should be constant functions. However, I cannot see why this must be the case from the theory. Any ideas?

---

Problem:

$$\begin{array}{ll} \text{maximize}_{{f,g}} & \displaystyle\int_{\mathbb{R}} g^u{f}^{1-u}\mathrm{d}\mu\\ \text{sumubject to} & \displaystyle\int_{\mathbb{R}} f \mathrm{d}\mu= 1,\quad \displaystyle\int_{\mathbb{R}} g\mathrm{d}\mu =1 \\ & f \geq {f_L} \\ & g \geq g_L \end{array}$$ where $u\in(0,1) $ and $$\int_{\mathbb{R}}f_L \mathrm{d}\mu< 1,\quad\int_{\mathbb{R}}g_L \mathrm{d}\mu< 1$$

Here, $f$ and $g$ are distinct density functions, $f_L,g_L$ are some known positive functions on $\mathbb{R}$ and $\mu$ is the Lebesgue measure.

---

Here is the Mathematica code (for the discrete case), if you would like to see:

u = 0.5;
n = 30
a = Table[0.8*PDF[BinomialDistribution[30, 0.4], i], {i, n}];
c = Table[0.8*PDF[BinomialDistribution[30, 0.6], i], {i, n}];
X = Array[x, n];
Y = Array[y, n];

FindMaximum[{X^(1 - u).Y^u, 
  Flatten[{ Total[X] == 1 , Total[Y] == 1, 
    MapThread[#1 <= #2  &, {a, X}], 
    MapThread[#1 <= #2  &, {c, Y}]}]}  , Flatten[Join[{X, Y}]]]

The result will give a vector with $60$ values, the first $30$ values are $f$ and the other $30$ values are for $g$. Then one can see that $g/f$ is always constant in two different domains, which correspond to the Lagrangians above. The same is true for the continuous case but I do not have a code for this case.