Octave - 2D Wave Equation

In this article I will give a step-by-step guide to implement the two dimensional wave equation in Octave.

Why use Octave?

The reasons for using Octave are:

  • By wrinting the simulation in Octave we reduce the amount of auxiliary code needed to perform the calculations.
  • Octave is optimized for doing large array (matrix) calculations on the CPU. This makes it a good choice for our simulation as we will work with big arrays.
  • We don't have to write a subprogram to visualize the results.
  • Octave is a free and open source software.

Now some people might say: "Why don't you use the GPU to do this calculation?" It is true that the algorithm used in this example would benefit enormously from being run in parallel on the GPU. Nevertheless writing a shader and the program needed to visualize the results is quite a tedious task to do. Additionally, compiled or complex code makes it harder to play around with the parameters of the simulation - which contradicts my main motivation of creating this simulation in the first place.

Installation of Octave

If you are already familiar with using Octave, you can skip this section and directly go to Structure. Octave is available to download in the software repositories of every major Linux distribution. You can install it by typing [Ubuntu/Debian]: sudo apt-get install octave into your terminal, or try this link: [AptUrl] On Windows you can download the installer from the project's download site or directly via their FTP server. After downloading, run the installer and follow the instructions.

This is how an Octave terminal looks like. Octave Window

Octave Syntax

Before we can start we have to get familiar with the syntax used by Octave to make array calculations. We can use

array(x, y);

