Algorithm for computing divergence of a vector field

I am using the 'divergence' command for calculating the divergence of a vector field in MATLAB. Does anyone know the method (any algorithm or paper) by which MATLAB does this? I need to mention the method in my report. Thanks!

Answer 1

[Roger Stafford]

Mathworks is unfortunately often reluctant to disclose the details of their algorithms, and I'm afraid this seems to be the case for the 'divergence' function. I see no mention of the method used in the documentation for this function.

However, with a little effort you can probably determine the mathematical nature of their algorithm. For the purposes of making the call

 div = divergence(X,Y,Z,U,V,W)

define X, Y, and Z with 'ndgrid' using monotone, non-uniformly spaced values for some reasonably small n by m by p size. Then first define U, V, and W linearly

 U = A*X + B*Y + C*Z + D
 V = E*X + F*Y + G*Z + H
 W = I*X + J*Y + K*Z + L

for some arbitrary coefficients A, B, C, ..., L. Then 'div' will surely be equal to the constant sum A+F+K throughout the entire field, including at the edges. (If not, you can request your money back.)

Next change U, V, W to that of quadratic functions:

 U = A*X.^2+B*Y.^2+C*Z.^2+D*Y.*Z+E*Z.*X+F*X.*Y+G*X+H*Y+I*J+K

and similarly with 20 more arbitrary coefficients in V and W. The exact divergence for these is of course

 2*A*X+E*Z+F*Y+G + corresponding terms from V and W

and can be compared with matlab's 'div'. If they are equal (to within a reasonable round-off error,) that means Mathworks has used at least a second order approximation for their partial derivatives. If not, they are content with a first order approximation. If the answer was yes and you are ambitious, you can go on to higher order functions of X, Y, and Z for U, V, W and eventually pin down exactly what order of approximations they are using. Be sure to check the values on the faces, edges, and corners.

Of course, if you are lucky perhaps one of Mathworks' representatives might eventually decide to reveal at least that much information and save you all that trouble.

Answer 2

[Stephan R. Kuberski]

Another way to get insight into the implemented method is to look at the source file:

edit divergence

From there you will find that computation makes solely use of the gradient function. Its (R2018b) documentation says: gradient calculates the central difference for interior data points. Although everything should be clear by that, you can also have a look at its implementation:

edit gradient