Numerical Study on Single and Multi-Dimensional Boundary Value Problems by the Method of Weighted Residual

Humaira Farzana^1*, Md. Abdul Alim²

¹Department of Arts & Sciences, Ahsanullah University of Science & Technology, Dhaka, Bangladesh.
²Department of Mathematics, Bangladesh University of Engineering & Technology, Dhaka, Bangladesh.

DOI: 10.4236/ajcm.2020.102012                

Abstract

The Galerkin method of weighted residual (MWR) for computing natural frequencies of some physical problems such as the Helmholtz equation and some second-order boundary value problems has been revised in this article. The use of one and two-dimensional characteristic Bernstein polynomials as the basis functions have been presented by the Galerkin MWR. The vibration of non-homogeneous membranes with Dirichlet type boundary conditions is also studied here. The useful properties of Bernstein polynomials, its derivatives and function approximations have also been illustrated. Besides, the efficiency and applicability of the proposed technique have been demonstrated through some numerical experiments.

Keywords

Eigenvalues, Heat Equation, Wave Equation, Bernstein Polynomials, Galerkin MWR

1. Introduction

Mathematical model derived from multi-dimensional differential equations play a crucial role in modeling a variety of scientific and engineering application problems. The study of vibration of membrane of different geometries has been the subject of interest to scientists and engineers since long back owing to their wide applications in every field of modern technology. Different shapes or geometries of membranes are the fundamental constituents in many scientific states. Many researchers concentrated on membrane structures and studied about the vibration theory of membrane employing numerous analytical or numerical techniques. However, to the best of our knowledge, none had attempted to provide approximate eigenvalues of boundary value problems employing Bernstein-Galerkin method of weighted residual (MWR).

Some physical problems arise in science and engineering are modeled by the elliptic, parabolic and hyperbolic partial differential equations. The process of solving certain linear evolution equations such as the heat or wave equations have received increased interest from several fields of science and engineering, either on bounded domains. Certain class of linear evolution equations for instance the heat or wave equations have a second order linear differential eigenvalue problem with two boundary conditions but no unique solution exists. A code for second order problems has been implemented and special attention has been paid to the approximation of the boundary conditions in the two-dimensional case.

A new numerical method based on Bernstein polynomials expansion is proposed for solving one dimensional (1-D) elliptic interface problems [1]. At present, formulations exploiting Galerkin method of weighted residual are constructed to determine the coefficient matrix derived from linear system of equations. Helmholtz equation appears in vibrating-string problems and in finding the temperature distribution in a rod. Calculations of eigenvalues of Helmholtz equation using boundary method are presented in [2]. Numerical technique utilizing Bernstein polynomials basis to give the approximate solution of a parabolic partial differential equation illustrated in [3]. The large deflection theory of membrane is applied [4] to derive the governing vibration equations of orthotropic membrane. The Bernstein polynomials are used to find numerical solution of elliptic boundary value problems with a point collocation method [5]. The authors solved the problem and obtained the power series formula of nonlinear vibration frequency of rectangular membrane with four edges fixed. Bernstein Ritz-Galerkin method for solving an initial boundary value problem for two-dimensional (2-D) wave equation has been studied by [6]. A spectral element model for the transverse vibration of a membrane has been developed by [7] using the boundary splitting method and the waveguide FEM-based spectral super element method (SSEM). Exact solutions are well known for the free vibration of a membrane with simple geometry such as elliptical [8], circular [9], general convex polygonal membranes [10], rectangular membranes [11] - [17]. Nevertheless, closed-form solutions rarely exist which leads many researchers to have numerical solution applying several techniques in view to minimize the error with less computational difficulties. Free vibrations of membranes have been studied in the literature by applying the method of superposition using wave-type base functions [10], the optimized Galerkin-Kantorovich method [11], classical Rayleigh-Ritz method [12], finite difference method [13], the Hybrid method [14], the discrete singular advance convolution method [15], the Kantorovich method [16]. Furthermore, the external and internal excitation methods [17], are also investigated. Two-dimensional (2-D) Helmholtz equation is solved for non-homogeneous membrane exploiting Sinc collocation method [18].

Double ultraspherical spectral Galerkin method with Chebychev polynomials is applied to solve elliptic partial differential equations (PDEs) [19]. Chebychev spectral quasi-inverse matrix diagonalization Galerkin method [20] is utilized to solve multidimensional PDE’s. Finite difference methods have been offered for the numerical solution of the one-dimensional parabolic equation in the study [21]. Eigenvalue analysis for second order boundary value problems arises in several engineering application areas has been investigated using analytical or numerical techniques [22]. Bernstein basis has been exploited to solve ordinary and partial differential to solve eigenvalues as well as boundary value problems employing numerous techniques are in [23] - [28].

In this study, we have presented Bernstein polynomials based Galerkin method of weighted residual (MWR) technique that offers accurate solutions, are put up with in terms of truncated series of smooth polynomial functions. The key advantages of applying the Galerkin MWR are its superior accuracy which can be achieved even with a few numbers of polynomials compared to the other techniques reported in the proposed work.

The significant benefit of the present method is that the trial Bernstein basis function can be modified effortlessly to any desired form so as to satisfy the Dirichlet and derivative boundary conditions for both one and two-dimensional problems.

Our goal is to construct an appropriate method with Dirichlet and derivative type boundary conditions and to reduce the errors as well as the computational difficulties. The main reason why the Galerkin MWR is chosen is its flexibility and simple implementation. The notable properties of Bernstein polynomials vanishing at its endpoints over the finite interval inspire us to execute it in the Galerkin technique. Furthermore, the smallest eigenvalue which characterizes potentially the most visual structures of the dynamical systems can be computed very accurately applying the said technique which converges to large significant digits and fairly close to the analytical results.

Some analytical methods involve a reduced amount of memory and arithmetic, but usually they are unable to fix the maximal dimension for the spectrum. Perturbation methods limit the extent of its application due to disadvantages in its theory in its existing form and all types of techniques required small parameter and consequently require some skill to implement. Finite difference methods (FDM) and finite elements methods (FEM) are well-known discretization techniques that comprise extensive calculation and huge memory storage in machine for storing matrices having larger dimensions. Additionally, the precision rapidly worsens for the upper eigenvalues. Furthermore, FDM gives value at particular points ineffective to evaluate the values at the desired points between two grid points and takes more computational cost for getting higher accuracy. The co-efficient matrix in Galerkin MWR is sparse and has symmetric banded matrix, which minimizes the computational effort and attains relatively smaller condition numbers. Although discretization methods attain high accuracy, in the end produce ill-conditioned system. On the other hand, the obtained matrices are full and non-symmetric which intensify the computational cost.

In section 2, we provide an outline of Bernstein polynomial basis in single and multi-dimensions. Approximation of functions to Bernstein basis functions and their derivatives required for numerical computations are also illustrated. Section 3 is offered to the construction of 2-D Helmholtz eigenvalue problems and formulation of one dimensional transient heat equation. In section 4, we give short description of hyperbolic equation, analytical solution of 2-D wave equation and Galerkin formulation of inhomogeneous rectangular membrane. Convergence of Galerkin method and error estimation is discussed in section 5. In section 6, we represent some numerical examples eigenvalue problems in one and two dimensions. In section 7, we illustrate some results with detailed discussions. Moreover, we validate our code by the offered method and eigenvalues are compared to the analytical ones. We have drawn conclusions of this study in section 8.

