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
The Wave equation in 1D reads
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.
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 wave equation is second order partial derivative because of this two
initial conditions
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.
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
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.
A simulation using the values
resulted in the wave shown below.
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
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 Uk has the dimension i x j, TL has dimension i x i and TR has dimension j x j. Their individual structures are given by
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:
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
A simulation using the values
resulted in the wave shown below.
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.
- 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