Numerical derivatives

It approximates the partial derivatives $\partial f/\partial x_k$ for $k=1,2,\ldots,N$ of a function $f$ of $N$ variables. If $x_k$ is in the interior, it uses the central difference: $$\frac{\partial f}{\partial x_k}(x_1,\dots,x_k,\ldots,x_N)\approx \frac{f(x_1,\ldots,x_k+h_k,\ldots,x_N)-f(x_1,\ldots,x_k-h_k,\ldots,x_N)}{2h_k}.$$ If $x_k$ is at the left boundary, it uses the one-sided difference: $$\frac{\partial f}{\partial x_k}(x_1,\ldots,x_k,\ldots,x_N)\approx \frac{f(x_1,\ldots,x_k+h_k,\ldots,x_N)-f(x_1,\ldots,x_k,\ldots,x_N)}{h_k}.$$ If $x_k$ is at the right boundary, it uses the one-sided difference: $$\frac{\partial f}{\partial x_k}(x_1,\ldots,x_k,\ldots,x_N)\approx \frac{f(x_1,\ldots,x_k,\ldots,x_N)-f(x_1,\ldots,x_k-h_k,\ldots,x_N)}{h_k}.$$ The above approximations will apply to the syntax forms below.

• F is an array of function values of $f$ evaluated on $N$ dimensional grid points.

• When F is a scalar or vector, it considers $f$ as a univariate function, and D1 contains the numerical derivatives $\partial f/\partial x_1$.
• When F is a matrix or has more dimensions, it considers $f$ as a multi-variate function, and D1 contains numerical derivatives $\partial f/\partial x_1$. The $x_1$ direction is the horizontal of F.
• h should be a scalar, specifying values of $h_1=h_2=\cdots=h_N$. It can be complex and is assumed 1 if h is not provided.

### [D1, D2, ..., DN] = gradient(F, h)

• F is an array of function values of $f$ evaluated on $N$ dimensional grid points.

• When F is a scalar or vector, it considers $f$ as a univariate function, and D1 contains the numerical derivatives $\partial f/\partial x_1$. An error will be raised if there are more than 1 output arguments.
• When F is a matrix or has more dimensions, it considers $f$ as a multi-variate function, and D1, D2, ..., DN contain the numerical derivatives $\partial f/\partial x_1,\partial f/\partial x_2,\ldots,\partial f/\partial x_N$, respectively. The $x_1$ and $x_2$ directions are the horizontal and vertical directions of F, respectively. An error will be raised if there are more than ndims(F) output arguments.
• h specifies values of $h_1=h_2=\cdots=h_N$. It can be complex, and is assumed 1 if h is not provided.

### [D1, D2, ..., Dk] = gradient(F, h1, h2, ..., hk)

• This is the same as [D1, D2, ..., Dk] = gradient(F) above, except that $h_1,h_2,\ldots,h_N$ are specified by h1, h1, ..., hN, respectively.

Example 1: Numerical derivatives of a linear function $f(x_1,x_2,x_3)=x_1+x_2+x_3$. The spacing is h = 0.6.

[X1,X2,X3]=ndgrid(-1:0.6:0.8);
F=X1+X2+X3;

D1(:, :, 1) =

1.000   1.000   1.000   1.000
1.000   1.000   1.000   1.000
1.000   1.000   1.000   1.000
1.000   1.000   1.000   1.000

D1(:, :, 2) =

1.000   1.000   1.000   1.000
1.000   1.000   1.000   1.000
1.000   1.000   1.000   1.000
1.000   1.000   1.000   1.000

D1(:, :, 3) =

1.000   1.000   1.000   1.000
1.000   1.000   1.000   1.000
1.000   1.000   1.000   1.000
1.000   1.000   1.000   1.000

D1(:, :, 4) =

1.000   1.000   1.000   1.000
1.000   1.000   1.000   1.000
1.000   1.000   1.000   1.000
1.000   1.000   1.000   1.000

D2(:, :, 1) =

1.000   1.000   1.000   1.000
1.000   1.000   1.000   1.000
1.000   1.000   1.000   1.000
1.000   1.000   1.000   1.000

D2(:, :, 2) =

1.000   1.000   1.000   1.000
1.000   1.000   1.000   1.000
1.000   1.000   1.000   1.000
1.000   1.000   1.000   1.000

D2(:, :, 3) =

1.000   1.000   1.000   1.000
1.000   1.000   1.000   1.000
1.000   1.000   1.000   1.000
1.000   1.000   1.000   1.000

D2(:, :, 4) =

1.000   1.000   1.000   1.000
1.000   1.000   1.000   1.000
1.000   1.000   1.000   1.000
1.000   1.000   1.000   1.000

D3(:, :, 1) =

1.000   1.000   1.000   1.000
1.000   1.000   1.000   1.000
1.000   1.000   1.000   1.000
1.000   1.000   1.000   1.000

D3(:, :, 2) =

1.000   1.000   1.000   1.000
1.000   1.000   1.000   1.000
1.000   1.000   1.000   1.000
1.000   1.000   1.000   1.000

D3(:, :, 3) =

1.000   1.000   1.000   1.000
1.000   1.000   1.000   1.000
1.000   1.000   1.000   1.000
1.000   1.000   1.000   1.000

D3(:, :, 4) =

1.000   1.000   1.000   1.000
1.000   1.000   1.000   1.000
1.000   1.000   1.000   1.000
1.000   1.000   1.000   1.000