A Conceptual Numerical Model of the Wave Equation Using the Complex Variable Boundary Element Method

Bryce D. Wilkins, Theodore V. Hromadka, Randy Boucher

Department of Mathematics, United States Military Academy, West Point, NY, USA.

DOI: 10.4236/am.2017.85057                

Abstract

In this work, a conceptual numerical solution of the two-dimensional wave partial differential equation (PDE) is developed by coupling the Complex Variable Boundary Element Method (CVBEM) and a generalized Fourier series. The technique described in this work is suitable for modeling initial-boundary value problems governed by the wave equation on a rectangular domain with Dirichlet boundary conditions and an initial condition that is equal on the boundary to the boundary conditions. The new numerical scheme is based on the standard approach of decomposing the global initial-boundary value problem into a steady-state component and a time-dependent component. The steady-state component is governed by the Laplace PDE and is modeled with the CVBEM. The time-dependent component is governed by the wave PDE and is modeled using a generalized Fourier series. The approximate global solution is the sum of the CVBEM and generalized Fourier series approximations. The boundary conditions of the steady-state component are specified as the boundary conditions from the global BVP. The boundary conditions of the time-dependent component are specified to be identically zero. The initial condition of the time-dependent component is calculated as the difference between the global initial condition and the CVBEM approximation of the steady-state solution. Additionally, the generalized Fourier series approximation of the time-dependent component is fitted so as to approximately satisfy the derivative of the initial condition. It is shown that the strong formulation of the wave PDE is satisfied by the superposed approximate solutions of the time-dependent and steady-state components.

Keywords

Complex Variable Boundary Element Method (CVBEM), Partial Differential Equations (PDEs), Numerical Solution Techniques, Laplace Equation, Wave Equation

1. Introduction

In this work, the Complex Variable Boundary Element Method (CVBEM) pro- cedure is extended to modeling applications of the two-dimensional wave partial differential equation (PDE). The wave PDE, given by $c^2\left(u_{xx}+u_{yy}\right)=u_{tt}$ , is com- monly used to model the propagation of waves in highly-stretched membranes that are uniformly dense with approximately constant tension. Consequently, the solution of the wave equation is important in disciplines such as physics, engineering, and other applied sciences that consider the propagation of such waves. The numerical methodology presented in this work uses a coupling of the CVBEM with a generalized Fourier series. The modeling outcome is novel since it is a function that satisfies the strong formulation of the wave PDE within the problem domain. The problem being considered here is

$\begin{array}{ll}\text{Domain:}\hfill & \Omega =\left\{\left(x,y\right)\in {ℝ}^2:0<x<L_1,0<y<L_2\right\}\hfill \\ \text{PDE:}\hfill & c^2\left(u_{xx}+u_{yy}\right)=u_{tt}\hfill \\ \text{Boundary}\text{Conditions:}\hfill & u\left(x\mathrm{,}y\mathrm{,}t\right)=f\left(x\mathrm{,}y\right)\forall \left(x\mathrm{,}y\right)\in \partial \Omega \hfill \\ \text{Initial}\text{Condition:}\hfill & u\left(x\mathrm{,}y\mathrm{,0}\right)=g\left(x\mathrm{,}y\right)\hfill \\ \hfill & u_t\left(x\mathrm{,}y\mathrm{,0}\right)=h\left(x\mathrm{,}y\right)\hfill \end{array}$

The proposed solution technique is based on the standard approach of de- composing the global initial-boundary value problem into two components; namely, a steady-state component and a time-dependent component. The steady- state problem is given as

$\begin{array}{ll}\text{Domain:}\hfill & \Omega =\left\{\left(x,y\right)\in {ℝ}^2:0<x<L_1,0<y<L_2\right\}\hfill \\ \text{PDE:}\hfill & \Delta u_1=0\hfill \\ \text{BoundaryConditions:}\hfill & u_1\left(x\mathrm{,}y\right)=f\left(x\mathrm{,}y\right)\forall \left(x\mathrm{,}y\right)\in \partial \Omega \hfill \end{array}$

The time-dependent problem is given as

$\begin{array}{ll}\text{Domain}:\hfill & \Omega =\left\{\left(x,y\right)\in {ℝ}^2:0<x<L_1,0<y<L_2\right\}\hfill \\ \text{PDE}:\hfill & c^2\Delta u_2=\frac{{\partial }^2{u}_2}{\partial t^2}\hfill \\ \text{Boundary}\text{Conditions}:\hfill & u_2\left(x,y,t\right)=0\forall \left(x,y\right)\in \partial \Omega \hfill \\ \text{Initial}\text{Condition}:\hfill & u_2\left(x,y,0\right)=g\left(x,y\right)-u_1\left(x,y\right)\hfill \\ \hfill & \frac{\partial }{\partial t}{u}_2\left(x,y,0\right)=h\left(x,y\right)\hfill \end{array}$

