#!/usr/bin/env python
# coding: utf-8
# # Hyperbolic PDE
#
#
# Copyright (C) 2010-2020 Luke Olson
# Copyright (C) 2020 Andreas Kloeckner
#
#
# MIT License
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
#
#
# -----
#
# One-way wave equation:
# $$ u_t + c u_x = 0$$
# $$ u(0,x) = f(x) \quad \text{on} \, [0,1] $$
#
# (Periodic boundary conditions)
#
# ## Solution
#
# $$ u(t,x) = f(x-ct) $$
# In[26]:
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
# ## Step 2: Set up grid
# - `dx`: grid spacing in the $x$-direction
# - `x`: grid coordinates
# In[27]:
nx = 100
x = np.linspace(0, 1, nx, endpoint=False)
dx = x[1] - x[0]
# ## Set time/speed parameters
# In[33]:
c = 1.0 # speed
T = 1.0 / c # end time
lmbda = 0.93
dt = dx * lmbda / c
nt = int(T//dt)
print('T = %g' % T)
print('tsteps = %d' % nt)
print(' dx = %g' % dx)
print(' dt = %g' % dt)
print('lambda = %g' % lmbda)
# ## Set Initial Condition
# (Square wave with amplitude 1)
# In[34]:
def f(x):
return ((x>0.4) & (x<0.6)).astype(np.float64)
# ## Plot solution
# In[42]:
import time
u = f(x)
fig = plt.figure(figsize=(8,8))
plt.title('u vs x')
line1, = plt.plot(x, u, lw=3, clip_on=False)
def timestepper(n):
uex = f((x - c * (n+1) * dt) % 1.0)
line1.set_data(x, uex)
return [line1]
from matplotlib import animation
from IPython.display import HTML
ani = animation.FuncAnimation(fig, timestepper, frames=nt, interval=20)
html = HTML(ani.to_jshtml())
plt.clf()
html
# In[ ]: