Numerically find a potential field from gradient I know that the gradient of a potential field/scalar field is a vector field, and given the function of the gradient I know how to integrate each component to get back the original scalar field. But is it possible to do this numerically? For example, if I know the value of the gradient as a vector at a finite number of points in 2 dimensions, how can I determine (or at least approximate) the value of the scalar potential function at each of those points?
 A: Let $G$ be a rectangular domain $[0, a] \times [0, b]$ with equally spaced grid $N+1 \times M+1$. Also let $h_x = \frac{a}{N}, h_y = \frac{b}{M}$. Also let the known graient be $\mathbf{f}_{i,j}, \; i = 0,\dots,N,\;j = 0,\dots,M$. We hope that $\mathbf{f} = (f_x,f_y) $ really is a conservative field.
Approach 1. Direct gradient approximation
For every point $(x_i, y_j)$ let's write
$$
\frac{\phi_{i+1,j}-\phi_{i,j}}{h_x} = \frac{(f_x)_{i,j}+(f_x)_{i+1,j}}{2},\quad i = 0,\dots,N-1,\;j = 0,\dots,M\\
\frac{\phi_{i,j+1}-\phi_{i,j}}{h_y} = \frac{(f_y)_{i,j}+(f_y)_{i,j+1}}{2},\quad i = 0,\dots,N,\;j = 0,\dots,M-1\\
$$
That's trapezoidal rule applied to
$$
\frac{\partial \phi}{\partial x} = f_x\\
\frac{\partial \phi}{\partial y} = f_y
$$
on every segment of your grid. There are too many equations ($2NM+N+M$) for $(N+1)(M+1)$ unknowns. That's because they are not independent and are consistent only if some condition imposed on $(f_x, f_y)$ is true. That condition is some difference approximation for conservativeness condition imposed on $\mathbf f$. Numerical integration with trapezoidal rule over any path consisting of grid segments should have the same value. If so, just drop extra equations and compute $\phi_{i,j}$ in following steps:
Since $\phi$ is defined up to an additive constant, fix $\phi$ at some point:
$$
\phi_{0,0} = 0.
$$
Compute $\phi_{i,0}$ for $i=1,\dots,N$ using the first difference equation running from $i = 1$ to $i = N$. For each $i$ now compute $\phi_{i,j}$ running from $j = 1$ to $j = M$ using the second difference equation.
The only drawback of this method is that you would obtain different values of $\phi_{i,j}$ if you choose another order of evaluation in the case when $\mathbf f$ fails to satisfy the difference conservative condition.
Approach 2. Poisson's equation
Taking divergence of $\nabla \phi = \mathbf f$ one obtain the Poisson's equation
$$
\Delta \phi = \operatorname{div} \mathbf f \text{ in } G.
$$
To make a well-possessed problem for Poisson's equation one need to impose some boundary conditions on $\partial G$. The natural conditions for that would be Neumann type conditions
$$
\frac{\partial \phi}{\partial n} = \mathbf n \mathbf f \text{ on } \partial G.
$$
This could be solved in many ways, but applying Galerkin numerical method gives us the
Approach 3. Weak formulation
Instead of solving $\nabla \phi = \mathbf f$ let's reformulate it in a weak way: find such $\phi$ that
$$
\int_G \nabla \phi \nabla v\, dx dy = \int_G \mathbf f \nabla v\, dx dy
$$
for every smooth function $v$ defined on $G$.
Let's take $\phi(x,y)$ as a continuous function, bilinear on every cell of the grid (that is $[x_i, x_{i+1}]\times [y_j, y_{j+1}]$). The values $\phi_{i,j}$ uniquelly define such function ($(N+1) \times (M+1)$ unknowns). Take $v(x,y)$ in exactly the same form. There are $(N+1)\times(M+1)$ linearly independend functions of this kind, thus giving us $(N+1)\times(M+1)$  equations to determine $\phi_{i,j}$. 
Let $\psi_{i,j}(x,y)$ be the function equal to $1$ at $(x_i, y_j)$ and equal to $0$ at every other grid point (but not at the every point, since $\psi_{i,j}(x,y)$ is continious. Note that $v_{i,j}$ support is $[x_{i-1},x_{i+1}]\times[y_{j-1},y_{j+1}]$, so integrating can be performed only on $[x_{i-1},x_{i+1}]\times[y_{j-1},y_{j+1}]$ instead of G. $\psi_{i,j}(x,y)$ is a finite element basis for $\phi(x,y)$ and $v(x,y)$, particulary
$$
\phi(x,y) = \sum_{i=0}^N\sum_{j=0}^M \phi_{i,j} \psi_{i,j}(x,y)
$$
Taking $\psi_{l,m}$ as different $v$ for weak formulation we obtain
$$
\sum_{i=0}^N\sum_{j=0}^M \phi_{i,j} \int_G \nabla\psi_{i,j} \nabla\psi_{l,m} dxdy =
\int_G \mathbf f \nabla\psi_{l,m} dxdy
$$
This is a huge system of linear equations for $\phi_{i,j}$. The matrix of the system, given by
$$
A_{(i,j),(l,m)} = \int_G \nabla\psi_{i,j} \nabla\psi_{l,m} dxdy
$$
has nonzero elements only if $(i,j)$ is not far from $(l,m)$, precisely $|i-l| \leq 1$ and $|j - m|\leq 1$. Thus the matrix is sparse. Elements of the matrix can be evaluated by hand or using some simple quadrature rule, precise for quadratic functions. Same is true for the right hand side.
The last to say is that the matrix is singular. That's because there is no unique solution to the problem, there are infinitely many. To resolve the problem one should fix some $\phi_{i,j}$, for example $\phi_{0,0} = 0$. That removes one unknown from the system and one equation at $(0,0)$. The remaining system is nonsingular.
