3
$\begingroup$

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)

cdfunc2

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.

cd2res

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)
 *)

cd2res2

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?

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

3 Answers 3

5
$\begingroup$

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

enter image description here

But who cares about this?

$\endgroup$
3
  • $\begingroup$ Please could you explain why you define the function like ct@D[expr_, x_] instead of just ct[expr_, x_]? Thank you $\endgroup$
    – Hefaestion
    May 18, 2023 at 13:44
  • 1
    $\begingroup$ @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 囧. $\endgroup$
    – xzczd
    May 18, 2023 at 13:51
  • 1
    $\begingroup$ @Hefaestion I've added a simplified implementation. Have a look. $\endgroup$
    – xzczd
    May 18, 2023 at 13:56
6
$\begingroup$

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, Δ}]

enter image description here

TraditionalForm @ %

enter image description here

TeXform @ %

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

$\endgroup$
1
  • 3
    $\begingroup$ 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] $\endgroup$ May 18, 2023 at 13:41
3
$\begingroup$

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.

$\endgroup$
2
  • $\begingroup$ I'm afraid this isn't what OP wants. The key point is to define a function that can generate difference formula nestedly. $\endgroup$
    – xzczd
    May 18, 2023 at 12:21
  • 2
    $\begingroup$ @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. $\endgroup$
    – Acus
    May 18, 2023 at 12:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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