Numerical Gradient
Written by Alejandro Pozas-Kerstjens
ICFO-The Institute of Photonic Sciences
alejandro.pozas@icfo.es
December 2016

Documentation

Description
The function NGrad[expr,{x1,...,xn},{(Subscript[x1, 0]),...,(Subscript[xn, 0])}] computes numerically the gradient of the function expr[x1,...,xn] evaluated at the point {(Subscript[x1, 0]), ...,(Subscript[xn, 0])} by repeated uses of ND[expr,x,(Subscript[x, 0])]. For one-variable functions, NGrad[expr,x,(Subscript[x, 0])] computes numerically the derivative by using ND[expr,x,(Subscript[x, 0])].

Algorithm
For every input variable (indexed by q), the function evaluates the expression in {(Subscript[x1, 0]),...,xq,...,(Subscript[xn, 0])} (the input point substituting the corresponding variable), computes the numerical derivative of the resulting 1-variable expression in (Subscript[xq, 0]) and stores the result in the corresponding position of the vector grad.

Warnings
The function needs the NumericalCalculus package to operate
When first running, Mathematica may output the following error

ND::shdw: Symbol ND appears in multiple contexts {NumericalCalculus`, Global`}; definitions in context NumericalCalculus` may shadow or be shadowed by other definitions. 

This is due to the symbol ND also being part of Global`. To solve this, re-evaluate the definition of NGrad and the function will be usable without any further issues

Function

NGrad[expr_,x_,x0_]:=(
<<NumericalCalculus`
Clear[vars];
Clear[grad];
Clear[q];
Clear[onevar];
If[Length[Flatten[{x}]]==1,ND[Function[{x},expr]@@{x},x,x0],
vars=Delete[Flatten[Position[x,_?(Not[NumericQ[#]]&)]],1];
grad=Table[0,{i,1,Length[vars]}];
For[q=1,q<Length[vars]+1,q++,
onevar=x0;
onevar[[vars[[q]]]]=x[[vars[[q]]]];
grad[[q]]=ND[Function[x,expr]@@onevar,x[[vars[[q]]]],x0[[vars[[q]]]]]];
grad]
)

Examples

NGrad[Exp[-(x^(2y))]Log[3z+(2^t)],{x,y,z,t},{1,2,1,4}]
NGrad[Exp[-(x^(2y))]Log[3z+(2^t)]/.z->1,{x,y,t},{1,2,4}]

f[x_,y_,z_,t_]:=Exp[-(x^(2y))]Log[3z+(2^t)]
NGrad[f[x,y,z,t],{x,y,z,t},{1,2,1,4}]
NGrad[f[x,y,1,t],{x,y,t},{1,2,4}]

NGrad[Exp[-(x^2)],x,1]

Same code, but with dummyuglyvariables to avoid any kind of interference with other codes when being implemented to other notebooks

NGrad[expr_,x_,x0_]:=(
<<NumericalCalculus`
Clear[dummyuglyvariable1];
Clear[dummyuglyvariable2];
Clear[dummyuglyvariable3];
Clear[dummyuglyvariable4];
If[Length[Flatten[{x}]]==1,ND[Function[{x},expr]@@{x},x,x0],
dummyuglyvariable1=Delete[Flatten[Position[x,_?(Not[NumericQ[#]]&)]],1];
dummyuglyvariable2=Table[0,{i,1,Length[dummyuglyvariable1]}];
For[dummyuglyvariable3=1,dummyuglyvariable3<Length[dummyuglyvariable1]+1,dummyuglyvariable3++,
dummyuglyvariable4=x0;
dummyuglyvariable4[[dummyuglyvariable1[[dummyuglyvariable3]]]]=x[[dummyuglyvariable1[[dummyuglyvariable3]]]];
dummyuglyvariable2[[dummyuglyvariable3]]=ND[Function[x,expr]@@dummyuglyvariable4,x[[dummyuglyvariable1[[dummyuglyvariable3]]]],x0[[dummyuglyvariable1[[dummyuglyvariable3]]]]]];
dummyuglyvariable2]
)