<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 <!-- File automatically generated using DocOnce (https
Step1:  The Python Function uN represents the expansion (2), and the printed
values represent $\boldsymbol{\hat{u}} = {\hat{u}k}{k=0}^{N-1}$. The expansion is fairly well resolved since
the highest values of ${\hat{u}k}{k=0}^{N-1}$ approach 0.
Note that the coefficients obtained are based on interpolation at
quadrature points and they do not agree completely with the coefficients truncated from an
infinite series $u(x) = \sum_{k=0}^{\infty} \hat{u}_k T_k$. The obtained series is
often denoted as $u_N(x) = I_N u(x)$, where $I_N$ is an interpolation operator.
Under the hood the coefficients are found by projection using quadrature for the integrals
Step2:  Differentiation
Let us now consider the $n$'th derivative of $u(x)$ instead, denoted here as $u^{(n)}$.
We want to find $u^{(n)}$ in the space $S_N$, which means that we want to obtain
the following expansion
$$
u_N^{(n)} = \sum_{k=0}^{N-1} \hat{u}^{(n)}_k T_k.
$$
First note that this expansion is not the same as the derivative of
the previously found $u_N$, which is
$$
(u_N)^{(n)} = \sum_{k=0}^{N-1} \hat{u}_k T^{(n)}_k,
$$
where $T^{(n)}_k$ is the $n$'th derivative of $T_k$, a polynomial of order $k-n$.
We again use projection to find $u_N^{(n)} \in S_N$ such that
$$
(u_N^{(n)}-u^{(n)}, v)_{\omega^{-1/2}} = 0, \quad \forall v \in S_N.
$$
Inserting for $u_N^{(n)}$ and $u^{(n)} = (u_N)^{(n)}$ we get
<!-- Equation labels as ordinary links -->
<div id="_auto1"></div>
$$
\begin{equation}
    \sum_{j=0}^{N-1}(T_j, T_k){\omega^{-1/2}} \hat{u}_j^{(n)} = (T_j^{(n)}, T_k){\omega^{-1/2}} \hat{u}_j, 
\label{_auto1} \tag{3}
\end{equation}
$$
<!-- Equation labels as ordinary links -->
<div id="_auto2"></div>
$$
\begin{equation}
    \sum_{j=0}^{N-1} d^{(0)}{kj} \hat{u}_j^{(n)}  = \sum{j=0}^{N-1} d^{(n)}_{kj} \hat{u}_j,
\label{_auto2} \tag{4}
\end{equation}
$$
where $d^{(n)}{kj} = (T_j^{(n)}, T_k){\omega^{-1/2}}$.
We compute $\hat{u}k^{(n)}$ by inverting the diagonal $d^{(0)}{kj}$
<!-- Equation labels as ordinary links -->
<div id="eq
Step3:  where uN1 $=u_N^{(1)}$ and uN2 $=u_N^{(2)}$.
Alternatively, doing all the work that goes on under the hood
Step4:  We can look at the sparsity patterns of $(d^{(1)}{kj})$ and $(d^{(2)}{kj})$
Step5:  just to see that they are upper triangular. We now ask is there a better and faster
way to get uN1 and uN2? A better approach would involve only sparse
matrices, like the diagonal $(d^{(0)}_{kj})$. But how do we get there?
Most textbooks on spectral methods use fast recursive methods to
find the coefficients ${\hat{u}_k^{(n)}}$. Here we will show a fast Galerkin approach.
It turns out that a simple change of test space/function will be sufficient.
Let us first replace the test space $S_N$ with
the Dirichlet space $D_N={v \in S_N | v(\pm 1) = 0}$ using basis
functions $v=T_k-T_{k+2}$ and see
what happens. Because of the two boundary conditions,
the number of degrees of freedom is reduced by two, and we need to use a
space with $N+2$ quadrature points in order to get a square matrix system.
The method now becomes classified as Chebyshev-Petrov-Galerkin, as we
wish to find $u_N^{(1)} \in S_N$ such that
$$
(u_N^{(1)}-u^{(1)}, v){\omega^{-1/2}} = 0, \quad \forall v \in D{N+2}.
$$
The implementation is straightforward
Step6:  and since uN11 = uN1 we see that we have achived the same result as in
the regular projection. However, the matrices in use now look like
Step7:  So $(d^{(0)}{kj})$ now contains two nonzero diagonals, whereas $(d^{(1)}{kj})$ is
a matrix with one single diagonal. There is no longer a full differentiation
matrix, and we can easily perform this projection for millions of degrees of freedom.
What about $(d^{(2)}_{kj})$? We can now use biharmonic test functions that
satisfy four boundary conditions in the space $B_N = {v \in S_N | v(\pm 1) = v'(\pm 1) =0}$,
and continue in a similar fashion
Step8:  We get that uN22 = uN2, so the Chebyshev-Petrov-Galerkin projection works. The matrices involved are now
Step9:  So there are now three nonzero diagonals in $(d^{(0)}{kj})$, whereas the differentiation matrix
$(d^{(2)}{kj})$ contains only one nonzero diagonal.
Why does this work so well? The Chebyshev polynomials and their derivatives satisfy the following orthogonality relation
<!-- Equation labels as ordinary links -->
<div id="eq
Step10:  which evidently is exactly the same result. It turns out that this is not strange. If we multiply
(11) by $\pi/2$, rearrange a little bit and move to a matrix form we get
<!-- Equation labels as ordinary links -->
<div id="_auto5"></div>
$$
\begin{equation}
c_k \pi/2 \hat{u}^{(1)}k -  \pi/2 \hat{u}^{(1)}{k+2} = (k+1)\pi \hat{u}_{k+1} 
\label{_auto5} \tag{12}
\end{equation}
$$
<!-- Equation labels as ordinary links -->
<div id="_auto6"></div>
$$
\begin{equation}
\underbrace{(c_k \pi/2 \delta_{kj} - \pi/2 \delta_{k, j-2})}{(D^0){kj}} \hat{u}^{(1)}j = \underbrace{(k+1)\pi \delta{k,j-1}}{(D^1){kj}} \hat{u}_{j} 
\label{_auto6} \tag{13}
\end{equation}
$$
<!-- Equation labels as ordinary links -->
<div id="_auto7"></div>
$$
\begin{equation}
D^0 \boldsymbol{\hat{u}} = D^1 \boldsymbol{\hat{u}} 
\label{_auto7} \tag{14}
\end{equation}
$$
<!-- Equation labels as ordinary links -->
<div id="_auto8"></div>
$$
\begin{equation}
\boldsymbol{\hat{u}^{(1)}} = (D^0)^{-1} D^1 \boldsymbol{\hat{u}}
\label{_auto8} \tag{15}
\end{equation}
$$
which is exactly how $\boldsymbol{\hat{u}^{(1)}}$ was computed above with the Chebyshev-Petrov-Galerkin approach
(compare with the code line uN11 = D0.solve(D1.matvec(uN, w0), uN11)). Not convinced? Check that the matrices
D0 and D1 are truly as stated above. The matrices below are printed as dictionaries with diagonal
number as key (main is 0, first upper is 1 etc) and diagonal values as values
<ASSISTANT_TASK:>
Python Code:
from shenfun import *
import sympy as sp
x = sp.Symbol('x')
ue = sp.sin(sp.pi*x)
N = 16
SN = FunctionSpace(N, 'C')
uN = Function(SN, buffer=ue)
uN
Explanation: <!-- File automatically generated using DocOnce (https://github.com/doconce/doconce/):
doconce format ipynb sparsity.do.txt  -->
Demo - Sparse Chebyshev-Petrov-Galerkin methods for differentiation
Mikael Mortensen (email: mikaem@math.uio.no), Department of Mathematics, University of Oslo.
Date: October 26, 2021
Summary. This demo explores how to use sparse Chebyshev-Petrov-Galerkin methods for finding Chebyshev coefficients of
the derivatives of smooth functions. We will compare the methods to the more commonly adopted
recursion methods that are found in most spectral textbooks.
Introduction
The Chebyshev polynomials of the first kind can be defined as
<!-- Equation labels as ordinary links -->
<div id="eq:chebTU"></div>
$$
\begin{equation}
\label{eq:chebTU} \tag{1}
    T_k(x) = \cos(k\theta),
\end{equation}
$$
where $\theta = \cos^{-1} x$, $k$ is a positive integer and $x \in [-1, 1]$.
The Chebyshev polynomials span the discrete space $S_N = \text{span}{T_k}_{k=0}^{N-1}$,
and a function $u(x)$ can be approximated in this space as
<!-- Equation labels as ordinary links -->
<div id="eq:uT"></div>
$$
\begin{equation}
u_N(x) = \sum_{k=0}^{N-1} \hat{u}_k T_k(x).
\label{eq:uT} \tag{2}
\end{equation}
$$
Consider the expansion of the function $u(x)=\sin(\pi x)$, created in shenfun as
End of explanation
SM = FunctionSpace(0, 'C')
uM = Function(SM, buffer=ue, abstol=1e-16, reltol=1e-16)
print(uM[:N] - uN[:N])
print(len(uM))
Explanation: The Python Function uN represents the expansion (2), and the printed
values represent $\boldsymbol{\hat{u}} = {\hat{u}k}{k=0}^{N-1}$. The expansion is fairly well resolved since
the highest values of ${\hat{u}k}{k=0}^{N-1}$ approach 0.
Note that the coefficients obtained are based on interpolation at
quadrature points and they do not agree completely with the coefficients truncated from an
infinite series $u(x) = \sum_{k=0}^{\infty} \hat{u}_k T_k$. The obtained series is
often denoted as $u_N(x) = I_N u(x)$, where $I_N$ is an interpolation operator.
Under the hood the coefficients are found by projection using quadrature for the integrals:
find $u_N \in S_N$ such that
$$
(u_N-u, v)_{\omega^{-1/2}} = 0, \quad \forall v \in S_N,
$$
where $\omega = (1-x^2)$ and the scalar product notation
$(a, b){\omega^{-1/2}} = \sum{j=0}^{N-1} a(x_j)b(x_j)\omega_j \approx \int_{-1}^{1} a(x)b(x) \omega(x)^{-1/2} dx$,
where ${\omega_j}{j=0}^{N-1}$ are the quadrature weights. The quadrature approach ensures
that $u(x_j) = u_N(x_j)$ for all quadrature points ${x_j}{j=0}^{N-1}$.
In shenfun we compute the following under the hood: insert for $u_N = \sum_{j=0}^{N-1} \hat{u}_j T_j$,
$u=\sin(\pi x)$ and $v = T_k$ to get
$$
\sum_{j=0}^{N-1}(T_j, T_k){\omega^{-1/2}} \hat{u}_j = (\sin(\pi x), T_k){\omega^{-1/2}},
$$
This has now become a linear algebra problem, and we recognise the matrix $d^{(0)}{kj} = (T_j, T_k){\omega^{-1/2}}=c_k \pi /2 \delta_{kj}$,
where $\delta_{kj}$ is the Kronecker delta function, and $c_0=2$ and $c_k=1$ for $k>0$.
The problem is solved trivially since $d^{(0)}_{kj}$ is diagonal,
and thus
$$
\hat{u}k = \frac{2}{c_k \pi} (\sin(\pi x), T_k){\omega^{-1/2}}, \quad \forall \, k\in I^N,
$$
where $I^N = {0, 1, \ldots, N-1}$.
We can compare this to the truncated coefficients, where the integral $(\sin(\pi x), T_k){\omega^{-1/2}}$
is computed with high precision. To this end we could use adaptive quadrature, or symbolic integration
with sympy, but it is sufficient to use a large enough number of polynomials to fully resolve the
function. Below we find this number to be 22 and we see that the absolute error in
the highest wavenumber $\hat{u}{15} \approx 10^{-11}$.
End of explanation
uN1 = project(Dx(uN, 0, 1), SN)
uN2 = project(Dx(uN, 0, 2), SN)
uN1
Explanation: Differentiation
Let us now consider the $n$'th derivative of $u(x)$ instead, denoted here as $u^{(n)}$.
We want to find $u^{(n)}$ in the space $S_N$, which means that we want to obtain
the following expansion
$$
u_N^{(n)} = \sum_{k=0}^{N-1} \hat{u}^{(n)}_k T_k.
$$
First note that this expansion is not the same as the derivative of
the previously found $u_N$, which is
$$
(u_N)^{(n)} = \sum_{k=0}^{N-1} \hat{u}_k T^{(n)}_k,
$$
where $T^{(n)}_k$ is the $n$'th derivative of $T_k$, a polynomial of order $k-n$.
We again use projection to find $u_N^{(n)} \in S_N$ such that
$$
(u_N^{(n)}-u^{(n)}, v)_{\omega^{-1/2}} = 0, \quad \forall v \in S_N.
$$
Inserting for $u_N^{(n)}$ and $u^{(n)} = (u_N)^{(n)}$ we get
<!-- Equation labels as ordinary links -->
<div id="_auto1"></div>
$$
\begin{equation}
    \sum_{j=0}^{N-1}(T_j, T_k){\omega^{-1/2}} \hat{u}_j^{(n)} = (T_j^{(n)}, T_k){\omega^{-1/2}} \hat{u}_j, 
\label{_auto1} \tag{3}
\end{equation}
$$
<!-- Equation labels as ordinary links -->
<div id="_auto2"></div>
$$
\begin{equation}
    \sum_{j=0}^{N-1} d^{(0)}{kj} \hat{u}_j^{(n)}  = \sum{j=0}^{N-1} d^{(n)}_{kj} \hat{u}_j,
\label{_auto2} \tag{4}
\end{equation}
$$
where $d^{(n)}{kj} = (T_j^{(n)}, T_k){\omega^{-1/2}}$.
We compute $\hat{u}k^{(n)}$ by inverting the diagonal $d^{(0)}{kj}$
<!-- Equation labels as ordinary links -->
<div id="eq:fhat"></div>
$$
\begin{equation}
\hat{u}k^{(n)} = \frac{2}{c_k \pi} \sum{j=0}^{N-1} d^{(n)}_{kj} \hat{u}_j, \quad \forall \, k \in I^{N}.
\label{eq:fhat} \tag{5}
\end{equation}
$$
The matrix $d^{(n)}{kj}$ is upper triangular, and the last $n$ rows are zero. Since $d^{(n)}{kj}$ is
dense the matrix vector product $\sum_{j=0}^{N-1} d^{(n)}_{kj} \hat{u}_j$ is costly
and also susceptible to roundoff errors if the structure of the matrix is
not taken advantage of. But computing it in shenfun
is straightforward, for $n=1$ and $2$:
End of explanation
u = TrialFunction(SN)
v = TestFunction(SN)
D0 = inner(u, v)
D1 = inner(Dx(u, 0, 1), v)
D2 = inner(Dx(u, 0, 2), v)
w0 = Function(SN) # work array
uN1 = Function(SN)
uN2 = Function(SN)
uN1 = D0.solve(D1.matvec(uN, w0), uN1)
uN2 = D0.solve(D2.matvec(uN, w0), uN2)
uN1
Explanation: where uN1 $=u_N^{(1)}$ and uN2 $=u_N^{(2)}$.
Alternatively, doing all the work that goes on under the hood
End of explanation
%matplotlib inline
import matplotlib.pyplot as plt
fig, (ax1, ax2) = plt.subplots(1, 2)
ax1.spy(D1.diags(), markersize=2, color='r')
ax2.spy(D2.diags(), markersize=2, color='b')
Explanation: We can look at the sparsity patterns of $(d^{(1)}{kj})$ and $(d^{(2)}{kj})$
End of explanation
DN = FunctionSpace(N+2, 'C', bc=(0, 0))
v = TestFunction(DN)
D0 = inner(u, v)
D1 = inner(Dx(u, 0, 1), v)
uN11 = Function(SN)
uN11 = D0.solve(D1.matvec(uN, w0), uN11)
print(uN11-uN1)
Explanation: just to see that they are upper triangular. We now ask is there a better and faster
way to get uN1 and uN2? A better approach would involve only sparse
matrices, like the diagonal $(d^{(0)}_{kj})$. But how do we get there?
Most textbooks on spectral methods use fast recursive methods to
find the coefficients ${\hat{u}_k^{(n)}}$. Here we will show a fast Galerkin approach.
It turns out that a simple change of test space/function will be sufficient.
Let us first replace the test space $S_N$ with
the Dirichlet space $D_N={v \in S_N | v(\pm 1) = 0}$ using basis
functions $v=T_k-T_{k+2}$ and see
what happens. Because of the two boundary conditions,
the number of degrees of freedom is reduced by two, and we need to use a
space with $N+2$ quadrature points in order to get a square matrix system.
The method now becomes classified as Chebyshev-Petrov-Galerkin, as we
wish to find $u_N^{(1)} \in S_N$ such that
$$
(u_N^{(1)}-u^{(1)}, v){\omega^{-1/2}} = 0, \quad \forall v \in D{N+2}.
$$
The implementation is straightforward
End of explanation
fig, (ax1, ax2) = plt.subplots(1, 2)
ax1.spy(D0.diags(), markersize=2, color='r')
ax2.spy(D1.diags(), markersize=2, color='b')
Explanation: and since uN11 = uN1 we see that we have achived the same result as in
the regular projection. However, the matrices in use now look like
End of explanation
BN = FunctionSpace(N+4, 'C', bc=(0, 0, 0, 0))
v = TestFunction(BN)
D0 = inner(u, v)
D2 = inner(Dx(u, 0, 2), v)
uN22 = Function(SN)
uN22 = D0.solve(D2.matvec(uN, w0), uN22)
print(uN22-uN2)
Explanation: So $(d^{(0)}{kj})$ now contains two nonzero diagonals, whereas $(d^{(1)}{kj})$ is
a matrix with one single diagonal. There is no longer a full differentiation
matrix, and we can easily perform this projection for millions of degrees of freedom.
What about $(d^{(2)}_{kj})$? We can now use biharmonic test functions that
satisfy four boundary conditions in the space $B_N = {v \in S_N | v(\pm 1) = v'(\pm 1) =0}$,
and continue in a similar fashion:
End of explanation
fig, (ax1, ax2) = plt.subplots(1, 2)
ax1.spy(D0.diags(), markersize=2, color='r')
ax2.spy(D2.diags(), markersize=2, color='b')
Explanation: We get that uN22 = uN2, so the Chebyshev-Petrov-Galerkin projection works. The matrices involved are now
End of explanation
f1 = np.zeros(N+1)
ck = np.ones(N); ck[0] = 2
for k in range(N-2, -1, -1):
    f1[k] = (f1[k+2]+2*(k+1)*uN[k+1])/ck[k]
print(f1[:-1]-uN1)
Explanation: So there are now three nonzero diagonals in $(d^{(0)}{kj})$, whereas the differentiation matrix
$(d^{(2)}{kj})$ contains only one nonzero diagonal.
Why does this work so well? The Chebyshev polynomials and their derivatives satisfy the following orthogonality relation
<!-- Equation labels as ordinary links -->
<div id="eq:orthon"></div>
$$
\begin{equation}
\label{eq:orthon} \tag{6}
  \int_{-1}^{1} T^{(n)}j T^{(n)}_k \omega^{n-1/2} dx = \alpha^{n}_k \delta{kj}, \quad \text{for}\, n \ge 0,
\end{equation}
$$
where
<!-- Equation labels as ordinary links -->
<div id="_auto3"></div>
$$
\begin{equation}
    \alpha^n_k = \frac{c_{k+n}\pi k (k+n-1)!}{2(k-n)!}.
\label{_auto3} \tag{7}
\end{equation}
$$
So when we choose a test function $\omega^n T^{(n)}_k$ and a trial function $T_j$,
we get the diagonal differentiation matrix
<!-- Equation labels as ordinary links -->
<div id="_auto4"></div>
$$
\begin{equation}
  d^{(n)}{kj} = \int{-1}^{1} T^{(n)}j (\omega^n T^{(n)}_k) \omega^{-1/2} dx = \alpha^{n}_k \delta{kj}, \quad \text{for}\, n \ge 0.
\label{_auto4} \tag{8}
\end{equation}
$$
The two chosen test functions above are both proportional to $\omega^n T^{(n)}k$. More precisely,
$T_k-T{k+2} = \frac{2}{k+1} \omega T^{(1)}{k+1}$ and the biharmonic test function is
$T_k-\frac{2(k+2)}{k+3}T{k+2} + \frac{k+1}{k+3}T_{k+4} = \frac{4 \omega^2T^{(2)}_{k+2}}{(k+2)(k+3)}$.
Using these very specific test functions correponds closely to using the Chebyshev
recursion formulas that are found in most textbooks. Here they are adapted to
a Chebyshev-Petrov-Galerkin method, where we simply choose test and trial functions and everything
else falls into place in a few lines of code.
Recursion
Let us for completion show how to
find $\hat{u}_N^{(1)}$ with a recursive approach. The Chebyshev polynomials
satisfy
<!-- Equation labels as ordinary links -->
<div id="eq:Trec1"></div>
$$
\begin{equation}
2T_k = \frac{1}{k+1}T'{k+1}- \frac{1}{k-1} T'{k-1}, \quad k \ge 1.
\label{eq:Trec1} \tag{9}
\end{equation}
$$
By using this and setting $u' = \sum_{k=0}^{\infty} \hat{u}^{(1)}k T_k = \sum{k=0}^{\infty} \hat{u}_k T'_k$
we get
<!-- Equation labels as ordinary links -->
<div id="eq:Trec2"></div>
$$
\begin{equation}
2k\hat{u}k = c{k-1}\hat{u}^{(1)}{k-1} - \hat{u}^{(1)}{k+1}, \quad k \ge 1.
\label{eq:Trec2} \tag{10}
\end{equation}
$$
Using this recursion on a discrete series, together with $\hat{u}^{(1)}{N} = \hat{u}^{(1)}{N-1} = 0$, we get
(see [canuto] Eq. (2.4.24))
<!-- Equation labels as ordinary links -->
<div id="eq:Trec3"></div>
$$
\begin{equation}
c_k \hat{u}^{(1)}k = \hat{u}^{(1)}{k+2} + 2(k+1)\hat{u}_{k+1}, \quad 0 \le k \le N-2,
\label{eq:Trec3} \tag{11}
\end{equation}
$$
which is easily implemented in a (slow) for-loop
End of explanation
import pprint
DN = FunctionSpace(N+2, 'C', bc=(0, 0))
v = TestFunction(DN)
D0 = inner(u, v)
D1 = inner(Dx(u, 0, 1), v)
pprint.pprint(dict(D0))
pprint.pprint(dict(D1))
Explanation: which evidently is exactly the same result. It turns out that this is not strange. If we multiply
(11) by $\pi/2$, rearrange a little bit and move to a matrix form we get
<!-- Equation labels as ordinary links -->
<div id="_auto5"></div>
$$
\begin{equation}
c_k \pi/2 \hat{u}^{(1)}k -  \pi/2 \hat{u}^{(1)}{k+2} = (k+1)\pi \hat{u}_{k+1} 
\label{_auto5} \tag{12}
\end{equation}
$$
<!-- Equation labels as ordinary links -->
<div id="_auto6"></div>
$$
\begin{equation}
\underbrace{(c_k \pi/2 \delta_{kj} - \pi/2 \delta_{k, j-2})}{(D^0){kj}} \hat{u}^{(1)}j = \underbrace{(k+1)\pi \delta{k,j-1}}{(D^1){kj}} \hat{u}_{j} 
\label{_auto6} \tag{13}
\end{equation}
$$
<!-- Equation labels as ordinary links -->
<div id="_auto7"></div>
$$
\begin{equation}
D^0 \boldsymbol{\hat{u}} = D^1 \boldsymbol{\hat{u}} 
\label{_auto7} \tag{14}
\end{equation}
$$
<!-- Equation labels as ordinary links -->
<div id="_auto8"></div>
$$
\begin{equation}
\boldsymbol{\hat{u}^{(1)}} = (D^0)^{-1} D^1 \boldsymbol{\hat{u}}
\label{_auto8} \tag{15}
\end{equation}
$$
which is exactly how $\boldsymbol{\hat{u}^{(1)}}$ was computed above with the Chebyshev-Petrov-Galerkin approach
(compare with the code line uN11 = D0.solve(D1.matvec(uN, w0), uN11)). Not convinced? Check that the matrices
D0 and D1 are truly as stated above. The matrices below are printed as dictionaries with diagonal
number as key (main is 0, first upper is 1 etc) and diagonal values as values:
End of explanation
<END_TASK>