<!-- dom:TITLE: Demo - Mixed bases for the Helmholtz problem -->
# Demo - Mixed bases for the Helmholtz problem
<!-- dom:AUTHOR: Mikael Mortensen Email:mikaem@math.uio.no at Department of Mathematics, University of Oslo. -->
<!-- Author: -->  
**Mikael Mortensen** (email: `mikaem@math.uio.no`), Department of Mathematics, University of Oslo.

Date: **March 22, 2021**

**Summary.** This demo shows how to solve the Helmholtz equation using different
bases for test and trial spaces. The use of different bases leads for
some optimal combinations to highly sparse and well-conditioned
coefficient matrices.












## The Helmholtz problem

We will consider Helmholtz equation with homogeneous Dirichlet boundary conditions

<!-- Equation labels as ordinary links -->
<a id="_auto1"></a>

$$
\begin{equation}
    \alpha u -  u^{''} = f \quad \text{in} \, {I}=(-1, 1), \quad u(\pm 1) = 0,
\label{_auto1} \tag{1}
\end{equation}
$$

where $\alpha \in \mathbb{R^+}$. The relevant function space for the Dirichlet problem is

<!-- Equation labels as ordinary links -->
<a id="_auto2"></a>

$$
\begin{equation}
    S_N=\text{span}\{T_k\}_{k=0}^{N-1}, \quad V_{N} = \{v \in {S}_N\,|\, v(\pm 1) = 0\},
\label{_auto2} \tag{2}
\end{equation}
$$

and the Chebyshev-Galerkin (CG) method is to find $u_N \in V_N$ such that

<!-- Equation labels as ordinary links -->
<a id="eq:dirGalerkin"></a>

$$
\begin{equation}
\label{eq:dirGalerkin} \tag{3}
    \alpha (u_N, v)_{\omega^{\sigma}} -(u^{''}_N, v)_{\omega^{\sigma}} = (f, v)_{\omega^{\sigma}}, \forall \, v \in V_N,
\end{equation}
$$

where $(u,v)_{\omega^{\sigma}}=\int_{{I}}uv\omega^{\sigma} dx$ is the scalar product in the weighted space $L^2_{\omega^{\sigma}}({I})$.

Shenfun has implemented three different Chebyshev Dirichlet
basis functions

<!-- Equation labels as ordinary links -->
<a id="eq:shen"></a>

$$
\label{eq:shen} \tag{4}
\phi_k = T_k-T_{k+2}, \quad k=0,1, \ldots, N-3,
$$

<!-- Equation labels as ordinary links -->
<a id="eq:heinrichs"></a>

$$
\label{eq:heinrichs} \tag{5}
\varphi_k = (1-x^2)T_k, \quad k=0,1, \ldots, N-3,
$$

<!-- Equation labels as ordinary links -->
<a id="eq:dirichletU"></a>

$$
\label{eq:dirichletU} \tag{6}
\psi_k = U_k-\frac{k+1}{k+3}U_{k+2}, \quad k=0,1, \ldots, N-3.
$$

These three bases are all linearly dependent and they are all bases
for $V_N$.

## Implementation

We can get all three function spaces as


```
from shenfun import *
N = 40
V0 = FunctionSpace(N, 'C', basis='ShenDirichlet')
V1 = FunctionSpace(N, 'C', basis='Heinrichs')
V2 = FunctionSpace(N, 'C', basis='DirichletU')
```

where $V0 = \text{span}\{\phi_k\}_{k=0}^{N-3}$,
$V1 = \text{span}\{\varphi_k\}_{k=0}^{N-3}$ and
$V2 = \text{span}\{\psi_k\}_{k=0}^{N-3}$. Now, to solve the Helmholtz problem we simply need to choose
test and trial bases. Shen's original method is using
`V0` for both. To assemble the stiffness and mass matrices
for this choice do


```
u = TrialFunction(V0)
v = TestFunction(V0)
A = inner(v, div(grad(u)))
B = inner(v, u)
```

A manufactured solution can be chosen using [Sympy](https://www.sympy.org)
We choose

<!-- Equation labels as ordinary links -->
<a id="_auto3"></a>

$$
\begin{equation}
u(x) = \sin \left( 2 \pi \cos \left( 2 \pi x \right) \right)
\label{_auto3} \tag{7}
\end{equation}
$$

implemented as


```
import sympy as sp
x = sp.Symbol('x', real=True)
ue = sp.sin(2*sp.pi*sp.cos(2*sp.pi*x))
```

The right hand side $f$ of Helmholtz equation is


```
alpha = 1
f = sp.simplify(alpha*ue-ue.diff(x, 2))
```

To solve the problem we can do


```
fj = Array(V0, buffer=f)  # Get f on quadrature mesh
f_hat = inner(v, fj)      # Compute right hand side
M = alpha*B - A           # Get coefficient matrix
u_hat = Function(V0)      # Container for the solution
sol = la.Solver(M)        # Solver
u_hat = sol(f_hat, u_hat) # Solve
```

Compare with exact solution.


```
uj = Array(V0, buffer=ue)
error = inner(1, (u_hat.backward()-uj)**2)
print('Error =', error)
```

Now that was the solution for test and trial bases from the same
basis ([4](#eq:shen)). Let us create a function that takes any
test and any trial basis, any manufactured solution and any $\alpha$
in the Helmholtz equation. We let the function return either
the L2-error norm, the condition number of the Helmholtz
coefficient matrix, or the matrix itself.


```
def main(N, test, trial, alpha=1, method=0, ue=sp.sin(2*sp.pi*sp.cos(2*sp.pi*x))):
    """Solve Helmholtz problem and return L2-error, condition number or matrix

    Parameters
    ----------
    N : int
        Number of quadrature points
    test, trial : int
        Test and trial functions.
        0 = :math:`T_k-T_{k+2}`
        1 = :math:`(1-x^2)T_k`
        2 = :math:`U_k-\frac{k+1}{k+3}U_{k+2}`
    alpha : Helmholtz parameter
    method : int
        0 = Return L2-error norm
        1 = Return condition number of matrix
        2 = Return Helmholtz matrix
    ue : Sympy function, optional
        The manufactured solution with homogeneous boundary conditions.

    Note
    ----
    Inhomogeneous boundary conditions require a small rewrite, but is
    not difficult.

    """
    bases = {0: 'ShenDirichlet', 1: 'Heinrichs', 2: 'DirichletU'}
    test = FunctionSpace(N, 'C', basis=bases[test])
    trial= FunctionSpace(N, 'C', basis=bases[trial], bc=(0, 0))
    # Check that boundary conditions are homogeneous
    assert abs(ue.subs(x, -1)) < 1e-8 and abs(ue.subs(x, 1)) < 1e-8
    u = TrialFunction(trial)
    v = TestFunction(test)
    f = sp.simplify(alpha*ue-ue.diff(x, 2))
    fj = Array(test, buffer=f) # Get f on quadrature mesh
    f_hat = inner(v, fj)      # Compute right hand side
    B = inner(v, u)
    A = inner(v, div(grad(u)))
    M = alpha*B-A
    if method == 1:
        return np.linalg.cond(M.diags('csr').toarray())
    if method == 2:
        return M

    u_hat = Function(trial)
    sol = la.Solver(M)
    u_hat = sol(f_hat, u_hat)
    uj = Array(trial, buffer=ue)
    error = np.sqrt(inner(1, (u_hat.backward()-uj)**2))
    return error
```

Let us first try basis ([4](#eq:shen)) as test function and
([5](#eq:heinrichs)) as trial function. Use otherwise
default parameters.


```
error = main(100, 0, 1)
print(error)
```

So the error is small in deed. Perhaps more interesting, let's
look at the sparsity pattern of the coefficient matrix


```
M = main(100, 0, 1, method=2)
import plotly.express as px
z = np.where(abs(M.diags().toarray()) > 1e-6, 0, 1).astype(bool)
fig = px.imshow(z, binary_string=True)
fig.show()
#plt.spy(M.diags(), markersize=0.2) # or use matplotlib
```

The coefficient matrix has 4 non-zero diagonals. You can now experiment
with different test and trial functions, but you will not get a better
result than that. Try basis ([5](#eq:heinrichs)) for both test and trial
function, and you'll get 5 nonzero diagonals.

To see the convergence rate call `main` for a range of
different mesh sizes


```
error = []
N = (2**4, 2**6, 2**8, 2**10)
for n in N:
    error.append(main(n, 0, 1))
fig = px.line(x=N, y=error, log_y=True)
fig.update_layout(yaxis=dict(showexponent='all', exponentformat='e'))
fig.show()
```

<!-- ======= Bibliography ======= -->

1900a9b49aded9d5d17233251696c407f81138e0
ipynb
Jupyter Notebook
content/mixingbases.ipynb
mikaem/shenfun-demos
c2ad13d62866e0812068673fdb6a7ef68ecfb7f2
BSD-2-Clause
content/mixingbases.ipynb
mikaem/shenfun-demos
c2ad13d62866e0812068673fdb6a7ef68ecfb7f2
BSD-2-Clause
2021-09-21T16:10:01.000Z
2021-09-21T16:10:01.000Z
content/mixingbases.ipynb
mikaem/shenfun-demos
c2ad13d62866e0812068673fdb6a7ef68ecfb7f2
BSD-2-Clause
Qwen/Qwen-72B
1. YES
2. YES