The steady-state problem is modeled using the CVBEM [1] [2] [3] [4] , which produces an approximation of $u_1$ denoted by ${\hat{u}}_1$ . The time-dependent problem is modeled by a generalized Fourier series, which produces an approximation of $u_2$ denoted by ${\hat{u}}_2$ . The global approximation function is the sum $\hat{u}\left(x,y,t\right)={\hat{u}}_1\left(x,y\right)+{\hat{u}}_2\left(x,y,t\right)$ .

A necessary condition for using this methodology is that the initial condition of the global initial-boundary value problem must be equal on the boundary to the boundary conditions of the global BVP. That is, $g\left(x,y\right)=f\left(x,y\right)$ for all $\left(x,y\right)\in \partial \Omega$ . This is a necessary condition because it ensures that the boundary conditions for the time-dependent problem will be identically zero. Therefore, a two-dimensional Fourier sine series in the spatial variables can be used to provide a particular solution to the time-dependent problem. It is noted that this methodology can still be used when this condition is not met, however, while the global approximation function might closely approximate the initial condition within the problem domain, it will not satisfy the initial condition on $\partial \Omega$ since the time-dependent problem would have non-homogeneous boundary conditions.

The CVBEM approximation function is both continuous and differentiable wherever the CVBEM basis functions are analytic. The approximate solution to the time-dependent component, given later in Equation (2), is also both continuous and differentiable. Therefore, the sum of the two approximate solutions is con- tinuous and differentiable wherever the CVBEM basis functions are analytic.

When the CVBEM basis functions are not analytic throughout the plane, special treatment of the resulting singularities or branch cuts is required in order to ensure that the CVBEM approximation function is at least analytic through- out the problem domain. The recent paper [5] examines several different types of CVBEM basis functions and includes a discussion about the treatment of basis functions that have branch cuts originating from nodes located exterior to the problem domain.

Suppose that the CVBEM basis functions are analytic throughout the problem domain (this can always be done by choosing entire functions, such as complex polynomials, as the CVBEM basis functions). Then, the global approximation function is continuous and differentiable throughout the problem domain. Thus, computational estimates can be developed throughout the problem domain, and vector fields representing the gradient of the global approximation function can be developed throughout the problem domain.

These properties of CVBEM approximation functions are not shared with the usual real-variable domain numerical schemes such as the finite element method (FEM) and finite difference method (FDM), which develop point estimates at nodal points defined by meshing the problem domain and require interpolation to infill approximations of the potential function at points in the domain that are not nodes of the computational mesh.

The differentiability within the problem domain of the coupled CVBEM and generalized Fourier series approximation functions distinguishes this methodology from commonly-used domain methods such as the Finite Difference Method (FDM), Finite Element Method (FEM), and Finite Volume Method (FVM), which only estimate values of the solution at a finite number of nodes defined by a discretization (mesh) of the problem domain.

The purpose of this research is to the extend the well-known CVBEM mo- deling technique to an important hyperbolic PDE. The coupled CVBEM pro- cedure that is presented in this work is significant because it applies some of the numerical advantages of solving Laplace’s equation with a boundary element method to a solution for a PDE that does not have a boundary integral repre- sentation. The result is a numerical technique that satisfies the strong formulation of the wave PDE. This is important because satisfying the strong formulation of a PDE is not generally a principle of the aforementioned domain discretization tech- niques, which only either satisfy the weak formulation of the wave equation or approximately solve the strong formulation of the wave equation.

2. Steady-State Component and the Complex Variable Boundary Element Method

As previously mentioned, the steady-state problem is governed by the Laplace equation. One numerical technique for approximating the solution to these problems is the Complex Variable Boundary Element Method, which has been the topic of numerous papers and books [6] [7] [8] [9] [10] . The CVBEM belongs to the general class of numerical methods known as boundary element methods (BEMs). As is the case with BEMs, the CVBEM relies on the use of a boundary integral equation in its formulation. The boundary integral equation for the CVBEM is the Cauchy integral formula.

