Demo - Poisson equation 1D
=======================

Solve Poisson equation in 1D with homogeneous Dirichlet bcs

$$    
\begin{align}
\nabla^2 u(x) &= f(x), \quad \forall \, x \in [-1, 1]\\
u(\pm 1) &= 0 
\end{align}
$$

Use either Chebyshev basis $P=\{T_k(x)\}_{k=0}^{N-1}$ or Legendre $P=\{L_k(x)\}_{k=0}^{N-1}$ and define Shen's composite basis as

$$
V^N = \{P_k - P_{k+2}\, | \, k=0, 1, \ldots, N-3\}.
$$

We also need the weighted inner product

$$
 (u, v)_w = \int_{-1}^1 u v w \, \text{dx},
$$

where $w(x)$ is a weight associated with the chosen basis. For Legendre the weight is simply $w(x)=1$, whereas for Chebyshev it is $w(x)=1/\sqrt{1-x^2}$. The weight ensures that inner products are orthogonal. We use quadrature to solve the integral

$$
\int_{-1}^1 u v w \, \text{dx} = \sum_{i=0}^{N-1} u(x_i) v(x_i) \omega_i,
$$

where $\{\omega_i\}_{i=0}^{N-1}$ are the quadrature weights associated with the chosen basis (C/L) and quadrature rule. For Chebyshev there are `Chebyshev-Gauss` or `Chebyshev-Gauss-Lobatto`, whereas for Legendre there are `Legendre-Gauss` or `Legendre-Gauss-Lobatto`. The quadrature points are denoted as $\{x_i\}_{i=0}^{N-1}$ and their value depends on the chosen basis and quadrature rule.

With the weighted product in place we can create variational problems from the original PDE by multiplying with a test function and integrating over the domain. For a Legendre basis we get: 

Find $u \in V^N$ such that

$$     (\nabla u, \nabla v) = -(f, v), \quad \forall \, v \in V^N.$$

For a Chebyshev basis the integration by parts is complicated due to the weight and the variational problem used is instead: 

Find $u \in V^N$ such that

$$     (\nabla^2 u, v)_w = (f, v)_w, \quad \forall \, v \in V^N.$$

We now break the problem down to linear algebra. With any choice of basis or quadrature rule we use $\phi_k(x)$ to represent the test function $v$ and thus

$$
\begin{align}
v(x) &= \phi_k(x) \\
u(x) &= \sum_{j=0}^{N-3} \hat{u}_j \phi_j(x)
\end{align}
$$

Insert into the variational problem for Legendre

$$
\begin{align}
(\nabla \sum_{j=0}^{N-3} \hat{u}_j \phi_j, \nabla \phi_k) &= -(f, \phi_k) \\
\sum_{j=0}^{N-3} \underbrace{(\nabla \phi_j, \nabla \phi_k)}_{a_{kj}} \hat{u}_j  &= -\underbrace{(f, \phi_k)}_{\tilde{f}_k} \\
A \boldsymbol{\hat{u}} &= -\boldsymbol{\tilde{f}}
\end{align}
$$

where $A = (a_{kj})_{0 \ge k, j \ge N-3}$ is the stiffness matrix, $\hat{\boldsymbol{u}} = \{\hat{u}_j\}_{j=0}^{N-3}$ and $\boldsymbol{\tilde{f}} = \{\tilde{f}_k\}_{k=0}^{N-3}$.



```python
from shenfun import *
import matplotlib.pyplot as plt

N = 100
V = Basis(N, 'Legendre', quad='LG', bc=(0, 0))
x, w = V.points_and_weights()
v = TestFunction(V)
u = TrialFunction(V)
A = inner(grad(u), grad(v))
```

Use a manufactured solution to create a right hand side $f(x)$


```python
import sympy
x = sympy.symbols('x')
u = (1-x**2)*(sympy.cos(4*x)*sympy.sin(6*x))
f = u.diff(x, 2)
```


```python
fl = sympy.lambdify(x, f, 'numpy')
fj = Array(V)
```


```python
fj[:] = fl(V.mesh())
```


```python
f_tilde = inner(v, -fj)
```


```python
u_hat = Function(V)
u_hat = A.solve(f_tilde, u=u_hat)
```


```python
uj = u_hat.backward()
```


```python
ul = sympy.lambdify(x, u, 'numpy')
ue = ul(V.mesh())
```


```python
T = Basis(N, 'L', quad='LG')
```


```python
du = project(Dx(u_hat, 0, 1), T)
```

5193080b970c5a529ec25a334109f25a00ba0542
ipynb
Jupyter Notebook
index.ipynb
mikaem/conda
b6e8285360226288e31336bcc45323aa63f05ce7
BSD-3-Clause
index.ipynb
mikaem/conda
b6e8285360226288e31336bcc45323aa63f05ce7
BSD-3-Clause
index.ipynb
mikaem/conda
b6e8285360226288e31336bcc45323aa63f05ce7
BSD-3-Clause
Qwen/Qwen-72B
1. YES
2. YES