Numerical Analysis: Wave PDE

GitHub
Source code is available in the GitHub repository.

Description

The project’s full title is:

Solving a Time–Dependent 1D Wave PDE
A Numerical Linear Algebraic Approach

Problem Statement

The governing PDE often used to model wave phenomena is the Wave equation: $$ \ddot{u} = c^2 \nabla^2 u $$ Our setting is one-dimensional $(\nabla = \frac{\partial}{\partial x} \widehat{\mathbf{i}})$ and we constrain $x \in [0,1]$. Thus, this PDE is equivalent to,

$$ \frac{\partial^2}{\partial t^2}u(x,t) - c^2 \frac{\partial^2}{\partial x^2}u(x,t) = 0 \, \bigg|_{x \in (0, 1) \textrm{ and } t>0} $$ Here the function $u(x,t)$ is defined as, $$ \begin{matrix} u: & [0,1] \times \R_{\geq 0} & \to & \R \\ & (x,t) & \mapsto & u(x,t) \end{matrix} $$ The PDE boundary conditions are, $$ u(0,t) = u(1,t) = 0 $$ And its initial conditions are, $$ u(x,0) = v(x) \textrm{ and } \frac{\partial u}{\partial t}(x, 0) = w(x) \, \bigg|_{x \in (0, 1)} $$ Here the functions $v, w: [0,1] \to \R$ are known.

This PDE could model, for instance, vertical oscillations over a length $1$.

Here $u(x, t)$ is thus the vertical position at point $x$ and time $t$. Likewise, $v$ and $w$ are the initial positions and initial velocities. $c \in \R^+$ is a coefficient (e.g., depends on the specific properties of the string).

For sufficiently regular functions $v$ and $w$, this PDE admits a unique solution exists of the form,

$$ u(x,t) = \sum_{k=1}^{\infty} \alpha_k \cos(k \pi c t + \phi_k) \sin (k \pi x) $$ Here $\alpha_k, \phi_k \in \R$ depend only on $v$ and $w$.

Project Objectives

The goal is to numerically solve this PDE, and thus find the numerical approximations to the exact solution $u(x, t)$ for any $x \in (0,1)$ and $t$.

This problem is discretized twice:

1. Spatial discretization in $x$ to construct a system of second-order time-dependent differential equations with appropriate conditions.
2. Temporal discretization in $t$ to construct an iterative solution scheme and its initial values.

Furthermore, a convergence analysis and MATLAB implementation are completed.

Space Discretization

Let the spatial step size be $h = \frac{1}{n + 1}$. Thus, the points $x_j = jh$ for $j \in \{1, \ldots , n\}$ define a uniform mesh on $[0, 1]$.

For a fixed $t$, $\frac{\partial^2}{\partial x^2}u(x_j, t)$ is approximated by the second-order finite difference, $$ \frac{u(x_{j-1}, t) - 2u(x_j, t) + u(x_{j+1}, t)}{h^2} $$

For $j \in \{0, \ldots, n\}$, we introduce the approximations $u_j(t)$ of the exact values $u(x_j, t)$. Thus, from the original continuous system construct a discrete approximate system $(P_h)$.

This is expressed as the following system of time-dependent second-order differential equations with boundary and initial conditions: for $j \in \{1, \ldots, n\}$ $$ (P_h) \begin{dcases} u''_j(t)+ c^2 \frac{u(x_{j-1}, t) - 2u(x_j, t) + u(x_{j+1}, t)}{h^2} = 0 \\ u_0(t) = u_{n+1}(t) = 0 \\ u_j(0) = v(x_j) \textrm{ and } u'_j(0) = w(x_j) \end{dcases} $$ This can be rewritten in matrix form. Proceed by defining the second-order finite difference matrix operator, $$ \mathrm{A} = \begin{bmatrix} 2 & -1 & & \\ -1 & \ddots & \ddots \\ & \ddots & \ddots & -1 \\ & & -1 & 2 \end{bmatrix} $$ and the vectors, $$ \mathbf{v} = \begin{bmatrix} v(x_1) \\ \vdots \\ v(x_n) \end{bmatrix}, \quad \mathbf{w} = \begin{bmatrix} w(x_1) \\ \vdots \\ w(x_n) \end{bmatrix}, \quad \mathbf{u} = \begin{bmatrix} u_1(t) \\ \vdots \\ u_n(t) \end{bmatrix}. $$ The problem is thus rewritten to its matrix form, $$ (P_h) \begin{dcases} \mathbf{u}''(t) + \frac{c^2}{h^2}\mathrm{A}\mathbf{u}(t), & \textrm{for } t>0 \\ \mathbf{u}(0) = \mathbf{v}, \;\; \mathbf{u}'(0) = \mathbf{w}. \end{dcases} $$

Time Discretization

Introduce the time step size $\Delta t$, thus the time-steps are $t_k = k \Delta t$ for $k \in \N$.

$u''_j(t_k)$ is approximated by the finite difference, $$ \frac{-2 u_{j}\left(t_{k}\right)+u_{j}\left(t_{k+1}\right)+u_{j}\left(t_{k-1}\right)}{\Delta t^{2}} $$ Using $u^k_j$ as an approximation of $u_j(t_k)$, deduce from $(P_h)$ for $j \in \{1, \ldots, n\}$ the equalities, $$ \frac{-2 u_{j}^{k}+u_{j}^{k+1}+u_{j}^{k-1}}{\Delta t^{2}}+c^{2} \frac{2 u_{j}^{k}-u_{j+1}^{k}-u_{j-1}^{k}}{h^{2}}=0 $$ Rewrite into matrix form using the notation $\mathbf{u}^k = \begin{bmatrix} u^k_1 & \cdots & u^k_n \end{bmatrix}^\mathsf{T}$, $$ \frac{-2 \mathbf{u}^{k}+\mathbf{u}^{k+1}+\mathbf{u}^{k-1}}{\Delta t^{2}}+\frac{c^{2}}{h^{2}} \mathrm{A} \mathbf{u}^{k}=0 $$ Deduce an explicit expression for $\mathbf{u}^{k+1}$, $$ \mathbf{u}^{k+1} = (2\mathrm{I} - \lambda \mathrm{A})\mathbf{u}^{k} - \mathbf{u}^{k-1}, $$ with $\lambda=\left(\frac{c \Delta t}{h}\right)^{2}$.

