cds-numerical-methods/Week 6/10 Hyperbolic PDEs.ipynb

CDS: Numerical Methods Assignments¶

• See lecture notes and documentation on Brightspace for Python and Jupyter basics. If you are stuck, try to google or get in touch via Discord.

• Solutions must be submitted via the Jupyter Hub.

• Make sure you fill in any place that says YOUR CODE HERE or "YOUR ANSWER HERE".

Submission¶

1. Name all team members in the the cell below
2. make sure everything runs as expected
3. restart the kernel (in the menubar, select Kernel$\rightarrow$Restart)
4. run all cells (in the menubar, select Cell$\rightarrow$Run All)
5. Check all outputs (Out[*]) for errors and resolve them if necessary
6. submit your solutions in time (before the deadline)

team_members = "Koen Vendrig, Kees van Kempen"

Dynamic Behavior from Hyperbolic PDEs¶

In the following you will derive and implement finite-difference methods to study generalized hyperbolic partial differential equations.

In [15]:

import numpy as np
from matplotlib import pyplot as plt

Task 1: finite-difference hyperbolic PDE solver¶

Our aim is to implement a Python function to find the solution of the following PDE:

\begin{align*} \frac{\partial^2}{\partial t^2} u(x,t) - \alpha^2 \frac{\partial^2}{\partial x^2} u(x,t) = 0, \qquad 0 \leq x \leq l, \qquad 0 \leq t \end{align*}

with the boundary conditions

\begin{align*} u(0,t) = u(l,t) &= 0, &&\text{for } t > 0 \\ u(x,0) &= f(x) \\ \frac{\partial}{\partial t} u(x,0) &= g(x) &&\text{for } 0 \leq x \leq l. \end{align*}

By approximating the partial derivatives with finite differences we can recast the problem into the following form

\begin{align*} \vec{w}_{j+1} = \mathbf{A} \vec{w}_{j} - \vec{w}_{j-1}. \end{align*}

Here $\vec{w}_j$ is a vector of length $m$ in the discretized spatial coordinate $x_i$ at time step $t_j$. The spatial coordinates are defined as $x_i = i h$, where $i=0,1,\dots,m-1$ and $h = l/(m-1)$. The time steps are defined as $t_j = j k$, where $j = 0, 1, \dots, n-1$.

The tri-diagonal matrix $\mathbf{A}$ has size $(m-2)\times(m-2)$ and is defined by

\begin{align*} \mathbf{A} = \left( \begin{array}{cccc} 2(1-\lambda^2) & \lambda^2 & & 0\\ \lambda^2 & \ddots & \ddots & \\ & \ddots & \ddots & \lambda^2 \\ 0 & & \lambda^2 & 2(1-\lambda^2) \end{array} \right), \end{align*}

where $\lambda = \alpha k / h$. This $(m-2)\times(m-2)$ structure accounts for the first set of boundary conditions. Note that the product $\mathbf{A} \vec{w}_{j}$ is thus only performed over the $m-2$ subset, i.e. $i=1,2,\dots,m-2$. The other boundary conditions are accounted for by initializing the first two time steps with

\begin{align*} w_{i,j=0} &= f(x_i) \\ w_{i,j=1} &= (1-\lambda^2) f(x_i) + \frac{\lambda^2}{2} f(x_{i+1}) + \frac{\lambda^2}{2} f(x_{i-1}) + kg(x_i). \end{align*}

Implement a Python function of the form $\text{pdeHyperbolic(a, x, t, f, g)}$, where $\text{a}$ represents the PDE parameter $\alpha$, $\text{x}$ and $\text{t}$ are the discretized spatial and time grids, and $\text{f}$ and $\text{g}$ are the functions defining the boundary conditions. This function should return a two-dimensional array $\text{w[:,:]}$, which stores the spatial vector $\vec{w}_j$ at each time coordinate $t_j$.

In [35]:

def pdeHyperbolic(a, x, t, f, g):
    # TODO: Docstring.
    
    n = len(t)
    m = len(x)
    
    # TODO: The stepsize is defined by fixed h and, but the input x and t
    #       could allow variable step size. What should it be?
    #h = x[1] - x[0]
    l = x[-1]
    h = l/(m - 1)
    k = t[1] - t[0]
    λ = a*k/h
    
    A = np.eye(m - 2)*2*(1 - λ**2) + ( np.eye(m - 2, m - 2, 1) + np.eye(m - 2, m - 2, -1) )*λ**2
    
    w = np.zeros((n, m))
    
    # Set initial values for w[i, 0] and w[i, 1]
    w[0] = f(x)
    w[1, 1:m - 1] = (1 - λ**2)*f(x[1:m - 1]) + λ**2/2*f(x[2:m]) - λ**2/2*f(x[0:m - 2]) + k*g(x[1:m - 1])
    
    for j in range(2, n):
        w[j, 1:m - 1] = np.dot(A, w[j, 1:m - 1]) - w[j - 1, 1:m - 1]
    
    return w

