I'm having trouble verifying
$$ \int_x^{x+\delta_x} \int_y^{y+\delta_y} f(u,v)\ dv \ dv \approx f(x,y)\delta_x \delta_y \ , $$
for small $\delta_x$ and $\delta_y$.
I'm not sure how to approach this problem. My first guess would be to use a 2 dimensional Taylor expansion in a way, but it seems that this approximation is made by an intuitive argument.
I find it reasonable to pick the value of $f$ at $(x,y)$, but how would one show the multiplication by $\delta_x \delta_y$
How would one derive this approximation being given that $\delta_x$ and $\delta_y$ are small?
Context.
I've come across this in studying continuous random variables, where $f(x,y)$ is the joint denisity function.
If $f$ is continuous at $r=(x,y)$ then for all $\epsilon > 0$ there exists a ball $B_r(\epsilon)$ so that for all $x\in B_r(\epsilon)$
$$ |f(x)-f(r)| < \epsilon $$
Now select $\delta_x,\delta_y$ so that $[x,x+\delta_x]\times[y,y+\delta_y]$ fit in that ball and you get $$ \int_x^{x+\delta_x}\int_y^{y+\delta_y} (f(u,v) - f(x,y) + f(x,y)) \,dudv = A + f(x,y)\delta_x\delta_y, $$ with $$ A = \int_x^{x+\delta_x}\int_y^{y+\delta_y} f(u,v) - f(x,y) \, du dv. $$ Hence we can estimate $|A|$ by $$ |A| \leq \int_x^{x+\delta_x}\int_y^{y+\delta_y} |f(u,v) - f(x,y)| \, du dv \leq \epsilon \delta_x \delta_y. $$ Hence the error can be made arbritrary small. You may want more control on the error on this equation if we assume all partial derivatives are uniformly bounded by C first differential then we can use with $e$ a unit vector of the correct direction and use $|(u,v)-(x,y)|=\delta$ $$ f(u,v) = f(x,y) + \int_0^\delta e \cdot \nabla f(r+et)\,dt $$
and
$$ |f(u,v) - f(x,y)| \leq 2 C \delta $$ Because $$ |e \cdot \nabla f|\leq \sum_i |e_i|C \leq 2 C $$ Then it's possible to plug this knowledge in the error term above and get a uniform bound on the error in the question.