Numerical integration of a force field I have the gradient of a function with three input variable: $f: \mathcal{R}_+^3 \rightarrow \mathcal{R}$. I want to calculate numerically $f$ from its gradient. The only boundary condition that I have is that it is 0 beyond $(a ,a, a)$. What is the best way to do it. It is a vector field in a three-dimensional domain. $\nabla f $ is numerically available. In one dimensional case, it is easy to do a simple integration using trapezoid rule, can I use it in this problem as well. Any resource would be great.
 A: You can apply the trapezoidal rule even in the multi-dimensional case. The key is the equation $$f(x) - f(y) = \int_0^1 \nabla f(tx+(1-t)y)\cdot(x-y)dt.$$ 
Here $f : \Omega \rightarrow \mathbb{R}$ is differentiable and $\nabla f : \Omega \rightarrow \mathbb{R}^n$ is continuous and $\Omega \subset \mathbb{R}^n$ is convex such that the line from $y$ to $x$ is contained in $\Omega$. With these assumptions, the function $\phi : [0,1] \rightarrow \Omega$ given by
$$ \phi(t) = f(tx + (1-t)y)$$
is differentiable with a continuous derivative given by the chain rule. We have
$$ \phi'(t) = \nabla f(tx+(1-t)y)\cdot(x-y),$$
where $\cdot$ denotes the Euclidian inner product. It follows that
$$ f(x) - f(y) = \phi(1) - \phi(0) = \int_0^1 \phi'(t)dt = \int_0^1 \nabla f(tx+(1-t)y)\cdot(x-y)dt.$$
With the ability to compute $\nabla f$ at any point along the line from $y$ to $x$ the trapezoidal rule can now be applied.
A: If you write down the values of $\nabla f$ as if they were obtained with a finite difference formula, you can relate the values of $f$ in a a grid point with values of $\nabla f$ and $f$ in surrounding points. So, you can start with a point in the boundary of the cube (where you know that $f=0$) and move backwards to the point where you need to compute $f$. For instance, knowing that $f(a, \frac{a}{10}, \frac{a}{10})=0$ and denoting $x_i = i \frac{a}{10}$, you can see that
\begin{cases}
f(x_{10}, \frac{a}{10}, \frac{a}{10})=0\\
f(x_i, \frac{a}{10}, \frac{a}{10})\approx f(x_{i+1}, \frac{a}{10}, \frac{a}{10})-\frac{a}{10} \frac{\partial f}{\partial x}(x_{i+1}, \frac{a}{10}, \frac{a}{10}), i=9,8, \cdots
\end{cases}
Repeating this procedure for every grid point in one face ($x=a$, $y=a$ or $z=a$) you can recover $f$ in the whole grid, and you don't even need all the components of the gradient. Using better finite difference approximations may require the solution of a linear system.
