WaveCam
Solving wave equations in realtime on the webcam buffer
WaveCam
An interactive wave-equation solver
This is a short guide on how to solve the wave equation using live images from the web camera as Dirichlet boundary conditions. We will use a finite difference scheme in both time and space, and the solver will be implemented in python. To easily access the web camera we will use the openCV (open computer vision) library and the calculations are optimized in numpy.
The reader should have some basic understanding of differential equations and the python programming language.
The solution scheme and the implementation is inspired by the work of Hans Petter Langtangen as presented in the course INF5620 - Numerical solution of PDEs.
The wave equation
We begin by stating the two dimensional wave equation with no damping, uniform velocity throughout the medium and a source term: \[(\frac{\partial}{\partial t})^2 u(x,y,t) + c^2 \nabla^2 u(x,y,t) + f(x,y,t) = 0\]
The source term $f(x,t)$ may be considered an outside force interacting with the medium, like a point oscillator. In acoustics, such a term would represent speakers and other sound sources.
The quantity $c$ is to be considered the wave velocity.
Solving the wave equation
To solve the this problem, we need to know how the initial wavefunction looks and its first derivative with respect to time. This initial conditions may be written as \[u(t_0) = u_0 = I(x,y)\] \[\frac{\partial}{\partial t} u_0 = V(x,y)\]
We need to pay special attention to what goes on at the boundaries of the domain. Two boundary conditions is commonly adressed:
Dirichlet: \(u\vert_b = b(x)\)
Neumann: \(\partial u\vert_b = c(x)\)
The Dirchlet condition forces the value of u to remain some constant distribution at the boundary. Consider for example a membrane on a drum, or a string on a guitar. The wave travelling in such a medium will never be able to change the values at the boundary, unless the system somehow breaks down.
The Neumann boundary will either cause positive,negative or zero netto influx at the boundaries. A perfectly reflecting boundary, such as a the walls of an imaginary cube half-filled with water, will have zero netto flux at the boundaries. A common way to implement this is by setting some virtual function value outside the domain to equal the one bordering on the inside. This is often called ghost cells.
Discretization of the domain
The first step in finding a numerical solution to the equation is to discretize the domain.
We let \[x \rightarrow x_i =i \Delta x\] \[y \rightarrow y_j =j \Delta y\] \[t \rightarrow t_k =k \Delta t\]
The discrete solution will be defined at all discrete coordinates of the domain. We may introduce a more convenient notation for the function value at these coordinate: \[u(x,y,t) \rightarrow u(x_i,y_j,t_k) = u_{ijk}\]
We will need to consider algebraic expressions with recursive patterns in the indexing. When doing this a further simplification of the notation may be introduced: \[u_{i+1,jk} = u_{i+1}\]
so that omitted indices imply that they are unshifted. (Or not of interest in the part we are working on.)
We are now ready to discretize the problem.
##Discretizing the problem
The Laplace-operator in two dimensions is \[\nabla^2 = (\frac{\partial}{\partial x})^2 + (\frac{\partial}{\partial y})^2\]
Discretizing the problem therefore means to find a numerical approximation to the double differentials in space and time. One simple and commonly employed such approximation is \[(\frac{\partial}{\partial x})^2 u(x) \approx \frac{u_{i-1} - 2u_i + u_{i+1}}{\Delta x^2}\]
By using this approximation in our equation we find that \[(\frac{\partial}{\partial t})^2 u(x,y,t) = c^2(\frac{u_{i-1} - 2u_i + u_{i+1}}{\Delta x^2} + \frac{u_{j-1} - 2u_j + u_{j+1}}{\Delta y^2}) + f(x,y,t)\]
In our case, we will set \[\Delta x = \Delta y\]
Notice also that \[u_i = u_j\]
This lets us rewrite the expression for the spacially double differential: \[(\frac{\partial}{\partial t})^2 u(x,y,t) = \frac{c^2}{\Delta x^2}(u_{i-1} + u_{i+1} + u_{j-1} + u_{j+1} - 4u_i) + f(x,y,t)\]
The expression in the parenthesis will need to be evaluated for each coordinate in the array. A visualization of this algorithmic pattern is plotted in the cell below.
#Visualizing the algorithmic kernel
from matplotlib.pyplot import *
figure(1)
hold("on")
plot([1,1],[0,3], color = "Black")
plot([2,2],[0,3], color = "Black")
plot([0,3],[1,1], color = "Black")
plot([0,3],[2,2], color = "Black")
text(1.1,1.4, "$-4u_i$", size = 30)
text(1.2,2.4, "$u_{j-1}$", size = 30)
text(1.2,0.4, "$u_{j+1}$", size = 30)
text(0.2,1.4, "$u_{i-1}$", size = 30)
text(2.2,1.4, "$u_{i+1}$", size = 30)
plt.axis('off')
hold("off")
show()
It is apparent from the pattern above why special attention is needed at the boundaries, since the function values outside the domain are undefined.
We are now done with the differentiation in space, and need to treat the time differentials. To simplify the expressions, we define the “right hand side”: \[RHS(u,f) = \frac{c^2}{\Delta x^2}(u_{i-1} + u_{i+1} + u_{j-1} + u_{j+1} - 4u_i) + f(x,y,t)\]
Ultimately, we seek a scheme that lets us evaluate our function at the next time step: \[u_{k+1} = F(u_{k}, u_{k-1})\]
By again approximating the double derivative, we find \[\frac{u_{k-1} - 2u_k + u_{k+1}}{\Delta t^2} = RHS(u,f)\]
By solving for $u_{k+1}$, we find \[u_{k+1} = RHS(u,f)\Delta t^2 + 2u_k - u_{k-1}\]
And we have arrived at an explicit scheme for the next time step.
Discretizing the initial conditions
At the first time step, we know $u_k$ but not $u_{k-1}$. However, we know V, and may therefore use a centered difference to get the simulation going: \[V_0 = \frac{u^{k+1}-u^{k-1}}{2\Delta t}\] \[\rightarrow u^{k-1} = u^{k+1}-2\Delta t V_0\]
If we insert this in the expression for the first step of the wave eqation, we find that \[u_{k+1} = \frac{1}{2} RHS(u,f)\Delta t^2 + u_k - \Delta t V_0\]
This is the modified first step, depending only upon the current k.
Implementing the solution scheme in python
A general implementation of the wave solver should take as input the initial conditions \(I(x,y)\) and \(V(x,y)\) (optionally in the form of arrays or functions), possibly a variable velocity \(c^2(x,y)\) (this will make RHS calculation only slightly more tricky), and the information needed to set up the mesh.
It is convenient to implement the solver as a class.
A possible “skeleton” for such a class is outlined below:
class fdwave():
def __init__(self, x,y,rho,b,q,I,V,f, dt):
"""
Initialize the solver with the following parameters
x [array containing the discretization in the x-direction] -> dx = x[1]-x[0] (for uniform mesh)
y [array containing the discretization in the y-direction]
rho [optional parameter if for damping]
b []
q [velocity squared, could be constant or a function of x,y]
I [Initial condition u_0]
V [Initial velocity]
f [Optional forcing term]
dt [time step]
Set up
u, u_p, u_pp
dx, dy
"""
def first_step(self):
"""
Take modified first step using V
"""
def advance(self):
"""
One ordinary step in time
"""
def rhs(self):
"""
Set up the right hand side (The spatial part)
"""
def solve(self, NT):
"""
Perform NT timesteps including a modified first one.
Advance the solution to NT*dt
- Calculate rhs
- Calculate u^k+1
- Impose boundary conditions
- Update arrays
"""
def impose_boundaries(self):
"""
Impose Dirichlet or Neumann conditions
"""
For a more detailed implementation, see wavecam.py.
Boundary conditions
A common choice for boundaries when doing finite difference based wave simulations, is to solve the equation for a rectangular or square mesh, using either a dirichlet conditions for constant 0 boundaries (perfect inverted reflection) or a Neumann condition for zero differential at the boundaries (a perfect non-inverted reflection).
The corresponding real-life analogy is respectively a rectangular 2D membrane fixed at the borders (for inverted reflections) or a liquid surface in a cubic tank (for non-inverted reflections).
Needless to say, these are exeptions from most everyday experiences with waves, reflected by irregular boundaries in the natural environment, and experiencing wave phenomenas through such simulations may create an impression that standing waves is a more common phenomena than it actually is. Another consequence of such boundaries is that the geometric distribution of waves is very regular and predictable, in contrast to the intermingled, complex structure of interfering waves.
A simple implementation of irregular boundaries with a high degree of interactivity (flexibility) would therefore have some pedagogical value as well as a very entertaining aspect, as the experimenter will be able to interact with the waves directly.
This could obviously easily be achieved by just stamping your feet in any nearby pond, but an even greater degree of insight might be gained by replicating the phenomena directly on your computer. “What it cannot create, I cannot understand”, as attributed to Richard Feynman.
We therefore seek
- A real-time (or live) wave solver
- A way of dynamically changing or adjusting the boundary conditions.
The approach we follow from here on is to utilize the webcamera, native to most modern laptops, to evaluate the boundaries.
The simplification of the geometry of the boundaries is mainly due to the often cumbersome way of creating irregular bounded arrays.
Open CV (Computer Vision) and webcam boundaries
Open CV is a library that simplifies the access to most computer’s webcam using python or C++. In just a minimum lines of code, we may access the image currently sampled by the webcam, manipulate it, and show it on the screen.
Open CV is released under a BSD license and hence it’s free for both academic and commercial use. Installers are available to most systems.
Below is an example script that continously streams the image from the webcam to a window on the screen, until “q” is pressed.
import cv2 #Import OpenCV
cap = cv2.VideoCapture(0) #Create an object with access to the webcam
while True:
ret, frame = cap.read() #read the current buffer from the webcam as array "frame"
#im = cv2.cvtColor(frame, cv2.COLOR_BGR2GRAY)/256.0 #Optional: Convert to grayscale and normalize as float (For manipulations)
cv2.imshow('frame',frame) #Display "frame" in a separate window
if cv2.waitKey(1) & 0xFF == ord('q'): #Option for ending the loop
break
cap.release() #Tell webcam that the session is over
cv2.destroyAllWindows() #Destroy windows created by OpenCV
(Note that the above code will not run inside notebook.)
A possible way of using the webcam as boundary conditions is therefore to intercept the loop above after webcam is read in the following manner
- Convert “frame” to grayscale
- Identify boundaries in “frame” as array values below a given threshold
- Impose these boundaries on the solution array $u$ from the solver by demanding that values at indices in $u$ corresponding to boundaries in “frame” is set to 0.
- Advance the solution one timestep.
- Superimpose $u$ and “frame” and show it as an on screen image.
The full implementation
A basic implementation in python may be downloaded here.
The Dirichlet boundary conditions is imposed inside the loop in the “stream_cam” function:
def stream_cam():
#load webcam and iteratively solve the equation
cap = cv2.VideoCapture(0)
Nx = 240
Ny = 320
ret = cap.set(3,Ny) #Scale down the webcam (needed to get a decent framerate)
ret = cap.set(4,Nx)
Lx = 48
Ly = 64
x = np.linspace(0,Lx,Nx) #Setting up the domain
y = np.linspace(0,Ly,Ny)
I = lambda x,y: np.zeros((len(x), len(y))) #np.sin(x[:,None] + y[:,None].T) #Initial condition
q = lambda x,y: np.ones((len(x), len(y))) #np.sin(x[:,None] + y[:,None].T) #Velocity field
V = lambda x,y: np.zeros((len(x),len(y)))#0*x[:,None]*y[:,None].T #Initial condition
f = lambda x,y,t: 10*np.sin(.5*t)*np.exp(-20*((x[:,None]-.5*Lx)**2 + (y[:,None].T-.5*Ly)**2)**2) #Source term
eq = fwave(x,y,1,0,q,I,V,f,0.01) #Initializing the solver
eq.first_step() #Perform the modified first step
while True:
eq.advance() #Advance solution one timestep
ret, frame = cap.read() #Read camera
im = cv2.cvtColor(frame, cv2.COLOR_BGR2GRAY)/256.0 #Normalize to float
#eq.q = im + .1 #Optional; use webcam as variable velocity
eq.u1[im>.4] = 0 #Impose Dirichlet boundary
eq.u2[im>.4] = 0 #Impose Dirichlet boundary
im += eq.u #Superimpose the image and the solution
cv2.imshow('frame',im) #Show superimposed array on screen
if cv2.waitKey(1) & 0xFF == ord('q'):
break
cap.release()
cv2.destroyAllWindows()
Future prospects
Some more interesting possibilities that arise from using the webcam this way has been explored in a separate code.
- Using the webcam array as variable wave velocity
- Using a given color as a parameter in the source term. (For example letting red colors be oscillators).
- Solving other equations/models with webcam boundaries. (Navier-Stokes, Lattice-Boltzmann, diffusion).
Other time-consuming ideas that have not yet been explored:
- Using GPU-acceleration for the calculations. (OpenCV supports this only in C++)
- Implementing Neumann boundaries. (Possibly by numerical gradients or convolution)
- General optimization of the code.
General information
Feel free to use the ideas and code discussed in this document as you please.
Questions, corrections and other topics is thankfully accepted by mail.