2. Bernstein Polynomials

2.1. Bernstein Polynomials in One Dimension [25] [26]

$B_{i,n}\left(x\right)=\left(\begin{array}{l}n\\ i\end{array}\right)x^i{\left(1-x\right)}^{n-i}$, $0\le i\le n$, (1)

where

$\left(\begin{array}{l}n\\ i\end{array}\right)=\frac{n!}{i!\left(n-i\right)!};n=0,1,2,\cdots$

For each positive integer n, there is a sequence of Bernstein polynomials over the finite interval $\left[0,1\right]$. Well convergence of Bernstein polynomials for some class of functions is found in [29]. Bernstein polynomials which approximates a functions of one real variable $f:\left[0,1\right]\to R$ is defined as

$B_{i,n}\left(f;x\right)={\displaystyle \underset{i=\text{0}}{\overset{n}{\sum }}f\left(\frac{i}{n}\right)}\left(\begin{array}{c}n\\ i\end{array}\right)x^i{\left(1-x\right)}^{n-i}$. (1a)

2.2. Two Dimensional Bernstein Polynomials

2-D Bernstein polynomials approximates a functions of two real variables $f:\left(\left[0,1\right]\times \left[0,1\right]\right)\to R$ is defined as

$B_{i,k}^{\left(n,m\right)}\left(f\left(x,y\right)\right)={\displaystyle \underset{i=0}{\overset{n}{\sum }}{\displaystyle \underset{k=0}{\overset{m}{\sum }}f}\left(\frac{i}{n}\right)}\left(\frac{j}{m}\right)B_{i,k}^{\left(n,m\right)}\left(x,y\right)$ (2)

where

$B_{i,k}^{\left(n,m\right)}\left(x,y\right)=\left(\begin{array}{c}n\\ i\end{array}\right)\left(\begin{array}{c}m\\ k\end{array}\right)x^i{\left(1-x\right)}^{n-i}{y}^k{\left(1-y\right)}^{m-k}$ (2a)

where $\left(\begin{array}{c}n\\ i\end{array}\right)=\frac{n!}{i!\left(n-i\right)!}$, $\left(\begin{array}{c}m\\ k\end{array}\right)=\frac{m!}{k!\left(m-k\right)!}$, $i=0,1,2,\cdots ,n$ ; $k=0,1,2,\cdots ,m$ with $\left(x,y\right)\in \left[0,1\right]$.

There are $\left(n+1\right)$ and $\left(m+1\right)$ $\left(n+m\right)$ -th degree Bernstein polynomials. For suitability, we usually set $B_{i,k}^{\left(n,m\right)}=0$, if $i<0$, $k<0$ or $i>n$, $k>m$.

Some interesting and useful properties of 2-D Bernstein polynomials are demonstrated in [30].

2.3. The Bernstein Approximation Theorem 1 [27]

Every continuous function f defined on [0, 1] can be uniformly approximated as closely as desired by a polynomial function. For any $\epsilon >0$, there exists a positive integer N such that for all $x\in \left[0,1\right]$, an integer $n\ge N$ we have

$\left|f\left(x\right)-B_n\left(f;x\right)\right|<\epsilon$. (3)

Hence given any power-form polynomial of degree N, it can be uniquely converted into a Bernstein polynomial of degree n for $n\ge N$.

Let $C:\left(\left[0,1\right]\times \left[0,1\right]\right)\to R$ is a continuous function space and is a linear norm space on $\left[0,1\right]\times \left[0,1\right]$ with the norm ${\Vert f\Vert }_{C\left(\left[0,1\right]\times \left[0,1\right]\right)}=\underset{\left(x,y\right)\in \left[0,1\right]\times \left[0,1\right]}{\max }\left|f\left(x,y\right)\right|$ and all partial derivatives of $f\left(x,y\right)$ exists, while $\left\{f_{n,m}\right\}$ is a sequence on the space $C\left(\left[0,1\right]\times \left[0,1\right]\right)$, $f\in C\left[0,1\right]\times \left[0,1\right]$.

${\Vert f_{n,m}-f\Vert }_{C\left(\left[0,1\right]\times \left[0,1\right]\right)}=0$.

which implies the uniform convergence to the function $f$. Let

$H=L^2\left\{\left[0,1\right]\times \left[0,1\right]\right\}$ and $\left\{B_{0,n}\left(x\right)B_{0,n}\left(y\right),B_{1,n}\left(x\right)B_{1,n}\left(y\right),B_{2,n}\left(x\right)B_{2,n}\left(y\right),\cdots ,B_{n,n}\left(x\right)B_{n,n}\left(y\right)\right\}\subset H$, $f\subset H$.

Let $V=\text{Span}\left\{B_{i,n}\left(x\right)B_{j,n}\left(y\right)\right\}$ ; $i,j=0,1,2,\cdots ,n$. Since V is a finite dimensional vector space, f has a unique best approximation from V say $v_0\in V$ $\exists v\in V$ such that ${\Vert f-v_0\Vert }_2\le {\Vert f-v\Vert }_2$.

$f\approx v_0={\displaystyle \underset{i=0}{\overset{n}{\sum }}{\displaystyle \underset{j=0}{\overset{n}{\sum }}{f}_{i,j}{B}_{i,n}\left(x\right)B_{j,n}\left(y\right)}}=C^{\text{T}}B\left(x,y\right)$. (4)

Lemma: Suppose that $f^{n+1}:\left[0,1\right]\times \left[0,1\right]\to R$ is a continuous differentiable function and all partial derivatives of $f\left(x,y\right)$ exist and are continuous. If $C^{\text{T}}B\left(x,y\right)$ is the best approximation for f out of V then the error bound is represented as

${\Vert f-C^{\text{T}}B\Vert }_2\le \frac{2l{\left(k_1+k_2\right)}^{n+2}}{\left(n+1\right)\text{!}\sqrt{\left(\text{2}n+\text{3}\right)\left(2n+4\right)}}$. (5)

as $k_1=\max \left\{1-x_0,x_0\right\}$, $k_2=\max \left\{1-y_0,y_0\right\}$ and the partial derivatives of f are bounded in magnitude by $l$. By applying the Taylor’s expansion in two variables we define

$\begin{array}{l}f\left(x,y\right)=f\left(x_0,y_0\right)+\left(x-x_0\right){\frac{\partial f}{\partial x}|}_{\left(x_0,y_0\right)}+\left(y-y_0\right){\frac{\partial f}{\partial y}|}_{\left(x_0,y_0\right)}\\ +\frac{1}{2!}[{\left(x-x_0\right)}^2{\frac{{\partial }^{2}f}{\partial x^2}|}_{\left(x_0,y_0\right)}+2\left(x-x_0\right)\left(y-y_0\right){\frac{{\partial }^{2}f}{\partial x\partial y}|}_{\left(x_0,y_0\right)}\\ +{\left(y-y_0\right)}^2{\frac{{\partial }^{2}f}{\partial y^2}|}_{\left(x_0,y_0\right)}]+\cdots +\frac{1}{n!}{\left[\left(x-x_0\right)\frac{\partial }{\partial x}+\left(y-y_0\right)\frac{\partial }{\partial y}\right]}^{n}f\left(x_0,y_0\right)\\ +\frac{1}{\left(n+1\right)!}{\left[\left(x-x_0\right)\frac{\partial }{\partial x}+\left(y-y_0\right)\frac{\partial }{\partial y}\right]}^{n+1}\end{array}$

