I have the following convex optimization problem: $$\begin{array}{ll} \text{maximize}_{{f,g}} & \displaystyle\int_{\mathbb{R}} g^u{f}^{1-u}\mathrm{d}\mu\\ \text{subject to} & \displaystyle\int_{\mathbb{R}} f \mathrm{d}\mu= 1,\quad \displaystyle\int_{\mathbb{R}} g\mathrm{d}\mu =1 \\ & f_L \leq {f} \leq f_U\\ & g_L \leq g \leq g_U\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$$
$$\int_{\mathbb{R}}f_U \mathrm{d}\mu> 1,\quad\int_{\mathbb{R}}g_U \mathrm{d}\mu> 1$$ Here, $f$ and $g$ are distinct density functions, $f_L,f_U,g_L,g_U$ are some known positive functions on $\mathbb{R}$ and $\mu$ is Lebesgue measure.
For illustration one can take:
$$f_L=0.8*f_{\mathcal{N}(-1,1)}$$ $$f_U=2*f_{\mathcal{N}(-1,1)}$$ $$g_L=0.8*f_{\mathcal{N}(1,1)}$$ $$g_U=2*f_{\mathcal{N}(1,1)}$$
where $f_{\mathcal{N}(\mu,\sigma^2)}$ is the Gaussian density with related parameters.
I want to solve this optimization problem but so far I have been unsuccesful. How can one solve this problem? Any programming approach is also highly appreciated. I am able to solve the discrete version of this problem in Mathematica for up to $35$ discrete values but I was nver able solve it for the continuous case. I tried to use Matlab, Mathematica as well as CVX tool for Matlab.
Analytically I attacked this question but there are some missing points. I think probably I need to write another question for the analytical case.
If you are okay with an approximation to the solution from a finite dimensional function space, then some of the software used to implement the finite element method (FEM) might be of use. This software is typically applied to solve PDEs in their weak form e.g find $u$ so that $a(u,v) = L(u)$ for all $v\in V$ where $V$ is some Hilbert space and $a$ is bilinear bounded coercive functional and $L$ is a bounded linear functional. If you can phrase your problem as equations on function spaces, then the FEM should be able to help you solve it.
In particular, I recommend Firedrake and Fenics. Both have python interfaces. To see an example of Lagrange multiplies implemented, check out this example. To summarize it, solving the pure neumann Poisson's equation requires you translate to it into a Lagrange multiplier problem and solve that using Fenics.
I will note that I see one possible flaw with this approach. I do not know how you solve a problem on the entire real line with FEM. Given your choice of example, I think that solving on a subset might require you to mess around a lot, but it could work.
Also, what exactly do you want from this? If you are looking for a precise formula, this is probably useless, but if you just want an approximation or a picture or some other quantity derived from $f$ and $g$ then this could be what you want.