Recent research related to the CVBEM has extended its use to modeling problems in three and higher spatial dimensions [11] as well as problems related to the diffusion equation [12] . Current efforts have been dedicated to improving the efficiency of the three-dimensional CVBEM as well as to extending the CVBEM to model applications of the time-dependent three-dimensional diffusion equation.

To solve the steady-state problem, the CVBEM is applied to the boundary conditions of the global BVP. The CVBEM formulation used in this particular work is based upon the use of complex polynomial basis functions, although any other family of analytic functions, including combinations of different families of analytic functions, could be used. It is noted that the choice of the analytic basis functions to be used in the CVBEM formulation may be motivated by the problem boundary geometry in order to simplify the modeling requirements and to better fit the boundary geometry.

The CVBEM approximation function is a linear combination of functions from the family of basis functions used in the CVBEM formulation. The coef- ficients are complex constants, and the basis functions are complex variable functions. Thus, the CVBEM approximation function has both a real and an imaginary part. The coefficients of the CVBEM linear combination are deter- mined by collocating the real part of the CVBEM approximation function with the Dirichlet boundary conditions.

Review of CVBEM Modeling Approach

A CVBEM approximation function, $\hat{\omega }$ , has the general form

$\hat{\omega }\left(z\right)={\displaystyle \underset{k=1}{\overset{p}{\sum }}}{c}_k{g}_k\left(z\right),$ (1)

where $c_k$ is the k^th complex coefficient, $g_k\left(z\right)$ is the k^th member of the family of CVBEM basis functions being used in the approximation, and p is the number of basis functions being used in the approximation.

The complex coefficients $c_k$ , which are composed of two real-valued con- stants, each result in two degrees of freedom to be determined as part of the CVBEM modeling process. Thus, there are 2p degrees of freedom in a CVBEM approximation function, so it is necessary to know at least 2p boundary conditions in order to uniquely determine the coefficients of the CVBEM approximation function when using collocation.

Various techniques for determining the complex coefficients of the CVBEM linear combination include Fredholm integral types of formulations, minimizing the departure between boundary condition values and the CVBEM boundary values in a Hilbert space with least squares, and fitting the CVBEM approximation function exactly to the specified boundary conditions by using collocation [13] . In this paper, collocation is used to determine the coefficient values.

Provided that the basis functions used in the CVBEM formulation are analytic within the problem domain and along the boundary, the result of modeling with the CVBEM is the development of an approximation function that is analytic over the entire problem domain and boundary whose real part solves the Laplace equation and can be fitted to satisfy the prescribed Dirichlet boundary con- ditions.

3. Time-Dependent Problem

The time-dependent component of the problem is modeled by a finite series of the form

$\begin{array}{c}{\hat{u}}_2\left(x,y,t\right)={\displaystyle \underset{m=1}{\overset{M}{\sum }}}{\displaystyle \underset{n=1}{\overset{N}{\sum }}}[a_{m,n}\sin \left(\frac{m\text{π}x}{L_1}\right)\sin \left(\frac{n\text{π}y}{L_2}\right)\cos \left(\text{π}\sqrt{\frac{m^2}{L_1^2}+\frac{n^2}{L_2^2}}t\right)\\ +b_{m,n}\sin \left(\frac{m\text{π}x}{L_1}\right)\sin \left(\frac{n\text{π}y}{L_2}\right)\sin \left(\text{π}\sqrt{\frac{m^2}{L_1^2}+\frac{n^2}{L_2^2}}t\right)]\end{array}$ (2)

Notice that these functions satisfy the strong formulation of the wave PDE. Sine functions are used specifically for the spatial variables in the approximation function so as to satisfy the homogeneous boundary conditions of the time- dependent problem.

In this work, the coefficients $a_{m\mathrm{,}n}$ are determined by collocation of the approximation function with points located in the interior of the problem do- main that are set to match the initial condition of the time-dependent problem. It is necessary to specify mn number of collocation points in order to calculate values for the $a_{m\mathrm{,}n}$ coefficients. These initial condition collocation points are typically located so as to be reasonably uniformly spaced throughout the pro- blem domain (future research will assess advantages of defining collocation point locations non-uniformly). Additionally, the coefficients $b_{m\mathrm{,}n}$ are deter- mined by collocation of the derivative of the approximation function with col- location points located in the problem domain that are set to match the de- rivative of the initial condition of the time-dependent problem. Again, it is necessary to specify mn number of collocation points in order to calculate values for the $b_{m\mathrm{,}n}$ coefficients. These collocation points are also located so as to be reasonably uniformly spaced throughout the problem domain.

