3
$\begingroup$

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]
$\endgroup$
2
  • $\begingroup$ 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. $\endgroup$ – xzczd May 18 '20 at 3:08
  • 1
    $\begingroup$ 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. $\endgroup$ – J. M.'s ennui May 18 '20 at 3:13
4
$\begingroup$

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"]

Mathematica graphics

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"]

Mathematica graphics

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

Mathematica graphics

Second order

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

Mathematica graphics

More points, gives better approximations

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

Mathematica graphics

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

Mathematica graphics

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

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

Mathematica graphics

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.

Mathematica graphics

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.

$\endgroup$
5
  • 1
    $\begingroup$ Ah, that FDFormula[] function is using Fornberg's method. I wrote a little about it here. $\endgroup$ – J. M.'s ennui May 18 '20 at 4:09
  • $\begingroup$ 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. $\endgroup$ – Nasser May 18 '20 at 4:47
  • $\begingroup$ Oops, sorry about that. I fixed the link in my first comment, please try it now. $\endgroup$ – J. M.'s ennui May 18 '20 at 4:51
  • $\begingroup$ 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. $\endgroup$ – zephyrus May 18 '20 at 19:13
  • 1
    $\begingroup$ @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. $\endgroup$ – J. M.'s ennui May 19 '20 at 3:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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