Apply central difference formula to a variable twice

Question

I am working with finite difference methods analytically and I would like to be able to perform operations on subscripted variables.

I would like to generate the following expression by applying a central difference function to a variable twice: $$\frac{y_{i+1} - 2 y_i + y_{i-1}}{\Delta x^2}$$

I have written the function:

CentralDiff[var_, i_] := (Subscript[var, i + 1] - Subscript[var, i - 1])/(
 2 Δx)

With this I can do the first derivative:

CentralDiff[y, 0]

Which returns: (-Subscript[y, -1] + Subscript[y, 1])/(2 Δx) The subscripts are correct.

If I try to get the 2nd derivative the :

CentralDiff[CentralDiff[y, 0], 0]

(*
(-Subscript[((-Subscript[y, -1] + Subscript[y, 1])/(2 Δx)), -1] + 
Subscript[(-Subscript[y, -1] + Subscript[y, 1])/(2 Δx), 1])/(2 Δx)
 *)

For the 2nd derivative using central differences the subscripts are not updated. I think I would have to define an operator for the subscripts perhaps? I don't know where to start with this. Perhaps there's an easier way?

Answer 1

You've made a bad design. The output of your CentralDiff is a complete difference formula, but the inputs of CentralDiff are just arguments of Subscript[…, …], this makes it hard to use the CentralDiff nestedly. Mixing up different coordinate system is also a bad idea in my view. (Notice $i-1$, $i$, $i+1$, … and $\Delta x$, $2 \Delta x$, $3 \Delta x$, … belong to 2 different coordinates! ) To achieve the goal, I'd recommend defining something like the following:

ClearAll[ct, delta]
SetAttributes[ct, HoldAll]

delta[a_ + b_] := delta@a + delta@b
delta[k_. delta[_]] := 0

ct@D[expr_, x_] := 
 Subtract @@ (expr /. {{x -> x + delta@x}, {x -> x - delta@x}})/(2 delta@x)

Now it can be used as follows (Notice the expected output you show in your question is wrong, unless you define your CentralDiff for half a grid):

ct@D[ct@D[f[x], x], x] // Simplify
(* (-2 f[x] + f[x - 2 delta[x]] + f[x + 2 delta[x]])/(4 delta[x]^2) *)

BTW I've used this method to generate difference formula in several posts:

Instability, Courant Condition and Robustness about solving 2D+1 PDE

Generate coefficient array from general formula of linear equation system

Free Convective Heat Transfer of Non-Newtonian Power Law Fluids from a Vertical Plate

3D FEM Vector Potential

---

Update

The ct above is coded several years ago, and nowadays it turns out to be unnecessarily advanced. The following is a simplified version:

ClearAll[ctD]
ctD[expr_, x_] := 
 With[{Δ = Δ@ToString@x}, ((expr /. x -> x + Δ) - (expr /. x -> x - Δ))/(2 Δ)]

ctD[ctD[f[x], x], x] // Simplify
(*
(-2 f[x] + f[x - 2 Δ["x"]] + f[x + 2 Δ["x"]])/(4 Δ["x"]^2)
 *)

Compared with ctD, the only advantage of ct I can think of is, we can use the 2D $\partial_x$ in notebook i.e. we can write something like

But who cares about this?

Answer 2

You can use DifferenceQuotient and Format as follows:

Format[f[y_], TraditionalForm] := 
    Subscript[f,  y /.  {x - Δ -> i - 1, x + Δ -> i + 1, x -> i}]

DifferenceQuotient[f[x - Δ], {x, 2, Δ}]

TraditionalForm @ %

TeXform @ %

$\huge \frac{-2 f_i+f_{i-1}+f_{i+1}}{\Delta ^2}$

Answer 3

Avoid sub/superscripts. Better define your own display/interpretation

MakeBoxes[CentralDiff[var_, i_], sf : StandardForm] := 
 TemplateBox[{MakeBoxes[var, sf], MakeBoxes[i, sf], 
   MakeBoxes[#, sf] & @@ {i + 1}, MakeBoxes[#, sf] & @@ {i - 1}, 
   MakeBoxes[2 \[CapitalDelta]x, sf]}, CentralDiff, 
  DisplayFunction :> (FractionBox[
      RowBox[{SubscriptBox[#1, #3], "-", 
        SubscriptBox[#1, #4]}], #5] &), 
  InterpretationFunction :> (RowBox[{"CentralDiff", "[", #1, ",", #2, 
       "]"}] &), SyntaxForm -> "fish", Tooltip -> None]

Then when you evaluate CentralDiff[y, 1], it will display nicely and you can interpret is as well.