4. Method Development

In this section, we do not consider the approximation functions ${\hat{u}}_1$ and ${\hat{u}}_2$ . Rather, we use the exact functions $u_1$ and $u_2$ to demonstrate that the global problem can be decomposed in the manner that we have suggested. Suppose that the initial and boundary conditions are specified as the continuous functions $g\left(x\mathrm{,}y\right)$ and $f\left(x,y\right)$ , respectively. We say that the initial condition is con- sistent with the boundary conditions if $g\left(x,y\right)=f\left(x,y\right)$ for each boundary coordinate $\left(x,y\right)\in \partial \Omega$ . Consequently, if there exists a boundary coordinate pair $\left(x,y\right)\in \partial \Omega$ such that $g\left(x,y\right)\ne f\left(x,y\right)$ , we say that the initial condition is inconsistent with the boundary conditions.

4.1. A Consistent Initial Condition

A function u satisfies the two-dimensional wave equation if $c\left(u_{xx}+u_{yy}\right)=u_{tt}$ , or equivalently, if $c\Delta u=u_{tt}$ . We can decompose the problem into a steady-state component and a time-dependent component. The solution to the steady-state problem is a function $u_1\left(x,y\right)$ such that $\Delta u_1=0$ . That is, the solution to the steady-state problem is a function that satisfies the well-known Laplace partial differential equation as well as the boundary conditions from the global BVP. The solution to the time-dependent component of the problem is a function

$u_2\left(x\mathrm{,}y\mathrm{,}t\right)$ that satisfies (1) the PDE $\Delta u_2=\frac{{\partial }^2{u}_2}{\partial t^2}$ , (2) the homogeneous boundary

conditions, and (3) the initial condition for the time-dependent component calculated as $g\left(x\mathrm{,}y\right)-u_1\left(x\mathrm{,}y\right)$ . The global solution is obtained by adding the solutions to the steady-state and time-dependent problems, $u=u_1+u_2$ .

We shall show that an arbitrary term of Equation (2) satisfies the PDE $\Delta u_2=\frac{{\partial }^2{u}_2}{\partial t^2}$ . This is sufficient to show that the linear combination in Equation (2) also satisfies the wave PDE.

$\begin{array}{c}\frac{{\partial }^2{u}_2}{\partial t^2}=-a_{m,n}{\text{π}}^2\left(\frac{m^2}{L_1^2}+\frac{n^2}{L_2^2}\right)\sin \left(\frac{m\text{π}x}{L_1}\right)\sin \left(\frac{n\text{π}y}{L_2}\right)\cos \left(\text{π}\sqrt{\frac{m^2}{L_1^2}+\frac{n^2}{L_2^2}}t\right)\\ -b_{m,n}{\text{π}}^2\left(\frac{m^2}{L_1^2}+\frac{n^2}{L_2^2}\right)\sin \left(\frac{m\text{π}x}{L_1}\right)\sin \left(\frac{n\text{π}y}{L_2}\right)\sin \left(\text{π}\sqrt{\frac{m^2}{L_1^2}+\frac{n^2}{L_2^2}}t\right)\\ =-a_{m,n}\left(\frac{m^2{\text{π}}^2}{L_1^2}\right)\sin \left(\frac{m\text{π}x}{L_1}\right)\sin \left(\frac{n\text{π}y}{L_2}\right)\cos \left(\text{π}\sqrt{\frac{m^2}{L_1^2}+\frac{n^2}{L_2^2}}t\right)\\ -b_{m,n}\left(\frac{m^2{\text{π}}^2}{L_1^2}\right)\sin \left(\frac{m\text{π}x}{L_1}\right)\sin \left(\frac{n\text{π}y}{L_2}\right)\sin \left(\text{π}\sqrt{\frac{m^2}{L_1^2}+\frac{n^2}{L_2^2}}t\right)\\ -a_{m,n}\left(\frac{n^2{\text{π}}^2}{L_2^2}\right)\sin \left(\frac{m\text{π}x}{L_1}\right)\sin \left(\frac{n\text{π}y}{L_2}\right)\cos \left(\text{π}\sqrt{\frac{m^2}{L_1^2}+\frac{n^2}{L_2^2}}t\right)\\ -b_{m,n}\left(\frac{n^2{\text{π}}^2}{L_2^2}\right)\sin \left(\frac{m\text{π}x}{L_1}\right)\sin \left(\frac{n\text{π}y}{L_2}\right)\sin \left(\text{π}\sqrt{\frac{m^2}{L_1^2}+\frac{n^2}{L_2^2}}t\right)\\ =\Delta u_2\end{array}$ (3)

