How to find Jacobian of a vector field given as a dataset?

I have a dataset given in 3x 1 x m x n x l. 3x1 vector at each point of m x n x l 3D tensor holds the velocity vector [u,v,w]. How do I estimate the Jacobian matrix at each point?

[John D'Errico]

A Jacobian is composed of derivatives of some function. It presumes differentiability.

What do you have? A list of discrete points. No more than that. While you may think of this as some sort of functional relationship that could be differentiated, it is only a list of distinct, discrete points. So talking about the Jacobian of a point set is meaningless.

My guess is you may be trying to estimate the coefficients of some model on that dataset, where a Jacobian would now be meaningful. That requires you have a model for the relationship. It must be differentiable as a function of the parameters. If so, then you don't really need to even compute the Jacobian, if you are using tools like lsqnonlin or lsqcurvefit. But this is just a conjecture.

[Mirlan Karimov]

Hi John,

Thank you for your reply. The velocity field data may be a result of a simulation or experiment. I do not understand what you mean when you say "talking about the Jacobian of a point set is meaningless". Of course, numerical Jacobian is defined at each point to understand how the velocity field maps from point to point.

What I have done initially is to calculate the Jacobian as follows:

%Calculate Gradients
[gradxVx,gradyVx,gradzVx] = gradient(vx,stepx,stepy,stepz);
[gradxVy,gradyVy,gradzVy] = gradient(vy,stepx,stepy,stepz);
[gradxVz,gradyVz,gradzVz] = gradient(vz,stepx,stepy,stepz);


gradV = zeros(3,3,length(yy),length(xx),length(zz));
gradVT = zeros(3,3,length(yy),length(xx),length(zz));

for i=1:length(yy)
    for j=1:length(xx)
        for k=1:length(zz)
        gradV(1,1,i,j,k) = gradxVx(i,j,k);
        gradV(2,1,i,j,k) = gradxVy(i,j,k);
        gradV(3,1,i,j,k) = gradxVz(i,j,k);
        
        gradV(1,2,i,j,k) = gradyVx(i,j,k);
        gradV(2,2,i,j,k) = gradyVy(i,j,k);
        gradV(3,2,i,j,k) = gradyVz(i,j,k);
        
        gradV(1,3,i,j,k) = gradzVx(i,j,k);
        gradV(2,3,i,j,k) = gradzVy(i,j,k);
        gradV(3,3,i,j,k) = gradzVz(i,j,k);
        end
    end
end

[J. Alex Lee]

It seems you're just looking for the gradient of your velocity vector, not a Jacobian. I agree with everything John said, it is nonsensical to speak of a Jacobian of a data set.

You haven't indicated that you have any position data to take the gradient. I speculate from your data format and your attempt that your data is in m x n x l format because it is structured, and your nodal indices imply spatial coordinates.

In that case, other than making sure your indices are what you want (in the storage loops), looks like you've achieved your goal with above.

Your storage method I think will be pretty slow becuase indexing arrays of dim>2 seems to be slow on matlab. Do you really need your gradient field in this format?

[Mirlan Karimov]

Well, it is slow if the resolution is very high, but generally not that bad. Yes, I need my grad field in this format for other calculations that need to be done at each spatial node.

Well, then it is my bad about the Jacobian, but I am certain that I have seen in some literatures where Jacobian of a velocity field is given as the velocity gradient (or deformation gradient).

Actually, after John's comment I got confused an started thinking that the deformation gradient and velocity gradient are not the same. Because in finding velocity gradient we want to know how the velocity field changes within the domain, but in finding the deformation gradient we want to know how each spatial point in the domain changes during the deformation of the domain. Then, Jacobian is nothing but the linear approximation to the flow map that maps the domain from the undeformed state to the deformed state.

[J. Alex Lee]

that is true, the "jacobian of mapping" is a thing.

velocity and deformation are also not to be confused...in general node positions need not be governed strictly by convection.

"deformation gradient" might also just be a generally confusing phrase, because "gradient" (to me, maybe this is not generally true) implies "real" spatial, i.e., must be computed based on the deformed coordinates...so really, what exactly ought it mean?

By the way, relatedly, if your matrix structure is holding "reference" coordinates of nodes, but as a consequence of your simulation you also have a deformed "real" position coordinates associated with each node, then velocity gradient computed by your method will not be correct; you'll want (I assume) to compute gradients based on your actual positions in real space.

[Mirlan Karimov]

Alex, the domain is fixed.

Yes, I found that velocity gradient as a term is quite confusing. As you noted, the gradient well suits the situations where there is deformation and we measure at each point how the point moves.

Anyway, my method works correctly as I benchmarked it with a data set generated from an analytical velocity field. I just wanted to know if there is a better way of doing this.

[J. Alex Lee]

Great, glad that you have a working method!

Just to note that I don't think "velocity gradient" is confusing at all...I think "deformation gradient" could potentially be confusing.

Back to the actual question at hand:

I think your method is as good as can be unless you want to try to compute with higher accuracy finite differences (TMW's gradient() function uses central differences; docs also state that edge gradients are approximated by one-sided differences, which are only 1st order accurate). Added effort for higher "accuracy" would probably go to waste unless your velocities were computed using similarly accurate FDs (if you method is based on FD).

If you're extracting the nodal solutions from an FEM-based code (I've seen it before on Answers), you might be better off computing the gradients within the FEM framework and exporting it, if you have control over that.