Skip to content

Commit

Permalink
feat: add quadratic example
Browse files Browse the repository at this point in the history
  • Loading branch information
alexfikl committed Oct 9, 2023
1 parent 6a1d8cd commit 6d46143
Show file tree
Hide file tree
Showing 3 changed files with 68 additions and 2 deletions.
66 changes: 66 additions & 0 deletions examples/caputo-derivative-quadratic.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,66 @@
# SPDX-FileCopyrightText: 2023 Alexandru Fikl <alexfikl@gmail.com>
# SPDX-License-Identifier: MIT

import numpy as np

from pycaputo.derivatives import CaputoDerivative, Side
from pycaputo.differentiation import CaputoSpectralMethod, diff
from pycaputo.grid import make_jacobi_gauss_lobatto_points
from pycaputo.utils import Array

try:
import matplotlib # noqa: F401
except ImportError as exc:
raise SystemExit(0) from exc
else:
from pycaputo.utils import figure, set_recommended_matplotlib

set_recommended_matplotlib()


def f(x: Array) -> Array:
return x**2


with figure("caputo-derivative-quadratic.pdf") as fig:
ax = fig.gca()

alphas = [
0.1,
0.2,
0.3,
0.4,
0.5,
0.6,
0.7,
0.8,
0.9,
1.1,
1.2,
1.3,
1.4,
1.5,
1.6,
1.7,
1.8,
1.9,
]

p = make_jacobi_gauss_lobatto_points(32, a=0.0, b=2.0)

for alpha in alphas:
d = CaputoDerivative(order=alpha, side=Side.Left)
method = CaputoSpectralMethod(d)

df_num = diff(method, f, p)
ax.plot(p.x, df_num, color="k", alpha=0.2)

ax.plot(p.x, f(p.x), label=r"$\alpha = 0$")
ax.plot(p.x, 2.0 * p.x, label=r"$\alpha = 1$")
ax.plot(p.x, np.full(p.x.shape, 2.0), label=r"$\alpha = 2$")

ax.set_xlabel("$x$")
ax.set_ylabel(r"$D^\alpha_C[f](x)$")
ax.legend()

# }}}
2 changes: 1 addition & 1 deletion pycaputo/grid.py
Original file line number Diff line number Diff line change
Expand Up @@ -59,7 +59,7 @@ def make_stretched_points(
where :math:`A` is the *strength* and :math:`x_c = 1/2` for the midpoint.
:arg n: number of points in :math:`[a, b]`.
:arg strength: a positive number that constrols the clustering of points at
:arg strength: a positive number that controls the clustering of points at
the midpoint, i.e. a larger number denotes more points.
"""
x = np.linspace(0, 1, n)
Expand Down
2 changes: 1 addition & 1 deletion pycaputo/interpolation.py
Original file line number Diff line number Diff line change
Expand Up @@ -77,7 +77,7 @@ def determine_truncation_error(offsets: Array, x: float, h: float = 1.0) -> floa
f(x) - \sum_{k \in \text{offsets}} \ell_k(x) f_{m + k} =
c \frac{\mathrm{d}^n f}{\mathrm{d} x^n}(\xi),
where the constant :math:`c` is approxiated as the truncation error. The
where the constant :math:`c` is approximated as the truncation error. The
Error itself also depends on the derivatives of :math:`f` at a undetermined
point :math:`\xi` in the interpolation interval.
Expand Down

0 comments on commit 6d46143

Please sign in to comment.