Since $u_1$ is a function of only x and y, it follows that $\frac{{\partial }^2{u}_1}{\partial t^2}=0$ . Therefore, $u_{tt}=\frac{{\partial }^2{u}_1}{\partial t^2}+\frac{{\partial }^2{u}_2}{\partial t^2}=\frac{{\partial }^2{u}_2}{\partial t^2}$ . To see that u satisfies the governing PDE, observe that

$\begin{array}{c}\Delta u=\Delta \left(u_1+u_2\right)=\Delta u_1+\Delta u_2\\ =0+\frac{{\partial }^2{u}_2}{\partial t^2}=\frac{{\partial }^2{u}_2}{\partial t^2}=u_{tt}\end{array}$ (4)

4.2. An Inconsistent Initial Condition

Now suppose $u\left(x\mathrm{,}y\mathrm{,0}\right)\ne f\left(x\mathrm{,}y\right)$ on $\Gamma$ . When the boundary and initial con- ditions are only specified at a discrete number of points, this case is equivalent to saying that there exists a point $\left(x_i\mathrm{,}{y}_i\right)$ at which $u\left(x_i\mathrm{,}{y}_i\mathrm{,0}\right)$ is specified, but $u\left(x_i\mathrm{,}{y}_i\mathrm{,0}\right)\ne u_1\left(x_i\mathrm{,}{y}_i\right)$ . In this case, the initial condition is not equal on the boundary to the global boundary conditions, so the initial condition is said to be inconsistent. Since the initial condition of the global problem is inconsistent with the boundary conditions, the difference between the initial condition of the global problem and the CVBEM approximation of the steady-state solution would not be approximately zero on the boundary. Thus, use of the two-dimen- sional Fourier sine series, as done in Section 4.1, to approximate the initial condition of the time-dependent component would not be appropriate since the boundary conditions of this problem would not be approximately zero on the boundary. Problems with inconsistent initial conditions can still be modeled using this procedure, however, it is noted that they outcome of this procedure will not satisfy the initial condition on $\partial \Omega$ .

5. Global Solution

The global solution is achieved by summing the solutions to the steady-state and time-dependent problems. Since both of the approximation functions from the two problems are continuous throughout the two-dimensional plane (except for at any branch points or cuts associated with the basis functions used in the Complex Variable Boundary Element Method formulation), their sum, which represents the solution to the global problem is also continuous throughout the plane except for at any point where the CVBEM basis functions are not analytic. Consequently, computational results can be developed continuously throughout the two-dimensional plane except where the CVBEM basis functions have branch points or branch cuts.

5.1. Satisfying the Boundary Conditions

Substituting $c_k={\alpha }_k+i{\beta }_k$ and $g_k\left(z\right)\mathrm{<mo>\; =\; </mo>}{\lambda }_k\left(x,y\right)+i{\mu }_k\left(x,y\right)$ in Equation (1), we can derive a CVBEM approximation function of the form

$\hat{\omega }\left(z\right)={\displaystyle \underset{k=1}{\overset{p}{\sum }}}{c}_k{g}_k\left(z\right)={\displaystyle \underset{k=1}{\overset{p}{\sum }}}\left({\alpha }_k+i{\beta }_k\right)\left({\lambda }_k\left(x,y\right)+i{\mu }_k\left(x,y\right)\right).$ (5)

The real part, ${\hat{u}}_1$ , of the CVBEM approximation function, which represents the potential function, is

$\begin{array}{c}{\hat{u}}_1\left(x,y\right)={\alpha }_1{\lambda }_1\left(x,y\right)-{\beta }_1{\mu }_1\left(x,y\right)+{\alpha }_2{\lambda }_2\left(x,y\right)-{\beta }_2{\mu }_2\left(x,y\right)\\ +\cdots +{\alpha }_p{\lambda }_p\left(x,y\right)-{\beta }_p{\mu }_p\left(x,y\right)\\ ={\displaystyle \underset{k=1}{\overset{p}{\sum }}}\left({\alpha }_k{\lambda }_k\left(x,y\right)-{\beta }_k{\mu }_k\left(x,y\right)\right).\end{array}$ (6)

The boundary conditions for the steady-state problem are Dirichlet and are specified as the same boundary conditions as the global initial-BVP. This implies a system of 2p equations given by

