Return to Answer

Initial answer

I have similar situation (but I numerically compute only first derivatives) and found that the most reliable and efficient solution is do not rely on ND but implement numerical derivatives by yourself using the definition

$$f'(a)=\lim_{h\to 0}\frac{f(a+h)-f(a)}{h}$$

Such approach requires at the first step finding optimal $h$ for each of the parameters which gives maximum precision. It cannot be too small because in this case the precision of the difference $f(a+h)-f(a)$ will vanish. In my case defining $h$ as

$$h=\frac{\log_2(|a|)}{2^n}$$

reduces the problem to finding optimal value of $n$ which is equal for all the parameters. For this purpose I compute numerical approximations for all partial derivatives with $n$ running from 15 to 60 (I compute with high WorkingPrecision) and see where the precision of these approximations becomes sufficient or maximal for all the partial derivatives.

UPDATE: efficient use of ND

Thanks to acl's great answer now it is clear that ND with Terms->1 uses exactly the same finite difference formula:

Needs["NumericalCalculus`"] ND[f[u], u, x, Terms -> 1, Scale -> h, WorkingPrecision -> Infinity] 

(-f[x] + f[h + x])/h

But one should be careful here: increasing the number of terms leads to division of the step $h$ by powers of 2 in successive higher order terms:

ND[f[u], u, x, Terms -> 2, Scale -> h, WorkingPrecision -> Infinity] // Simplify 

-((3 f[x] - 4 f[h/2 + x] + f[h + x])/h)

ND[f[u], u, x, Terms -> 3, Scale -> h, WorkingPrecision -> Infinity] // Simplify 

(-21 f[x] + 32 f[h/4 + x] - 12 f[h/2 + x] + f[h + x])/(3 h)

It means that for getting the best possible approximation with ND one should estimate the step $h$ by the method described in the initial section of this answer, then multiply it by $2^{t-1}$ (where $t$ is the value of the Terms option of ND) and use as value for the Scale option.

These considerations can be compiled into a function with more intuitive and consistent behavior:

Options[myND] = {Terms -> 1, WorkingPrecision -> MachinePrecision}; myND[expr_, x_, x0_, h_, opts : OptionsPattern[]] := ND[expr, x, x0, Scale -> h*2^(OptionValue[Terms]-1), opts] 

Initial answer

I have similar situation (but I numerically compute only first derivatives) and found that the most reliable and efficient solution is do not rely on ND but implement numerical derivatives by yourself using the definition

$$f'(a)=\lim_{h\to 0}\frac{f(a+h)-f(a)}{h}$$

Such approach requires at the first step finding optimal $h$ for each of the parameters which gives maximum precision. It cannot be too small because in this case the precision of the difference $f(a+h)-f(a)$ will vanish. In my case defining $h$ as

$$h=\frac{\log_2(|a|)}{2^n}$$

reduces the problem to finding optimal value of $n$ which is equal for all the parameters. For this purpose I compute numerical approximations for all partial derivatives with $n$ running from 15 to 60 (I compute with high WorkingPrecision) and see where the precision of these approximations becomes sufficient or maximal for all the partial derivatives.

UPDATE: efficient use of ND

Thanks to acl's great answer now it is clear that ND with Terms->1 uses exactly the same finite difference formula:

Needs["NumericalCalculus`"] ND[f[u], u, x, Terms -> 1, Scale -> h, WorkingPrecision -> Infinity] 

(-f[x] + f[h + x])/h

But one should be careful here: increasing the number of terms leads to division of the step $h$ by powers of 2 in successive higher order terms:

ND[f[u], u, x, Terms -> 2, Scale -> h, WorkingPrecision -> Infinity] // Simplify 

-((3 f[x] - 4 f[h/2 + x] + f[h + x])/h)

ND[f[u], u, x, Terms -> 3, Scale -> h, WorkingPrecision -> Infinity] // Simplify 

(-21 f[x] + 32 f[h/4 + x] - 12 f[h/2 + x] + f[h + x])/(3 h)

