wave_mpi, a C code which solves the 1D wave equation in parallel, using MPI.

This code solves the 1D wave equation of the form:

        Utt = c^2 Uxx
over the spatial interval [X1,X2] and time interval [T1,T2], with initial conditions:
        U(X,T1)  = U_T1(X),
        Ut(X,T1) = UT_T1(X),
and boundary conditions of Dirichlet type:
        U(X1,T) = U_X1(T),
U(X2,T) = U_X2(T).
The value C represents the propagation speed of waves.

The code uses the finite difference method, and marches forward in time, solving for all the values of U at the next time step by using the values known at the previous two time steps.

Central differences may be used to approximate both the time and space derivatives in the original differential equation.

Thus, assuming we have available the approximated values of U at the current and previous times, we may write a discretized version of the wave equation as follows:

        Uxx(X,T) = ( U(X+dX,T   ) - 2 U(X,T) + U(X-dX,T,  ) ) / dX^2
        Utt(X,T) = ( U(X,   T+dT) - 2 U(X,T) + U(X,   T-dt) ) / dT^2
If we multiply the first term by C^2 and solve for the single unknown value U(X,T+dt), we have:
        U(X,T+dT) =        (     C^2 * dT^2 / dX^2 ) * U(X+dX,T)
                    +  2 * ( 1 - C^2 * dT^2 / dX^2 ) * U(X,   T)
                    +      (     C^2 * dT^2 / dX^2 ) * U(X-dX,T)
                    -                                  U(X,   T-dT)
(Equation to advance from time T to time T+dT, except for FIRST step!)

However, on the very first step, we only have the values of U for the initial time, but not for the previous time step. In that case, we use the initial condition information for dUdT which can be approximated by a central difference that involves U(X,T+dT) and U(X,T-dT):

        dU/dT(X,T) = ( U(X,T+dT) - U(X,T-dT) ) / ( 2 * dT )
and so we can estimate U(X,T-dT) as
        U(X,T-dT) = U(X,T+dT) - 2 * dT * dU/dT(X,T)
If we replace the "missing" value of U(X,T-dT) by the known values on the right hand side, we now have U(X,T+dT) on both sides of the equation, so we have to rearrange to get the formula we use for just the first time step:
        U(X,T+dT) =   1/2 * (     C^2 * dT^2 / dX^2 ) * U(X+dX,T)
                    +       ( 1 - C^2 * dT^2 / dX^2 ) * U(X,   T)
                    + 1/2 * (     C^2 * dT^2 / dX^2 ) * U(X-dX,T)
                    +  dT *                         dU/dT(X,   T   )
(Equation to advance from time T to time T+dT for FIRST step.)

It should be clear now that the quantity ALPHA = C * dT / dX will affect the stability of the calculation. If it is greater than 1, then the middle coefficient (1-C^2 dT^2 / dX^2) is negative, and the sum of the magnitudes of the three coefficients becomes unbounded.

We wish to use MPI in order to accelerate this computation. We use the method of domain decomposition - that is, we assume we have P MPI processes, we divide the original interval into P subintervals, and we expect each process to update the data associated with its subinterval.

However, to compute the estimated solution U(X,T+dT) at the next time step requires information about U(X-dX,T) and U(X+dX,T). When process ID tries to make these estimates, it will need one value from process ID-1, and one from process ID+1, before it can make all the updates. MPI allows the processes to communicate this information using messages.

At the end of the complete calculation, we wish to print a table of the solution at the final time. To do this in an organized fashion, we want each process to send its final result to the master process (with ID = 0). Once all the data has been collected, the master process prints it.


The computer code and data files described and made available on this web page are distributed under the MIT license


wave_mpi is available in a C version and a C++ version and a FORTRAN90 version.

Related Data and codes:

COMMUNICATOR_MPI, a C code which creates new communicators involving a subset of initial set of MPI processes in the default communicator MPI_COMM_WORLD.

HEAT_MPI, a C code which solves the 1D Time Dependent Heat Equation using MPI.

HELLO_MPI, a C code which prints out "Hello, world!" using the MPI parallel codeming environment.

LAPLACE_MPI, a C code which solves Laplace's equation on a rectangle, using MPI for parallel execution.

mpi_test, C codes which illustrate the use of the MPI application code interface for carrying out parallel computations in a distributed memory environment.

MULTITASK_MPI, a C code which demonstrates how to "multitask", that is, to execute several unrelated and distinct tasks simultaneously, using MPI for parallel execution.

POISSON_MPI, a C code which computes a solution to the Poisson equation in a rectangle, using the Jacobi iteration to solve the linear system, and MPI to carry out the Jacobi iteration in parallel.

PRIME_MPI, a C code which counts the number of primes between 1 and N, using MPI for parallel execution.

QUAD_MPI, a C code which approximates an integral using a quadrature rule, and carries out the computation in parallel using MPI.

RANDOM_MPI, a C code which demonstrates one way to generate the same sequence of random numbers for both sequential execution and parallel execution under MPI.

RING_MPI, a C code which uses the MPI parallel codeming environment, and measures the time necessary to copy a set of data around a ring of processes.

SATISFY_MPI, a C code which demonstrates, for a particular circuit, an exhaustive search for solutions of the circuit satisfiability problem, using MPI to carry out the calculation in parallel.

SEARCH_MPI, a C code which searches integers between A and B for a value J such that F(J) = C, using MPI.



  1. Geoffrey Fox, Mark Johnson, Gregory Lyzenga, Steve Otto, John Salmon, David Walker,
    Solving problems on concurrent Processors,
    Volume 1: General Techniques and Regular Problems,
    Prentice Hall, 1988,
    ISBN: 0-13-8230226,
    LC: QA76.5.F627.

Source Code:

Last revised on 26 June 2020.