$\begin{array}{l}\cdot f\left(x_0+\theta \left(x-x_0\right)y_0+\theta \left(y-y_0\right)\right),0<\theta <1\\ =f_1\left(x,y\right)+\frac{1}{\left(n+1\right)!}{\left[\left(x-x_0\right)\frac{\partial }{\partial x}+\left(y-y_0\right)\frac{\partial }{\partial y}\right]}^{n+1}\\ \cdot f\left(x_0+\theta \left(x-x_0\right),y_0+\theta \left(y-y_0\right)\right)\end{array}$

$\left|f\left(x,y\right)-f_1\left(x,y\right)\right|=\frac{1}{\left(n+1\right)!}\left|{\left[\left(x-x_0\right)\frac{\partial }{\partial x}+\left(y-y_0\right)\frac{\partial }{\partial y}\right]}^{n+1}f\left(\alpha ,\beta \right)\right|$

where $\alpha =x_0+\theta \left(x-x_0\right)$, $\beta =y_0+\theta \left(y-y_0\right)$

$\begin{array}{l}{\Vert f-C^{\text{T}}B\Vert }_2^2\le {\Vert f-f_1\Vert }_2^2\\ ={\displaystyle \underset{0}{\overset{1}{\int }}{\displaystyle \underset{0}{\overset{1}{\int }}\frac{1}{\left(n+1\right){!}^2}}}{\left|{\left[\left(x-x_0\right)\frac{\partial }{\partial x}+\left(y-y_0\right)\frac{\partial }{\partial y}|\right]}^{n+1}f\left(\alpha ,\beta \right)\right|}^2\text{d}x\text{d}y\\ ={\displaystyle \underset{0}{\overset{1}{\int }}{\displaystyle \underset{0}{\overset{1}{\int }}\frac{1}{\left(n+1\right){!}^2}}}{\left|{\displaystyle \underset{m=0}{\overset{n+1}{\sum }}\left(\begin{array}{l}n+1\\ m\end{array}\right){\left(x-x_0\right)}^m{\left(y-y_0\right)}^{n+1-m}\frac{{\partial }^{n+1}}{\partial x^m\partial y^{n+1-m}}f\left(\alpha ,\beta \right)}\right|}^2\text{d}x\text{d}y\\ \le \frac{l^2}{\left(n+1\right){!}^2}{\displaystyle \underset{0}{\overset{1}{\int }}{\displaystyle \underset{0}{\overset{1}{\int }}{\left|{\displaystyle \underset{m=0}{\overset{n+1}{\sum }}\left(\begin{array}{l}n+1\\ m\end{array}\right){\left(x-x_0\right)}^m{\left(y-y_0\right)}^{n+1-m}}\right|}^2\text{d}x\text{d}y}}\\ =\frac{l^2}{\left(n+1\right){!}^2}{\displaystyle \underset{0}{\overset{1}{\int }}{\displaystyle \underset{0}{\overset{1}{\int }}{\left|{\left[\left(x-x_0\right)+\left(y-y_0\right)\right]}^{n+1}\right|}^2\text{d}x\text{d}y}}\le \frac{4l^2{\left(k_1+k_2\right)}^{2n+4}}{\left(n+1\right){!}^2\left(2n+3\right)\left(2n+4\right)}\end{array}$

${\Vert f-C^{\text{T}}B\Vert }_2^2\le \frac{2l{\left(k_1+k_2\right)}^{n+2}}{\left(n+1\right)!\sqrt{\left(2n+3\right)\left(2n+4\right)}}$. (6)

$\Vert f-C^{\text{T}}B\Vert \to 0$ as $n\to \infty$.

2.4. Derivatives of Multivariable Bernstein Polynomials

$B_{i,j}\left(x,y\right)={\displaystyle \underset{i=0}{\overset{n}{\sum }}{\displaystyle \underset{j=0}{\overset{m}{\sum }}{c}_{i,j}{B}_{i,n}\left(x\right)B_{j,m}\left(y\right)}}$. (7)

We can define p-th order partial derivative in x-axis direction in a following way

$\frac{{\partial }^p{B}_{i,j}\left(x,y\right)}{\partial x^p}={\displaystyle \underset{i=0}{\overset{n}{\sum }}{\displaystyle \underset{j=0}{\overset{m}{\sum }}{c}_{i,j}{D}^p{B}_{i,n}\left(x\right)B_{j,m}\left(y\right)}}$. (8)

Also q-th order partial derivative in y-axis direction in a following way

$\frac{{\partial }^q{B}_{i,j}\left(x,y\right)}{\partial x^q}={\displaystyle \underset{i=0}{\overset{n}{\sum }}{\displaystyle \underset{j=0}{\overset{m}{\sum }}{c}_{i,j}{D}^q{B}_{i,n}\left(x\right)B_{j,m}\left(y\right)}}$. (9)

Finally derivative of order $p+q$ is defined as follows

$\frac{{\partial }^{p+q}{B}_{i,j}\left(x,y\right)}{\partial x^p\partial x^q}={\displaystyle \underset{i=0}{\overset{n}{\sum }}{\displaystyle \underset{j=0}{\overset{m}{\sum }}{c}_{i,j}{D}^{p+q}{B}_{i,n}\left(x\right)B_{j,m}\left(y\right)}}$. (10)

3. Description of 2-D Helmholtz Problem

We consider the governing transient heat transfer on a two dimensional region $\Omega$ [31] [32].

$c\frac{{\partial }^{2}u}{\partial t^2}-\frac{\partial }{\partial x}\left(a_{11}\frac{\partial u}{\partial x}\right)-\frac{\partial }{\partial y}\left(a_{22}\frac{\partial u}{\partial y}\right)+a_{0}u=f\left(x,y,t\right)$. (11)

with the boundary conditions (b. c’s)

$u=\hat{u},q_n=a_{11}\frac{\partial u}{\partial x}{n}_x+a_{22}\frac{\partial u}{\partial y}{n}_y={\hat{q}}_n$ on $\gamma$. (11a)

and the initial condition

$u\left(x,y,0\right)=u_0\left(x,y\right)$ in $\Omega$. (11b)

where t specifies for the time and c, $a_{11}$, $a_{22}$, $a_0$, $u_0$, f, $\hat{u}$, ${\hat{q}}_n$ are regarded as functions of positions and/or time.

The corresponding homogeneous form of the above Equation (11) is

$c\frac{{\partial }^{2}u}{\partial t^2}-\frac{\partial }{\partial x}\left(a_{11}\frac{\partial u}{\partial x}\right)-\frac{\partial }{\partial y}\left(a_{22}\frac{\partial u}{\partial y}\right)+a_{0}u=0,u=0,q_n=0$ on $\gamma$. (12)

Approximate solution of Equation (12) is given as

$u\left(x,y,t\right)=U\left(x,y\right)e^{-\alpha t}$. (12a)

Thus Equation (12) becomes

$\frac{\partial }{\partial x}\left(a_{11}\frac{\partial u}{\partial x}\right)+\frac{\partial }{\partial y}\left(a_{22}\frac{\partial u}{\partial y}\right)+\lambda u=0.$ (13)

where

$\lambda =c\alpha +a_0$ and boundary conditions take the form $u=0$, $\frac{\partial u}{\partial n}=0$. (13a)

