MATLAB code to numerically solve a 1D Wave Equation and compare to exact solution

The applications of the Wave Equation are vast, from image processing to system vibrations. This code is directly applied to a perpendicularly-vibrating system linearly defined along the x-axis with some length L (so solving in the x-dimension between bounds 0 and L) for arbitrary time since t = 0s (in the solution, time bounds 0 and 0.5s are chosen). This is a driven system due to an external environmental applied force on the system. The code is automated, so you can change the parameters to observe the solution for a longer time interval or with finer step sizes.

The solution to the 1D Wave Equation presented here is very versatile: it can be applied to a plethora of perpendicularly-vibrating systems, from strings or cables to bridges! Also, the code used to solve the 1D Wave Equation is highly dynamic: the parameters can be set to control the x-dimension length L of the system, the time interval studied in the solution, and even the initial displacement of the system at t = 0s due to external stimuli!

The system is:

• fixed at both ends
• initially at rest
• initially displaced to a general position f(x) = sin(πx)

Note that the initial displacement f(x) of the system can be set to a different function to observe the system response to a different initial stimulus.

The general n-dimension Wave Equation is

In 1 dimension, the 1D Wave Equation is reduced to

Moreover, central difference approximations can be used to represent the second-order partial derivatives:

Dropping the error terms and substituting the difference equations into the 1D Wave Equation, the following computational molecule is obtained, where the Courant Number r = aΔt/Δx

Numerically, given the desired x and t step-sizes (Δx = h and Δt = k respectively), the x-t plane can be partitioned into a mesh, with mesh points of coordinates (ih, jk) for i = {0, 1, ..., xmax/h} and j = {0, 1, ..., tmax/k}. The numerical solution to the 1D Wave Equation is a set of u-values associated with each mesh point which forms the solution function u(x,t). In the mesh array, j gives the row and i gives the column index.

The numerical solution in MATLAB uses the computational molecule applied to 2 subsequent rows (j and j-1 rows) to find the next row (j+1 row). This is an iterative explicit mesh method. Row 0 and the first & last columns of the mesh grid are given by boundary conditions. Also, using the boundary condition that the first derivative of u partial to time is zero, with a backward difference approximation at t=0, a row j=-1 is defined where each mesh point is equal to the mesh point of same i in the j=0 row above. These two rows (-1 and 0) given by the boundary conditions are the first iteration used to find row 1 with the computational molecule. Then rows 0 and 1 are used, etc ... No boundary condition is given for an upper j row, so this iterative method is effective as it can be used for any time length from 0s to arbitrary t.

The analytical solution to the 1D Wave Equation that I calculated with Separation of Variables is:

where

For the given boundary condition where f(x) = sin(πx), bn = 1 for n=1 but bn = 0 for integer n>1 as then bn is the integral of orthogonal functions. Hence, here the exact solution is u(x,t) = sin(πx)cos(πt). This is compared to the numerical solution to show the effectiveness of this numerical method.