Finite differences#

A FiniteDifference operator acts on a time-dependent quantity to produce a finite-difference discretization. For example, the second-order Adams-Bashforth discretization of \(u(\textbf{x}, t)\)

\[u(\textbf{x}, t^n)\approx \tfrac{3}{2}u^n - \tfrac{1}{2}u^{n-1}\]

is written using the AB2 operator applied to u as

AB2(u)

which is equivalent to manually writing out

1.5 * u[0] - 0.5 * u[-1]

Finite differences can also be applied argument-wise to an expression of time-dependent quantities. A FiniteDifferenceArgwise operator can be composed from individual FiniteDifference operators using @. For example, the discretization of \(f(a, b)\) with second-order Adams-Bashforth applied to \(a(\textbf{x}, t)\) and Crank-Nicolson applied to \(b(\textbf{x}, t)\) is written as

(AB2 @ CN)(f)

which is equivalent to manually writing out

(1.5 * a[0] - 0.5 * a[-1]) * (0.5 * b[1] + 0.5 * b[0])

Finite differences from the Adams-Bashforth, Adams-Moulton, backward-differential formulae families are provided ready to use, and the FiniteDifference base class permits further custom finite difference operators to be defined.

operating on a time-dependent quantity $\(\mathcal{D}(u(\textbf{x},t^n))=\sum_{j\leq 1}\alpha_ju^{n+j}(\textbf{x})\)\( \)\(\mathcal{D}(c(t^n))=\sum_{j\leq 1}\alpha_jc^{n+j}\)$

operating on a time-independent quantity $\(\mathcal{D}(u(\textbf{x}))=u(\textbf{x})\)\( \)\(\mathcal{D}(c)=c\)$

from lucifex.fdm import FunctionSeries, ExprSeries, FiniteDifference, DT, AB, AM, CN
from lucifex.fem import Constant
from lucifex.mesh import interval_mesh


mesh = interval_mesh(1.0, 10)
order = 3

a = FunctionSeries((mesh, 'P', 1), 'a', order=order)
for t in (0.0, 0.1, 0.2):
    a.update(2 * t)
    a.forward(t) 

print('\n Time derivative:')
dt = 0.001
print(DT(a, dt))
print(DT(a))
dt = Constant(mesh, 0.001, name='dt')
print(DT(a, dt))

print('\n Adams-Bashforth explicit discretizations:')
for n in range(1, order + 1):
    print(AB(n)(a))

print('\n Adams-Moulton implicit discretizations:')
for n in range(1, order + 1):
    print(AM(n)(a))
    print(AM(n)(a, False))

print('\n Custom discretization:')
fd_custom = FiniteDifference({0: 2.0, -1: -1.0}, init=AB(1))
print(fd_custom(a))
 Time derivative:
(v_1 + -1.0 * a⁰) / 0.001
v_1 + -1.0 * a⁰
(v_1 + -1.0 * a⁰) / dt

 Adams-Bashforth explicit discretizations:
a⁰
1.5 * a⁰ + -0.5 * a⁻¹
1.9166666666666667 * a⁰ + -1.3333333333333333 * a⁻¹ + 0.4166666666666667 * a⁻²

 Adams-Moulton implicit discretizations:
v_1
a⁺¹
0.5 * v_1 + 0.5 * a⁰
0.5 * a⁺¹ + 0.5 * a⁰
0.4166666666666667 * v_1 + 0.6666666666666666 * a⁰ + -0.08333333333333333 * a⁻¹
0.4166666666666667 * a⁺¹ + 0.6666666666666666 * a⁰ + -0.08333333333333333 * a⁻¹

 Custom discretization:
2.0 * a⁰ + -1.0 * a⁻¹
c = Constant(mesh, 0.1)

for n in range(1, order + 1):
    print(AM(n)(c))

Argument-wise application of finite differences#

\[(\mathcal{D}_0 \circ \mathcal{D}_1 \circ \dots)(f(u_0, u_1, \dots)) = f(\mathcal{D}_0(u_0), \mathcal{D}_1(u_1), \dots)\]
b = FunctionSeries((mesh, 'P', 1), 'b', order=order)
for t in (0.0, 0.1, 0.2):
    b.update(-2 * t)
    b.forward(t) 

f = ExprSeries.from_args(a, b, lambda a, b: a * b, name='f')

print(f)
'ExprSeries(a⁻² * b⁻², a⁻¹ * b⁻¹; Unsolved; Unsolved)'
fd = AB(2) @ CN

print(fd(f))
print(fd(f, trial=b))
(1.5 * a⁰ + -0.5 * a⁻¹) * (0.5 * b⁺¹ + 0.5 * b⁰)
(1.5 * a⁰ + -0.5 * a⁻¹) * (0.5 * v_1 + 0.5 * b⁰)