Here the eigenvalues $\lambda$ ’s and eigenfunctions U’s satisfy the time independent form of Equation (13). Then Equation (13) is known as Helmholtz equation which can be solved using separation of variables technique. Now the approximate solution of Equation (13) is

$u_k\left(x,y,t\right)={\displaystyle \underset{k=1}{\overset{n}{\sum }}{a}_k}\left(t\right)U_k\left(x,y\right)$. (14)

where $U_k\left(x,y\right)$ denotes spatial variable and parameter $a_k\left(t\right)$ defines the function of time only.

Galerkin Formulation for One Dimensional Heat Equation

Here we applied a modified weighted residual technique which involves only first derivative instead of second derivative terms in the residual equation. This can be performed applying integration by parts in which the operational matrices for integration and the product is utilized to convert the given differential eigenvalue problems to a system of algebraic equations.

One dimensional parabolic partial differential equation

$\rho cA\frac{\partial u}{\partial t}-\frac{\partial }{\partial x}\left(kA\frac{\partial u}{\partial x}\right)=q\left(x,t\right),a<t<b,t\ge 0$ (15)

with the boundary and initial conditions

$u\left(a,t\right)=f_1(\; t\; )$

$u\left(b,t\right)=f_2(\; t\; )$

$u\left(x,0\right)=f\left(x\right)$.

which governs the transient heat transfer in one-dimensional system (i.e., a plane wall), here u denotes the temperature, k the thermal conductivity, $\rho$ the density, A the cross-sectional area, c the specific heat, q the heat generation per unit length. Here $u\left(x,t\right)$ represents the temperature of a rod for both space and time and $U_k\left(x\right)$ is trial or basis function. If non-homogeneous heat Equation (15) with a lateral heat loss reduces to homogeneous equation (q = 0).

$\rho cA\frac{\partial u}{\partial t}-\frac{\partial }{\partial x}\left(kA\frac{\partial u}{\partial x}\right)=0$. (16)

Let us assume the approximate solution of Equation (16)

The solution of Equation (15) can be written as

${\tilde{u}}_k\left(x,t\right)={\displaystyle \underset{k=1}{\overset{n}{\sum }}{a}_k}\left(t\right)U_k(\; x\; )$

The Galerkin weighted residual equation can be written as

${\displaystyle \underset{a}{\overset{b}{\int }}\left[\rho cA\frac{\partial \tilde{u}\left(x,t\right)}{\partial t}-\frac{\partial }{\partial x}\left\{kA\frac{\partial \tilde{u}\left(x,t\right)}{\partial x}\right\}\right]U_j\left(x\right)\text{d}x}=0,j=1,2,3,\cdots ,n$ (17)

${\tilde{u}}_k\left(x,t\right)={\displaystyle \underset{k=1}{\overset{n}{\sum }}{F}_{j,k}\frac{\text{d}{a}_k}{\text{d}t}}+{\displaystyle \underset{k=1}{\overset{n}{\sum }}{G}_{j,k}}{a}_k\left(t\right)$.

On simplification the by using the conditions the above residual Equation (17) takes the matrix form

$\left[F\right]\left\{\frac{\text{d}{a}_k\left(t\right)}{\text{d}t}\right\}-\left[G\right]\left\{a_k\left(t\right)\right\}=\left\{0\right\}$. (18)

where $\frac{\text{d}a}{\text{d}t}=\hat{u}$ refers to the time derivative.

Hence the matrix form of Equation (18)

where $\hat{u}=\frac{\text{d}c}{\text{d}t}$ refers to the time derivative.

$\left[G\right]=G_{j,k}={\displaystyle \underset{a}{\overset{b}{\int }}kA\frac{\text{d}{U}_j}{\text{d}x}}\frac{\text{d}{U}_k}{\text{d}x}\text{d}x$.

$\left[F\right]=F_{j,k}={\displaystyle \underset{a}{\overset{b}{\int }}\rho cA{U}_j}\left(x\right)U_k\left(x\right)\text{d}x$.

The solution of Equation (18) can be given as

$\tilde{u}\left(x,t\right)=\hat{u}\left(x\right){\text{e}}^{-\lambda t}$. (19)

Hence the matrix form of (18)

$\left[G\right]\left\{\hat{u}\right\}-\lambda \left[F\right]\left\{\hat{u}\right\}=\left\{0\right\}$. (20)

Solving (19) we can determine $\lambda$ and nonzero $\hat{u}\left(x\right)$ such that (19) holds and the boundary conditions of the problem met. Equation (19) describe an eigenvalue problem; $\lambda$ being the eigenvalues and $\hat{u}\left(x\right)$ be the eigenfunction.

4. Brief Description of Hyperbolic Equations

The transverse motion of a membrane, is governed by a partial differential equation of the form as given in (14), where, $u\left(x,y,t\right)$ denotes the transverse deflection, c is the material density of the membrane, $a_{11}$ and $a_{22}$ are the tensions in the x and y directions of the membrane, $a_0$ is the modulus of elastic foundation on which the membrane is stretched (often, $a_0=0$ ) i.e., there is no foundation), and $f\left(x,y,t\right)$ is the transversely distributed force. The above equation is known as wave equation, and is classified mathematically as a hyperbolic equation. The function $u\left(x,y,t\right)$ must be determined such that it satisfies Equation (19) in a region $\Omega$ together with specified boundary and initial conditions. The problem of finding the solution as given by Equation (14) holds for homogeneous boundary and initial conditions and $f=0$ is called an eigenvalue problem.

4.1. Solution of Two-Dimensional Wave Equation in Rectangular Membrane

We consider the equation for the vibration of a tightly stretched membrane as the membrane of a drum. We assume that the membrane is of uniform tension and the tension per unit length is same in all directions at every point and m be the mass of per unit of area of the membrane.

$\frac{{\partial }^{2}u}{\partial t^2}=c^2\left[\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}\right]$

where

$c^2=\frac{T}{m}$. (21)

The fact that we are keeping the edges of the membrane fixed is expressed by the boundary conditions

$u\left(x,0,t\right)=u\left(x,b,t\right)$ (21a)

$u\left(0,y,t\right)=u\left(a,y,t\right)$ (21b)

We must also specify how the membrane is initially deformed and set into motion. This is done via the initial conditions

$u\left(x,y,0\right)=0$ (21c)

The solution of Equation (21) applying the separation of variables gives the particular values of l and n.

$u\left(x,y,t\right)={\displaystyle \underset{l=1}{\overset{\infty }{\sum }}{\displaystyle \underset{n=1}{\overset{\infty }{\sum }}\sin \frac{l\pi x}{a}\sin \frac{n\pi y}{b}}}$.

${\alpha }_{ln}=\pi c\sqrt{\left[{\left(\frac{l}{a}\right)}^2+{\left(\frac{n}{b}\right)}^2\right]}$. (22)

These are the solutions of the wave Equation (22) which are zero on the boundary of the rectangular membrane. These functions are known as eigenfunctions and the numbers are known eigenvalues.

The general solution is

$u=A_{ln}\cos \left({\alpha }_{ln}t\right)+B_{ln}\sin \left({\alpha }_{ln}t\right)$. (23)

4.2. Vibration of Distributed Inhomogeneous System in a Rectangular Domain

