:

## Centered Divided Difference approximations

Maple

This code enables building Centered Divided Difference (CDD) approximations of derivatives of a univariate function.
Depending on the stencil we choose we can get arbitrary high order approximations.

The extension to bivariate functions is based upon what is often named tensorization in numerical analysis: for instance diff(f(x, y), [x, y] is obtained this way (the description here is purely notional)

1. Let CDD_x the CDD approximation of diff( f(x), x) ) .
CDD_x is a linear combination of shifted replicates of f(x)
2. Let s one of this shifted replicates
Let CDD_y(s) the CDD approximation of diff( s(y), y) ) .
3. Replace in CDD_x all shifted replicates by their corresponding expression CDD_y(s)

REMARKS:

• When I write for instance "approximation of diff(f(x), x)", this must be intended as a short for "approximation of diff(f(x), x) at point x=a"
• I agree that a notation like, for instance, diff(f(a), a) is not rigourous and that something like a Liebnitz notation would be better. Unfortunately I don't know how to get it when I use mtaylor.

 > restart:

CDDF stands for Cendered Divided Difference Formula

 > CDDF := proc(f, A, H, n, stencil)   description "f = target function,\nA = point where the derivatives are approximated,\nH = step,\nn = order of the derivative,\nstencil = list of points for the divided differenceCDDF\n";   local tay, p, T, sol, unknown, Unknown, expr:   tay := (s, m) -> convert(                      eval(                        convert(                          taylor(op(0, f)(op(1, f)), op(1, f)=A, m),                          Diff                        ),                        op(1, f)=A+s*H),                      polynom                    ):   p   := numelems(stencil):   T   := add(alpha[i]*tay(i, p+1), i in stencil):   T   := convert(%, diff):   if p > n+1 then     sol := solve([seq(coeff(T, h, i)=0, i in subsop(n+1=NULL, [\$0..p]))], [seq(alpha[i], i in stencil)])[];   else     sol := solve([seq(coeff(T, H, i)=0, i in subsop(n+1=NULL, [\$0..n]))], [seq(alpha[i], i in stencil)])[];   end if:   if `and`(is~(rhs~(sol)=~0)[]) then     WARNING("no solution found"):     return   else     unknown := `union`(indets~(rhs~(sol))[])[];     Unknown := simplify(solve(eval(T, sol) = Diff(op(0, f)(A), A\$n), unknown));     sol     := lhs~(sol) =~ eval(rhs~(sol), unknown=Unknown);     expr    := normal(eval(add(alpha[i]*op(0, f)(A+i*H), i in stencil), sol));   end if:   return expr end proc:
 > Describe(CDDF)
 # f = target function, # A = point where the derivatives are approximated, # H = # step, # n = order of the derivative, # stencil = list of points for the divided # differenceCDDF # CDDF( f, A, H, n, stencil )
 > # 2-point approximation of diff(f(x), x) | x=a CDDF(f(x), a, h, 1, [-1, 1]); convert(simplify(mtaylor(%, h=0, 4)), Diff);
 (1)
 > # 3-point approximation of diff(f(x), x\$2) | x=a CDDF(f(x), a, h, 2, [-1, 0, 1]); convert(simplify(mtaylor(%, h=0)), Diff);
 (2)
 > # 5-point pproximation of diff(f(x), x\$2) | x=a CDDF(f(x), a, h, 2, [\$-2..2]); convert(simplify(mtaylor(%, h=0, 8)), Diff);
 (3)
 > # 7-point approximation of diff(f(x), x\$2) | x=a CDDF(f(x), a, h, 2, [\$-3..3]); # simplify(taylor(%, h=0, 10)); convert(simplify(mtaylor(%, h=0, 10)), Diff);
 (4)
 > # 4-point staggered approximation of diff(f(x), x\$3) | x=a CDDF(f(x), a, h, 3, [seq(-3/2..3/2, 1)]); convert(simplify(mtaylor(%, h=0, 6)), Diff);
 (5)
 > # 6-point staggered approximation of diff(f(x), x\$3) | x=a CDDF(f(x), a, h, 3, [seq(-5/2..5/2, 1)]); # simplify(taylor(%, h=0, 8)); convert(simplify(mtaylor(%, h=0, 8)), Diff);
 (6)
 > # 5-point approximation of diff(f(x), x\$4) | x=a CDDF(f(x), a, h, 4, [\$-2..2]); convert(simplify(mtaylor(%, h=0, 8)), Diff);
 (7)
 > # 7-point approximation of diff(f(x), x\$4) | x=a CDDF(f(x), a, h, 4, [\$-3..3]); convert(simplify(mtaylor(%, h=0, 10)), Diff);
 (8)

