5
$\begingroup$

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.

$\endgroup$

2 Answers 2

6
$\begingroup$

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

enter image description here

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

enter image description here

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

$\endgroup$
2
  • $\begingroup$ Very nice, thanks! I'll also have another look at pdetoode/pdetoae to see what else I forgot about. $\endgroup$
    – Chris K
    Oct 1, 2023 at 13:17
  • $\begingroup$ Also, nice inclusion of the variable-specific step-size! $\endgroup$
    – Chris K
    Oct 2, 2023 at 2:08
3
$\begingroup$

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

Mathematica graphics

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.

$\endgroup$
1
  • $\begingroup$ Not exactly what I need, but why don't you leave it in case this works for someone else in the future? $\endgroup$
    – Chris K
    Oct 1, 2023 at 13:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.