We consider free vibrations of a uniformly stretched inhomogeneous membrane of rectangular shape with clamped edges in an elastic medium. We consider the following eigenvalue problem

$\left[\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}\right]+\left[\lambda \rho \left(x,y\right)-q\left(x,y\right)\right]u=0$. (24)

where the boundary conditions are

$u\left(0,y\right)=u\left(a,y\right)=0$ and $u\left(x,0\right)=u\left(x,b\right)=0$. (24a)

Let $\tilde{u}\left(x,y\right)$ denotes an approximation to the exact solution.

$\tilde{u}\left(x,y\right)={\theta }_0\left(x,y\right)+{\displaystyle \underset{i=1}{\overset{n}{\sum }}{a}_i}{\theta }_i\left(x,y\right)$ and ${\theta }_0\left(x,y\right)=0$ , is specified by the Dirichlet boundary conditions, of Equation (24a).

$R\left(x,y\right)=\left[\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}\right]+\left[\lambda \rho \left(x,y\right)-q\left(x,y\right)\right]u=0$. (24b)

where $R\left(x,y\right)=0$ is the residual function.

The Green’s theorem states that if A is a region in the xy plane bounded by a closed curve $\Gamma$ then we suppose:

$U\left(x,y\right)=F\left(x,y\right)i+G\left(x,y\right)j$ is a continuous vector field defined on a region A in Ω. Moreover, suppose F and G have continuous partial derivatives and that the boundary $\Gamma$ is a simple closed curve C with positive orientation. Then

${\displaystyle \underset{A}{\iint }\left(\frac{\partial G}{\partial x}-\frac{\partial F}{\partial y}\right)\text{d}x\text{d}y}={\displaystyle \underset{\Gamma }{\oint }F\text{d}x+G\text{d}y}$. (25)

We can write

${\displaystyle \underset{A}{\iint }\left(\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}\right)\text{d}x\text{d}y}={\displaystyle \underset{\Gamma }{\oint }F\text{d}x+G\text{d}y}$. (26)

where

$F=-\frac{\partial u}{\partial y}{\theta }_i$ and $G=\frac{\partial u}{\partial x}{\theta }_i$.

${\displaystyle \underset{A}{\iint }R\left(x,y\right){\theta }_i}\text{d}x\text{d}y={\displaystyle \underset{\Gamma }{\oint }F\text{d}x+G\text{d}y}.$ (26a)

${\displaystyle \underset{A}{\iint }R\left(x,y\right){\theta }_i\text{d}x\text{d}y}={\displaystyle \underset{\Gamma }{\oint }\left\{\frac{{\partial }^2\tilde{u}}{\partial x^2}+\frac{{\partial }^2\tilde{u}}{\partial y^2}+\left[\lambda \rho \left(x,y\right)-q\left(x,y\right)\right]\tilde{u}\right\}{\theta }_i}\text{d}x\text{d}y$ (26b)

${\displaystyle \underset{\Gamma }{\oint }F\text{d}x+G\text{d}y}={\displaystyle \underset{{\rm T}}{\oint }{\theta }_i\left\{\frac{\partial u}{\partial y}\text{d}x+\frac{\partial u}{\partial x}\text{d}y\right\}}$ (27)

Using Equations (26a), (27) in Equation (26b)

$\begin{array}{l}-{\displaystyle \underset{A}{\iint }\left\{\frac{\partial {\theta }_i}{\partial x}\frac{\partial \tilde{u}}{\partial x}+\frac{\partial {\theta }_i}{\partial y}\frac{\partial \tilde{u}}{\partial y}-q\left(x,y\right){\theta }_i\tilde{u}\right\}\text{d}x\text{d}y}\\ -\lambda {\displaystyle \underset{A}{\iint }\rho \left(x,y\right){\theta }_i\tilde{u}\text{d}x\text{d}y}+{\displaystyle \underset{{\rm T}}{\oint }{\theta }_i\left\{\frac{\partial \tilde{u}}{\partial y}\text{d}x+\frac{\partial \tilde{u}}{\partial x}\text{d}y\right\}}=0\end{array}$

Since

${\displaystyle \underset{{\rm T}}{\oint }{\theta }_i\left\{\frac{\partial \tilde{u}}{\partial y}\text{d}x+\frac{\partial \tilde{u}}{\partial x}\text{d}y\right\}}=0$ (28)

We finally obtain

$-{\displaystyle \underset{A}{\iint }\left\{\frac{\partial {\theta }_i}{\partial x}\frac{\partial \tilde{u}}{\partial x}+\frac{\partial {\theta }_i}{\partial y}\frac{\partial \tilde{u}}{\partial y}-q\left(x,y\right){\theta }_i\tilde{u}\right\}\text{d}x\text{d}y}=\lambda {\displaystyle \underset{A}{\iint }\rho \left(x,y\right){\theta }_i\tilde{u}\text{d}x\text{d}y}$ (29)

5. Error and Convergence of Galerkin MWR

Let $\left\{{\phi }_n|n=1,2,3,\cdots \right\}$ be the set of normalized eigenfunctions of the Sturm-Liouville system. If f be a function in Hilbert space H, then

$\underset{n\to \infty }{\lim }\Vert f-{\displaystyle \underset{k=1}{\overset{n}{\sum }}{c}_k}{\varphi }_k\Vert =0$. (30a)

where

$c_k={\displaystyle {\int }_0^1{\bar{\varphi }}_{k}f\left(x\right)r\left(x\right)\text{d}x}$. (30b)

We consider the linear operator L with a field of definition $F_L$ , that is, Lu is defined for $u\in F_L$. The inner product is

${\left(u,v\right)}_L=\left(u,Lv\right)={\displaystyle \underset{v}{\int }uLv\text{d}x}$.

The operator is symmetric if for elements $u,v\in F_L$

$\left(u,Lv\right)=\left(v,Lu\right)$.

The operator is positive definite if for any function in $F_{\text{}L}$, not identically zero:

$\left(u,Lu\right)\ge 0$ and is positive bounded below if for any $u\in F_L$ ;

$\left(u,Lu\right)\ge \beta \left(u,u\right)$, for $\beta >0$.

Uniform Convergence: For any $\epsilon >0$, we can find n such that

$\left|u\left(x\right)-{\tilde{u}}_n\left(x\right)\right|<\epsilon$ (31)

Convergence in energy requires $\Vert u-{\tilde{u}}_n\Vert <\epsilon$ where energy is defined as

$\left|u\right|=\left(u,Lu\right)$.

The corresponding error bounds are point wise error in Equation (31). A sequence $u_n$, converges weakly to an element u of a space if $\underset{n\to \infty }{Lt}\left(u_n,\varphi \right)=\left(u,\varphi \right)$ holds for all $\varphi$ in the space. Sequences produced by the Galerkin method sometimes converge weakly to a generalized solution.

A set of trial functions form a complete set of functions and they are complete in a space if any function in the space can be expanded in terms of the set of functions, for sufficiently large n [33].

$\Vert u-{\displaystyle \underset{k=1}{\overset{\infty }{\sum }}{c}_k{\varphi }_k}\Vert <\epsilon$. (32)

Since the powers of Bernstein polynomials $x^i{\left(1-x\right)}^{n-i}$ ( $n=0,1,2,3,\cdots$ ) have sufficient continuity and they form an infinite sequence of functions which is complete. Also the convergence of the eigenvalues by Galerkin WRM method is measured by the relative error.