A FEW 2D EXTENSIONS

 > # diff(f(x, y), [x, y]) approximation over a (2 by 2)-point stencil stencil := [-1, 1]: # step 1: approximate diff(f(x, y), x) over stencil < stencil > fx  := CDDF(f(x), a, h, 1, stencil): fx  := eval(% , f=(u -> f[u](y))): ix  := [indets(fx, function)[]]: # step 2: approximate diff(g(y), y) over stencil < stencil > where #         g represents any function in fx. fxy := add(map(u -> CDDF(u, b, k, 1, stencil)*coeff(fx, u), ix)): # step 3: rewrite fxy in a more convenient form [seq(u=op([0, 0], u)(op([0, 1], u), op(1, u)), u in indets(fxy, function))]: fxy := simplify( eval(fxy, %) ); convert(mtaylor(fxy, [h=0, k=0]), Diff)
 (9)
 > # Approximation of diff(f(x, y), [x, x, y, y] a (3 by 3)-point stencil stencil := [-1, 0, 1]: # step 1: approximate diff(f(x, y), x) over stencil < stencil > fx  := CDDF(f(x), a, h, 2, stencil): fx  := eval(% , f=(u -> f[u](y))): ix  := [indets(fx, function)[]]: # step 2: approximate diff(g(y), y) over stencil < stencil > where #         g represents any function in fx. fxy := add(map(u -> CDDF(u, b, k, 2, stencil)*coeff(fx, u), ix)): # step 3: rewrite fxy in a more convenient form [seq(u=op([0, 0], u)(op([0, 1], u), op(1, u)), u in indets(fxy, function))]: fxy := simplify( eval(fxy, %) ); convert(mtaylor(fxy, [h=0, k=0], 8), Diff)
 (10)
 > # Approximation of diff(f(x, y), [x, x, y] a (3 by 2)-point stencil stencil_x := [-1, 0, 1]: stencil_y := [-1, 1]: # step 1: approximate diff(f(x, y), x) over stencil < stencil > fx  := CDDF(f(x), a, h, 2, stencil_x): fx  := eval(% , f=(u -> f[u](y))): ix  := [indets(fx, function)[]]: # step 2: approximate diff(g(y), y) over stencil < stencil > where #         g represents any function in fx. fxy := add(map(u -> CDDF(u, b, k, 1, stencil_y)*coeff(fx, u), ix)): # step 3: rewrite fxy in a more convenient form [seq(u=op([0, 0], u)(op([0, 1], u), op(1, u)), u in indets(fxy, function))]: fxy := simplify( eval(fxy, %) ); convert(mtaylor(fxy, [h=0, k=0], 6), Diff)
 (11)
 > # Approximation of the laplacian of f(x, y) stencil := [-1, 0, 1]: # step 1: approximate diff(f(x, y), x) over stencil < stencil > fx  := CDDF(f(x), a, h, 2, stencil): fy  := CDDF(f(y), b, k, 2, stencil): fxy := simplify( eval(fx, f=(u -> f(u, b))) + eval(fy, f=(u -> f(a, u))) ); convert(mtaylor(fxy, [h=0, k=0], 6), Diff)
 (12)
 >