Good vibes

• Topic: finite difference method, wave motion

• Task A:

  1. Implement the method to advance time in the 1D wave simulation.

  2. Run the simulation with different time steps and see how it behaves.

  3. Try changing the initial shape of the wave.

• Task B:

  1. Implement the method to advance time in the 2D wave simulation.

  2. Run the simulation and check that you get realistic waves.

• Task C:

  1. Implement the boundary conditions for freely moving bounds in the 1D simulation.

  2. Run the simulation with different boundary conditions and check how the wave behaves now.

• Template:

  • waves.py

  • waves_2d.py

• Further reading:

  • https://en.wikipedia.org/wiki/Wave_equation

  • https://en.wikipedia.org/wiki/Finite_difference_method

waves.py

waves.add_gaussian(string, a=1, center=0.5)

Initialize the wave to a gaussian initial shape.

The function adjusts the shape of the wave by adding a narrow gaussian function to the wave function. Since this operation is additive, it can be called several times to add several pulses.

Parameters:

• string (array) – discretized wave function to be modified

• a (float) – pulse height

• center (float) – fractional position of the peak, should be between 0 and 1

waves.add_sinusoidal(string, n, a=1)

Initialize the wave to a sinusoidal initial shape.

The function adjusts the shape of the wave by adding a sinus function to the wave function. The sinus function will go to zero at the edges of the system. Since this operation is additive, it can be called several times to build the wave as a superposition of sinus functions.

Parameters:

• string (array) – discretized wave function to be modified

• n (int) – number of antinodes

• a (float) – amplitude

waves.add_triangle(string, a=1, peak=0.5)

Initialize the wave to a triangular initial shape.

The function adjusts the shape of the wave by adding a triangular function to the wave function. The triangle will go to zero at the edges of the system. Since this operation is additive, it can be called several times to add several triangles.

Parameters:

• string (array) – discretized wave function to be modified

• a (float) – amplitude

• peak (float) – fractional position of the peak, should be between 0 and 1

waves.animate(wave, L, left_boundary_free=False, right_boundary_free=False)

Animates the simulation.

Parameters:

• wave (list) – list of waves as a time series

• L (float) – length of simulated region \(L\)

waves.draw(frame, wave, x)

Draws the wave.

Used for animation.

Parameters:

• frame (int) – index of the frame to be drawn

• wave (list) – list of waves as a time series

• x (array) – x coordinates

waves.expand_discretization(wave, left_boundary_free=False, right_boundary_free=False)

Adds zeros to the ends of the wave vector, if necessary.

If the system has fixed boundaries, the wave function will always be set to zero there. It is therefore unnecesary to include it in the simulation proper. For these boundaries, the wave is simulated ignoring the zero end points, and this function reinserts the zeros in place after the simulation.

Parameters:

• wave (array) – discretized wave function

• left_boundary_free (bool) – If True, free medium boundary condition is applied at the left bound

• right_boundary_free (bool) – If True, free medium boundary condition is applied at the right bound

waves.first_step(initial_wave, matrix, initial_velocity=0, dt=0)

Performs the first time step.

Since the wave equation is a second degree partial differential equation with respect to time, one needs to know the wave function at two time steps in order to simulate the wave dynamics.

Often one knows instead the initial shape of the wave, \(u(0)\), and the initial velocity of the medium, \(\partial_t u(0)\). Especially if the medium is initially at rest, \(\partial_t u(0) = 0\).

This function calculates the shape of the wave after the first as

\[u(\Delta t) = \frac{1}{2} M u(0) + \partial_t u(0) \Delta t.\]

Parameters:

• initial_wave (array) – discretized wave function in the beginning, \(u(0)\)

• matrix (array) – matrix \(M\) - should be pre-calculated using propagator()

• initial_velocity (array) – velocity of the wave medium in the beginning, \(\partial_t u(0)\)

• dt (float) – time step \(\Delta t\)

Returns:

discretized wave function after the first time step, \(u(\Delta t)\)

Return type:

array

waves.main()

The main program. Simulates a 1D wave.

waves.print_progress(step, total)

Prints a progress bar.

Parameters:

• step (int) – progress counter

• total (int) – counter at completion

waves.propagator(n_nodes, v, dt, dx, left_boundary_free=False, right_boundary_free=False)

Calculates the matrix describing the time evolution of the system.

The discretized wave equation can be written in matrix form as

\[u(t + \Delta t) = M u(t) - u(t - \Delta t).\]

Here \(M\) is a sparse matrix which has the diagonal elements

\[M_{i,i} = 2 - 2 v^2 \left( \frac{\Delta t}{\Delta x} \right)^2\]

and off-diagonals

\[M_{i,i\pm 1} = v^2 \left( \frac{\Delta t}{\Delta x} \right)^2.\]

If the medium is fixed to \(u = 0\) at the boundaries, this is true for all off-diagonals. If the medium is free to move at the boundaries, the first or last row off-diagonals must be adjusted to

\[M_{0,1} = M_{N-1,N-2} = 2 v^2 \left( \frac{\Delta t}{\Delta x} \right)^2.\]

Note

This function is incomplete!

Parameters:

• n_nodes (int) – number of grid points \(N\) (i.e., the number of elements in the vector \(u\))

• v (float) – wave speed

• dt (float) – time step, \(\Delta t\)

• dx (float) – distance between grid points, \(\Delta x\)

• left_boundary_free (bool) – If True, free medium boundary condition is applied at the left bound

• right_boundary_free (bool) – If True, free medium boundary condition is applied at the right bound

Returns:

the matrix \(M\)

Return type:

array

waves.run_simulation(initial_wave, next_wave, matrix, dt, time, recording_dt=0)

Run a dynamic wave simulation.

The discretized wave equation can be written in matrix form as

\[u(t + \Delta t) = M u(t) - u(t - \Delta t).\]

This function iterates this step for the required time.

The function returns a list containing the wave function at different times.

Note

This function is incomplete!

Parameters:

• initial_wave (array) – discretized wave function at the start, \(u(0)\)

• next_wave (array) – discretized wave function after first step, \(u(\Delta t)\)

• matrix (array) – matrix \(M\) - should be pre-calculated using propagator()

• dt (float) – time step \(\Delta t\)

• time (float) – total simulation time

• recording_dt (float) – time interval between saved wavefunctions

Returns:

time evolution of the wave function

Return type:

list

waves_2d.py

waves_2d.add_gaussian(film, a, center_x, center_y)

Initialize the wave to a gaussian initial shape.

The function adjusts the shape of the wave by adding a narrow gaussian function to the wave function. Since this operation is additive, it can be called several times to add several pulses.

Parameters:

• film (array) – discretized wave function to be modified

• a (float) – pulse height

• center_x (float) – fractional position of the peak in x direction, should be between 0 and 1

• center_y (float) – fractional position of the peak in y direction, should be between 0 and 1

waves_2d.animate(wave)

Animates the simulation as a colormap.

Parameters:

wave (list) – list of waves as a time series

waves_2d.animate_3d(wave)

Animates the simulation as a surface.

Parameters:

wave (list) – list of waves as a time series

waves_2d.draw(frame, wave)

Draws the wave as a colormap.

Used for animation.

Parameters:

• frame (int) – index of the frame to be drawn

• wave (list) – list of waves as a time series

waves_2d.draw_3d(frame, wave, xx, yy)

Draws the wave as a surface.

Used for animation.

Parameters:

• frame (int) – index of the frame to be drawn

• wave (list) – list of waves as a time series

• xx (array) – meshgrid with x coordinates

• yy (array) – meshgrid with y coordinates

waves_2d.first_step(initial_wave, matrix, initial_velocity=0, dt=0)

Performs the first time step.

Since the wave equation is a second degree partial differential equation with respect to time, one needs to know the wave function at two time steps in order to simulate the wave dynamics.

Often one knows instead the initial shape of the wave, \(u(0)\), and the initial velocity of the medium, \(\partial_t u(0)\). Especially if the medium is initially at rest, \(\partial_t u(0) = 0\).

This function calculates the shape of the wave after the first as

\[u(\Delta t) = \frac{1}{2} M u(0) + \partial_t u(0) \Delta t.\]