to address a single value in this array ("array" can be any arbitrary name as long as it doesn't conflict with any built in function). But the big advantage of using Octave over languages like Java or C for this purpose starts when you want to address multiple elements at a time. Normally you would need two nested for loops to change the values of an 2D array, but with Octave you can do this with just one command:

array(from : to, from : to) = newvalue;

So if you want to set the values of an array in the range between 12 and 30 in x-direction and between 25 and 42 in y-direction to zero you just have to write:

array(12:30, 25:42) = 0;

To speed up the program it is recomended to allocate memory before the calculation begins. I normally use the "zeros"-command for this purpose:

array = zeros(sizex, sizey);

Where sizex and sizey is the number of elements in x and y direction. Now we can start to write the actual code.

Navigating through directories in Octave is similar to the way it is done in a terminal or cmd. You can change your working directory by using cd <newdirectory> , where <newdirectory> has to be replaced by the path to the directory you want to go to. For example if you saved your filename.m script in your Documents folder, you use cd /home/<username>/Documents on Linux and Unix machines or cd C:\Users\<username>\Documents on Windows machines, to get there.

To run a a script, cd to the folder where you saved the wave.m file and enter:

source("wave.m")

On Linux you can also run a script by entering

octave wave.m

in a terminal. But you have to keep in mind that this will only work if the script is in your current working directory.

Structure

The program consists of the following parts:

  1. Variable declaration
  2. Setup of initial conditions
  3. Main calculation which runs inside a for loop until a stop time is reached
  4. The visualisation also runs inside this for loop. It is done by plotting the values in 3D using an Octave plot command.

To start this project open your favourite text editor and create a new file called wave.m. The file extension .m indicates an Octave/MATLAB script file.

Setup

The variable declaration is implemented using the following lines of code:

s = 128;              % Array size (spatial resolution of the simulation)
% Create Arrays
p = zeros(s,s);       % Past
c = zeros(s,s);       % Current
f = zeros(s,s);       % Future
speed = 10;           % Propagation speed
dt = 0.0001;          % Timestep (time resolution of the simulation)
dx = 0.01;            % Distance between elements
r = speed * dt / dx; 
stopTime = 10;        % Stop time
n = 300;              % Frequency (of dynamic sources)

This sets up three arrays, which will contain the calculated values from three iterations. As you may have realized the % sign comments out text in the Octave language. In case you want to change the value of the propagation speed or the timestep, it is important to keep the value of r smaller than 0.5. Otherwise the algorithm used to simulate the wave equation will increase the errors produced in every iteration, resulting in a blow up of the numeric values.

Initial conditions

Next we create a for loop to set up the initial conditions:

% Initial Loop
for i = s/2 - s/16 : s/2 + s/16
    % Initial Conditions
    h = s/2 - s/16 : s/2 + s/16;
    c(i,j) = - 2 * cos(0.5 * 2 * pi / (s/8) * j) * cos(0.5 * 2 * pi / (s/8) * i);
    p(i,1:s) = c(i,1:s);
end

This code segment creates initial conditions which are similar to the pattern created shortly after a water drop hits the surface of a lake. The variable h is a vector which helps us circumvent a second for loop.

Main loop

The main implementation of the wave equation starts with the formula we have derived in another article:

$$ f(x,y) = 2c(x,y) - p(x,y)+r^{2} \lbrack c(x-\Delta x,y)+c(x+\Delta x,y)+c(x,y-\Delta y)+c(x,y+\Delta y)-4c(x,y)\rbrack $$

As we can see this equation includes the variables we have defined earlier. Before we can start to translate the equation to Octave code we have to start the main for loop.

% Main Loop
for t = 0:dt:stopTime

To specify a certain increment or decrement of the count variable in every iteration, we can just put value between colons. At this point we are ready to implement the core feature - the 2D wave equation!

%wave equation
f(2:s-1, 2:s-1) = 2 * c(2:s-1, 2:s-1) - p(2:s-1, 2:s-1) + r^2 * ( c(1:s-2, 2:s-1) + c(3:s, 2:s-1) + c(2:s-1, 1:s-2) + c(2:s-1, 3:s) - 4*c(2:s-1, 2:s-1) );

Don't be overwhelmed with the number of values and operators inside the brackets. Octave starts indexing arrays with 1. Thus our first column ( f(2:s-1, 2:s-1) ) addresses all elements in the 2D array called f exept the ones located at the edges. To prevent errors during code execution we have to make sure that the number of elements of the array we want to assign matches the number of elements of the array we are assigning it to. While working with single values we could just use x-1 to address the element "left" to the element at our current position x. Because we are working with arrays we have to shift all indices of the array one element to the desired location. So as an example consider the case where we want to address the elements "left" to the elements c(2:s-1, 2:s-1). To do this we have to change the x indices to look like the following: c(1:s-2, 2:s-1).

After this short digression into indexing we now come back to our simulation. The next step is to update the arrays:

% Update Values
p(1:s, 1:s) = c(1:s, 1:s);
c(1:s, 1:s) = f(1:s, 1:s);

Plotting

At that point the array c contains the calculated waveform of the current timestep. To visualize the results we have to keep an eye on our computer's resources: plotting in three dimensions needs a lot of computing power. To keep the simulation relatively smooth on moderate hardware, we are going to plot every 10th simulation step:

if mod(t/dt, 10) == 0
   % Lines of code to be executed
end

t/dt will return a natural number which represents the number of timesteps since the simulation has started. The modulo of this number and 10 will return the remainder of the division of these two numbers. Consequently comparing the modulo with 0 will only return true every tenth time the count variable t was incremented - which excerts exactly the behavior we were looking for. We are going to visualize the data using a surface plot:

if mod(t/dt, 10) == 0
        b = surf(c);                % Surface Plot
        axis([0 s 0 s -2 2])        % Axis Scale
        caxis([-1 1])               % Color Axis Scale
        pause ( .0001 )
end

The function b = surf(c) creates a 3D surface plot of the values in c and assigns it to the variable b. To get a better view of the results we are going to fix the scale of the axis with the command: axis([0 s 0 s -2 2]) . This sets the scale of the x and the y-axis from 0 to s and it just shows the z values from -2 to 2. Next we fix the color scale to a range between -1 and 1. This is done using caxis([-1 1]). The last step is to add a delay to give Octave some time to draw the plot.

Here you can see the waves produced by a central water drop. Drop 1 Drop 2 Drop Reflection

Transparent boundaries and corners

To create more advanced simulations, we have to introduce new boundary conditions. At the moment our code assigns the value 0 to all boundary elements, which means that every wave which reaches the edge of the simulated area gets reflected inwards. If we continously add a signal (for example the displacement behind a double slit) the maximum displacements of the waves will increase inside our simulation area. When a certain threshold is exceeded, the numbers produces by the simulation start to blow up. One way of avoiding this is to introduce so-called "transparent boundaries". They modify the borders to act as if the simulated area extends indefinitely. The original source of the algorithm used in our example can found here. Unfortunately the authors did not provide the derivation of the boundary conditions. To implement the function we need to add the following lines to our main loop:

f(2:s-1,1) = (2 * c(2:s-1,1) + (r-1) * p(2:s-1,1) + 2*r*r*(c(2:s-1,2  ) + c(3:s,1) + c(1:s-2,1) - 3 * c(2:s-1,1)))/(1+r); % Y:1
f(2:s-1,s) = (2 * c(2:s-1,s) + (r-1) * p(2:s-1,s) + 2*r*r*(c(2:s-1,s-1) + c(3:s,s) + c(1:s-2,s) - 3 * c(2:s-1,s)))/(1+r); % Y:s
f(1,2:s-1) = (2 * c(1,2:s-1) + (r-1) * p(1,2:s-1) + 2*r*r*(c(2,2:s-1  ) + c(1,3:s) + c(1,1:s-2) - 3 * c(1,2:s-1)))/(1+r); % X:1
f(s,2:s-1) = (2 * c(s,2:s-1) + (r-1) * p(s,2:s-1) + 2*r*r*(c(s-1,2:s-1) + c(s,3:s) + c(s,1:s-2) - 3 * c(s,2:s-1)))/(1+r); % Y:s

Each of these lines can handle one edge of our simulated field. As you may have already figured out a "normal" simulation cell in the center of our simulated field has four neighboured cells. As we move to the edge of this field it gets reduced to just three neighbours. But there are still four corners which have only two neighboured cells- so how can we simulate them to function as "transparent corners"?

f(1,1) = ( 2 * c(1,1) + (r-1) * p(1,1) + 2*r*r* ( c(2,1  ) + c(1,2  ) - 2*c(1,1))) / (1+r); % X:1 ; Y:1
f(s,s) = ( 2 * c(s,s) + (r-1) * p(s,s) + 2*r*r* ( c(s-1,s) + c(s,s-1) - 2*c(s,s))) / (1+r); % X:s ; Y:s
f(1,s) = ( 2 * c(1,s) + (r-1) * p(1,s) + 2*r*r* ( c(2,s  ) + c(1,s-1) - 2*c(1,s))) / (1+r); % X:1 ; Y:s
f(s,1) = ( 2 * c(s,1) + (r-1) * p(s,1) + 2*r*r* ( c(s-1,1) + c(s,2  ) - 2*c(s,1))) / (1+r); % X:s ; Y:1

The structure is similar to the one used for the transparent boundaries, exept now with one element removed. To compensate the missing element in the difference calculation we substract two times the current displacement instead of three times. Edge (three adjacent cells):

c(2:s-1,2) + c(3:s,1) + c(1:s-2,1) - 3 * c(2:s-1,1)

Corner (two adjacent cells):

c(2,1) + c(1,2) - 2*c(1,1)

Dynamic sources

The last part which has to be implemented to simulate the double slit experiment is a dynamic source. A dynamic source continuousely changes the values of a specified area in the simulation field. If we want to simulate the double slit experiment we have to add two sources which represent the displacement behind the two slits caused by a continous sinusoidal signal. The change in time can be related with the time which has elapsed since the start of the simulation. The general formulation for a sinusoidal signal is:

y(t) = A * sin(n * t);

A represents the amplitude of the signal, n ist the frequency and t is the time. With that knowledge we can now create two sources in specific locations, chosen to represent the double slit:

f(s/2+s/4-1 : s/2+s/4+1 , 1:2) = 1.5 * sin(t*n);
f(s/2-s/4-1 : s/2-s/4+1 , 1:2) = 1.5 * sin(t*n);

These lines of code have to be inserted in the main loop below the code of the wave equation. When executed two spots on one side of the field will start to periodically move in the vertical direction. Both of them have the same phase, so the appearance of the interference pattern can be modified by changing the frequency of the signal (value of the variable n) as well as the distance between the sources.

This is the wave pattern produced by the double slit experiment: Doubleslit Experiment

Performance

In case the simulation does not run fluently, I recommend changing the value of s - which is the variable to set up the field size (resulting in s^2 cells). So if the value of s is halved the number of elements will be quartered. This results in a tremendous saving of computational ressources. With s=32 I even got the program running on my Raspberry Pi 2!

Download

Here you can download the file. Just unzip it in a directory of your choice and run it with Octave!

Download