Advanced Search
All Posts : Applications, Examples and Libraries


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/2)*(f(a-h)-f(a+h))/h

  

Diff(f(a), a)+(1/6)*(Diff(Diff(Diff(f(a), a), a), a))*h^2

 (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);

-(-f(a-h)+2*f(a)-f(a+h))/h^2

  

Diff(Diff(f(a), a), a)+(1/12)*(Diff(Diff(Diff(Diff(f(a), a), a), a), a))*h^2

 (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);

-(1/12)*(f(a-2*h)-16*f(a-h)+30*f(a)-16*f(a+h)+f(a+2*h))/h^2

  

Diff(Diff(f(a), a), a)-(1/90)*(Diff(Diff(Diff(Diff(Diff(Diff(f(a), a), a), a), a), a), a))*h^4

 (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);

-(1/180)*(-2*f(a-3*h)+27*f(a-2*h)-270*f(a-h)+490*f(a)-270*f(a+h)+27*f(a+2*h)-2*f(a+3*h))/h^2

  

Diff(Diff(f(a), a), a)+(1/560)*(Diff(Diff(Diff(Diff(Diff(Diff(Diff(Diff(f(a), a), a), a), a), a), a), a), a))*h^6

 (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);

-(f(a-(3/2)*h)-3*f(a-(1/2)*h)+3*f(a+(1/2)*h)-f(a+(3/2)*h))/h^3

  

Diff(Diff(Diff(f(a), a), a), a)+(1/8)*(Diff(Diff(Diff(Diff(Diff(f(a), a), a), a), a), a))*h^2

 (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);

(1/8)*(f(a-(5/2)*h)-13*f(a-(3/2)*h)+34*f(a-(1/2)*h)-34*f(a+(1/2)*h)+13*f(a+(3/2)*h)-f(a+(5/2)*h))/h^3

  

Diff(Diff(Diff(f(a), a), a), a)-(37/1920)*(Diff(Diff(Diff(Diff(Diff(Diff(Diff(f(a), a), a), a), a), a), a), a))*h^4

 (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);

(f(a-2*h)-4*f(a-h)+6*f(a)-4*f(a+h)+f(a+2*h))/h^4

  

Diff(Diff(Diff(Diff(f(a), a), a), a), a)+(1/6)*(Diff(Diff(Diff(Diff(Diff(Diff(f(a), a), a), a), a), a), a))*h^2

 (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);

(1/6)*(-f(a-3*h)+12*f(a-2*h)-39*f(a-h)+56*f(a)-39*f(a+h)+12*f(a+2*h)-f(a+3*h))/h^4

  

Diff(Diff(Diff(Diff(f(a), a), a), a), a)-(7/240)*(Diff(Diff(Diff(Diff(Diff(Diff(Diff(Diff(f(a), a), a), a), a), a), a), a), a))*h^4

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

(1/4)*(f(a-h, b-k)-f(a-h, b+k)-f(a+h, b-k)+f(a+h, b+k))/(k*h)

  

Diff(Diff(f(a, b), a), b)+(1/6)*(Diff(Diff(Diff(Diff(f(a, b), a), a), a), b))*h^2+(1/6)*(Diff(Diff(Diff(Diff(f(a, b), a), b), b), b))*k^2

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

-(2*f(a, b-k)-4*f(a, b)+2*f(a, b+k)-f(a-h, b-k)+2*f(a-h, b)-f(a-h, b+k)-f(a+h, b-k)+2*f(a+h, b)-f(a+h, b+k))/(h^2*k^2)

  

Diff(Diff(Diff(Diff(f(a, b), a), a), b), b)+(1/12)*(Diff(Diff(Diff(Diff(Diff(Diff(f(a, b), a), a), a), a), b), b))*h^2+(1/12)*(Diff(Diff(Diff(Diff(Diff(Diff(f(a, b), a), a), b), b), b), b))*k^2

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

(1/2)*(2*f(a, b-k)-2*f(a, b+k)-f(a-h, b-k)+f(a-h, b+k)-f(a+h, b-k)+f(a+h, b+k))/(h^2*k)

  

Diff(Diff(Diff(f(a, b), a), a), b)+(1/12)*(Diff(Diff(Diff(Diff(Diff(f(a, b), a), a), a), a), b))*h^2+(1/6)*(Diff(Diff(Diff(Diff(Diff(f(a, b), a), a), b), b), b))*k^2

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

(f(a-h, b)-2*f(a, b)+f(a+h, b))/h^2+(f(a, b-k)-2*f(a, b)+f(a, b+k))/k^2

  

Diff(Diff(f(a, b), a), a)+Diff(Diff(f(a, b), b), b)+(1/12)*(Diff(Diff(Diff(Diff(f(a, b), a), a), a), a))*h^2+(1/12)*(Diff(Diff(Diff(Diff(f(a, b), b), b), b), b))*k^2

 (12)
> 

 


 

Download CDD.mw
derivative approximation mtaylor

Please Wait...