${\hat{u}}_1\left(x_i,y_i\right)={\displaystyle \underset{k=1}{\overset{p}{\sum }}}\left({\alpha }_k{\lambda }_k\left(x_i,y_i\right)-{\beta }_k{\mu }_k\left(x_i,y_i\right)\right)=u_1\left(x_i,y_i\right)$ (7)

for $i=1,2,\cdots ,2p$ . Equivalently, in matrix form

$\left[A_1\right]\left\{c\right\}=\left\{u_1\right\},$ (8)

where $\left\{u_1\right\}$ is a vector containing the specified potential boundary conditions from the global BVP, $\left[A_1\right]$ is the matrix obtained from evaluating the CVBEM basis functions at each of the collocation points, and $\left\{c\right\}$ is a vector containing the unknown coefficients ${\alpha }_1\mathrm{,}{\beta }_1\mathrm{,}\cdots \mathrm{,}{\alpha }_p\mathrm{,}{\beta }_p$ . Once the coefficients are determined by solving the linear system in (8), they can be substituted back into Equation (6), which is the CVBEM approximation of the solution to the steady-state problem. The resulting function can be evaluated throughout the two-dimen- sional plane wherever the CVBEM basis functions are analytic.

In both two-dimensional and higher-dimensional problems, the use of a ran- dom number generator is sometimes employed to locate the collocation points approximately uniformly spatially on the problem boundary. However, they may also be located exactly uniformly spatially on the problem boundary, and the specific problem may dictate which approach is desirable.

Computational issues may arise depending on the choice of basis functions used in the CVBEM approximation function. For example, basis functions in- volving branch cuts, such as complex logarithms or reciprocals of complex monomials, among other types of functions, include considerations of positioning branch cuts to lie outside of the problem domain and boundary in order to enlarge the area of applicability of the CVBEM approximation. Procedures for handling these branch cuts have been well-documented in several papers and books including, [3] and the most recent book [4] and so are not repeated here.

5.2. Satisfying the Initial Condition

Since the exact solution to the steady-state problem is unknown, we need to use the CVBEM approximation of $u_1$ to determine the initial condition for the time-dependent problem. Thus, the initial condition of the time-dependent pro- blem is given by $u_2\left(x,y,0\right)=g\left(x,y\right)-{\hat{u}}_1\left(x,y\right)$ . The coefficients $a_{m\mathrm{,}n}$ are se- lected so as to satisfy

${\hat{u}}_2\left(x_i,y_i,0\right)={\displaystyle \underset{m=1}{\overset{M}{\sum }}}{\displaystyle \underset{n=1}{\overset{N}{\sum }}}{a}_{m,n}\sin \left(\frac{m\text{π}{x}_i}{L_1}\right)\sin \left(\frac{n\text{π}{y}_i}{L_2}\right)=g\left(x_i,y_i\right)-{\hat{u}}_1\left(x_i,y_i\right)$ (9)

at each of the initial condition collocation points. Equivalently, in matrix form

$\left[A_2\right]\left\{a\right\}=\left\{u_2\right\},$ (10)

where $\left\{u_2\right\}$ is a vector containing the specified initial conditions from the time-dependent problem, $\left[A_2\right]$ is the $mn\times mn$ matrix obtained from eva-

luating $sin\left(\frac{m\text{π}x}{L_1}\right)sin\left(\frac{n\text{π}y}{L_2}\right)$ at each of the collocation points, and $\left\{a\right\}$ is a vector containing the unknown coefficients $a_{m\mathrm{,}n}$ .

5.3. Satisfying the Derivative of the Initial Condition

The derivative of the initial condition is specified as $h\left(x\mathrm{,}y\right)$ . At $t=0$ , the time derivative of the finite series approximation of the time-dependent component is

$\frac{\partial }{\partial t}{\hat{u}}_2\left(x,y,0\right)={\displaystyle \underset{m=1}{\overset{M}{\sum }}}{\displaystyle \underset{n=1}{\overset{N}{\sum }}}{b}_{m,n}\left(\text{π}\sqrt{\frac{m^2}{L_1^2}+\frac{n^2}{L_2^2}}\right)\sin \left(\frac{m\text{π}x}{L_1}\right)\sin \left(\frac{n\text{π}y}{L_2}\right).$ (11)

The coefficients $b_{m\mathrm{,}n}$ are selected so as to satisfy