Initialization

Initialize the iterations by choosing values for $\mathbf{u}^{0}$ and $\mathbf{u}^{1}$, the approximations of $\mathbf{u}$ at $t = 0$ and $t = \Delta t$.

Choose $\mathbf{u}^{0} = \mathbf{v}$ since $\mathbf{u}(0) = \mathbf{v}$.

To choose a good value for $\mathbf{u}^{1}$, exploit Taylor’s theorem,

$$ u_j (\Delta t) = u_j (0) + \Delta t \, u'_j (0) + R_1(\Delta t) $$ Express the remainder $R_1$ explicitly in Lagrange mean-value form.
$\exists \xi \in (0, \Delta t) :$ $$ R_{k}(\Delta t)=\frac{u_j^{(k+1)}\left(\xi \right)}{(k+1) !}(\Delta t - 0)^{k+1} $$ For $k = 1$, obtain, $\exists \xi \in (0, \Delta t) :$ $$ u_j (\Delta t) = u_j (0) + \Delta u'_j (0) + \frac{\Delta t^2}{2}u''_j(\xi) $$ Approximate $u''_j(\xi)$ with $u''_j(0)$, which is evaluated as in $(P_h)$, thus obtain, $$ u_{j}(\Delta t) \approx u_{j}(0)+\Delta t u_{j}^{\prime}(0)+\frac{c^{2} \Delta t^{2}}{2h^{2}}\big[u_{j-1}(0)-2 u_{j}(0)+u_{j+1}(0)\big] $$ Using this approximation, choose $\mathbf{u}^1$ to be, $$ \mathbf{u}^1 = \mathbf{v} + \Delta t \, \mathbf{w} - \frac{\lambda}{2}\mathrm{A}\mathbf{v} = (\mathrm{I}-\tfrac{1}{2}\lambda \mathrm{A})\mathbf{v} + \Delta t \, \mathbf{w} $$ The iterative scheme is thus, for $k \geq 1$, $$ \begin{dcases} \mathbf{u}^{0} &= \mathbf{v}, \\ \mathbf{u}^{1} &= \big(\mathrm{I}-\tfrac{1}{2}\lambda \mathrm{A}\big)\mathbf{v} + \Delta t \, \mathbf{w}, \\ \mathbf{u}^{k+1} &= (2\mathrm{I} - \lambda \mathrm{A})\mathbf{u}^{k} - \mathbf{u}^{k-1}, \end{dcases} $$ with $\lambda=\left(\frac{c \Delta t}{h}\right)^{2}$.

Convergence

The sequence $(\mathbf{u}^k)$ must not diverge as $k$ increases. The iterative solution above is said to be stable if the sequence $(\mathbf{u}^k)$ verifies: $$ \exists C>0 : \left\Vert\mathbf{u}^k\right\Vert_\infty \leq C \quad \forall k \in \N $$ To assure stability, we verify on the time and space steps ($\Delta t$ and $h$) the Courant–Friedrichs–Lewy (CFL) convergence condition, $$ \frac{\Delta t}{h} \leq \frac{1}{c} $$ When the CFL condition is satisfied, the iterative scheme is proven to converge quadratically.

This means that, for a fixed time $T > 0$, by approaching $u(x, T)$ over $[0, 1]$ using time step $ \Delta t=\frac{T}{M} \leq \frac{h}{c} $ then the error is bound as follows, $$ \max_{1 \leq j \leq n}\left|u_{j}^{M}-u\left(x_{j}, T\right)\right| \leq K h^{2} $$ Here $K$ is a constant independent from $M$ and $n$.

MATLAB Implementation

We numerically check the CFL condition and the convergence order in the particular case $c = 10$ with the initial conditions, $$ v_\mathrm{ex}(x)=\sin (4 \pi x), \qquad w_\mathrm{ex}(x)=-40 \pi \sin (4 \pi x) $$ Here, the exact solution $u_\mathrm{ex}$ is thus, $$ u_\mathrm{ex}(x,t) = \sqrt{2} \cos \left(40 \pi t+\frac{\pi}{4}\right) \sin (4 \pi x) $$ This is periodic in time $(\frac{1}{20})$ and in space $(\frac{1}{2})$.

Implemented functions

The following functions have been implemented in MATLAB.

$\mathtt{prodmatvec(x,n)}$

$\mathtt{exinitpos(y)}$

$\mathtt{exinitvel(y)}$

$\mathtt{solex(y,t)}$

$\mathtt{sketchsolex(np,t)}$

$\mathtt{solverope(c, n, initpos,initvel,M,T)}$

Example usage:

solverope(10, 4, @exinitpos, @exinitvel, 5, 0.02)
>> ans = 
    [-0.9427  1.5253  -1.5253  0.9427]

$\mathtt{study(T)}$

Solution

The obtained solution is animated below by varying the parameter $t$ and observing changes across $x$.

References

Trefethen & Bau, Numerical Linear Algebra, 1997.

Strang, Computational Science & Engineering, 2007.