Parameters:

• initial_wave (array) – discretized wave function in the beginning, \(u(0)\)

• matrix (array) – matrix \(M\) - should be pre-calculated using propagator()

• initial_velocity (array) – velocity of the wave medium in the beginning, \(\partial_t u(0)\)

• dt (float) – time step \(\Delta t\)

Returns:

discretized wave function after the first time step, \(u(\Delta t)\)

Return type:

array

waves_2d.main()

The main program. Simulates a 2D wave.

waves_2d.matrix_to_vector(matrix)

Reshapes a square matrix to a vector.

The wave functions are treated as \(N \times N\) matrices e.g. during plotting, but as \(1 \times N^2\) vectors during calculation.

This function reshapes the matrix into a vector.

Parameters:

matrix (array) – the matrix to reshape

waves_2d.print_progress(step, total)

Prints a progress bar.

Parameters:

• step (int) – progress counter

• total (int) – counter at completion

waves_2d.propagator(n, v, dt, dx, top_boundary_free=False, bottom_boundary_free=False, left_boundary_free=False, right_boundary_free=False)

Calculates the matrix describing the time evolution of the system.

The discretized wave equation can be written in matrix form as

\[u(t + \Delta t) = M u(t) - u(t - \Delta t).\]

Here \(M\) is a sparse matrix which is built out of \(N \times N\) blocks on its diagonal.

Each block controls wave dynamics on a horizontal line of grid points in the 2D system. They contain the diagonal elements

\[M_{i,i} = 2 - 4 v^2 \left( \frac{\Delta t}{\Delta x} \right)^2\]

and off-diagonals

\[M_{i,i\pm 1} = v^2 \left( \frac{\Delta t}{\Delta x} \right)^2\]

except for the first and last row of the block where the element outisde the block remains zero.

In addition, adjacent blocks are also coupled by the elements

\[M_{i,i\pm N} = v^2 \left( \frac{\Delta t}{\Delta x} \right)^2.\]

These elements represent the coupling of grid points on a vertical line in the 2D system.

If the medium is free to move at the boundaries, the first or last row off-diagonals of each block must be adjusted to \(2 v^2 \left( \frac{\Delta t}{\Delta x} \right)^2\). This is also true for top and bottom boundaries and the top and bottom off-blocks.

Parameters:

• n_nodes (int) – number of grid points \(N\) (i.e., the number of elements in the vector \(u\))

• v (float) – wave speed

• dt (float) – time step, \(\Delta t\)

• dx (float) – distance between grid points, \(\Delta x\)

• top_boundary_free (bool) – If True, free medium boundary condition is applied at the top bound

• bottom_boundary_free (bool) – If True, free medium boundary condition is applied at the bottom bound

• left_boundary_free (bool) – If True, free medium boundary condition is applied at the left bound

• right_boundary_free (bool) – If True, free medium boundary condition is applied at the right bound

Returns:

the matrix \(M\)

Return type:

array

waves_2d.run_simulation(initial_wave, next_wave, matrix, dt, time, recording_dt=0)

Run a dynamic wave simulation.

The discretized wave equation can be written in matrix form as

\[u(t + \Delta t) = M u(t) - u(t - \Delta t).\]

This function iterates this step for the required time.

The function returns a list containing the wave function at different times.

Note

This function is incomplete!

Parameters:

• initial_wave (array) – discretized wave function at the start, \(u(0)\)

• next_wave (array) – discretized wave function after first step, \(u(\Delta t)\)

• matrix (array) – matrix \(M\) - should be pre-calculated using propagator()

• dt (float) – time step \(\Delta t\)

• time (float) – total simulation time

• recording_dt (float) – time interval between saved wavefunctions

Returns:

time evolution of the wave function

Return type:

list

waves_2d.vector_to_matrix(vector)

Reshapes a vector to a square matrix.

The wave functions are treated as \(N \times N\) matrices e.g. during plotting, but as \(1 \times N^2\) vectors during calculation.

This function reshapes the vector into a matrix.

Note: the length of the vector must be a square number.

Parameters:

vector (array) – the vector to reshape