VR Musical Instruments


This is a 2D wave simulation based on a differential equation describing classical waves in two dimensions. The simulation is based on approximations on the current and previous states thus allowing real-time interaction with adding waves by clicking and adding walls for the waves to reflect off of. The application is creating in C++ using Raylib. The web online demo is compiled using Emscripten. The desktop application uses a thread pool to take full advantage of the hardware to allow for faster calculations. You can continue reading to understand the math behind the simulation.

The Wave Equation

This is the two-dimensional wave equation which describes wave propagation. 2ut2=c2(2ux2+2uy2)\frac{\part^2 u}{\part t^2} = c^2 \left( \frac{\part^2 u}{\part x^2} + \frac{\part^2 u}{\part y^2} \right)
  • uu represents the wave function u(x,y,t)u(x,y,t). It can represent various quantities, such as pressure in a sound wave or the height of water waves given a particular position and time.
  • tt is time.
  • xx and yy are the spatial coordinates.
  • cc is the speed of wave propagation.

Starting on the left 2ut2\frac{\part^2 u}{\part t^2}, this is the second time derivative of the wave function which represents the acceleration of the wave in time. The right side of the equation 2ux2+2uy2\frac{\part^2 u}{\part x^2} + \frac{\part^2 u}{\part y^2} is the second spatial derivative of the wave function for each dimension. This represents the "curvature" of the wave in space. This whole equation thus shows that the acceleration of the wave in time is governed by the curvature of the wave in space multiplied by some constant wave speed cc.

Time Derivative

Since we are dealing with discrete values for both space (a grid) and time (individual frames/time steps), we need a way to approximate these second derivatives. Essentially, given a current grid state and previous states, we want to determine what the next state/frame should be in time. For this, we will use the finite difference method which allows us to approximate derivatives given a certain number of states in either time or space.

The second order finite difference formula for one dimension is f(x)f(x+h)2f(x)+f(xh)h2 f''(x) \approx \frac{f(x+h)-2f(x)+f(x-h)}{h^2}
  • xx is the position in space.
  • hh is the distance between discrete points in space.

Since we want to find the derivative of time and not space, we can adjust this formula. Instead of looking at discrete distances left and right of a position, we can look forward and backward in time. Now hh represents the delta timestep between states rather than between spatial grid points. To represent this better, we can say f(x+h)f(x+h) is the state in the future ffuture(x)f_\text{future}(x), f(x)f(x) is the current state fpresent(x)f_\text{present}(x), and f(xh)f(x-h) is the state in the past fpast(x)f_\text{past}(x). We can also say hh is the delta time between states or Δt\Delta t.

Now we get f(x)ffuture(x)2fpresent(x)+fpast(x)Δt2f''(x) \approx \frac{f_\text{future}(x)-2f_\text{present}(x)+f_\text{past}(x)}{\Delta t^2}
  • xx is the position in space.
  • hh is the distance between discrete points in space.

However, this assumes we have one spatial dimension xx, so we can substitute xx in the functions for (x,y)(x,y). f(x,y)ffuture(x,y)2fpresent(x,y)+fpast(x,y)Δt2f''(x,y) \approx \frac{f_\text{future}(x,y)-2f_\text{present}(x,y)+f_\text{past}(x,y)}{\Delta t^2}
Our goal is to determine the state of the function in the future at a specific (x,y)(x,y) coordinate. This means we need to solve for ffuture(x,y)f_\text{future}(x,y) ffuture(x,y)=f(x,y)Δt2fpast(x,y)+2fpresent(x,y)f_\text{future}(x,y) = f''(x,y)\Delta t^2 - f_\text{past}(x,y) + 2f_\text{present}(x,y)
Going back to our original wave propagation equation, 2ut2\frac{\part^2 u}{\part t^2} is represented as f(x,y)f''(x,y) in our equation. They are the same, just difference notation and also a partial derivative because the wave function u(x,y,t)u(x,y,t) is a function of both space and time and thus taking the derivative of time makes it a partial derivative instead of a complete derivative. ffuture(x,y)=2ut2Δt2fpast(x,y)+2fpresent(x,y)f_\text{future}(x,y) = \frac{\part^2 u}{\part t^2}\Delta t^2 - f_\text{past}(x,y) + 2f_\text{present}(x,y)
In our simulation, we can store the present fpresent(x,y)f_\text{present}(x,y) and past fpast(x,y)f_\text{past}(x,y) states of our simulation in separate buffers and also keep track of the timestep Δt\Delta t, however, in order to fully calculate the future state, we need the second partial derivative of the wave function with respect to time 2ut2\frac{\part^2 u}{\part t^2}. Using the original wave propagation function, we can substitute the right side of the equation in. ffuture(x,y)=c2(2ux2+2uy2)Δt2fpast(x,y)+2fpresent(x,y)f_\text{future}(x,y) = c^2\left( \frac{\part^2 u}{\part x^2} + \frac{\part^2 u}{\part y^2} \right)\Delta t^2 - f_\text{past}(x,y) + 2f_\text{present}(x,y)

Spatial Derivative

We now need to approximate the second spatial derivative of the wave function 2ux2+2uy2\frac{\part^2 u}{\part x^2} + \frac{\part^2 u}{\part y^2}. We can use the same method we used for time; the finite difference method. However, unlike time where there is a single dimension, we have two spatial dimensions thus requiring us to use the two-dimensional second-order finite difference formula. f(x,y)f(xh,y)+f(x+h,y)4f(x,y)+f(x,yh)+f(x,y+h)Δg2f''(x,y) \approx \frac{f(x-h,y)+f(x+h,y)-4f(x,y)+f(x,y-h)+f(x,y+h)}{\Delta g^2}
  • xx and yy are coordinates in space.
  • Δg\Delta g is the distance between grid points.

Even though this equation looks pretty complicated, it can actually be visualized nicely through a stencil which can be imagined as an overlay applied on top of a grid where the number is a multiplier.
stencil grid with multipliers

We can substitute this into our equation for predicting the future. ffuture(x,y)=c2fpresent(x,y)Δt2fpast(x,y)+2fpresent(x,y)f_\text{future}(x,y) = c^2f''_\text{present}(x,y)\Delta t^2 - f_\text{past}(x,y) + 2f_\text{present}(x,y)
Now we have everything necessary to solve for ffuture(x,y)f_\text{future}(x,y). We simply need to provide neighboring grid cells surrounding (x,y)(x,y), a past value of (x,y)(x,y) and some constants such as wave speed cc, timestep Δt\Delta t, and grid spacing Δg\Delta g.

Embed Demo

<embed width="100%" height="100%" src="https://static.pixeled.site/wave-simulation/index.html"></embed>

© 2023 Matthew Oros