$\frac{\partial }{\partial t}{\hat{u}}_2\left(x_i,y_i,0\right)={\displaystyle \underset{m=1}{\overset{M}{\sum }}}{\displaystyle \underset{n=1}{\overset{N}{\sum }}}{b}_{m,n}\left(\text{π}\sqrt{\frac{m^2}{L_1^2}+\frac{n^2}{L_2^2}}\right)\sin \left(\frac{m\text{π}{x}_i}{L_1}\right)\sin \left(\frac{n\text{π}{y}_i}{L_2}\right)=h\left(x_i,y_i\right).$ (12)

at each of the collocation points for the derivative of the initial condition. Equivalently, in matrix form

$\left[A_3\right]\left\{b\right\}=\left\{h\right\}\mathrm{,}$ (13)

where $\left\{h\right\}$ is a vector containing the value of the derivative of the initial condition at each of the collocation points, $\left[A_3\right]$ is the $mn\times mn$ matrix

obtained from evaluating the expression $\left(\text{π}\sqrt{\frac{m^2}{L_1^2}+\frac{n^2}{L_2^2}}\right)\sin \left(\frac{m\text{π}x}{L_1}\right)\sin \left(\frac{n\text{π}y}{L_2}\right)$ at each of the collocation points, and $\left\{b\right\}$ is a vector containing the unknown coefficients $b_{m\mathrm{,}n}$ .

The coefficients of the CVBEM approximation function are determined so that the CVBEM outcome approximately satisfies the specified Dirichlet boun- dary conditions of the global initial-boundary value problem. Further, the boun- dary conditions for the time-dependent component are identically zero on the boundary due to the use of the two-dimensional Fourier sine series for the spatial variables. Consequently, upon summing the approximate solutions from the two components, the value of the global approximation function on the boundary is exactly the value of the CVBEM outcome on on the boundary, which approximately satisfies the global boundary conditions.

The initial condition of the time-dependent component is calculated as the difference between the initial condition of the global initial-boundary value problem and the CVBEM approximation of the steady-state solution. Therefore, upon superposing the approximate solutions of the steady-state and time- dependent problems, the global initial condition is approximately satisfied. As indicated by the generalized Fourier series in Equation (2), there are two sets of coefficients to be determined; namely, the $a_{m\mathrm{,}n}$ and $b_{m\mathrm{,}n}$ coefficients. When evaluated at $t=0$ , the only terms of the time-dependent approximation func- tion that are nonzero are those corresponding to the $a_{m\mathrm{,}n}$ coefficients. Hence, the $a_{m\mathrm{,}n}$ coefficients are determined so as to approximately satisfy the calcu- lated initial condition of the time-dependent component. Further, when evalu- ated at $t=0$ , the only terms of the time derivative of the time-dependent appro- ximation function that are nonzero are those corresponding to the $b_{m\mathrm{,}n}$ co- efficients. Thus, the $b_{m\mathrm{,}n}$ coefficients are determined so as to approximately satisfy the specified derivative of the initial condition.

6. Limitations

To ensure that the boundary conditions of the global problem are satisfied, the steady-state solution is fitted to the global boundary conditions. Since the global boundary conditions are satisfied exclusively by the steady-state solution, it is necessary to prescribe homogeneous boundary conditions for the time-depen- dent problem. Therefore, only sine functions are used in the Fourier series approximation of the time-dependent problem. By using the Fourier sine series, the solution to the time-dependent problem is zero along the boundary of a rectangular domain of dimensions $L_1\times L_2$ . However, this is a property of the rectangular geometry of the domain and would not necessarily be true for other geometries. Thus, this technique is currently generally limited to problems with rectangular domain geometries or unions of rectangular subdomains. Additionally, this methodology is limited in that it will not satisfy the initial condition on the boundary if $g\left(x\mathrm{,}y\right)\ne h\left(x\mathrm{,}y\right)$ for all $\left(x,y\right)\in \partial \Omega$ since the boundary con- ditions for the transient component would no longer be homogeneous, which means that the two-dimensional Fourier sine series would not satisfy the boundary conditions of the transient problem.

7. Conclusions

The Complex Variable Boundary Element Method approach to solving partial differential equations is well-documented in other publications and the reader is referred to [4] for details regarding the mathematical foundations of the mo- deling approach. In this work, the Complex Variable Boundary Element Method is extended to modeling applications of the two-dimensional wave equations. The technique is based on resolving the global problem into two components, namely steady-state and time-dependent components, that are modeled sepa- rately. The solutions to each component are then summed, which yields a func- tion that satisfies the strong formulation of the wave PDE and approximately satisfies the boundary and initial conditions.