Task 2¶

Use your implementation to solve the following problems. Compare in the first problem your numerical solution to the analytic one and use it to debug your code.

Problem 1:¶

\begin{align*} \frac{\partial^2}{\partial t^2} u(x,t) - \frac{\partial^2}{\partial x^2} u(x,t) &= 0, &&\text{for }0 \leq x \leq 1 \text{ and } 0 \leq t \leq 1\\ u(0,t) = u(l,t) &= 0, &&\text{for } t > 0 \\ u(x,0) &= \operatorname{sin}\left(2 \pi x \right), \\ \frac{\partial}{\partial t} u(x,0) &= 2 \pi \operatorname{sin}\left(2 \pi x \right) &&\text{for } 0 \leq x \leq 1. \end{align*}

Compare your numerical solution to the analytic one and use it to debug your code. The corresponding analytic solution is

\begin{equation} u(x,t) = \operatorname{sin}(2 \pi x) (\operatorname{cos}(2 \pi t) + \operatorname{sin}(2 \pi t)). \end{equation}

Then write a unit test for your function based on this problem.

In [39]:

# There is an 𝑙 in the boundary conditions we assume should be a 1. (-1 pt)
a = 1
x = np.linspace(0, 1, 300)
t = np.linspace(0, 1, 300)
f = lambda x: np.sin(2*np.pi*x)
g = lambda x: 2*np.pi*np.sin(2*np.pi*x)

w = pdeHyperbolic(a, x, t, f, g)
u = lambda x, t: np.sin(2*np.pi*x)*(np.cos(2*np.pi*t) + np.sin(2*np.pi*t))

#plt.plot(t[1:], w[1:, ])
sub = np.arange(1, 300)
#plt.plot(t[sub], w[sub, sub])
#plt.plot(t, u(x, t))
print(np.max(w))

# Now we hope that w[j, i] \approx u(x[i], t[j]).
jlist = np.arange(1, 300)
ilist = np.arange(1, 300)

plt.plot(t[jlist], u(x[ilist], t[jlist]))
plt.plot(t[jlist], w[jlist, ilist])

0.9999862004036565

Out[39]:

[<matplotlib.lines.Line2D at 0x14e5e1970610>]

In [ ]:

def test_pdeHyperbolic():
    # YOUR CODE HERE
    raise NotImplementedError()
    
test_pdeHyperbolic()

Problem 2:¶

\begin{align*} \frac{\partial^2}{\partial t^2} u(x,t) - \frac{\partial^2}{\partial x^2} u(x,t) &= 0, &&\text{for }0 \leq x \leq 1 \text{ and } 0 \leq t \leq 2\\ u(0,t) = u(l,t) &= 0, &&\text{for } t > 0 \\ u(x,0) &= \left\{ \begin{array}{cc} +1 & \text{for } x<0.5 \\ -1 & \text{for } x \geq 0.5 \end{array} \right., \\ \frac{\partial}{\partial t} u(x,0) &= 0 &&\text{for } 0 \leq x \leq 1. \end{align*}

Use $m=200$ and $n=400$ to discretize the spatial and time grids, respectively.

In [ ]:

# YOUR CODE HERE
raise NotImplementedError()

Task 3¶

Animate your solutions! To this end you can use the following code:

# use matplotlib's animation package
import matplotlib.pylab as plt
import matplotlib
import matplotlib.animation as animation
# set the animation style to "jshtml" (for the use in Jupyter)
matplotlib.rcParams['animation.html'] = 'jshtml'

# create a figure for the animation
fig = plt.figure()
plt.grid(True)
plt.xlim( ... )     # fix x limits
plt.ylim( ... )     # fix y limits

# Create an empty plot object and prevent its showing (we will fill it each frame)
myPlot, = plt.plot([0], [0])
plt.close()

# This function is called each frame to generate the animation (f is the frame number)
def animate(f):                   
    myPlot.set_data( ... )  # update plot

# Show the animation
frames = np.arange(1, np.size(t))  # t is the time grid here
myAnimation = animation.FuncAnimation(fig, animate, frames, interval = 20)
myAnimation

In [ ]:

# Animate problem 1 here ...

# YOUR CODE HERE
raise NotImplementedError()

In [ ]:

# Animate problem 2 here ...

# YOUR CODE HERE
raise NotImplementedError()