$e_k=\left|\frac{{\lambda }^{Exact}-{\lambda }^{\left(Approx.\right)}}{{\lambda }^{Exact}}\right|<\delta$.

$\delta \le {10}^{-10}$ depends upon the problems.

We consider a two dimensional boundary value problem.

$-{\displaystyle \underset{i,j=1}{\overset{3}{\sum }}{B}_{i,j}\left(x\right)\frac{\partial }{\partial x_i}\left(\frac{\partial u}{\partial x_j}\right)}+{\displaystyle \underset{i=1}{\overset{3}{\sum }}{C}_i\left(x\right)}\frac{\partial u}{\partial x_j}+F\left(x\right)u=0$ (33)

is non-self-adjoint unless $C_i=0$. For convergence of the Galerkin method, we assume there exists a constant ${\varsigma }_0$ such that for any point in V and real numbers $l_i$,

$-{\displaystyle \underset{i,j=1}{\overset{3}{\sum }}{B}_{i,j}{l}_i{l}_j}\ge {\displaystyle \underset{i=1}{\overset{3}{\sum }}{\zeta }_o{l}_i^2}$. (34)

The coefficients $B_{i,j}$, and their first derivatives are continuous in $\Omega$ and the coefficients C and F are continuous in $\bar{\Omega }$. We assume the problem has a unique solution. The proof of the theorem is illustrated in and the references therein [33]. Here $\Omega$ is a finite bounded domain and $\bar{\Omega }=S\cup \Omega$.

$L_{0}u=-{\displaystyle \underset{i,j=1}{\overset{3}{\sum }}\frac{\partial }{\partial x_i}\left(B_{i,j}\frac{\partial u}{\partial x_j}\right)}$. (35)

where the operator $L_0$ is positive bounded below for the set of functions which vanishes on S. The convergence proofs are done in terms of the operator $L_0$ and the remaining terms in the differential equation are bounded in terms of this operator and its inverse.

Theorem2: Assume the problem (11) is unique for the boundary conditions $u=0$ on $S_1$ , $\frac{\partial u}{\partial n}=0$ on $S_2$, $\frac{\partial u}{\partial n}+\mu n=0$ on $S_3$.

Exploiting the Galerkin method using trial functions which are complete in energy of $L_0$ , are linearly independent, and satisfy $u=0$ on $S_1$. Then the approximate solution ${\tilde{u}}_n$ , converges in energy of $L_0$ to the solution.

Let us consider the two-dimensional the following eigenvalue problem

$-{\displaystyle \underset{i,j=1}{\overset{3}{\sum }}{B}_{i,j}\left(x\right)\frac{\partial }{\partial x_i}\left(\frac{\partial u}{\partial x_j}\right)}+{\displaystyle \underset{i=1}{\overset{3}{\sum }}{C}_i\left(x\right)}\frac{\partial u}{\partial x_j}+\lambda F\left(x\right)u=0$. (36)

$u=0$ on S where $B_{i,j}$, satisfies Equation (34).

Theorem 3: If the Galerkin method is applied to Equation (36) using a system of functions which is complete in energy and linearly independent, then the approximate eigenvalues converge to the true eigenvalue. This proves convergence of the Galerkin method in section 5 when applied to Equation (36).

6. Numerical Experiments

In this section we apply Galerkin MWR exploiting Bernstein polynomial only to the eigenvalue problems frequently encountered in science and engineering in one and two dimensions. For validation of our proposed technique formulated in section 4, we demonstrate three linear problems in one dimension and two linear problems in two dimensions.

In this section we apply the current techhnique exploiting Bernstein polynomial only to the eigenvalue problems frequently encountered in science and engineering in one and two dimensions. For validation of our proposed technique formulated in section 4, we demonstrate three linear problems in one dimension and two linear problems in two dimensions.

Example 1: To determine the boundaries of parametric instability regions for the following modified Mathieu equation using the method of accelerated convergence which is described in details in [22].

$\frac{{\text{d}}^{2}u}{\text{d}{t}^2}+\left[\lambda -q\left(t,e\right)\right]u=0$. (37)

$\left|\lambda \right|<\infty$, $q\left(t,e\right)=\frac{e\cos \left(2\pi t\right)}{1-e\cos \left(2\pi t\right)}$, $\left|e\right|<1$.

For $\left|e\right|\ll 1$, Equation (37) differs from the classical Mathieu equation.

To compute the eigenvalues depending on parameter e which determines the boundaries of the parametric instability [22], we need to construct periodic solutions of that:

$u\left(-1\right)=u\left(1\right)=0$, $u^{\prime }\left(-1\right)=u^{\prime }\left(1\right)=0$. (37a)

Since the function $q\left(t,e\right)$ is even with respect to t, $q\left(-t,e\right)=q\left(t,e\right)$ in (the conditions of perioiodicity in (37a) are equivalent to the conditions of the first kind given as

$u\left(0\right)=u\left(1\right)=0$ ; (37b)

and the second kind given as

$u^{\prime }\left(0\right)=u^{\prime }\left(1\right)=0$ ; (37c)

here the vibration of elastic Crankshaft with concentrated load yield the following equation with periodic co-efficient and modulation depth e illustrated in [22] is considered.

$\frac{{\text{d}}^{2}u}{\text{d}{x}^2}+\frac{\lambda }{1-e\cos \left(2\pi x\right)}u=0$, $0\le \left|e\right|<1$. (38)

Example 2: Let us consider one dimensional Helmholtz equation [2].

$\{\begin{array}{l}\frac{{\text{d}}^{2}y}{\text{d}{x}^2}={\lambda }^{2}u\\ u\left(0\right)=u\left(1\right)=0\end{array}$ (39)

Example 3:

We consider a plane wall, initially at a uniform temperature [32].

$-\frac{{\text{d}}^{2}u}{\text{d}{x}^2}=-\lambda u$, (40)

Set 1: $u\left(0\right)=0$, $u\left(1\right)=0$, (40a)

Set 2: $u\left(0\right)=0$, ${\left(\frac{\text{d}u}{\text{d}x}+u\right)|}_{x=1}=0$. (40b)

Table 2 displays the first ten eigenvalues for n = 15, exploiting Bernstein polynomials.

Example 4: The well-known transient heat conduction [31] [32] in a square region with boundary and initial conditions is considered.

$\frac{\partial u}{\partial t}=\left[\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}\right]=f\left(x,y,t\right)$, $\Omega =\left\{\left(x,y\right):0<x<1,0<y<1\right\}$. (41)

Subject to the boundary conditions

The boundary conditions are:

$u\left(x,1,t\right)=0$, $u\left(1,y,t\right)=0$. (41a)

$\frac{\partial u}{\partial x}\left(0,y,t\right)=0$ and $\frac{\partial u}{\partial y}\left(x,0,t\right)=0$. (41b)

and the initial condition is

$u\left(x,y,0\right)=0$. (41c)

Homogeneous form of (41) for $f=0$ is

$\frac{\partial u}{\partial t}=\left[\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}\right]$. (41d)

The derivative boundary conditions will be applied integrand of the residual equation. Hence for executing only the Dirichlet boundary conditions in Equation (41a), we modify the polynomial basis as

