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:
|
# 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) |
|