Wave equation simulator

The Python code in this repository solves a 1D wave equation on a grid. Additionally, a 2D wave equation simulator was added. The topics which will be discussed are:

• Discretization of the wave equation in 1D and 2D
• Choosing and implementation of a boundary condition
• Implementation of initial conditions
• Stability analysis
• Results
• Benchmark comparing 1D and 2D simulation
• Sources

Wave equation simulator 1D

Discretizing the wave equation

The Wave equation in 1D reads

In the next step, this partial differential equation is discretized by using a finite difference method

The above equation is a central finite difference. Using this expression simplifies the wave equation to

Now this does not look particularly simple. But setting fully discretizing space and time results in the following equation for updating the equation from one time step to another

The above equation is visualized in the image below, where the amplitude of the i-th point at the (j+1)-th time step is given the point at the time step and its neighbors. Further it relies on the state of the amplitudes at the (j-1)-ths time step.

The state of N amplitudes at time step j can be imagined as a vector

This transforms the update equation to a linear equation involving an update matrix T, the j-th state and the (j-1)-th state

The frame of zeros in the matrix T correspond to the boundary condition of fixed end points which was already quietly implemented without mentioning. Further only the diagonal and the off-diagonal elements of T are populated. The discussion regarding the boundary condition will follow now.

Boundary condition

As stated in the last chapter the chosen boundary condition are fixed endpoints for all time steps. For a 1D grid consisting of N amplitudes and T time steps this corresponds to

The frame of zeros in the matrix T makes sure that the boundary condition is always fulfilled.

Initial conditions

The wave equation is second order partial derivative because of this two initial conditions

These equations have to be discretized now. The first condition is fulfilled by setting

This result can now be inserted in the basic update equation to calculate the first time step. However, doing this causes a problem since the it involves the amplitudes at time step -1 which are not defined at this point. But, the amplitudes at this time step can be calculated using the second initial condition. To do this the derivative is discretized via a central finite difference

This equation can now be rearranged to get an expression for the amplitudes at the -1 time step

Using the derived expressions results in the following expression for the amplitude of the first time step

Using the matrix T results in a rather handy expression for the amplitudes of the first time step

This equation reflects the impact of the initial conditions on the solution of the wave equation.

Stability analysis

In this section the stability of the discrete scheme will be examined using the Von-Neumann-stability formalism. This theory examines the stability of a discrete scheme by looking at the propagation of errors. The total error of the i-th amplitude at the j-th time step is decomposed into a Fourier series , where A is the so called amplification factor. The amplification factor is assumed to behave according to a power law in time

A scheme is unstable if the modulus of the amplification factor is greater than 1, |A| > 1. Further, the scheme diverges if the amplification factor of one Fourier component is greater than 1. The errors propagate in the scheme via the update equation. Plugging one Fourier component into the update equation from the j-th to the (j+1)-th time step results in

The possible solutions for A read

For a better understanding of this section let me introduce the Courant number. This is the central variable defining the stability of a scheme, which has to fulfill the condition

As long as this condition is fulfilled the scheme is guaranteed to be stable.

Results

A simulation using the values

and the initial conditions

resulted in the wave shown below.

Wave equation simulator 2D

Discretization of the 2D wave equation

Deriving the discretized wave equation in 2D involves exactly the same procedure as its 1D counterpart. Because of this only the resulting equation which describes the transition of a grid point at the j-th position in the i-th row from the k-th to the k+1-th time step is stated

Note that the grid spacings in x- and y-direction are the same. Otherwise there would be two different courant numbers. For the first time step the initial conditions are now 2D and the following time step is given by the equation

Both of these equations can be interpreted as matrix equations where the amplitudes on the grid are represented by the entries of a matrix U. The equation allowing a transition from the k-th to the k+1-th time step is given by

F and G are matrices given by the initial conditions evaluated on the grid points. Both matrices in the above equations have the structure as the matrix in the 1D case. The only differences to the 1D case are their dimensions. If the matrix U_k has the dimension i x j, T_L has dimension i x i and T_R has dimension j x j. Their individual structures are given by

Boundary conditions

In the 2D simulation three different boundary conditions are implemented: fixed edges, loose edges and a cyclical boundary condition. These boundary conditions are implemented via the time step matrices:

• Fixed edges:
  
  

  

• Loose edges:
  
  

  

• Cyclical:
  
  

  

Stability analysis

The calculation of the stability condition is largely the same as in 1D. The only difference is regarding the error at the l-th grid point in x, the m-th grid point in y-direction and the n-th time step, which is now given by a 2D Fourier series. Again it suffices to only observe one Fourier component

Plugging this into the update equation and simplifying the result leads to the stability condition

Results

A simulation using the values

and the initial conditions

resulted in the wave shown below.

Benchmark

This is a bench mark for different number of time steps and different dimensions of the problem for random initial amplitudes. The initial velocity of the amplitudes are set to zero resulting in the following bench mark.

It appears that the 2D simulation is rather slow. If I have time I will maybe speed it up in the future using Cython or something like this.

Sources

• Numpy
• Matplotlib
• Discretization of the 1D wave equation and stability analysis, accessed: 21.02.21, 22:45 MEZ
• Stability analysis of the 2D wave equation, accessed: 27.03.21, 16:43 MEZ