Standard Library function that computes the Jacobian of a function represented by n m-dimensional arrays

f—An n-element list of m-dimensional arrays all of the same dimensions d; f represents a n-valued function of m variables.

j—An n-element list of m-element lists of m-dimensional arrays: (j(p))(q) is the m-dimensional array (of dimensions d) which represents the derivitive of the pth dependent variable with respect to the qth independent variable.

x—An m-element list of vectors defining the independent variables; by default, x(i) = findgen( d(i) ).

An n-valued function of m variables can be represented by an n-element list of m-dimensional arrays all of the same dimensions di (i = 0 ... m-1). The mn partial derivatives are computed by the JACOBIAN function and returned in an n-element list of m-element lists of m-dimensional arrays of dimensions di. Each array of the input function represents a dependent variable that has been sampled on a grid in the independent variable space. Each output array represents one of the partial derivatives sampled on the same grid. The default grid is the uniform square-cell grid of array indices, while non-uniform rectangular-cell grids can be accommodated with the X keyword. General curvilinear grids can be accommodated via the chain rule of differential calculus. Consider for example a curvilinear grid of polar coordinates:

a = interpol([-1,1],500) 

r = sqrt(tensor_add(a*a,a*a))                    ; radius 

t = cprod(list(a,a)) 

t = reform(atan(t(*,1),t(*,0)),500,500)          ; angle

contour, r, a, a, nlev=9, po=[50,50,450,450], /devi 

contour, t, a, a, nlev=9, po=[50,50,450,450], /devi, /noer

And suppose the data:

z = hanning(500,500) 

shade_surf, z

has been sampled on the polar grid defined by arrays r and t. Derivatives with respect to the polar coordinates are obtained via the chain rule:

dz/dr = (dz/dx)(dx/dr) + (dz/dy)(dy/dr)   

dz/dt = (dz/dx)(dx/dt) + (dz/dy)(dy/dt) 

where x and y are most conveniently chosen to be array index coordinates.

dzdxy = jacobian(list(z)) 

drtdxy = jacobian(list(r,t)) 

dxydrt = inverse(drtdxy) 

dzdr = dzdxy(0,0)*dxydrt(0,0) + dzdxy(0,1)*dxydrt(1,0) 

dzdt = dzdxy(0,0)*dxydrt(0,1) + dzdxy(0,1)*dxydrt(1,1) 

tvscl, dzdr 

tvscl, dzdt

The key to this example is that even though z is sampled on a curvilinear grid defined by r and t, all three variables are related to the same set of array indices. The indices then serve as convenient parametric coordinates for the chain rule. This example was limited to two dimensions, but the same technique applies to higher dimensions: Derivatives with respect to curvilinear coordinates can always be computed by using the array indices as parametric coordinates and applying the chain rule.

f = LIST( RANDOMU(s,10,20), RANDOMU(s,10,20), RANDOMU(s,10,20) )

j = JACOBIAN( f )

FOR p=0L,2 DO FOR q=0L,1 DO PM, SAME( (j(p))(q), DERIVN(f(p),q) )

; PV-WAVE prints:

;   1

;   1

;   1

;   1

;   1

;   1

CURVATURES, DERIVN, EUCLIDEAN, NORMALS