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?