Note index \(j\) in \(\ell_j(x)\), the \(j\)’th basis function.
Numerator contains the product of all differences \((x-x^m)\) for all \(m\) except \(m=j\)
Denominator contains the product of all differences \((x^j-x^m)\) for all \(m\) except \(m=j\)
We can write
\[
\ell_j(x) = \prod_{\substack{0 \le m \le k \\ m \ne j}} \frac{x-x^m}{x^j-x^m}.
\]
Lagrange basis in Sympy
A basis is a collection of basis functions
\[
\{\ell_j(x)\}_{j=0}^k
\]
def Lagrangebasis(xj, x=x):"""Construct Lagrange basis for points in xj Parameters ---------- xj : array Interpolation points (nodes) x : Sympy Symbol Returns ------- Lagrange basis as a list of Sympy functions """from sympy import Mul n =len(xj) ell = [] numert = Mul(*[x - xj[i] for i inrange(n)])for i inrange(n): numer = numert/(x - xj[i]) denom = Mul(*[(xj[i] - xj[j]) for j inrange(n) if i != j]) ell.append(numer/denom)return ell
The basis for linear interpolation of our first example is
ell = Lagrangebasis(xj[3:5], x=x)ell
[4.0 - 5.0*x, 5.0*x - 3.0]
We can plot these two linear functions between \(x_3\) and \(x_4\)
With the basis we can now compute the Lagrange interpolation polynomial
\[
L(x) = \sum_{j=0}^k u^j \ell_j(x)
\]
def Lagrangefunction(u, basis):"""Return Lagrange polynomial Parameters ---------- u : array Mesh function values basis : tuple of Lagrange basis functions Output from Lagrangebasis """ f =0for j, uj inenumerate(u): f += basis[j]*ujreturn f
L = Lagrangefunction(u[3:5], ell)xf = np.linspace(xj[3]-0.1, xj[4]+0.1, 10)plt.figure(figsize=(6, 3))plt.plot(xj, u, 'bo')plt.plot(xf, sp.lambdify(x, L)(xf), 'g') plt.plot(0.75, L.subs(x, 0.75), 'ko')
We can choose any, and any number of basis functions
All Lagrange basis functions are such that on the chosen mesh points \(\{x^{i}\}_{i=0}^{k}\) we have
\[
\ell_j(x^i) = \delta_{ij} = \begin{cases} 1 \text{ for } i=j \\
0 \text{ for } i\ne j
\end{cases}
\]
Note
The Lagrange basis functions do not depend on the mesh function values, only on the mesh points. The mesh points do not need to be uniform, any mesh will do as long as all points are different.
Interpolation of derivatives
The solution has been obtained in mesh points \(\boldsymbol{u}=(u(x_j))_{j=0}^N\). How can we compute \(u'(x)\) for any given \(x\)?
Can we use the derivative matrix \(D^{(1)}\) and thus
\[
\boldsymbol{f} = D^{(1)} \boldsymbol{u}
\]
such that the mesh function \(\boldsymbol{f}=(u'(x_j))_{j=0}^N\)?
Yes, but you then need to interpolate \(\boldsymbol{f}\)!
Is there a better way?
Yes! Just take the derivative of the Lagrange function
For the numerical solution \(u\) and the exact solution \(u^{e}\) the \(L^2\) error norm can in general, for any domain \(\Omega\) and number of dimensions, be defined as
which is very small because the Lagrange interpolator happens to be close to exact exact for all the midpoints for this function. It is not zero in general.
Numerical integration can also, for example, use the trapezoidal or Simpson methods
We use an initial condition \(u(x, y, 0) = I(x, y)\) and \(\frac{\partial u}{\partial t}(x, y, 0)=0\) and homogeneous Dirichlet boundary conditions \(u(x, y, t)=0\) for the entire boundary.
For 2D waves we have frequencies in two directions
Any solution to the wave equation can be written as combinations of waves
\[
u(x, y, t) = g(k_x x + k_y y - |\boldsymbol{k}|ct)
\]
Note
Try to validate by inserting for \(u(x, y, t)= g(k_x x + k_y y - |\boldsymbol{k}|ct)\) into the wave equation. Use for example a variable \(\alpha = k_x x + k_y y - |\boldsymbol{k}|ct\) and the chain rule, such that \(\frac{\partial g}{\partial t} = \frac{\partial g}{\partial \alpha} \frac{\partial \alpha}{\partial t}\) etc.
Discretization of the wave equation in 2D
For all internal points we have the second order accurate \[
\frac{u^{n+1}_{i,j} - 2u^n_{i,j} + u^{n-1}_{i, j}}{\Delta t^2} =
c^2 \left(\frac{u^n_{i+1,j} - 2u^n_{i,j} + u^n_{i-1, j}}{\Delta x^2} + \frac{u^n_{i,j+1} - 2u^n_{i,j} + u^n_{i, j-1}}{\Delta y^2}\right)
\]
or on vectorized form using \(U^n = \Big( u(x_i, y_j, n \Delta t)\Big)_{i,j=0}^{N, N}\) and \(D_x^{(2)}\) and \(D_y^{(2)}\) as second derivative matrices for \(x\) and \(y\) directions, respectively
The initial condition \(\frac{\partial u}{\partial t}(x, y, 0) = 0\) can be implemented like we did in lecture 5 for the wave equation with one spatial dimension. We use a ghost node:
where \(\omega = |\boldsymbol{k}| c\) and \(\hat{\imath}=\sqrt{-1}\).
Note
\(\boldsymbol{x} = x \boldsymbol{i} + y \boldsymbol{j}\) and \(\boldsymbol{k} = k_x \boldsymbol{i} + k_y \boldsymbol{j}\), such that \(\boldsymbol{k} \cdot \boldsymbol{x} = k_x x + k_y y\).
A mesh solution is then
\[
\begin{align}
u^n_{ij} &= u(x_i, y_j, t_n) = e^{\hat{\imath} (i k_x \Delta x + j k_y \Delta y + \omega n \Delta t)} \\
&=(e^{\hat{\imath} \omega \Delta t})^n e^{\hat{\imath} (i k_x \Delta x + j k_y \Delta y)} \\
&= A^n E(i, j)
\end{align}
\]
with amplification factor \(A=e^{\hat{\imath} \omega \Delta t}\) and \(E(i, j) = e^{\hat{\imath} (i k_x \Delta x + j k_y \Delta y)}\).
To calculate the stability of this scheme we need to insert for \(u^n_{ij} = A^n E(i, j)\) and compute the amplification factor \(A\). For stability we require
\[
|A| \le 1
\]
This is the same stability criterion as used in the 1D case.
Insert for \(\small u^n_{ij} = A^n E(i, j)\) in the discretized wave equation
