A Problem of Mixed Numerical Derivatives

Michael Monagan, CECM@SFU

Given f(x), a function of x, the Maple diff command can often compute
a formula for f'(x).  E.g. for the Bessel function J_v(x), we have

> diff(BesselJ(v,x),x);
                                          v BesselJ(v, x)
                     -BesselJ(v + 1, x) + ---------------
                                                 x

But Maple cannot compute the derivative of J_v(x) wrt the parameter v.
A consequence of this is that Maple cannot graph this derivative.
In such cases numerical derivatives are necessary.
I'll describe how to compute f'(x) and f''(x) accurate to n digits
of precision based on the formulae 

     f'(x) = (f(x+h)-f(x-h))/2/h + O(h^2)

    f''(x) = (f(x+h) - 2 f(x) + f(x-h))/h^2 + O(h^2).

The problem, which I do not know how to solve, is how to do this for
mixed derivatives.  E.g. how to estimate the partial derivative
diff(f(x,y),x,y) at a point x=a,y=b.  I.e. D[1,2](f)(a,b) in Maple.

This work was motivated by professor Yury Brychkov's work where he
gave formulae for derivatives of special functions w.r.t. parameters 
but restricted to special cases of the parameter.