The steady-state component of the problem is solved by application of the CVBEM, and the time-dependent component of the problem is modeled by a truncated series of functions that would be used to solve the PDE analytically. The steady-state solution is fitted so as to satisfy the boundary conditions of the global problem. The boundary conditions of the time-dependent problem are specified to be zero continuously so that when the solutions to the steady-state and time-dependent components are added together, the value of the global solution at each boundary collocation point is as specified by the boundary con- ditions of the global problem. Although the current work focuses on the wave equation, the modeling approach presented applies to other PDE formulations as well as to problems in higher dimensions.

An alternative to the CVBEM is the recently-developed strong-form colloca- tion method known as the Singular Boundary Method (SBM), which uses a linear combination of the fundamental solution of the governing partial dif- ferential equation to approximate the field variables. In that sense, the SBM is a coupling of the boundary element method and the method of fundamental solu- tions. The SBM can easily be applied to problems with complex geometries and can be extended to model problems in three spatial dimensions. The reader is referred to [14] [15] [16] for more details.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Hromadka, T. and Guymon, G. (1984) The Complex Variable Boundary Element Method. International Journal of Numerical Methods Engineering. https://doi.org/10.1007/978-3-642-82361-9
[2]	Hromadka, T. and Lai, C. (1987) The Complex Variable Boundary Element Method. Springer-Verlag, New York.
[3]	Hromadka, T.V. and Whitley, R.J. (1998) Advances in the Complex Variable Boundary Element Method. Springer, New York. https://doi.org/10.1007/978-1-4471-3611-8
[4]	Hromadka, T. and Whitley, R. (2014) Foundations of the Complex Variable Boundary Element Method. Springer, New York. https://doi.org/10.1007/978-3-319-05954-9
[5]	Wilkins, B.D., Hromadka, T., Whitley, R., Johnson, A.N., Boucher, R., Outing, D.A., Phillips, M.D., McInvale, H.D. and Horton, S. (2016) Assessment of Complex Variable Basis Functions in the Approximation of Ideal Fluid Flow Problems. International Journal of Computational Methods and Experimental Measurements.
[6]	Hromadka, T. and Whitley, R. (2005) Approximating Three-Dimensional Steady-State Potential Flow Problems Using Two-Dimensional Complex Polynomials. Engineering Analysis with Boundary Elements, 29, 190-194.
[7]	Dean, T. and Hromadka, T. (2009) A Collocation Cvbem Using Program Mathematica. Engineering Analysis with Boundary Elements, 34, 417-422.
[8]	Hromadka, T. and Whitley, R. (1996) A New Formulation for Developing CVBEM Approximation Functions. Engineering Analysis with Boundary Elements, 18, 39-41.
[9]	Dumir, P. and Kumar, R. (1993) Complex Variable Boundary Element Method for Torsion of Anisotropic Bars. Applied Mathematical Modelling, 17, 80-88.
[10]	Hromadka, T. and Pardoen, G. (1985) Application of the CVBEM to Non-Uniform st. Venant Torsion. Computer Methods in Applied Mechanics and Engineering, 53, 149-161.
[11]	Hromadka, T. (2002) A Multi-Dimensional Complex Variable Boundary Element Method. Topics in Engineering, Vol. 40, WIT Press, Billerica.
[12]	Wilkins, B.D., Greenberg, J., Redmond, B., Baily, A., Flowerday, N., Kratch, A., Hromadka, T.V., Boucher, R., McInvale, H.D. and Horton, S. (2016) An Unsteady Two-Dimensional Complex Variable Boundary Element Method (cvbem). International Journal of Computational Methods and Experimental Measurements.
[13]	Bohannon, A. and Hromadka, T. (2009) The Complex Polynomial Method with a Least-Squares Fit to Boundary Conditions. Engineering Analysis with Boundary Elements, 33, 1100-1102.
[14]	Chen, W. and Wang, F. (2016) Singular Boundary Method Using Time-Dependent Fundamental Solution for Transient Diffusion Problems. Engineering Analysis with Boundary Elements, 68, 115-123.
[15]	Li, J., Fu, Z. and Chen, W. (2016) Numerical Investigation on the Obliquely Incident Water Wave Passing through the Submerged Breakwater by Singular Boundary Method. Computers and Mathematics with Applications, 71, 381-390.
[16]	Qu, W., Chen, W. and Zheng, C. (2017) Diagonal form Fast Multipole Singular Boundary Method Applied to the Solution of High-Frequency Acoustic Radiation and Scattering. International Journal for Numerical Methods in Engineering.