# Apply central difference formula to a variable twice

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?

• perhaps DifferenceQuotient[f[x - \[CapitalDelta]], {x, 2, \[CapitalDelta]}]?
– kglr
May 18 at 13:06

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

• Please could you explain why you define the function like ct@D[expr_, x_] instead of just ct[expr_, x_]? Thank you May 18 at 13:44
• @Hefaestion It's pretty OK to simply define ct[expr_, x_] as you've suggested, then the line SetAttributes[ct, HoldAll] can be removed. (I'd suggest defining it as ctD so one can modifying D in a easier manner, though. ) I initially designed it to be ct@D because I intuitively thought it would be better to keep the separate D. Subsequent practices suggest this doesn't seem to be necessary, but I'm too lazy to modify the code 囧. May 18 at 13:51
• @Hefaestion I've added a simplified implementation. Have a look. May 18 at 13:56

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}$$

• Also defining a subvalue makes it easy to extend to higher derivatives like diffQuot[h_][s_] := DifferenceQuotient[s, {x, h}] then Composition[diffQuot[h], diffQuot[-h], diffQuot[h], diffQuot[-h]]@f[x] May 18 at 13:41

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.

• I'm afraid this isn't what OP wants. The key point is to define a function that can generate difference formula nestedly. May 18 at 12:21
• @xzczd Agree, a bad design here is the main problem. The question formulation was somehow a misleading, since I thought he simply wants to manipulate notation.
– Acus
May 18 at 12:37