It means that for getting the best possible approximation with ND one should estimate the step $h$ by the method described in the initial section of this answer, then multiply it by $2^{t-1}$ (where $t$ is the value of the Terms option of ND) and use as value for the Scale option.

These considerations can be compiled into a function with more intuitive and consistent behavior:

Options[myND] = {Terms -> 1, WorkingPrecision -> MachinePrecision}; myND[expr_, x_, x0_, h_, opts : OptionsPattern[]] := ND[expr, x, x0, Scale -> h*2^(OptionValue[Terms]-1), opts] 

Initial answer

I have similar situation (but I numerically compute only first derivatives) and found that the most reliable and efficient solution is do not rely on ND but implement numerical derivatives by yourself using the definition

$$f'(a)=\lim_{h\to 0}\frac{f(a+h)-f(a)}{h}$$

Such approach requires at the first step finding optimal $h$ for each of the parameters which gives maximum precision. It cannot be too small because in this case the precision of the difference $f(a+h)-f(a)$ will vanish. In my case defining $h$ as

$$h=\frac{\log_2(|a|)}{2^n}$$

reduces the problem to finding optimal value of $n$ which is equal for all the parameters. For this purpose I compute numerical approximations for all partial derivatives with $n$ running from 15 to 60 (I compute with high WorkingPrecision) and see where the precision of these approximations becomes sufficient or maximal for all the partial derivatives.

UPDATE: efficient use of ND

Thanks to acl's great answer now it is clear that ND with Terms->1 uses exactly the same finite difference formula:

ND[f[u], u, x, Terms -> 1, Scale -> h, WorkingPrecision -> Infinity] 

(-f[x] + f[h + x])/h

But one should be careful here: increasing the number of terms leads to division of the step $h$ by powers of 2 in successive higher order terms:

ND[f[u], u, x, Terms -> 2, Scale -> h, WorkingPrecision -> Infinity] // Simplify 

-((3 f[x] - 4 f[h/2 + x] + f[h + x])/h)

ND[f[u], u, x, Terms -> 3, Scale -> h, WorkingPrecision -> Infinity] // Simplify 

(-21 f[x] + 32 f[h/4 + x] - 12 f[h/2 + x] + f[h + x])/(3 h)

It means that for getting the best possible approximation with ND one should estimate the step $h$ by the method described in the initial section of this answer, then multiply it by $2^{t-1}$ (where $t$ is the value of the Terms option of ND) and use as value for the Scale option.

These considerations can be compiled into a function with more intuitive and consistent behavior:

Options[myND] = {Terms -> 1, WorkingPrecision -> MachinePrecision}; myND[expr_, x_, x0_, h_, opts : OptionsPattern[]] := ND[expr, x, x0, Scale -> h*2^(OptionValue[Terms]-1), opts] 

Initial answer

I have similar situation (but I numerically compute only first derivatives) and found that the most reliable and efficient solution is do not rely on ND but implement numerical derivatives by yourself using the definition

$$f'(a)=\lim_{h\to 0}\frac{f(a+h)-f(a)}{h}$$

Such approach requires at the first step finding optimal $h$ for each of the parameters which gives maximum precision. It cannot be too small because in this case the precision of the difference $f(a+h)-f(a)$ will vanish. In my case defining $h$ as

$$h=\frac{\log_2(|a|)}{2^n}$$

reduces the problem to finding optimal value of $n$ which is equal for all the parameters. For this purpose I compute numerical approximations for all partial derivatives with $n$ running from 15 to 60 (I compute with high WorkingPrecision) and see where the precision of these approximations becomes sufficient or maximal for all the partial derivatives.

UPDATE: efficient use of ND

Thanks to acl's great answer now it is clear that ND with Terms->1 uses exactly the same finite difference formula:

Needs["NumericalCalculus`"] ND[f[u], u, x, Terms -> 1, Scale -> h, WorkingPrecision -> Infinity] 

(-f[x] + f[h + x])/h

But one should be careful here: increasing the number of terms leads to division of the step $h$ by powers of 2 in successive higher order terms:

