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

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!

Answer 1

[Roger Stafford]

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.

[chi shing]

Many thanks again! So the non-zero numbers basically come from the round-off error instead of from the wrong syntax, right?

[Roger Stafford]

I would presume so.