MATLAB: How to compute DIV of CURL of vector ? Taking load wind as example

curldivdivergence of the curlMATLAB

I wanna to see whether the divergence of the curl of vector field is zero or not.
However, the result isn't zero. My script is below.
load wind; [cx,cy,cz] = curl(x,y,z,u,v,w]; ans = divergence(x,y,z,cx,cy,cz);
but ans is not zero. Suppose all vectors should be zero after this operation. Is there any mistake?
Many thanks!

Best Answer

That the divergence of the curl of a vector field is identically zero follows easily from the famous identity in calculus that
d^2F(x,y)/dx/dy = d^2F(x,y)/dy/dx
for a function F(x,y) provided these two second partial derivatives are continuous functions of x and y.
In your attempt to verify the above numerically, matlab is necessarily using discrete finite differences as approximations for these first and second derivatives. There is a limit, however, to how accurately matlab can approximate a second derivative. That is because the computation of a finite second difference involves the subtraction of large numbers which are close to one another in value so that the effect of round-off errors becomes relatively large.
As an example of this phenomenon, here are some second finite differences approximating the second derivative of sin(x) at x = pi/4. The exact second derivative as we all know is -sin(pi/4).
format long
[(sin(pi/4+.1)-2*sin(pi/4)+sin(pi/4-.1))/(.1)^2;
(sin(pi/4+.01)-2*sin(pi/4)+sin(pi/4-.01))/(.01)^2;
(sin(pi/4+.001)-2*sin(pi/4)+sin(pi/4-.001))/(.001)^2;
(sin(pi/4+.0001)-2*sin(pi/4)+sin(pi/4-.0001))/(.0001)^2;
(sin(pi/4+.00001)-2*sin(pi/4)+sin(pi/4-.00001))/(.00001)^2;
(sin(pi/4+.000001)-2*sin(pi/4)+sin(pi/4-.000001))/(.000001)^2;
(sin(pi/4+.0000001)-2*sin(pi/4)+sin(pi/4-.0000001))/(.0000001)^2]
-sin(pi/4)
ans =
-0.70651772191905
-0.70710088865056
-0.70710672239738
-0.70710679533903
-0.70710770572191
-0.70732308898869
-0.72164496600635
ans =
-0.70710678118655
Notice that the most accurate finite difference occurs for dx = .0001 . Larger values of dx incur errors because too large an increment is used, while smaller values show the increasing effect of round-off errors. For no increment do we get more than about seven-place accuracy, whereas the values in the sine function itself are good to some fifteen decimal places.