ND[f[u], u, x, Terms -> 2, Scale -> h, WorkingPrecision -> Infinity] // Simplify 

-((3 f[x] - 4 f[h/2 + x] + f[h + x])/h)

ND[f[u], u, x, Terms -> 3, Scale -> h, WorkingPrecision -> Infinity] // Simplify 

(-21 f[x] + 32 f[h/4 + x] - 12 f[h/2 + x] + f[h + x])/(3 h)

It means that for getting the best possible approximation with ND one should estimate the step $h$ by the method described in the initial section of this answer, then multiply it by $2^{t-1}$ (where $t$ is the value of the Terms option of ND) and use as value for the Scale option.

These considerations can be compiled into a function with more intuitive and consistent behavior:

Options[myND] = {Terms -> 1, WorkingPrecision -> MachinePrecision}; myND[expr_, x_, x0_, h_, opts : OptionsPattern[]] := ND[expr, x, x0, Scale -> h*2^(OptionValue[Terms]-1), opts] 

I have similar situation (but I numerically compute only first derivatives) and found that the most reliable and efficient solution is do not rely on ND but implement numerical derivatives by yourself using the definition

$$f'(a)=\lim_{h\to 0}\frac{f(a+h)-f(a)}{h}$$

Such approach requires at the first step finding optimal $h$ for each of the parameters which gives maximum precision. It cannot be too small because in this case the precision of the difference $f(a+h)-f(a)$ will vanish. In my case defining $h$ as

$$h=\frac{\log_2(|a|)}{2^n}$$

reduces the problem to finding optimal value of $n$ which is equal for all the parameters. For this purpose I compute numerical approximations for all partial derivatives with $n$ running from 15 to 60 (I compute with high WorkingPrecision) and see where the precision of these approximations becomes sufficient or maximal for all the partial derivatives.

Initial answer

I have similar situation (but I numerically compute only first derivatives) and found that the most reliable and efficient solution is do not rely on ND but implement numerical derivatives by yourself using the definition

$$f'(a)=\lim_{h\to 0}\frac{f(a+h)-f(a)}{h}$$

Such approach requires at the first step finding optimal $h$ for each of the parameters which gives maximum precision. It cannot be too small because in this case the precision of the difference $f(a+h)-f(a)$ will vanish. In my case defining $h$ as

$$h=\frac{\log_2(|a|)}{2^n}$$

reduces the problem to finding optimal value of $n$ which is equal for all the parameters. For this purpose I compute numerical approximations for all partial derivatives with $n$ running from 15 to 60 (I compute with high WorkingPrecision) and see where the precision of these approximations becomes sufficient or maximal for all the partial derivatives.

UPDATE: efficient use of ND

Thanks to acl's great answer now it is clear that ND with Terms->1 uses exactly the same finite difference formula:

ND[f[u], u, x, Terms -> 1, Scale -> h, WorkingPrecision -> Infinity] 

(-f[x] + f[h + x])/h

But one should be careful here: increasing the number of terms leads to division of the step $h$ by powers of 2 in successive higher order terms:

ND[f[u], u, x, Terms -> 2, Scale -> h, WorkingPrecision -> Infinity] // Simplify 

-((3 f[x] - 4 f[h/2 + x] + f[h + x])/h)

ND[f[u], u, x, Terms -> 3, Scale -> h, WorkingPrecision -> Infinity] // Simplify 

(-21 f[x] + 32 f[h/4 + x] - 12 f[h/2 + x] + f[h + x])/(3 h)

It means that for getting the best possible approximation with ND one should estimate the step $h$ by the method described in the initial section of this answer, then multiply it by $2^{t-1}$ (where $t$ is the value of the Terms option of ND) and use as value for the Scale option.

These considerations can be compiled into a function with more intuitive and consistent behavior:

Options[myND] = {Terms -> 1, WorkingPrecision -> MachinePrecision}; myND[expr_, x_, x0_, h_, opts : OptionsPattern[]] := ND[expr, x, x0, Scale -> h*2^(OptionValue[Terms]-1), opts]