Fix a point $(p_1, p_2) \in (0,1)^2$ and let $\epsilon \in (0,1)$.
Consider now the following region
\begin{equation} \mathcal{D} = \{(x,y)\ :\ (x,y) \in (0,1)^2,\ n_1 x^i + n_2 y^i - (n_1 p_1^i + n_2 p_2^i) \leq \epsilon, \ \text{for all }i \in \mathbb{N}\} \end{equation}
Question 1
I would like to find the area of $\mathcal{D}$ as a function of $n_1, n_2, p_1, p_2, \epsilon$.
Question 2
Let \begin{equation} \mathcal{D}^* = \{(x,y)\ :\ (x,y) \in (0,1)^2,\ n_1 x^{2^i} + n_2 y^{2^i} - (n_1 p_1^{2^i} + n_2 p_2^{2^i}) \leq \epsilon, \ \text{for all } i \in \mathbb{N}\} \end{equation} Since the constraints change slower for higher values of $i$, $\mathcal{D}^*$ seems to be a good approximation of $\mathcal{D}$. Obviously $\mathcal{D} \subset \mathcal{D^*}$. Is there a way to bound the difference \begin{equation} \mathrm{Area}(\mathcal{D}^*) - \mathrm{Area}(\mathcal{D})\ ? \end{equation}
For example let $n_1 = 100,\ n_2 = 95,\ p_1 = 0.8,\ p_2 = 0.86,\ \epsilon = 0.6$.
In the following plot we see the first $20$ constraints.
This is the region $\mathcal{D}$.
These are the constraints of $\mathcal{D^*}$ for the previous example.
And this is $\mathcal{D^*}$.
I would like to outline a finite (finite I that while for example there are not necessarily nice solutions maybe to required points of intersection, they can be found with newton's method for example) steps to answer question 1.
First to define the boundary.
$n_1x^i+n_2y^i-(n_1p_1^i+n_2p_2^i)=\epsilon$
Solve for $y$.
$y=\left( \frac{\epsilon +n_1p_1^i+n_2p_2^i-n_1x^i}{n_2} \right)^{\frac{1}{i}}$
Now the intersection of $i$th and $i+1$th curve we plug in $y$ above to the $i+1$ equation.
$n_1x^{i+1}+n_2\left( \frac{\epsilon +n_1p_1^i+n_2p_2^i-n_1x^i}{n_2} \right)^{\frac{i+1}{i}}-(n_1p_1^{i+1}+n_2p_2^{i+1})=\epsilon$
Now using assumption, that at some point the intersections will be beyond the region and these will not effect the area. Substitute $x=0$.
$n_2\left( \frac{\epsilon +n_1p_1^i+n_2p_2^i}{n_2} \right)^{\frac{i+1}{i}}-(n_1p_1^{i+1}+n_2p_2^{i+1})=\epsilon$
This equation is to be solved yielding a function $I(n_1,n_2,p_1,p_2)=i_{\Bbb R}$ which gives the exponent yielding intersections at the boundary. All $i\in \Bbb N \le i_{\Bbb R}$ must be integrated.
The area is then
$\sum_{i=2}^{i\lt i_{\Bbb R}}\int_{a_i}^{b_i}\left(\frac{\epsilon+n_1p_1^i+n_2p_2^i-n_1x^i}{n_2}\right)^{\frac{1}{i}}dx$
Where $a_i,b_i$ are the two corresponding intersection points of the $i$th and $i+1$th curve and the $i+1$th and $i+2$th curve. The most inner section and two outer sections can be handled in a slightly ad hoc manner compared to above.