Replace derivative terms with finite difference symbols

Question

I would like to replace the derivative terms in an expression with their finite difference approximations. Note, I don't want to actually evaluate the approximations at a given point, just display the appropriate symbols. Explicitly, if I have a function:

$f[x, y]$ and an expression with terms like $f_{xx}$, how can I replace the differentials with the central difference representations $f_{xx} \rightarrow \frac{f[i + 1, j] - f[i - 1, j]}{(\Delta x)^2}$

I have tried to use Replace but even when I give an explicit rule, this does not seem to replace any terms.

ϕ = {ϕ1[x, y], ϕ2[x, y], ϕ3[x, y]};
F = D[ϕ, x, x];
replacements = {
   {D[ϕ1, x, x] -> (f[i + 1, j] - f[i - 1, j])/(δx)}
   };
Replace[F, replacements]

Answer 1

There used to be a nice FDFormula at WRI site, but it is gone now. But I used it before. Here is the result.

I'll show some examples then the code at end

 getFormula[1, {-1, 0, 1}, "centered"]

The first argument to getFormula is derivative order. So 1 for first order, 2 for second order. The second argument is list of points to generate the difference approximation on. The last argument is the type you want. Either centered, forward or backward.

The function returns the difference formula and also the error in the approximation (the big O).

Here are more examples

  getFormula[1, {-1, 0, 1}, "forward"]

  getFormula[1, {-1, 0, 1}, "backward"]

Second order

   getFormula[2, {-1, 0, 1}, "centered"]

More points, gives better approximations

   getFormula[2, {-2, -1, 0, 1, 2}, "centered"]

   getFormula[2, {-1, 0, 1}, "backward"]

4th order. Need to supply more grid points in this case, otherwise will get an error.

    getFormula[4, {-2, -1, 0, 1, 2}, "centered"]

Code

    (*FDFormula from 
    http://reference.wolfram.com/mathematica/tutorial/NDSolvePDE.html*)

FDFormula[(m_Integer)?Positive, (n_Integer)?Positive, (s_Integer)?
   NonNegative] := 
   Module[{do, F}, F = Table[f[Subscript[x, i + k]], {k, -s, n - s}]; 
      W = 
   PadRight[
    CoefficientList[Normal[Series[x^s*Log[x]^m, {x, 1, n}]/h^m], x], 
    Length[F], 0]; 
      Wfact = 1/PolynomialGCD @@ W; W = Simplify[W*Wfact]; 
      taylor[h_] = 
   Normal[Series[f[Subscript[x, i] + h], {h, 0, n + 2}]]; 
      error = Drop[CoefficientList[
     Expand[Table[taylor[h*k], {k, -s, n - s}] . W/Wfact], h], 1]; 
      do = Position[error, e_ /; e != 0][[1, 1]]; error = error[[do]]; 
      error = error /. (f_)[Subscript[x, i]] -> f; error = h^do*error; 
      {Derivative[m][f][Subscript[x, i]] \[TildeEqual] F . W/Wfact, 
   error}]

This uses the above function

getFormula[order_, gridPoints_, type_String] := Module[{s},
  s = Which[type == "centered", (Length[gridPoints] - 1)/2,
    type == "forward", 0,
    True, Length[gridPoints] - 1];
  Print[s];
  FDFormula[order, Length[gridPoints] - 1, s]
  ]

I've used this in past to make a detailed Manipulate. But I never send it to Wolfram demo site.

This Demonstration illustrates the effect of numerical errors on the approximation of derivatives when using the finite-difference scheme with different step sizes and different orders of accuracy. You can select to approximate up to the fourth derivative, the desired local truncation accuracy order O(h^n), and the finite difference scheme to use (centered, forward, or backward).

Also @xzczd has a finite difference formula generator function on this site. I do not have the link right now. That might also be something to look at. I've seen him use it to answer many questions.

ps. if you want to download the full Manipulate shown above, you can go to this page and search for "difference" and you'll find it near the top of the page there. One day I might submit to WRI demo site when I clean it a little more.