Reference Request: Variational Problem I want to solve approximately the following variational problem: 
Given a function $c:[-1,1]^2\rightarrow [0,1]$, constants $p_1...p_n\in \mathbb{R}^+$, $\alpha_1...\alpha_n\in \mathbb{R}$, and $\beta_1...\beta_n\in \mathbb{R}$ let 
$$V =\sup_{f_i,g_i}\left\{\sum_{i=1}^n \;\;p_i\cdot \int_0^1\int_0^1  f_i(x,y)\cdot g_i(x,y)\; dx\; dy \;\;: \sum_{i}\alpha_i\cdot f_i = c = \sum_i \beta_i\cdot g_i \right\}$$ 
where the supremum is taken over all square integrable functions 
$f_i,g_i:[-1,1]^2\rightarrow [0,1]$. 
Observe that since all functions have domain $[-1,1]^2$ and range 
$[0,1]$, the value of each integral $\int_{0}^1\int_{0}^1 f(x,y)\cdot g(x,y) \;dx\; dy$ lies between $0$ and $1$. Therefore $V$ lies between $0$ and $\sum_{i} p_i$. 
Assumptions:
All functions $f_i$ and $g_i$ and $c$ can be square integrable.
If it makes the problem simpler, $c:[-1,1]^2\rightarrow [0,1]$ can be assumed to be $c(x,y)=1$ if $x=y$ and $c(x,y)=0$ otherwise.  
Questions: 


*

*What numerical methods are available to approximate the value of $V$ up to an $\varepsilon$ additive factor? In other words, I want to find a $V'$ such that $|V-V'|\leq \varepsilon$, where $\varepsilon$ is a given precision parameter.

*What is the rate of convergence of such methods?  
 A: A bit long for a comment.
Let's clean up the formulation a bit:


*

*First, the domain $[-1,1]^2$ of definition does play any role, and hence, we assume that all respective quantities are functions on some $\Omega$ (and keep in mind that $\Omega=[-1,1]^2$ is the desired case). 

*Then we see that we are supremizing over $f = (f_1,\dots,f_n)\in L^2(\Omega)^n$ and $g=(g_1,\dots,g_n)\in L^2(\Omega)^n$, so in total over $(f,g)\in L^2(\Omega)^{2n}$ which is a nice Hilbert space.

*Now we introduce new variables $\tilde f_i = p_i f_i$ and observe that the objective is now $$\sum_i p_i\langle f_i,g_i\rangle = \sum_i \langle p_if_i,g_i\rangle = \sum_i \langle \tilde f_i,g_i\rangle = \langle \tilde f,g\rangle_{L^2(\Omega)^n}.$$

*Let's turn to the constraints: You have the bounds $0\leq f_i,g_i\leq 1$ which turn into $0\leq g_i\leq 1$ and $0\leq \tilde f_i\leq p_i$. The equality constraints are of the form $Af=c$ with
$$ A:L^2(\Omega)^n\to L^2(\Omega),\quad Af = \sum_i \alpha_i f_i$$
$$ B:L^2(\Omega)^n\to L^2(\Omega),\quad Bg = \sum_i \beta_i g_i.$$
This turns into
$$Af = \sum_i \alpha_if_i = \sum_i \tfrac{\alpha_i}{p_i}p_if_i = \sum_i \tfrac{\alpha_i}{p_i}\tilde f_i = \tilde A\tilde f.$$
With all this, the problem is
$$
\begin{array}{lrl}
\sup\limits_{(\tilde f,g)\in L^2(\Omega)^{2n}} \langle \tilde f,g\rangle_{L^2(\Omega)^n} & 0&\leq \tilde f_i\leq p_i,\\
 & 0&\leq g_i\leq 1,\\
& \tilde A\tilde f &= c,\\
& Bg &= c.
\end{array}
$$
The objective is not a convex function in $\tilde f$ and $g$ jointly but it's smooth. Hence, your problem classifies as a non-linear optimization problem with linear constraints.
