# Derivatives


```python
import numpy as np
from itertools import product
import matplotlib.pylab as plt
import scipy.special as spec
```

# Finite difference dispersion

There is a one to one corespondence between a finite step derivative and the dispersion of the corresponding momentum operator on the lattice.

E.g. for a Laplace finite step derivative of the form

$$
  \frac{\partial^2}{\partial^2x}f(x) \mapsto
  \frac{1}{\epsilon^2}\sum_{n=-N_s}^{N_s} c_n^{N_s} f(\epsilon n_x + \epsilon n) \, ,
  x = \epsilon n_x \, ,
$$
where $\epsilon$ is the lattice spacing, corresponds to the momentum dispersion
$$
    p^2 \mapsto D^{N_s}(p) = - \frac{1}{\epsilon^2}\sum_{n=-N_s}^{N_s} c_n^{N_s} \cos(n p \epsilon) \, .
$$

In general, one wants to identify the coefficients $c_n^{N_s}$ such that $D^{N_s}(p) = p^2 + \mathcal O \left( (\epsilon p)^{2 N_s} \right)$. 


If one now expands 
$$
    \cos(n p \epsilon) = \sum_{m=0}^\infty \frac{(-)^m}{(2m)!} (n p \epsilon)^{2m}\, ,
$$
one finds that the coefficients $c_n$ are determined by the liner equation
\begin{align}
    v_m &= \sum_{n=0}^{N_s} A_{m n} \gamma_n
    \, , & 
    A_{mn} &= \frac{(-)^m}{(2m)!} n^{2m} \, ,
    \\
    v_1 &= 1 
    \, , & 
    v_m &= 0 \, \forall \, m \neq 1
    \\
    c_0 &= - \gamma_0 
    \, , &
    c_{\pm n} &= - \frac{\gamma_n}{2} \, \forall \, n > 0
\end{align}
where $n$ and $m$ run from $0$ to $N_s$.



```python
def derivative_coeffs(Nstep: int):
    """Computes the derivative coefficient for the lapace operator up to step range Nstep.
    
    The dispersion of the related momentum squared operator is equal to p**2 up to order P**(2*NStep).
    
    **Arguments**
        Nstep: int
            The number of momentum steps in each direction.
    """
    v = np.zeros(Nstep + 1)
    v[1] = 1

    nn, mm = np.meshgrid(*[np.arange(Nstep + 1)] * 2)
    A = 1 / spec.gamma(2 * mm + 1) * (-1) ** mm * nn ** (2 * mm)
    
    gamma = np.linalg.solve(A, v)

    coeffs = {}
    for nstep, coeff in enumerate(gamma):
        if nstep == 0:
            coeffs[nstep] = -coeff
        else:
            coeffs[+nstep] = -coeff / 2
            coeffs[-nstep] = -coeff / 2

    return coeffs

```


```python
fig, ax = plt.subplots(figsize=(3, 2), dpi=250)

p = np.linspace(-1, 1, 1000)

for nstep in range(2, 10, 2):
    coeffs = derivative_coeffs(nstep)
    Dp = np.sum([-cn * np.cos(n * p) for n, cn in coeffs.items()], axis=0)
    ax.plot(p, p ** 2 - Dp, label=f"$D^{{({nstep})}}(p)$", lw=1)

ax.legend(fontsize=6, frameon=False)
ax.set_xlabel("$p$")
ax.set_ylabel("$p^2 - D(p)$")
ax.set_yscale("log")
ax.set_ylim(1.0e-12, 1)

plt.show(fig)
```


```python
fig, ax = plt.subplots(figsize=(3, 2), dpi=250)

p = np.linspace(-1, 1, 1000)

for nstep in range(10, 30, 5):
    coeffs = derivative_coeffs(nstep)
    nrange = np.array(list(coeffs.keys()))
    crange = np.array(list(coeffs.values()))
    inds = np.argsort(nrange)
    ax.plot(nrange[inds], np.abs(crange[inds]), label=f"$N_s = {nstep}$", ls="-", marker=".", ms=1, lw=1)

ax.legend(fontsize=6, frameon=False)
ax.set_xlabel("$p$")
ax.set_ylabel("$|c_n|$")
ax.set_yscale("log")

plt.show(fig)
```


```python

```

35aabe97af7cd4d2e36c333f12a67992d40ed9c7
ipynb
Jupyter Notebook
notebooks/derivations/Derivatives.ipynb
ckoerber/luescher-nd
d1bc6bff0c6ee9f4dc0d1d0bb4bcfa842c44cceb
BSD-3-Clause
2021-04-12T22:19:38.000Z
2021-07-26T14:06:49.000Z
notebooks/derivations/Derivatives.ipynb
ckoerber/luescher-nd
d1bc6bff0c6ee9f4dc0d1d0bb4bcfa842c44cceb
BSD-3-Clause
2019-12-16T19:49:19.000Z
2021-06-02T00:50:31.000Z
notebooks/derivations/Derivatives.ipynb
ckoerber/luescher-nd
d1bc6bff0c6ee9f4dc0d1d0bb4bcfa842c44cceb
BSD-3-Clause
Qwen/Qwen-72B
1. YES
2. YES