Crafting a replacement rule to convert derivatives to finite differences

Question

I need to construct a replacement rule to replace (first) derivatives with centered finite differences. What I've got so far is

DtoFD := f_'[var_] -> (f[var + ϵ/2] - f[var - ϵ/2])/ϵ

which works on simple cases like

f'[x] /. DtoFD
(* (-f[x - ϵ/2] + f[x + ϵ/2])/ϵ *)

Unfortunately I didn't consider the possibility of multivariate functions like

D[f[x, y], x] /. DtoFD
(* f^(1,0)[x, y] *)

where it fails.

Could someone come up with a rule that works in the multivariate case? First differences would be enough, but I'll also take arbitrary order if you've got it!

EDIT: I should add my motivation is to automatically set up a Jacobian matrix for a system of ordinary difference / differential equations, where some of the right hand sides might be numerically defined "black box" functions. So this substitution rule serves as a backup in case some derivatives can't be found with D first.

Answer 1

Aha, this question is for me.

diff = Simplify@*
   ReplaceAll[
    Derivative[o__][u_][x__] :> 
     d[u[x], ConstantArray @@@ ({{x}, {o}}\[Transpose]) // Flatten]]; 

d[expr_, {}] := expr
d[expr_, {var_, var2___}] := 
 1/ϵ[ToString@var] Subtract @@ (d[expr, {var2}] /. 
                                {{var -> var + ϵ[ToString@var]/2}, 
                                 {var -> var - ϵ[ToString@var]/2}})

D[f[x, y], y, y] // diff

D[f[x, y], x] + D[f[x, y], y] // diff

And as you know, pdetoode/pdetoae, which is a more complete discretizer, can be found here :) .

Answer 2

a rule that works in the multivariate case?

Why not use $\text{`NDSolve`FiniteDifferenceDerivative}$ ?

h = .5;
grid = h* Range[0, 5];
result = NDSolve`FiniteDifferenceDerivative[Derivative[1], grid,f[#, y] & /@ grid];
(Normal[result]) // Column

Gives

The above is the finite difference formula for derivative of $f(x,y)$ over $x$, along the $x$ direction at each grid point. There are 6 grid points in this example.

You could do the same also for 2D grid if you want, so that everything is done for you.

The above does not give a formula as you wanted, but the actual value of the approximation to the derivative on each grid point.

If this is not what you want, just let me know and will delete this.