$B_{i,k}^{\left(n,m\right)}\left(x,y\right)=\left[\left(\begin{array}{c}n\\ i\end{array}\right)\left(\begin{array}{c}m\\ k\end{array}\right)x^i{\left(1-x\right)}^{n-i}{y}^k{\left(1-y\right)}^{m-k}\right]\times \left[\left(1-x\right)\left(1-y\right)\right]$ (42)

The exact solution of (42) is $u\left(x,t\right)=T_0\hat{u}\left(x,y\right){\text{e}}^{-\lambda t}$, where $\lambda$ refers to the required eigenvalue which corresponds to the eigenfunctions $\stackrel{⌢}{u}\left(x,y\right)$. The analytical eigenvalue ${\lambda }_{m,n}={\pi }^2\left(m^2+n^2\right)/4$, $m,n=1,3,5,\cdots$

Example 5: The homogeneous membrane of sides 4 ft by 2 ft, fixed on all its four edges i.e., $u=0$ on $\Gamma$ is considered [32]. We assume that the tension of the membrane is $T=\text{12}.\text{5}\text{lb}/{\text{ft}}^{-1}$ (i.e. $a_{11}=a_{22}=12.5$ ). The exact natural frequencies of the rectangular membrane are obtained by considering the full domain

$\frac{{\partial }^{2}u}{\partial t^2}=5\left[\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}\right]$ (43)

The fact that we are keeping the edges of the membrane fixed is expressed by the boundary conditions

$u\left(0,y,t\right)=u\left(4,y,t\right)=0$ (43a)

$u\left(x,0,t\right)=u\left(x,2,t\right)=0$ (43b)

and the initial conditions

$u\left(x,y,0\right)=0$ (43c)

For implementing boundary conditions in Equations (43a), (43b) we modify the polynomial basis as

$B_{i,k}^{\left(n,m\right)}\left(x,y\right)=\left[\left(\begin{array}{c}n\\ i\end{array}\right)\left(\begin{array}{c}m\\ k\end{array}\right)x^i{\left(4-x\right)}^{n-i}{y}^k{\left(2-y\right)}^{m-k}\right]\times \left[x\left(4-x\right)y\left(2-y\right)\right]$ (44)

7. Result and Discussions

In Table 1, we computed even and odd eigen solutions for different modulation depth e ranges from 0.1 to 0.9. Analytical eigen-solutions using perturbation method are listed in [22]. Our computed relative errors for the smallest modes are less than 5% which shows that the present technique produced significantly small errors. From Table 2, the smallest eigenvalues attains the accuracy up to 10⁻¹⁶ and error increases rapidly for higher eigenvalues than the lower values which is better than boundary method. As we increase the grid points or nodes from n = 15 to n = 30, the error decays fast for all the eigenvalues and consistent accuracy is obtained up to 10⁻¹³.

We observed that increasing of nodes reveal the stable behaviour of all the eigenvalues for n = 30. It is also observed that our present approach accomplishes accurate results and is compatible to the existence new boundary approach for one dimensional Helmholtz equation. Comparison of relative errors obtained for Gal. Relative errors between our proposed Galerkin MWR and FEM [32] are compared which are displayed in Figure 1. In Table 3(a) and Table 3(b) we computed seven eigen frequencies for the two sets of boundary conditions exploiting current formulation and comparison to the exact and finite element results are shown with relative errors. From these tables, we observed that the smallest error decreases up to ${10}^{-7}$ which is much smaller than that of FEM [32]. Eigen frequencies using different degree of polynomial are displayed in Table 4 of example 4. We noticed that the smallest order of the accuracies for the largest eigenvalues are 10⁻², 10⁻³, 10⁻⁴ respectively with 4, 8, 10, 12 Bernstein polynomials and the first three eigenvalues converge to the exact results for n = 13 illustrated in Table 4. This confirms that the accuracy increases with the increased degree of polynomials. Our current method produces small errors in percentages (%) than those attained by FEM [31].

Plot of the relative errors between exact and approximate eigenvalues are depicted in Figure 2 for different values of n. Also, the relative errors obtained by Galerkin MWR for different degree of Bernstein polynomials are depicted in Figure 3. Table 5 shows that the maximum relative error for the smallest and the largest eigenvalues using present method are 10⁻⁶ and 10⁻⁷ respectively, whereas applying FEM [32], the errors obtained as 10⁻³, respectively. From the results in tabular form we conclude that the Galerkin MWR is much accurate, efficient and compatible with other techniques.

Table 1. Comparison of eigenvalues for vibration of crankshaft of Example 1.

Table 2. Comparison of relative errors obtained by Galerkin MWR with for Example 2.

(a) (b)

Table 3. (a) Table of calculated eigenvalues and their relative errors computed by the current method for example 3; (b). List of eigenvalues with their relative errors attained by present method and FEM for Example 3.

Table 4. List of eigenvalues with their relative errors attained by present method and FEM for Example 4.

Table 5. Comparison of natural frequencies and relative errors using Gal. WRM with FEM for Example 5.

Figure 1. Comparison of relative errors obtained for Gal. MWR and FEM.

Figure 2. comparison of eigenvalues using GAL MWR and exact result.

Figure 3. Relative errors attained by GAL MWR for increasing polynomial.

8. Conclusions

For the vibration of crankshaft computed results utilizing Galerkin MWR are very close to the analytical results [22] and error increased as the modulation depth closer to one. In the case of 1-D heat equation for the plane wall [32] and 2-D heat equation [31], the smallest errors by present method attains superior accuracy than finite element methods using linear and quadrilateral elements worked out in those studies. Eigenvalues achieved from 1D Helmholtz equation for smaller eigenvalues converge faster as the degree of Bernstein polynomials increases and attains many accurate results. For homogeneous rectangular membrane, relative errors in the present method are smaller than that of finite element method [32]. This proves that our present method is much more efficient than various numerical methods available in the literature and well suited which has much applicability in physical and engineering models.

Acknowledgements

