<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 function $u(x)$ can now be created with all N coefficients
equal to zero as
Step2:  When using Chebyshev polynomials the computational domain is always
$[-1, 1]$. However, we can still use a different physical domain,
like
Step3:  and under the hood shenfun will then map this domain to the reference
domain through
$$
u(x) = \sum_{k=0}^{N-1} \hat{u}_k \psi_k(2(x-0.5))
$$
Approximating analytical functions
The u function above was created with only zero
valued coefficients, which is the default. Alternatively,
a Function may be initialized using a constant
value
Step4:  but that is not very useful. A third method to initialize
a Function is to interpolate using an analytical
Sympy function.
Step5:  Here the analytical Sympy function will first be evaluated
on the entire quadrature mesh of the T function space,
and then forward transformed to get the coefficients. This
corresponds to a projection to T. The projection is
Find $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 the
notation $(\cdot, \cdot)_w$ represents a weighted inner product.
In this projection $u_h$ is the Function, $u$ is the sympy function and we use sympy
to 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 the
basis 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 unknown
expansion 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 size
The number of basis functions can also be left open during creation
of the function space, through
Step6:  This is useful if you want to approximate a function and
are 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 using
Step7:  We see that $N=45$ is required to resolve this function. This agrees
well with what is reported also by Chebfun.
Note that in this process a new FunctionSpace() has been
created under the hood. The function space of u can be
extracted using
Step8:  To further show that shenfun is compatible with Chebfun we can also
approximate the Bessel function
Step9:  which gives 83 basis functions, in close agreement with Chebfun (89).
The difference lies only in the cut-off criteria. We cut frequencies
with a relative tolerance of 1e-12 by default, but if we make this criteria
a little bit stronger, then we will also arrive at a slightly higher number
Step10:  Plotting the function on its quadrature points looks
a bit ragged, though
Step11:  To improve the quality of this plot we can instead evaluate the
function on more points
Step12:  Alternatively, we can refine the function, which simply
pads zeros to $\hat{u}$
Step13:  The padded expansion coefficients are now given as
Step14:  More features
Since we have used a regular Chebyshev basis above, there
are many more features that could be explored simply by going through
Numpy's Chebyshev module.
For example, we can create a Chebyshev series like
Step15:  The Chebyshev series in Numpy has a wide range of possibilities,
see here.
However, we may also work directly with the Chebyshev
coefficients already in u. To find the roots of the
polynomial that approximates the Bessel function on
domain $[0, 100]$, we can do
Step16:  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
Step17:  Here np.where returns the indices where the condition is true. The condition
is that the imaginary part is zero, whereas the real part is within the
true domain $[0, 100]$.
Notice.
Using directly cheb.chebroots(c) does not seem to work (even though the
series has been generated with the non-standard domain) because
Numpy only looks for roots in the reference domain $[-1, 1]$.
We could also use a function space with boundary conditions built
in, like
Step18:  As we can see this leads to a function space of dimension
very similar to the orthogonal space.
The major advantages of working with a space with boundary conditions
built in only comes to life when solving differential equations. As
long as we are only interested in approximating functions, we may just
as well stick to the orthogonal spaces. This goes for Legendre as
well as Chebyshev.
Multidimensional functions
Multidimensional tensor product spaces are created
by taking the tensor products of one-dimensional function spaces.
For example
Step19:  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 the
1D case. We can here interpolate 2D Sympy functions
Step20:  Like for 1D the coefficients are computed through projection,
where the exact function is evaluated on all quadrature points
in the mesh.
The Cartesian mesh represents the quadrature points of the
two function spaces, and can be visualized as follows
Step21:  We may alternatively plot on a uniform mesh
Step22:  Curvilinear coordinates
With shenfun it is possible to use curvilinear coordinates,
and not necessarily with orthogonal basis vectors. With
curvilinear coordinates the computational coordinates are
always straight lines, rectangles and cubes. But the physical
coordinates can be very complex.
Consider the unit disc with polar coordinates. Here
the 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)
Step23:  Note that we are using a Fourier space for the azimuthal
direction, since the solution here needs to be periodic.
We can now create functions on the space using an
analytical function in computational coordinates
Step24:  However, when this is plotted it may not be what you expect
Step25:  We see that the function has been plotted in computational coordinates,
and not on the disc, as you probably expected. To plot on
the disc we need the physical mesh, and not the computational
Step26:  Notice.
The periodic plot does not wrap all around the circle. This is
not wrong, we have simply not used the same point twice, but it
does not look very good. To overcome this problem we can wrap the
grid all the way around and re-plot.
Step27:  Adaptive functions in multiple dimensions
If you want to find a good resolution for a function in multiple
dimensions, the procedure is exactly like in 1D. First create function
spaces with 0 quadrature points, and then call Function
<ASSISTANT_TASK:>
Python Code:
from shenfun import *
N = 8
T = FunctionSpace(N, 'Chebyshev', domain=(-1, 1))
Explanation: <!-- File automatically generated using DocOnce (https://github.com/doconce/doconce/):
doconce format ipynb functions.do.txt  -->
Demo - Working with Functions
Mikael Mortensen (email: mikaem@math.uio.no), Department of Mathematics, University of Oslo.
Date: August 7, 2020
Summary. This is a demonstration of how the Python module shenfun can be used to work with
global spectral functions in one and several dimensions.
Construction
A global spectral function $u(x)$ is represented on a one-dimensional
domain (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 a
position inside the domain. ${\hat{u}k}{k=0}^{N-1}$ are the
expansion coefficient for the series, and often referred to as the
degrees 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 the
function space is the domain, which is
specified when a function space is created. To create a function space
$T=\text{span}{T_k}_{k=0}^{N-1}$ for
the first N Chebyshev polynomials of the first kind on the default domain $[-1, 1]$,
do
End of explanation
u = Function(T)
Explanation: The function $u(x)$ can now be created with all N coefficients
equal to zero as
End of explanation
T = FunctionSpace(N, 'Chebyshev', domain=(0, 1))
Explanation: When using Chebyshev polynomials the computational domain is always
$[-1, 1]$. However, we can still use a different physical domain,
like
End of explanation
T = FunctionSpace(N, 'Chebyshev', domain=(-1, 1))
u = Function(T, val=1)
Explanation: and under the hood shenfun will then map this domain to the reference
domain through
$$
u(x) = \sum_{k=0}^{N-1} \hat{u}_k \psi_k(2(x-0.5))
$$
Approximating analytical functions
The u function above was created with only zero
valued coefficients, which is the default. Alternatively,
a Function may be initialized using a constant
value
End of explanation
import sympy as sp
x = sp.Symbol('x', real=True)
u = Function(T, buffer=4*x**3-3*x)
print(u)
Explanation: but that is not very useful. A third method to initialize
a Function is to interpolate using an analytical
Sympy function.
End of explanation
T = FunctionSpace(0, 'Chebyshev', domain=(-1, 1))
Explanation: Here the analytical Sympy function will first be evaluated
on the entire quadrature mesh of the T function space,
and then forward transformed to get the coefficients. This
corresponds to a projection to T. The projection is
Find $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 the
notation $(\cdot, \cdot)_w$ represents a weighted inner product.
In this projection $u_h$ is the Function, $u$ is the sympy function and we use sympy
to 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 the
basis 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 unknown
expansion 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 size
The number of basis functions can also be left open during creation
of the function space, through
End of explanation
u = Function(T, buffer=sp.cos(20*x))
print(len(u))
Explanation: This is useful if you want to approximate a function and
are 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 using
End of explanation
Tu = u.function_space()
print(Tu.N)
Explanation: We see that $N=45$ is required to resolve this function. This agrees
well with what is reported also by Chebfun.
Note that in this process a new FunctionSpace() has been
created under the hood. The function space of u can be
extracted using
End of explanation
T1 = FunctionSpace(0, 'Chebyshev', domain=(0, 100))
u = Function(T1, buffer=sp.besselj(0, x))
print(len(u))
Explanation: To further show that shenfun is compatible with Chebfun we can also
approximate the Bessel function
End of explanation
u = Function(T1, buffer=sp.besselj(0, x), reltol=1e-14)
print(len(u))
Explanation: which gives 83 basis functions, in close agreement with Chebfun (89).
The difference lies only in the cut-off criteria. We cut frequencies
with a relative tolerance of 1e-12 by default, but if we make this criteria
a little bit stronger, then we will also arrive at a slightly higher number:
End of explanation
%matplotlib inline
import matplotlib.pyplot as plt
Tu = u.function_space()
plt.plot(Tu.mesh(), u.backward())
Explanation: Plotting the function on its quadrature points looks
a bit ragged, though:
End of explanation
xj = np.linspace(0, 100, 1000)
plt.plot(xj, u(xj))
Explanation: To improve the quality of this plot we can instead evaluate the
function on more points
End of explanation
up = u.refine(200)
Tp = up.function_space()
plt.plot(Tp.mesh(), up.backward())
Explanation: Alternatively, we can refine the function, which simply
pads zeros to $\hat{u}$
End of explanation
print(up)
Explanation: The padded expansion coefficients are now given as
End of explanation
import numpy.polynomial.chebyshev as cheb
c = cheb.Chebyshev(u, domain=(0, 100))
Explanation: More features
Since we have used a regular Chebyshev basis above, there
are many more features that could be explored simply by going through
Numpy's Chebyshev module.
For example, we can create a Chebyshev series like
End of explanation
z = Tu.map_true_domain(cheb.chebroots(u))
Explanation: The Chebyshev series in Numpy has a wide range of possibilities,
see here.
However, we may also work directly with the Chebyshev
coefficients already in u. To find the roots of the
polynomial that approximates the Bessel function on
domain $[0, 100]$, we can do
End of explanation
z2 = z[np.where((z.imag == 0)*(z.real > 0)*(z.real < 100))].real
print(z2[:5])
Explanation: 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
End of explanation
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))
Explanation: Here np.where returns the indices where the condition is true. The condition
is that the imaginary part is zero, whereas the real part is within the
true domain $[0, 100]$.
Notice.
Using directly cheb.chebroots(c) does not seem to work (even though the
series has been generated with the non-standard domain) because
Numpy only looks for roots in the reference domain $[-1, 1]$.
We could also use a function space with boundary conditions built
in, like
End of explanation
C0 = FunctionSpace(20, 'C')
C1 = FunctionSpace(20, 'C')
T = TensorProductSpace(comm, (C0, C1))
u = Function(T)
Explanation: As we can see this leads to a function space of dimension
very similar to the orthogonal space.
The major advantages of working with a space with boundary conditions
built in only comes to life when solving differential equations. As
long as we are only interested in approximating functions, we may just
as well stick to the orthogonal spaces. This goes for Legendre as
well as Chebyshev.
Multidimensional functions
Multidimensional tensor product spaces are created
by taking the tensor products of one-dimensional function spaces.
For example
End of explanation
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())
Explanation: 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 the
1D case. We can here interpolate 2D Sympy functions
End of explanation
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')
Explanation: Like for 1D the coefficients are computed through projection,
where the exact function is evaluated on all quadrature points
in the mesh.
The Cartesian mesh represents the quadrature points of the
two function spaces, and can be visualized as follows
End of explanation
X = T.local_mesh(broadcast=True, uniform=True)
plt.contourf(X[0], X[1], u.backward(kind={'chebyshev': 'uniform'}))
Explanation: We may alternatively plot on a uniform mesh
End of explanation
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))
Explanation: Curvilinear coordinates
With shenfun it is possible to use curvilinear coordinates,
and not necessarily with orthogonal basis vectors. With
curvilinear coordinates the computational coordinates are
always straight lines, rectangles and cubes. But the physical
coordinates can be very complex.
Consider the unit disc with polar coordinates. Here
the 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
End of explanation
u = Function(T, buffer=(1-r)*r*sp.sin(sp.cos(theta)))
Explanation: Note that we are using a Fourier space for the azimuthal
direction, since the solution here needs to be periodic.
We can now create functions on the space using an
analytical function in computational coordinates
End of explanation
X = T.local_mesh(True)
plt.contourf(X[0], X[1], u.backward(), 100)
Explanation: However, when this is plotted it may not be what you expect
End of explanation
X = T.local_cartesian_mesh()
plt.contourf(X[0], X[1], u.backward(), 100)
Explanation: We see that the function has been plotted in computational coordinates,
and not on the disc, as you probably expected. To plot on
the disc we need the physical mesh, and not the computational
End of explanation
up = u.backward()
xp, yp, up = wrap_periodic([X[0], X[1], up], axes=[1])
plt.contourf(xp, yp, up, 100)
Explanation: Notice.
The periodic plot does not wrap all around the circle. This is
not wrong, we have simply not used the same point twice, but it
does not look very good. To overcome this problem we can wrap the
grid all the way around and re-plot.
End of explanation
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)
Explanation: Adaptive functions in multiple dimensions
If you want to find a good resolution for a function in multiple
dimensions, the procedure is exactly like in 1D. First create function
spaces with 0 quadrature points, and then call Function
End of explanation
<END_TASK>