# Replace derivative terms with finite difference symbols

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]

• 1. D[\[Phi]1, x, x] should be D[\[Phi]1[x, y], x, x]. 2. Replace[F, replacements] should be Replace[F, replacements, 1], check the document of Replace for more info. Alternatively, use ReplaceAll (/.) instead. 3. \[Delta]x should be \[Delta]x^2. – xzczd May 18 '20 at 3:08
• The following is a bit fragile; someone else will have to modify this to make it a proper answer: D[{ϕ1[x, y], ϕ2[x, y], ϕ3[x, y]}, {x, 2}] /. Derivative[ords__][f_][args__] :> Fold[(DifferenceQuotient[f[args], {#[[2]], #[[1]] - 1, h}] - DifferenceQuotient[f[args], {#[[2]], #[[1]] - 1, -h}])/h &, Transpose[{{ords}, {args}}]] // Simplify. – J. M.'s ennui May 18 '20 at 3:13

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 =
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.

• Ah, that FDFormula[] function is using Fornberg's method. I wrote a little about it here. – J. M.'s ennui May 18 '20 at 4:09
• thanks. Nice article. The link to Bengt Fornberg on your page there is gone also ;). This is why if I find something useful on the net, I download, because I know it might not be there next time I look. – Nasser May 18 '20 at 4:47
• Oops, sorry about that. I fixed the link in my first comment, please try it now. – J. M.'s ennui May 18 '20 at 4:51
• Thank you! This seems like a good start, but I would need to extend this in nontrivial ways to account for the fact my function is a fn of 2 variables $f(x,y)$ and need to create FD terms in both. – zephyrus May 18 '20 at 19:13
• @zephyrus, you could just apply it separately for each variable: apply Nasser's function to, say, x, and then use the result of that with Nasser's function again, but y is now the variable. – J. M.'s ennui May 19 '20 at 3:17