The first author is very much thankful to Muhammad Sajjad Hossain, Assistant Professor of Mathematics, Department of Arts & Sciences, Ahsanullah University of Science and Technology, Dhaka, for his knowledgeable input, inexorable assistance and for always being keen to help in editing this article.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Liu J., Zheng, Z. and Xu, Q. (2012) Bernstein-Polynomials-Based Highly Accurate Methods for One-Dimensional Interface Problems. Journal of Applied Mathematics, 2012, 1-11. https://doi.org/10.1155/2012/859315
[2]	Reutskiy, S.Y. (2006) The Method of Fundamental Solutions for Helmholtz Eigenvalues Problems in Simply and Multiply Connected Domains. Engineering Analysis with Boundary Elements, 30, 150-159. https://doi.org/10.1016/j.enganabound.2005.08.011
[3]	Yousefi, S.A., Behroozifar, M. and Dehghan, M. (2011) The Operational Matrices of Bernstein Polynomials for Solving the Parabolic Equation Subject to Specification of Mass. Journal of Computational and Applied Mathematics, 235, 5272-5283. https://doi.org/10.1016/j.cam.2011.05.038
[4]	Zhou-Lian, Z., Chang-Jiang L.H., Xiao, T. and Chen, S-L. (2009) Free Vibration Analysis of Rectangular Orthotropic Membranes in Large Deflection. Mathematical Problems in Engineering, 2009, 1-9. https://doi.org/10.1155/2009/634362
[5]	Mirkov, N. and Rasuo, B. (2013) Bernstein Polynomial Collocation Method for Elliptic Boundary Value Problems. Proceedings in Applied Mathematics and Mechanics, 13, 421-422. https://doi.org/10.1002/pamm.201310206
[6]	Yousefi, S.A. and Barikbiz, M.Z. (2010) Bernstein Ritz-Galerkin Method for Solving a Initial Boundary Value Problem that Combines Neumann and Integral Condition for the Wave Equation. Numerical Methods for Partial Differential Equations, 26, 1236-1246. https://doi.org/10.1002/num.20521
[7]	Park, J., Park, I. and Lee, U. (2014) Transverse Vibration and Waves in a Membrane: Frequency Do Main Spectral Element Modeling and Analysis. Mathematical Problems in Engineering, 2014, 1-14. https://doi.org/10.1155/2014/642782
[8]	Sato, K. (1974) Free Vibration Analysis of a Composite Elliptical Membrane Consisting of Confocal Elliptical Parts. Journal of Sound and Vibration, 34, 161-171. https://doi.org/10.1016/S0022-460X(74)80301-9
[9]	Chen, S.-L. and Zheng, Z.-L. (2003) Large Deformation of Circular Membrane under the Concentrated Force. Applied Mathematics and Mechanics (English Edition), 24, 28-31. https://doi.org/10.1007/BF02439374
[10]	Kim, Y.Y. and Kang, J.H. (1996) Free Vibration Analysis of Membranes Using Wave-Type Base Functions. The Journal of the Acoustical Society of America, 99, 2938-2946. https://doi.org/10.1121/1.414822
[11]	Laura, P.A.A., Rossi, R.E. and Gutierrez, R.H. (1997) The Fundamental Frequency of Non-Homogeneous Rectangular Membrane. Journal of Sound and Vibration, 204, 229-306. https://doi.org/10.1006/jsvi.1996.0931
[12]	Vega, D.A., Vera, S.A. and Laura, P.A.A. (1999) Fundamental Frequency of Vibration of Rectangular Membranes with an Internal Oblique Support. Journal of Sound and Vibration, 224, 780-783. https://doi.org/10.1006/jsvi.1999.2219
[13]	Masad, A. (1996) Free Vibrations of a Non-Homogeneous Rectangular Membrane. Journal of Sound and Vibration, 195, 674-678. https://doi.org/10.1006/jsvi.1996.0454
[14]	Ho, S.H. and Chen, C.K. (2000) Free Vibration Analysis of Non-Homogeneous Rectangular Membranes Using a Hybrid Methods. Journal of Sound and Vibration, 233, 547-555. https://doi.org/10.1006/jsvi.1999.2808
[15]	Ersoy, H., Civalekb, O. and Ozpolatb, L. (2010) Free Vibration Analysis of Rectangular Membranes with Variable Density Using the Discrete Singular Convolution Approach. Asian Journal of Civil Engineering, 11, 83-94.
[16]	Cortinez, V.H. and Laura, P.A.A. (1992) Vibrations of Non-Homogeneous Rectangular Membranes. Journal of Sound and Vibration, 156, 217-225. https://doi.org/10.1016/0022-460X(92)90694-S
[17]	Reutskiy, S.Y. (2007) The Methods of External and Internal Excitation for Problems of Free Vibrations of Non-Homogeneous Membranes. Engineering Analysis with Boundary Elements, 31, 906-918. https://doi.org/10.1016/j.enganabound.2007.04.003
[18]	Amore, P. (2009) A New Method for Studying the Vibration of Non-Homogeneous Membranes. Journal of Sound and Vibration, 321, 104-114. https://doi.org/10.1016/j.jsv.2008.09.024
[19]	Doha, E.H. and Abd-Elhameed, W.M. (2005) Accurate Spectral Solutions for the Parabolic and Elliptic Partial Differential Equations by the Ultra-Spherical Tau Method. Journal of Computational and Applied Mathematics , 181, 24-45. https://doi.org/10.1016/j.cam.2004.11.015
[20]	Sahuck, O.H. (2019) An Efficient Spectral Method to Solve Multi-Dimensional Linear Partial Differential Equations Using Chebyshev Polynomials. Mathematics, 7, 1-21. https://doi.org/10.3390/math7010090
[21]	Dehghan, M. (2006) A Computational Study of the One-Dimensional Parabolic Equation Subject to Non-Classical Boundary Specifications. Numerical Methods for Partial Differential Equations, 22, 220-257. https://doi.org/10.1002/num.20071
[22]	Leonid, A.D. and Nesterov, S.I. (2006) High Precision Method for Eigenvalue Problem. Chapman Hall/CRC, Boca Raton, FL.
[23]	Alshbool, M.H.T., Bataineh, A.S., Hashim, I. and Isik, O.R. (2015) Approximate Solutions of Singular Differential Equations with Estimation Error by Using Bernstein Polynomials. International Journal of Pure and Applied Mathematics, 100, 109-125. https://doi.org/10.12732/ijpam.v100i1.10
[24]	Zhou, J. and Li, H. (2015) A Ritz-Galerkin Approximation to the Solution of Parabolic Equation with Moving Boundaries. Boundary Value Problems, 236, 1-17. https://doi.org/10.1186/s13661-015-0503-5
[25]	Doha, E.H., Bhrawy, A.H. and Saker, M.A. (2011) On the Derivatives of Bernstein Polynomials, an Application for the Solution of High Even-Order Differential Equations. Boundary Value Problems, 2011, Article ID: 829543. https://doi.org/10.1155/2011/829543
[26]	Farzana, H., Islam, M.S. and Bhowmik, S.K. (2015) Computation of Eigenvalues of the Fourth Order Sturm-Liouville BVP by Galerkin Weighted Residual Method. British Journal of Mathematics and Computer Science, 9, 73-85. https://doi.org/10.9734/BJMCS/2015/15370
[27]	Islam, M.S., Farzana, H. and Bhowmik, S.K. (2017) Numerical Solutions of Sixth Order Eigenvalue Problems Using Galerkin Weighted Residual Method. Differential Equations and Dynamical Systems, 25, 187-205. https://doi.org/10.1007/s12591-016-0323-9
[28]	Farzana, H. and Bhowmik, S.K. (2019) Comparative Study on Sixth Order Boundary Value Problems with Application to Linear Hydrodynamic Stability Problem and Benard Layer Eigenvalue Problem. Differential Equations and Dynamical Systems. https://doi.org/10.1007/s12591-019-00509-4
[29]	Cheng, F. (1983) The Rate of Convergence of Bernstein Polynomials of Functions of Bounded Variation. Journal of Approximation Theory, 39, 259-274. https://doi.org/10.1016/0021-9045(83)90098-9
[30]	Martinez, F.L. (1989) Some Properties of Two-Dimensional Bernstein Polynomials. Journal of Approximation Theory, 59, 300-306. https://doi.org/10.1016/0021-9045(89)90095-6
[31]	Thowhida, A., Islam, M.S. and Ahmad, S.U. (2010) Eigenvalue Analysis of 2D Helmholtz Equation on Quadrilateral Elements. Dhaka University Journal of Science, 58, 141-142.
[32]	Reddy, J.N. (1993) An Introduction to Finite Element Method. In: Mechanical Engineering, McGraw-Hill, New York.
[33]	Finlayson, A.B. (1972) The Method of Weighted Residual and Variational Principles with Application in Fluid Mechanics, Heat and Mass Transfer. Academic Press, New York.