<jupyter_start><jupyter_text>Demo - Working with Functions  **<NAME>** (email: `<EMAIL>`), Department of Mathematics, University of Oslo.Date: **August 7, 2020****Summary.** This is a demonstration of how the Python module [shenfun](https://github.com/spectralDNS/shenfun) can be used to work withglobal spectral functions in one and several dimensions. ConstructionA global spectral function $u(x)$ is represented on a one-dimensionaldomain (a line) as $$u(x) = \sum_{k=0}^{N-1} \hat{u}_k \psi_k(x)$$ where $\psi_k(x)$ is the $k$'th basis function and $x$ is aposition inside the domain. $\{\hat{u}_k\}_{k=0}^{N-1}$ are theexpansion coefficient for the series, and often referred to as thedegrees of freedom. There is one degree of freedom per basis function.We can use any number of basis functions,and the span of the chosen basis is then a function space. Also part of thefunction space is the domain, which isspecified when a function space is created. To create a function space$T=\text{span}\{T_k\}_{k=0}^{N-1}$ forthe first N Chebyshev polynomials of the first kind on the default domain $[-1, 1]$,do<jupyter_code>from shenfun import *
N = 8
T = FunctionSpace(N, 'Chebyshev', domain=(-1, 1))<jupyter_output><empty_output><jupyter_text>The function $u(x)$ can now be created with all N coefficientsequal to zero as<jupyter_code>u = Function(T)<jupyter_output><empty_output><jupyter_text>When using Chebyshev polynomials the computational domain is always$[-1, 1]$. However, we can still use a different physical domain,like<jupyter_code>T = FunctionSpace(N, 'Chebyshev', domain=(0, 1))<jupyter_output><empty_output><jupyter_text>and under the hood shenfun will then map this domain to the referencedomain through $$u(x) = \sum_{k=0}^{N-1} \hat{u}_k \psi_k(2(x-0.5))$$  Approximating analytical functionsThe `u` function above was created with only zerovalued coefficients, which is the default. Alternatively,a [Function](https://shenfun.readthedocs.io/en/latest/shenfun.forms.htmlshenfun.forms.arguments.Function) may be initialized using a constantvalue<jupyter_code>T = FunctionSpace(N, 'Chebyshev', domain=(-1, 1))
u = Function(T, val=1)<jupyter_output><empty_output><jupyter_text>but that is not very useful. A third method to initializea [Function](https://shenfun.readthedocs.io/en/latest/shenfun.forms.htmlshenfun.forms.arguments.Function) is to interpolate using an analyticalSympy function.<jupyter_code>import sympy as sp
x = sp.Symbol('x', real=True)
u = Function(T, buffer=4*x**3-3*x)
print(u)<jupyter_output>[-6.93889390e-17 -5.47405982e-17 -1.82237827e-18  1.00000000e+00
  1.96261557e-17 -3.33066907e-16  1.49457239e-16 -2.93881977e-16]<jupyter_text>Here the analytical Sympy function will first be evaluatedon the entire quadrature mesh of the `T` function space,and then forward transformed to get the coefficients. Thiscorresponds to a projection to `T`. The projection isFind $u_h \in T$, such that $$(u_h - u, v)_w = 0 \quad \forall v \in T,$$ where $v \in \{T_j\}_{j=0}^{N-1}$ is a test function,$u_h=\sum_{k=0}^{N-1} \hat{u}_k T_k$ is a trial function and thenotation $(\cdot, \cdot)_w$ represents a weighted inner product.In this projection $u_h$ is the [Function](https://shenfun.readthedocs.io/en/latest/shenfun.forms.htmlshenfun.forms.arguments.Function), $u$ is the sympy function and we use sympyto exactly evaluate $u$ on all quadrature points$\{x_j\}_{j=0}^{N-1}$. With quadrature we then have $$(u, v)_w = \sum_{j\in\mathcal{I}^N} u(x_j) v(x_j) w_j \forall v \in T,$$ where $\mathcal{I}^N = (0, 1, \ldots, N-1)$ and $\{w_j\}_{j\in \mathcal{I}^N}$are the quadrature weights. The left hand side of the projection is $$(u_h, v)_w = \sum_{j\in\mathcal{I}^N} u_h(x_j) v(x_j) w_j \forall v \in T.$$ A linear system of equations arise when inserting for thebasis functions $$\left(u, T_i\right)_w = \tilde{u}_i \forall i \in \mathcal{I}^N,$$ and $$\begin{align*}\left(u_h, T_i \right)_w &= (\sum_{k\in \mathcal{I}^N} \hat{u}_k T_k , T_i)_w \\   &= \sum_{k\in \mathcal{I}^N} \left( T_k, T_i\right)_w \hat{u}_k\end{align*}$$ with the mass matrix $$a_{ik} = \left( T_k, T_i\right)_w   \forall (i, k) \in \mathcal{I}^N \times \mathcal{I}^N,$$ we can now solve to get the unknownexpansion coefficients. In matrix notation $$\hat{u} = A^{-1} \tilde{u},$$ where $\hat{u}=\{\hat{u}_i\}_{i\in \mathcal{I}^N}$,$\tilde{u}=\{\tilde{u}_i\}_{i \in \mathcal{I}^N}$ and$A=\{a_{ki}\}_{(i,k) \in \mathcal{I}^N \times \mathcal{I}^N}$ Adaptive function sizeThe number of basis functions can also be left open during creationof the function space, through<jupyter_code>T = FunctionSpace(0, 'Chebyshev', domain=(-1, 1))<jupyter_output><empty_output><jupyter_text>This is useful if you want to approximate a function andare uncertain how many basis functions that are required.For example, you may want to approximate the function $\cos(20 x)$.You can then find the required [Function](https://shenfun.readthedocs.io/en/latest/shenfun.forms.htmlshenfun.forms.arguments.Function) using<jupyter_code>u = Function(T, buffer=sp.cos(20*x))
print(len(u))<jupyter_output>45<jupyter_text>We see that $N=45$ is required to resolve this function. This agreeswell with what is reported also by [Chebfun](https://www.chebfun.org/docs/guide/guide01.html).Note that in this process a new [FunctionSpace()](https://shenfun.readthedocs.io/en/latest/shenfun.forms.htmlshenfun.forms.arguments.FunctionSpace) has beencreated under the hood. The function space of `u` can beextracted using<jupyter_code>Tu = u.function_space()
print(Tu.N)<jupyter_output>45<jupyter_text>To further show that shenfun is compatible with Chebfun we can alsoapproximate the Bessel function<jupyter_code>T1 = FunctionSpace(0, 'Chebyshev', domain=(0, 100))
u = Function(T1, buffer=sp.besselj(0, x))
print(len(u))<jupyter_output>83<jupyter_text>which gives 83 basis functions, in close agreement with Chebfun (89).The difference lies only in the cut-off criteria. We cut frequencieswith a relative tolerance of 1e-12 by default, but if we make this criteriaa little bit stronger, then we will also arrive at a slightly higher number:<jupyter_code>u = Function(T1, buffer=sp.besselj(0, x), reltol=1e-14)
print(len(u))<jupyter_output>87<jupyter_text>Plotting the function on its quadrature points looksa bit ragged, though:<jupyter_code>%matplotlib inline

import matplotlib.pyplot as plt
Tu = u.function_space()
plt.plot(Tu.mesh(), u.backward())<jupyter_output><empty_output><jupyter_text>To improve the quality of this plot we can instead evaluate thefunction on more points<jupyter_code>xj = np.linspace(0, 100, 1000)
plt.plot(xj, u(xj))<jupyter_output><empty_output><jupyter_text>Alternatively, we can refine the function, which simplypads zeros to $\hat{u}$<jupyter_code>up = u.refine(200)
Tp = up.function_space()
plt.plot(Tp.mesh(), up.backward())<jupyter_output><empty_output><jupyter_text>The padded expansion coefficients are now given as<jupyter_code>print(up)<jupyter_output>[ 6.43655435e-02 -1.37262942e-01  1.26593632e-01 -1.31630372e-01
  1.19810309e-01 -1.19499576e-01  1.07387826e-01 -9.95196201e-02
  8.81197795e-02 -7.07299943e-02  6.13974837e-02 -3.40222966e-02
  2.84578300e-02  6.05886310e-03 -6.23468323e-03  4.01956708e-02
 -3.44202241e-02  5.58353154e-02 -4.58848577e-02  4.33923041e-02
 -3.39778387e-02  6.26445963e-03 -3.37831264e-03 -3.27304442e-02
  2.63630973e-02 -4.12313235e-02  3.07429560e-02 -7.73458857e-03
  4.10441868e-03  3.29922472e-02 -2.47812478e-02  2.71654360e-02
 -1.81421444e-02 -2.02112290e-02  1.53660102e-02 -3.02943723e-02
  1.97324599e-02  1.94051084e-02 -1.41460820e-02  2.47321146e-02
 -1.49668436e-02 -3.02763148e-02  1.98156190e-02 -2.08517347e-03
 -6.91643640e-05  3.13799203e-02 -1.81924272e-02 -3.77110141e-02
  2.23051599e-02  2.82766986e-02 -1.67226896e-02 -1.59965877e-02
  9.40673776e-03  7.34392211e-03 -4.28398348e-03 -2.84497604e-03
  1.64428109e-03  9.53197714e-04 -5.45450869e-04 -2.80998822e-04
  1.59136941e-04  7.38240[...]<jupyter_text>More featuresSince we have used a regular Chebyshev basis above, thereare many more features that could be explored simply by going through[Numpy's Chebyshev module](https://numpy.org/doc/stable/reference/routines.polynomials.chebyshev.html).For example, we can create a Chebyshev series like<jupyter_code>import numpy.polynomial.chebyshev as cheb
c = cheb.Chebyshev(u, domain=(0, 100))<jupyter_output><empty_output><jupyter_text>The Chebyshev series in Numpy has a wide range of possibilities,see [here](https://numpy.org/doc/stable/reference/generated/numpy.polynomial.chebyshev.Chebyshev.htmlnumpy.polynomial.chebyshev.Chebyshev).However, we may also work directly with the Chebyshevcoefficients already in `u`. To find the roots of thepolynomial that approximates the Bessel function ondomain $[0, 100]$, we can do<jupyter_code>z = Tu.map_true_domain(cheb.chebroots(u))<jupyter_output><empty_output><jupyter_text>Note that the roots are found on the reference domain $[-1, 1]$and as such we need to move the result to the physical domain using`map_true_domain`. The resulting roots `z` are both real and imaginary,so to extract the real roots we need to filter a little bit<jupyter_code>z2 = z[np.where((z.imag == 0)*(z.real > 0)*(z.real < 100))].real
print(z2[:5])<jupyter_output>[ 2.40482556  5.52007811  8.65372791 11.79153444 14.93091771]<jupyter_text>Here `np.where` returns the indices where the condition is true. The conditionis that the imaginary part is zero, whereas the real part is within thetrue domain $[0, 100]$.**Notice.**Using directly `cheb.chebroots(c)` does not seem to work (even though theseries has been generated with the non-standard domain) becauseNumpy only looks for roots in the reference domain $[-1, 1]$.We could also use a function space with boundary conditions builtin, like<jupyter_code>Td = FunctionSpace(0, 'C', bc=(sp.besselj(0, 0), sp.besselj(0, 100)), domain=(0, 100))
ud = Function(Td, buffer=sp.besselj(0, x))
print(len(ud))<jupyter_output>82<jupyter_text>As we can see this leads to a function space of dimensionvery similar to the orthogonal space.The major advantages of working with a space with boundary conditionsbuilt in only comes to life when solving differential equations. Aslong as we are only interested in approximating functions, we may justas well stick to the orthogonal spaces. This goes for Legendre aswell as Chebyshev. Multidimensional functionsMultidimensional tensor product spaces are createdby taking the tensor products of one-dimensional function spaces.For example<jupyter_code>C0 = FunctionSpace(20, 'C')
C1 = FunctionSpace(20, 'C')
T = TensorProductSpace(comm, (C0, C1))
u = Function(T)<jupyter_output><empty_output><jupyter_text>Here $\text{T} = \text{C0} \otimes \text{C1}$, the basis function is$T_i(x) T_j(y)$ and the Function `u` is $$u(x, y) = \sum_{i=0}^{N-1} \sum_{j=0}^{N-1} \hat{u}_{ij} T_i(x) T_j(y).$$ The multidimensional Functions work more or less exactly like for the1D case. We can here interpolate 2D Sympy functions<jupyter_code>y = sp.Symbol('y', real=True)
u = Function(T, buffer=sp.cos(10*x)*sp.cos(10*y))
X = T.local_mesh(True)
plt.contourf(X[0], X[1], u.backward())<jupyter_output><empty_output><jupyter_text>Like for 1D the coefficients are computed through projection,where the exact function is evaluated on all quadrature pointsin the mesh.The Cartesian mesh represents the quadrature points of thetwo function spaces, and can be visualized as follows<jupyter_code>X = T.mesh()
for xj in X[0]:
    for yj in X[1]:
        plt.plot((xj, xj), (X[1][0, 0], X[1][0, -1]), 'k')
        plt.plot((X[0][0], X[0][-1]), (yj, yj), 'k')<jupyter_output><empty_output><jupyter_text>We may alternatively plot on a uniform mesh<jupyter_code>X = T.local_mesh(broadcast=True, uniform=True)
plt.contourf(X[0], X[1], u.backward(kind='uniform'))<jupyter_output><empty_output><jupyter_text>Curvilinear coordinatesWith shenfun it is possible to use curvilinear coordinates,and not necessarily with orthogonal basis vectors. Withcurvilinear coordinates the computational coordinates arealways straight lines, rectangles and cubes. But the physicalcoordinates can be very complex.Consider the unit disc with polar coordinates. Herethe position vector $\mathbf{r}$ is given by $$\mathbf{r} = r\cos \theta \mathbf{i} + r\sin \theta \mathbf{j}$$ The physical domain is $\Omega = \{(x, y): x^2 + y^2 < 1\}$,whereas the computational domain is the Cartesian product$D = \{(r, \theta) \in [0, 1] \times [0, 2 \pi]\}$.We create this domain in shenfun through<jupyter_code>r, theta = psi = sp.symbols('x,y', real=True, positive=True)
rv = (r*sp.cos(theta), r*sp.sin(theta))
B0 = FunctionSpace(20, 'C', domain=(0, 1))
F0 = FunctionSpace(20, 'F')
T = TensorProductSpace(comm, (B0, F0), coordinates=(psi, rv))<jupyter_output><empty_output><jupyter_text>Note that we are using a Fourier space for the azimuthaldirection, since the solution here needs to be periodic.We can now create functions on the space using ananalytical function in computational coordinates<jupyter_code>u = Function(T, buffer=(1-r)*r*sp.sin(sp.cos(theta)))<jupyter_output><empty_output><jupyter_text>However, when this is plotted it may not be what you expect<jupyter_code>X = T.local_mesh(True)
plt.contourf(X[0], X[1], u.backward(), 100)<jupyter_output><empty_output><jupyter_text>We see that the function has been plotted in computational coordinates,and not on the disc, as you probably expected. To plot onthe disc we need the physical mesh, and not the computational<jupyter_code>X = T.local_cartesian_mesh()
plt.contourf(X[0], X[1], u.backward(), 100)<jupyter_output><empty_output><jupyter_text>**Notice.**The periodic plot does not wrap all around the circle. This isnot wrong, we have simply not used the same point twice, but itdoes not look very good. To overcome this problem we can wrap thegrid all the way around and re-plot.<jupyter_code>up = u.backward()
xp, yp, up = wrap_periodic([X[0], X[1], up], axes=[1])
plt.contourf(xp, yp, up, 100)<jupyter_output><empty_output><jupyter_text>Adaptive functions in multiple dimensionsIf you want to find a good resolution for a function in multipledimensions, the procedure is exactly like in 1D. First create functionspaces with 0 quadrature points, and then call [Function](https://shenfun.readthedocs.io/en/latest/shenfun.forms.htmlshenfun.forms.arguments.Function)<jupyter_code>B0 = FunctionSpace(0, 'C', domain=(0, 1))
F0 = FunctionSpace(0, 'F')
T = TensorProductSpace(comm, (B0, F0), coordinates=(psi, rv))
u = Function(T, buffer=((1-r)*r)**2*sp.sin(sp.cos(theta)))
print(u.shape)<jupyter_output>(5, 12)
{'fraction_non_alphanumeric': 0.1304109213143994, 'fraction_numerical': 0.07052205529258421, 'mean_word_length': 5.608340888485947, 'pattern_counts': {'":': 0, '<': 108, '<?xml version=': 0, '>': 107, 'https://': 10, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '28', 'dir': 'jupyter-structured-clean-dedup', 'id': '24933', 'license': 'BSD-2-Clause', 'n_tokens_mistral': 5795, 'n_tokens_neox': 5114, 'n_words': 1496, 'path': 'docs/source/functions.ipynb', 'repo_name': 'liqihao2000/shenfun', 'size': '14598'}