# On Solving One-Dimensional Partial Differential Equations With Spatially Dependent Variables Using the Wavelet-Galerkin Method

Simon Jones, Mathias Legrand

► To cite this version:

Simon Jones, Mathias Legrand. On Solving One-Dimensional Partial Differential Equations With Spatially Dependent Variables Using the Wavelet-Galerkin Method. Journal of Applied Mechanics, 2013, 10.1115/1.4023637. hal-01330602

HAL Id: hal-01330602

https://hal.science/hal-01330602

Submitted on 17 Jun 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# On Solving One-Dimensional Partial Differential Equations With Spatially Dependent Variables Using The Wavelet-Galerkin Method

Simon Jones<sup>(1)</sup>, Mathias Legrand<sup>(2)</sup>

**Abstract** The discrete orthogonal wavelet-Galerkin method is illustrated as an effective method for solving partial differential equations (PDE's) with spatially varying parameters on a bounded interval. Daubechies scaling functions provide a concise but adaptable set of basis functions and allow for implementation of varied loading and boundary conditions. These basis functions can also effectively describe $C^0$ continuous parameter spatial dependence on bounded domains. Doing so allows the PDE to be discretized as a set of linear equations composed of known inner products which can be stored for efficient parametric analyses. Solution schemes for both free and forced PDE's are developed; natural frequencies, mode shapes, and frequency response functions for an Euler-Bernoulli beam with piecewise varying thickness are calculated. The wavelet-Galerkin approach is shown to converge to the first four natural frequencies at a rate greater than that of the linear finite element approach; mode shapes and frequency response functions converge similarly.

**Keywords** Wavelets — Wavelet-Galerkin method — Partial Differential Equation — Boundary Value Problem — Euler-Bernoulli beam

## 1. Introduction

Mathematical modeling of physical systems often requires the use of a governing partial differential equation (PDE). Phenomena such as elastodynamics of a beam, heat conduction through a plate, or fluid flow through a pipe can all be simplified to respective PDE's with known analytic solutions. When the system of interest has greater complexity it is still possible to describe the system dynamics using a PDE, though complex geometries can make analytical solutions intractable. In such cases solutions can be approximated using numerical methods, such as the finite element method, the finite difference method, the boundary element method, method of lines, general Galerkin methods, or multigrid methods [1, 2, 3].

Relatively recently a new method of solving PDE's has emerged, known as the wavelet-Galerkin method [4, 5]. Wavelets are well localized, oscillatory functions which have been shown to provide a basis on unbounded domains $L^2(\mathbb{R})$ and bounded domains $L^2[a, b]$ [6]. Unlike the Fourier basis, a wavelet basis can define a sparse representation of piecewise signals, including transients and singularities [4]. It is hypothesized that these localized properties of the wavelets make them useful test and weighting functions for use in bounded domain Galerkin approaches where problem parameters or field variables feature high gradients or even $C^0$ continuity.

A number of wavelet functions have been developed to accentuate desirable characteristics, such as orthogonality, compact support, symmetry, and number of vanishing moments. References such as Mallat [4], Strang and Nguyen [5], and Daubechies [7] provide excellent detail on the construction and theory behind a wide range of popular wavelets.

The plethora of wavelets from which to select has resulted in extensive literature regarding the wavelet-Galerkin method. For instance, Williams and Amaratunga provide a review of orthogonal wavelet use in engineering [8] and specifically to solutions of linear boundary value problems [9]. Ala et al. [10] use semi-orthogonal spline wavelets to solve the linear electric field integral, while Joly et al. [11] use wavelet packets to define adaptive methods of solving linear time dependent PDE's. The wavelet-Galerkin method has also been used for nonlinear PDE's, for example by Beylkin and Keiser [12] to efficiently capture shock-like responses in equations described by the semigroup approach. Chen et al. [13] provide a seminal derivation of the discrete orthogonal wavelet-Galerkin method on a bounded interval while investigating the solution to Burger's equation, and Lakestani et al. [14, 15] consider Duffing's equation with various boundary

conditions using B-spline wavelets. Restrepo and Leaf [16] look specifically at periodic solutions using orthogonal wavelets. Periodic wavelets are not generally scale-invariant due to “wrap-around” thus special consideration is required when changing the wavelet-scale during the investigation (scale-invariance requires wavelets at any scale to be a pure dilation of the mother-wavelet). Pernot and Lamarque [17] also investigate time-periodic systems, focusing on transient vibrations and stability analysis. Much work has also been done regarding wavelet-Galerkin solutions on bounded domains; the standard wavelet-Galerkin formulation requires the use of wavelets which span outside the domain boundaries thus amendments must be made to the functions to prevent boundary issues. Beylkin [18], Chen et al. [13], and Romine and Peyton [19] all investigate the computation of inner products and other operators of discrete orthogonal wavelets on bounded domains, while Plonka et al. [20] account for domain bounds using the Chebyshev transform.

In all the work referenced above, the differential equations being solved involved constant parameters. Recently Hein and Feklistova [21] consider a wavelet-Galerkin approach to solving free-vibration problems involving beams with smoothly varying properties using Haar wavelets. Solutions to bounded PDE's with spatially dependent variables can be difficult to approximate accurately using standard Fourier-Galerkin methods; if parameters differ at the domain boundaries this will act as a jump discontinuity for the periodic Fourier functions and slow solution convergence. Other non-periodic polynomial basis (i.e. Legendre, Chebyshev) can be used effectively for such problems [2, 22], but if parameter variation is only $C^0$ continuous, a large number of polynomials will be required to capture the rapid gradient changes which again slows convergence. The problem could also be spatially discretized using a finite element approach, but the mesh density must be carefully selected to accurately capture any property variation. In such cases, it is hypothesized that wavelets may allow for a more efficient representation of parameter variation.

In the current investigation, a method for solving one-dimensional partial differential equations with $C^0$ continuous, spatially dependent variables is illustrated using Daubechies wavelets [7]. Derivations for various boundary and loading conditions are given, and solution schemes for both homogeneous and forced equations are detailed.

The paper is broken into three sections. In Section 2 the notation and relevant theory for orthogonal wavelets is reviewed. In Section 3 the discrete orthogonal wavelet-Galerkin method is de-

<footnote>(1) Department of Mechanical Engineering, Rose-Hulman Institute of Technology, Terre Haute, Indiana, USA 47803</footnote>

<footnote>(2) Department of Mechanical Engineering, McGill University, Montreal, Quebec, Canada H3C 1P4</footnote>
tailed using transverse vibrations of a Euler-Bernoulli beam as an illustrative example. Finally, in Section 4 the method is employed for a clamped-pinned beam with a $C^0$ continuous, spatially varying cross-section; the results for both free and forced vibration are compared to a finite-element approximation to assess the effectiveness of the method.

## 2. Notation of Discrete Orthogonal Wavelets

The family of compactly supported orthogonal wavelets were first reported by Daubechies [7], and have since been utilized in a number of areas including signal analysis, image compression, and numerical analysis [5]. As notation for orthogonal wavelets varies in the literature, a brief review is included below for clarity.

The discrete orthogonal wavelet family is defined by a set of L filter coefficients {$p_ℓ : ℓ = 0, 1, ..., L-1$}, where L is an even integer. The fundamental two-scale equations in wavelet theory are defined as

$$ \phi(x) = \sum_{\ell=0}^{L-1} p_{\ell} \phi(2x - \ell) \quad (1) $$

$$ \psi(x) = \sum_{\ell=0}^{L-1} (-1)^{\ell} p_{1-\ell} \phi(2x - \ell) \quad (2) $$

where $\phi(x)$ is the scaling function and $\psi(x)$ is the wavelet function, with fundamental support over the finite intervals $[0, L-1]$ and $[1-L/2, L/2]$, respectively. These equations can be used to determine the value of the scaling and wavelet function at dyadic points $x = n/2^J$, $n = 0, 1, ...$ using the algorithm provided by Chen et al. [13].

The wavelet filter coefficients $p_ℓ$ were derived by Daubechies to produce scaling and wavelet functions with specific properties [7, 13], some of which include:

* the area under the scaling function is unity

$$ \int_{-\infty}^{\infty} \phi(x) dx = 1 \quad (3) $$

* the coefficients sum to two

$$ \sum_{\ell=0}^{L-1} p_{\ell} = 2 \quad (4) $$

* the scaling function and its translates are orthogonal

$$ \int_{-\infty}^{\infty} \phi(x) \phi(x-k) dx = \delta_{0,k}, \quad k \in \mathbb{Z} \quad (5) $$

* the scaling and wavelet functions are orthogonal

$$ \int_{-\infty}^{\infty} \phi(x) \psi(x-k) dx = 0, \quad k \in \mathbb{Z} \quad (6) $$

* the wavelet function has $L/2$ vanishing moments

$$ \int_{-\infty}^{\infty} x^m \psi(x) dx = 0, \quad m = 0, 1, ..., L/2. \quad (7) $$

A basis of $L^2(\mathbb{R})$ can be formed from the linear spans of the scaling and wavelet functions at resolution level J, respectively [6]

$$ \phi_{J,k}(x) = 2^{J/2} \phi(2^J x - k) \quad k \in \mathbb{Z} \quad (8) $$

$$ \psi_{J,k}(x) = 2^{J/2} \psi(2^J x - k) \quad k \in \mathbb{Z}. \quad (9) $$

The span of scaling functions at level J is commonly denoted $V_J$, while the wavelet span is denoted $W_J$. The characteristics of the wavelet basis allows for a multiresolution analysis of a bounded domain on $L^2[a, b]$ by decomposition of the space into a chain of closed subspaces with care taken to properly account for functions at the boundaries [16]

$$ V_0 \subset V_1 \subset V_2 \subset V_3 \subset \cdots \subset L^2[a, b] \quad (10) $$

such that

$$ \lim_{J \to \infty} V_J = L^2[a, b]. \quad (11) $$

Noting that $W_J$ is the orthogonal complement to $V_J$ in $V_{J+1}$

$$ V_{J+1} = V_J \oplus W_J \quad (12) $$

Equation (10) implies

$$ V_{J+1} = V_0 \oplus W_0 \oplus W_1 \oplus \cdots \oplus W_J. \quad (13) $$

## 3. Discrete Orthogonal Wavelet-Galerkin Method

One arguable drawback to orthogonal wavelets is the two-scale functional form of the fundamental equations (1) and (2). The corresponding scaling and wavelet functions are highly non-smooth and fractal in nature: as one increases the resolution of either function the shape does not converge but rather continues to increase in complexity [23]. This makes accurately estimating inner products of scaling or wavelet functions using numerical integration or quadrature error prone. Interestingly, the exact solution to many orthogonal wavelet operators have been derived using the two-scale relation allowing accurate implementation of the wavelet-Galerkin scheme [18, 13, 19, 24]. The equation of motion for a Euler-Bernoulli beam is used below as an illustrative example of such an implementation.

### 3.1 Transverse Vibrations of a Euler-Bernoulli Beam

The equation of motion for transverse, damped vibrations of a clamped-pinned Euler-Bernoulli beam assuming no body forces is [25]

$$ (E(x)I(x)y_{,xx}(x,t))_{,xx} + \rho(x)A(x)\ddot{y}(x,t) + \eta(x)A(x)\dot{y}(x,t) = 0 \quad (14) $$

where $y(0,t) = y_{,x}(0,t) = y(L_b, t) = 0$ and $EIy_{,xx}(L_b, t) = M(t)$. It is assumed the loading and displacements are harmonic

$$ y(x,t) = u(x)e^{i\omega t}. \quad (15) $$

Here $y$ is the transverse beam displacement, $E$ is the elastic modulus, $I$ is the second moment of area, $\rho$ is the mass per unit volume, $A$ is the cross-sectional area, $L_b$ is the length of the beam, $\eta$ is the damping factor, $\omega$ is the frequency of vibration and $i = \sqrt{-1}$.

The beam in this investigation is assumed rectangular in cross-section with a base dimension $b$ and height $h$. As shown in Figure 1, the variation in height with respect to the $x$-direction is only $C^0$ continuous, due to its piece-wise construction. This implies both $I = 1/12bh^3$ and $A = bh$ are also $C^0$ continuous functions of $x$; all other parameters are assumed to remain constant. Accounting for the assumed harmonic response and the spatial dependency of $u(x)$, $I(x)$ and $A(x)$ results in the following equation of motion

$$ I(x)u_{,xxxx}(x) + 2I_{,x}(x)u_{,xxx}(x) + I_{,xx}(x)u_{,xx}(x) - \beta(\omega)A(x)u(x) = 0 \quad (16) $$
<figure><img src="image_2.png" /><figcaption>Figure 1. Thickness variation h(x) of the beam</figcaption></figure>

where

$$ \beta(\omega) = \frac{\rho\omega^2 - i\omega\eta}{E}. \qquad (17) $$

The total length of the beam is $L_b = 2m$. Let the domain $x \in [0, L_b]$ be discretized by $2^J$ points, hence the distance between discrete points $\Delta x$ is given by $\Delta x = L_b/2^J$. Furthermore, it is convenient in the following sections to consider a domain of $[0, 1]$, thus a change of variable is introduced

$$ x = x/L_b. \qquad (18) $$

## 3.2 Galerkin Approximation

In accordance with the Galerkin method [2], an approximate trial solution for $u(x)$ is introduced using scaling functions at level $J$ as the test functions

$$ u(x) \approx \sum_{k=2-L}^{2^J-1} u_k \phi_{J,k}(x) = 2^{J/2} \sum_{k=2-L}^{2^J-1} u_k \phi(2^J x - k) \quad (19) $$

where $J > 0$. The limits of summation are selected such that an approximation of $u(x)$ can be made for the whole domain $x \in [0, 1]$. When $x = 0$, the summation in Equation (19) involves scaling functions evaluated at integer values $-k$. The scaling function $\phi(z)$ has only $L-2$ non-zero integer values corresponding to $z = 1, 2, \dots, L-2$ [13]. Therefore, the only non-zero $\phi_{J,k}(x)$ correspond to $k = 2-L, 3-L, \dots, -1$. Similarly when $x = 1$ the only non-zero values correspond to $k = 2^J - L + 2, 2^J - L + 3, \dots, 2^J - 1$. Therefore the required summation spans the values $k = 2-L, 3-L, \dots, 2^J - 1$.

It is important to stress that Equation (19) is a reduced basis approximation of the function $u(x)$ because no wavelets are included. According to Equation (13) the exact function can be represented as

$$ u(x) = \sum_{k} u_k \phi_{J,k}(x) + \sum_{j=J}^{N} \sum_{\ell} w_{j,\ell} \psi_{j,\ell}(x) \qquad (20) $$

where $N$ is the maximum resolution given by the sampling rate of the function $u(x)$ [4]. As discrete orthogonal wavelet functions are analogous to high-pass filters [5], neglecting the wavelet terms in Equation (20) is akin to removing portions of the high frequency content of $u(x)$, resulting in a smoothed approximation of the true function. As shown in Equation (10) however, the scaling function approximation approaches $L^2[a, b]$ as the level $J$ is increased. Therefore to improve the accuracy of the approximation the wavelet functions at level $J$ can be included, or more simply the level of the scaling functions can be increased to $J+1$ as detailed in Equation (12). Restricting the test functions to scaling functions at a single level simplifies the algorithms while also providing a reduced basis which is adaptable to the desired level of accuracy.

Another consideration before continuing with the Galerkin approach is how to deal with the spatially dependent terms $I(x)$ and $A(x)$ from Equation (16). Since the height function $h(x)$ is piecewise linear, directly substituting this dependence into

the governing equation would result in different equations to be satisfied depending on the $x$ location. If one considers a more general case where the height dependence is not piecewise linear, the Galerkin formulation would require finding the integral of a triple product of two scaling functions and an arbitrary function of $x$ (e.g. the integral of the weighting function, test function, and the $I(x)$ function). There are a limited set of functions where this triple product integral can be found exactly [13]; in general this integral would have to be approximated using numerical integration which is known to be error prone [23].

The novelty of the proposed approach lies in the decomposition of the spatially varying functions, $A(x)$ and $I(x)$, using discrete orthogonal scaling functions such that

$$ A(x) \approx \sum_{j=2-L}^{2^J-1} A_j \phi_{J,j}(x) = 2^{J/2} \sum_{j=2-L}^{2^J-1} A_j \phi(2^J x - j) \quad (21) $$

$$ I(x) \approx \sum_{j=2-L}^{2^J-1} I_j \phi_{J,j}(x) = 2^{J/2} \sum_{j=2-L}^{2^J-1} I_j \phi(2^J x - j). \quad (22) $$

The scaling function approximations provide relatively accurate representations of the piecewise functions due to the vanishing moment properties in Equation (7) allowing perfect reproduction of polynomials up to order $(L/2 - 1)$ [13]. This approach results in a single governing equation, regardless of $x$ location, and the resulting triple product integrals can be found exactly (see below for full explanation).

The coefficients $A_j$ and $I_j$ are calculated using the inner product [4]

$$ A_j = \int_0^1 A(x) \phi_{J,j}(x) \, dx \quad \text{and} \quad I_j = \int_0^1 I(x) \phi_{J,j}(x) \, dx. \quad (23) $$

This inner product can be performed using the discrete wavelet transform (DWT) or fast wavelet transform (FWT) [4, 5]. Attention must be paid to the calculation of the coefficients near the boundary $x=0$. Standard decomposition requires the scaling functions to extend into the negative $x$ domain by $L-2$ points, thus a decision on how to represent this boundary must be made. A number of options exist including scaling function extrapolation [9], symmetric extension [4], or the use of boundary scaling functions [18, 13, 19]. In the current investigation it is assumed a symmetric extension is appropriate (i.e. $A(-x) = A(x)$); due to the beam geometry in this problem all the methods listed would actually provide the same results.

The Galerkin method defines the residual $R$ as the substitution of Equations (19), (21) and (22) into Equation (16)

$$ R = \sum_{j=2-L}^{2^J-1} \sum_{k=2-L}^{2^J-1} \left[ I_j u_k (\phi_{J,j}^{(0)}(x) \phi_{J,k}^{(4)}(x) + 2\phi_{J,j}^{(1)}(x) \phi_{J,k}^{(3)}(x) + \phi_{J,j}^{(2)}(x) \phi_{J,k}^{(2)}(x)) - \beta(\omega) A_j u_k (\phi_{J,j}^{(0)}(x) \phi_{J,k}^{(0)}(x)) \right] \quad (24) $$

where the superscript (n) refers to the order of differentiation. The n-th derivative of the scaling function $\phi_{J,k}^{(n)}(x)$ is defined as

$$ \phi_{J,k}^{(0)}(x) = \phi_{J,k}(x) \qquad (25) $$

$$ \phi_{J,k}^{(n)}(x) = \frac{d^n \phi_{J,k}}{dx^n}(x) = \frac{d}{dx} \phi_{J,k}^{(n-1)}(x) \qquad (26) $$

By amalgamating Equations (8), (18) and (26), and accounting for the vanishing moment condition in Equation (7), it is possible to write

$$ \phi_{J,k}^{(n)}(x) = \frac{2^{J/2}}{\Delta x^n} \phi^{(n)}(2^J x - k), \quad n = 0, 1, \dots, L/2 - 1 \quad (27) $$
It is important to note the limit of $n = L/2 - 1$ due to the vanishing moment condition; in the current investigation where a fourth-order derivative is required this condition implies a minimum scaling function order of $L = 10$. A method for determining values of the scaling function and its derivatives at dyadic points $x = k/2^J$ is given by Chen et al. [13].

Continuing with the Galerkin method, scaling functions of level $J$ are selected as the weighting functions and the inner product of the residual and the weighting functions is set to zero [2]

$$ \int_{0}^{1} R \phi_{J,\ell}(x) \, dx = 0, \quad \ell = 2-L, 3-L, \dots, 2^J - 1 \quad (28) $$

which, using Equation (24), results in

$$ \sum_{j=2-L}^{2^J-1} \sum_{k=2-L}^{2^J-1} \left[ I_j u_k (a_{\ell,j,k} + 2b_{\ell,j,k} + c_{\ell,j,k}) - \beta(\omega) A_j u_k (d_{\ell,j,k}) \right] = 0 \quad (29) $$

for $\ell = 2-L, 3-L, \dots, 2^J - 1$, where

$$ a_{\ell,j,k} = \int_{0}^{1} \phi_{J,\ell}(x) \phi_{J,j}^{(0)}(x) \phi_{J,k}^{(4)}(x) \, dx \quad (30) $$

$$ b_{\ell,j,k} = \int_{0}^{1} \phi_{J,\ell}(x) \phi_{J,j}^{(1)}(x) \phi_{J,k}^{(3)}(x) \, dx \quad (31) $$

$$ c_{\ell,j,k} = \int_{0}^{1} \phi_{J,\ell}(x) \phi_{J,j}^{(2)}(x) \phi_{J,k}^{(2)}(x) \, dx \quad (32) $$

$$ d_{\ell,j,k} = \int_{0}^{1} \phi_{J,\ell}(x) \phi_{J,j}^{(0)}(x) \phi_{J,k}^{(0)}(x) \, dx. \quad (33) $$

Chen et al. refer to triple product integrals of this form as three-term connection coefficients [13] and show that the exact solution to these proper integrals can be found. Reference [26] suggests improved algorithms for calculation of these connection coefficients. The standard notation for the three-term connection coefficients is

$$ \Omega_{j,k}^{m,n}(x) = \int_{0}^{x} \phi(y) \phi^{(m)}(y-j) \phi^{(n)}(y-k) \, dy \quad (34) $$

which requires Equations (30) to (33) to be written as

$$ a_{\ell,j,k} = \frac{2^{J/2}}{\Delta x^4} (\Omega_{j-\ell,k-\ell}^{0,4}(2^J - \ell) - \Omega_{j-\ell,k-\ell}^{0,4}(-\ell)) \quad (35) $$

$$ b_{\ell,j,k} = \frac{2^{J/2}}{\Delta x^4} (\Omega_{j-\ell,k-\ell}^{1,3}(2^J - \ell) - \Omega_{j-\ell,k-\ell}^{1,3}(-\ell)) \quad (36) $$

$$ c_{\ell,j,k} = \frac{2^{J/2}}{\Delta x^4} (\Omega_{j-\ell,k-\ell}^{2,2}(2^J - \ell) - \Omega_{j-\ell,k-\ell}^{2,2}(-\ell)) \quad (37) $$

$$ d_{\ell,j,k} = 2^{J/2} (\Omega_{j-\ell,k-\ell}^{0,0}(2^J - \ell) - \Omega_{j-\ell,k-\ell}^{0,0}(-\ell)). \quad (38) $$

Finally, Equation (29) can be written in matrix form as

$$ [\mathbf{G} - \beta(\omega)\mathbf{H}] \mathbf{u}_k = \mathbf{0} \quad (39) $$

where

$$ G_{\ell,k} = \frac{2^{J/2}}{\Delta x^4} \sum_{j=2-L}^{2^J-1} I_j \left[ \begin{aligned} & (\Omega_{j-\ell,k-\ell}^{0,4}(2^J - \ell) - \Omega_{j-\ell,k-\ell}^{0,4}(-\ell)) \\ & + 2(\Omega_{j-\ell,k-\ell}^{1,3}(2^J - \ell) - \Omega_{j-\ell,k-\ell}^{1,3}(-\ell)) \\ & + (\Omega_{j-\ell,k-\ell}^{2,2}(2^J - \ell) - \Omega_{j-\ell,k-\ell}^{2,2}(-\ell)) \end{aligned} \right] \quad (40) $$

and

$$ H_{\ell,k} = 2^{J/2} \sum_{j=2-L}^{2^J-1} A_j (\Omega_{j-\ell,k-\ell}^{0,0}(2^J - \ell) - \Omega_{j-\ell,k-\ell}^{0,0}(-\ell)) $$

for $k, \ell = 2-L, 3-L, \dots, 2^J - 1$.

It should be noted that the three-term connection coefficients are parameter independent, thus they can be stored in matrix form and reused efficiently for parametric analyses or even different PDE's requiring the same scaling function derivative combinations. This makes the wavelet-Galerkin method an efficient tool overall.

### 3.3 Imposing Boundary Conditions

The final step before a solution can be found involves prescription of the appropriate boundary conditions on the basis functions. As shown in Figure 1, the beam in this investigation is clamped at $x = 0$ and pinned at $x = 1$, which is equivalent to $u(0) = 0$, $u_{,x}(0) = 0$, $u(1) = 0$, and $E\hat{I}u_{,xx}(1) = 0$ respectively, where $\hat{I}$ is the moment of inertia at $x = 1$. Consider the clamped boundary conditions at $x = 0$; using Equation (19) these can be written in the scaling function domain as

$$ u(0) = 2^{J/2} \sum_{k=2-L}^{2^J-1} u_k \phi(-k) = 0 \quad (41) $$

$$ u_{,x}(0) = \frac{2^{J/2}}{\Delta x} \sum_{k=2-L}^{2^J-1} u_k \phi^{(1)}(-k) = 0 \quad (42) $$

Noting that $\phi(x)$ and $\phi^{(1)}(x)$ have non-zero values for integers $x = 1, 2, \dots, L - 2$ only [13], these equations simplify to

$$ \sum_{k=2-L}^{-1} u_k \phi(-k) = 0 \quad (43) $$

$$ \sum_{k=2-L}^{-1} u_k \phi^{(1)}(-k) = 0. \quad (44) $$

Now consider the pinned boundary condition at $x = 1$. If the end is free of any applied moments, the BC's can be written

$$ u(1) = 2^{J/2} \sum_{k=2-L}^{2^J-1} u_k \phi(2^J - k) = 0 \quad (45) $$

$$ E \hat{I} u_{,xx}(1) = \frac{2^{J/2}}{\Delta x^2} E \hat{I} \sum_{k=2-L}^{2^J-1} u_k \phi^{(2)}(2^J - k) = 0. \quad (46) $$

Accounting again for the non-zero values of $\phi(x)$ and $\phi^{(2)}(x)$, these can be simplified to

$$ \sum_{k=2^J-L+2}^{2^J-1} u_k \phi(2^J - k) = 0 \quad (47) $$

$$ \sum_{k=2^J-L+2}^{2^J-1} u_k \phi^{(2)}(2^J - k) = 0. \quad (48) $$

If instead there is a moment $M(t) = Me^{i\omega t}$ applied at $x = 1$, as depicted in Figure 1, the respective BC becomes

$$ \sum_{k=2^J-L+2}^{2^J-1} u_k \phi^{(2)}(2^J - k) = \frac{M \Delta x^2}{2^{J/2} E \hat{I}}. \quad (49) $$
### 3.4 Free Vibrations Solution Scheme

The natural frequencies of the beam can be found by assessing the homogeneous solution to Equation (39). This requires setting the forcing moment $M = 0$ and the damping factor $\eta = 0$.

Equation (39) is a generalized eigenvalue problem in $\beta(\omega)$: the eigenvalues can be used to determine the natural frequencies of the system $\omega_n$ using Equation (17); the corresponding eigenvectors can be used to compute the modeshapes using Equation (19). Thus the four boundary conditions (43), (44), (47) and (48) must be imposed in a manner so as not to disturb the generalized eigenvalue form.

One efficient approach is to constrain four participation factors using the boundary conditions. To accomplish this it is convenient to first reorder Equation (39) so that

$$ (\hat{\mathbf{G}} - \beta(\omega)\hat{\mathbf{H}}) \begin{Bmatrix} \mathbf{u}_k^m \\ \mathbf{u}_k^s \end{Bmatrix} = \mathbf{0} \qquad (50) $$

where $\mathbf{u}_k^m$ are the master scaling function participation factors and $\mathbf{u}_k^s$ are the slave factors which will be constrained. $\hat{\mathbf{G}}$ and $\hat{\mathbf{H}}$ are the reordered $\mathbf{G}$ and $\mathbf{H}$ matrices, respectively. The size of $\mathbf{u}_k^m$ is $(2^J + L - 6) \times 1$ and $\mathbf{u}_k^s$ is $4 \times 1$ to account for the four boundary conditions.

Boundary conditions (43), (44), (47) and (48) can be used to form the constraint equations

$$ \begin{Bmatrix} \mathbf{u}_k^m \\ \mathbf{u}_k^s \end{Bmatrix} = \hat{\mathbf{C}}\mathbf{u}_k^m \qquad (51) $$

where $\hat{\mathbf{C}}^\top = [\mathbf{I}^\top, \mathbf{C}^\top]$ is the constraint matrix. $\mathbf{I}$ is an identity matrix of order $(2^J + L - 6)$ and $\mathbf{C}$ is a $4 \times (2^J + L - 6)$ matrix derived from Equations (43), (44), (47) and (48). The slave participation factors should be selected such that the condition number of $\mathbf{C}$ is minimized; for the current boundary conditions the slave factors are $u_{-1}$, $u_{-2}$, $u_{2^J - 2}$, $u_{2^J - 1}$. A sample calculation of $\mathbf{C}$ is given in the Appendix.

Substituting Equation (51) into (50) and pre-multiplying by the transpose of $\hat{\mathbf{C}}$ results in the constrained system of equations

$$ \hat{\mathbf{C}}^\top(\hat{\mathbf{G}} - \beta(\omega)\hat{\mathbf{H}})\hat{\mathbf{C}}\mathbf{u}_k^m = 0 \qquad (52) $$

thus the boundary conditions have been imposed without disturbing the generalized eigenvalue form of the equation. The eigenvector coefficients for the slave participation factors can be determined using Equation (51).

### 3.5 Forced Vibrations Solution Scheme

As an example of how loading can be incorporated in the model, a harmonic moment loading condition at the right end of the beam is considered as shown in Figure 1. It should be noted that more general loading conditions, such as distributed loads, can also be included using the wavelet-Galerkin method.

Similar to the free response case, the system of equations (39) is constrained using the appropriate boundary condition (43), (44), (47) and (49). Due to the applied moment in Equation (49) the constraint equation must be amended to

$$ \begin{Bmatrix} \mathbf{u}_k^m \\ \mathbf{u}_k^s \end{Bmatrix} = \hat{\mathbf{C}}\mathbf{u}_k^m + \hat{\mathbf{F}}. \qquad (53) $$

$\hat{\mathbf{F}}$ is of size $(2^J + L - 2) \times 1$ and all entries are zero except for the row associated with Equation (49) in $\mathbf{C}$ where the entry is $(M\Delta x^2)/(2^{J/2}E\hat{\mathbf{I}})$. The resulting constrained system of equations is

$$ \hat{\mathbf{C}}^\top(\hat{\mathbf{G}} - \beta(\omega)\hat{\mathbf{H}})\hat{\mathbf{C}}\mathbf{u}_k^m = -\hat{\mathbf{C}}^\top(\hat{\mathbf{G}} - \beta(\omega)\hat{\mathbf{H}})\hat{\mathbf{F}} \qquad (54) $$

This system can be solved using standard matrix inversion for participation factors $\mathbf{u}_k^m$. The solution $y(x, t)$ for the beam can be found by solving for $u(x)$ using Equation (19), finding $y(x, t) = u(x)e^{i\omega t}$, and finally transforming back to the $x$ domain using $x = L_b x$.

## 4. Analysis and Discussion

The discrete orthogonal wavelet-Galerkin method is employed to approximate the transverse vibrations of an Euler-Bernoulli beam as described by Equation (14). The beam, whose properties are given in Table 1, varies in thickness as shown in Figure 1. The 10 coefficient Daubechies scaling function (DB10) is used for all investigations; the level $J$ is adjusted to determine the effect on prediction accuracy.

<table><caption>Table 1. Beam properties</caption><thead><tr><th>Parameter</th><th>Free Vibration</th><th>Forced Vibration</th></tr></thead><tbody><tr><td><i>E</i> (GPa)</td><td>70</td><td>70</td></tr><tr><td><i>ρ</i> (kg/m<sup>3</sup>)</td><td>2700</td><td>2700</td></tr><tr><td><i>b</i> (m)</td><td>0.1</td><td>0.1</td></tr><tr><td><i>L</i><sub>b</sub> (m)</td><td>2.0</td><td>2.0</td></tr><tr><td><i>η</i> (N s/m<sup>4</sup>)</td><td>0</td><td>5×10<sup>4</sup></td></tr><tr><td><i>M</i> (N m)</td><td>0</td><td>10</td></tr></tbody></table>

A finite element model (FEM) comprised of linear beam elements [27] is also constructed using the beam parameters listed in Table 1 to act as a comparison model.

### 4.1 Free Vibration

Figure 2 shows the convergence of the first four natural frequencies using the finite element model (FEM) and the wavelet-Galerkin model (W-G); the number of unknowns shown equate to $J = 5, 6, 7, 8$.

<figure><img src="image_3.png" /><figcaption>Figure 2. Convergence of the first four natural frequencies using the finite element method (FEM) and wavelet-Galerkin method (W-G)</figcaption></figure>

As shown in the figure, the wavelet-Galerkin approach converges to a solution more rapidly than the finite element method. This is attributed to the scaling functions' ability to better represent the piecewise variation in cross-sectional area and second moment of inertia than can the linear finite elements. The FEM discretization results in a step-wise approximation of the thickness variation thus a fine mesh is required to converge to a solution. The vanishing moment properties of the Daubechies scaling functions allow perfect reconstruction of the linear variations in area and cubic variation in moment of inertia with only localized error at the piecewise joints. These joint locations have
discontinuous slope and require higher level wavelets to perfectly reconstruct, hence the convergence as *J* is increased.

For this example, computational speed is essentially proportional to the number of unknowns. The element matrices can be pre-computed and stored for the FEA method in a similar way to the connection coefficients for the W-G method. The matrices of interest are relatively small thus the time required for the pre- and post-multiplication in Equation (52) and (54) is negligible. The computational time required for the FEA or W-G method is thus governed by the solution to the eigenvalue or linear algebra problem, which is the same for both methods. It is currently unclear how the wavelet-Galerkin method would scale up to much larger problems; this is an active area of research.

The associated normalized displacement mode shapes predicted by FEM and W-G at $J = 7$ are shown in Figure 3. The predicted mode shapes from both models overlay each other (< 0.1% relative error rms), attesting to the validity of the proposed wavelet-Galerkin method.

<figure><img src="image_4.png" /><figcaption>**Figure 3.** First four mode shapes predicted by the finite element method (FEM) and wavelet-Galerkin method (W-G) using 128 degrees of freedom ($J = 7$)</figcaption></figure>

## 4.2 Forced Vibration

The frequency response function for the displacement of the midspan of the beam ($x = 1$ m) due to a 10 N m harmonic moment at $x = 2$ m is provided in Figure 4. Included in the figure is the FEM prediction for 128 degrees of freedom and the W-G prediction at levels $J = 4, 5, 6$. The peaks of the curve correspond to the natural frequencies shown in Figure 2; these peaks are bounded due to the damping included for this portion of the investigation (see Table 1).

The plot shows how the wavelet-Galerkin solution converges rapidly as the level $J$ is increased. As in the free vibration case, the values of the natural frequencies are over predicted at smaller $J$, most noticeably for higher modes. This is attributed to an inadequate number of scaling functions necessary to precisely reconstruct the increasingly rapid changes in slope of the higher modes, as well as the thickness variation of the beam.

## 5. Conclusion

The discrete orthogonal wavelet-Galerkin method is introduced as an efficient method of solving partial differential equation (PDE) with spatially varying parameters on a bounded domain. The Daubechies scaling functions provide a concise but adaptable set of basis functions: increasing the level of the scaling functions results in a more encompassing subspace of $L^2[a, b]$ and thus

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Frequency response function at midspan ($x = 1$ m) under a 10 N m harmonic moment acting at $x = 2$ m</figcaption></figure>

a more accurate approximation to the solution of the PDE. The spatial dependence of the parameters can be described using the same scaling functions, which allows the PDE to be written as a set of linear equations composed of three-term connection coefficients. The connection coefficient matrices are parameter and even problem independent thus they can be stored and reused for efficient parametric analyses of various PDE's.

An example problem of transverse vibrations of a pinned-clamped Euler-Bernoulli beam with piecewise varying thickness is used to describe the method. Solution schemes for both homogeneous and forced PDE's are developed for a number of boundary conditions. The natural frequencies computed using the wavelet-Galerkin method are shown to converge more rapidly than do those predicted by a linear finite element model with equivalent degrees of freedom. The associated modeshape are also shown to match between the two models. Finally, the frequency response function predicted by the wavelet-Galerkin method is shown to converge as the level of scaling function is increased as expected. These results suggest the wavelet-Galerkin method is well suited to approximating the solution of PDE's with spatially dependent parameters.

## Acknowledgements

This work was supported in part by the Tomlinson Postdoctoral Fellowship.

## Appendix

Here a sample constraint matrix $C$ from Section 3.4 is constructed using the boundary conditions given in Equations (43), (44), (47) and (48) for $L = 6$ and $J = 4$. Substituting these parameters into the respective boundary condition equations results in

$$ \begin{bmatrix} \phi^{(0)}(4) & \phi^{(0)}(3) \\ \phi^{(1)}(4) & \phi^{(1)}(3) \end{bmatrix} \begin{Bmatrix} u_{-4} \\ u_{-3} \end{Bmatrix} + \begin{bmatrix} \phi^{(0)}(2) & \phi^{(0)}(1) \\ \phi^{(1)}(2) & \phi^{(1)}(1) \end{bmatrix} \begin{Bmatrix} u_{-2} \\ u_{-1} \end{Bmatrix} = \begin{Bmatrix} 0 \\ 0 \end{Bmatrix} \quad (55) $$

and

$$ \begin{bmatrix} \phi^{(0)}(4) & \phi^{(0)}(3) \\ \phi^{(2)}(4) & \phi^{(2)}(3) \end{bmatrix} \begin{Bmatrix} u_{12} \\ u_{13} \end{Bmatrix} + \begin{bmatrix} \phi^{(0)}(2) & \phi^{(0)}(1) \\ \phi^{(2)}(2) & \phi^{(2)}(1) \end{bmatrix} \begin{Bmatrix} u_{14} \\ u_{15} \end{Bmatrix} = \begin{Bmatrix} 0 \\ 0 \end{Bmatrix}. \quad (56) $$

Solving for the slave participation factors gives

$$ \begin{Bmatrix} u_{-2} \\ u_{-1} \end{Bmatrix} = \begin{bmatrix} a_{11} & a_{12} \\ a_{21} & a_{22} \end{bmatrix} \begin{Bmatrix} u_{-4} \\ u_{-3} \end{Bmatrix} \quad \text{and} \quad \begin{Bmatrix} u_{14} \\ u_{15} \end{Bmatrix} = \begin{bmatrix} b_{11} & b_{12} \\ b_{21} & b_{22} \end{bmatrix} \begin{Bmatrix} u_{12} \\ u_{13} \end{Bmatrix} \quad (57) $$
where

$$a_{11} = (\phi^{(0)}(1)\phi^{(1)}(4) - \phi^{(0)}(4)\phi^{(1)}(1))/\Delta_a \quad (58)$$

$$a_{12} = (\phi^{(0)}(1)\phi^{(1)}(3) - \phi^{(0)}(3)\phi^{(1)}(1))/\Delta_a \quad (59)$$

$$a_{21} = (\phi^{(0)}(4)\phi^{(1)}(2) - \phi^{(0)}(2)\phi^{(1)}(4))/\Delta_a \quad (60)$$

$$a_{22} = (\phi^{(0)}(3)\phi^{(1)}(2) - \phi^{(0)}(2)\phi^{(1)}(3))/\Delta_a \quad (61)$$

$$b_{11} = (\phi^{(0)}(1)\phi^{(2)}(4) - \phi^{(0)}(4)\phi^{(2)}(1))/\Delta_b \quad (62)$$

$$b_{12} = (\phi^{(0)}(1)\phi^{(2)}(3) - \phi^{(0)}(3)\phi^{(2)}(1))/\Delta_b \quad (63)$$

$$b_{21} = (\phi^{(0)}(4)\phi^{(2)}(2) - \phi^{(0)}(2)\phi^{(2)}(4))/\Delta_b \quad (64)$$

$$b_{22} = (\phi^{(0)}(3)\phi^{(2)}(2) - \phi^{(0)}(2)\phi^{(2)}(3))/\Delta_b \quad (65)$$

and

$$\Delta_a = \phi^{(0)}(2)\phi^{(1)}(1) - \phi^{(0)}(1)\phi^{(1)}(2) \quad (66)$$

$$\Delta_b = \phi^{(0)}(2)\phi^{(2)}(1) - \phi^{(0)}(1)\phi^{(2)}(2). \quad (67)$$

Combining these gives

$$\left\{ \begin{matrix} u_{-2} \\ u_{-1} \\ u_{14} \\ u_{15} \end{matrix} \right\} = \left[ \begin{matrix} a_{11} & a_{12} & 0 & \cdots & 0 \\ a_{21} & a_{22} & 0 & \cdots & 0 \\ 0 & \cdots & 0 & b_{11} & b_{12} \\ 0 & \cdots & 0 & b_{21} & b_{22} \end{matrix} \right] \left\{ \begin{matrix} u_{-4} \\ u_{-3} \\ u_0 \\ u_1 \\ \vdots \\ u_{11} \\ u_{12} \\ u_{13} \end{matrix} \right\} \quad (68)$$

or $\mathbf{u}_k^s = \mathbf{Cu}_k^m$.

## References

[1] Mark Gockenbach. *Partial Differential Equations: Analytical and Numerical Methods*. 2nd edition. Society for Industrial and Applied Mathematics, 2011. ISBN: 9780898719352.

[2] Clive Fletcher. *Computational Galerkin Methods*. 1st edition. Springer-Verlag Berlin Heidelberg, 1984. ISBN: 9783642859519. DOI: 10.1007/978-3-642-85949-6.

[3] Leon Lapidus and George Pinder. *Numerical Solution of Partial Differential Equations in Science and Engineering*. John Wiley and Sons, 1999. ISBN: 9780471359449.

[4] Stéphane Mallat. *A wavelet tour of signal processing*. 3rd edition. Academic Press, 2009. ISBN: 9780123743701.

[5] Gilbert Strang and Truong Nguyen. *Wavelets and Filter Banks*. 2nd edition. Wellesley-Cambridge Press, 1996. ISBN: 9780961408879.

[6] Albert Cohen, Ingrid Daubechies, and Pierre Vial. "Wavelets on the interval and fast wavelet transforms". *Applied and Computational Harmonic Analysis* 1.1 (1993), pages 54–81. [hal-01311753].

[7] Ingrid Daubechies. *Ten lectures on wavelets*. Society for Industrial and Applied Mathematics, 1992. DOI: 10.1137/1.9781611970104.

[8] John Williams and Kevin Amaratunga. "Introduction to wavelets in engineering". *International Journal for Numerical Methods in Engineering* 37.14 (1994), pages 2365–2388. [hal-01311772].

[9] Kevin Amaratunga and John Williams. "Wavelet-Galerkin solution of boundary value problems". *Archives of Computational Methods in Engineering* 4.3 (1997), pages 243–285. [hal-01311725].

[10] Guido Ala, Maria Luisa Di Silvestre, Elisa Francomano, and Adele Tortorici. "An advanced numerical model in solving thin-wire integral equations by using semi-orthogonal compactly supported spline wavelets". *IEEE Transactions on Electromagnetic Compatibility* 45.2 (2003), pages 218–228. [hal-01321546].

[11] Pascal Joly, Yvon Maday, and Valérie Perrier. "Towards a method for solving partial differential equations by using wavelet packet bases". *Computer Methods in Applied Mechanics and Engineering* 116.1-4 (1994), pages 301–307. [hal-01321537].

[12] Gregory Beylkin and James Keiser. "On the adaptive numerical solution of nonlinear partial differential equations in wavelet bases". *journal of computational physics* 132.2 (1997). [hal-01322927].

[13] Ming-Quayer Chen, Chyi Hwang, and Yen-Ping Shih. "The computation of wavelet-Galerkin approximation on a bounded interval". *International Journal for Numerical Methods in Engineering* 39 (1996), pages 2921–2944. [hal-01330583].

[14] Mehrard Lakestani, Mohsen Razzaghi, and Mehdi Dehghan. "Numerical solution of the controlled Duffing oscillator by semi-orthogonal spline wavelets". *Physica Scripta* 74.3 (2006), pages 362–366. [hal-01330591].

[15] Mehrard Lakestani and Mehdi Dehghan. "The solution of a second-order nonlinear differential equation with Neumann boundary conditions using semi-orthogonal B-spline wavelets". *International Journal of Computer Mathematics* 83 (2006), pages 685–694. [hal-01330587].

[16] Juan Mario Restrepo and Gary Leaf. "Inner product computations using periodized Daubechies wavelets". *International Journal for Numerical Methods in Engineering* 40.19 (1997), pages 3557–3578. [hal-01323007].

[17] Stéphane Pernot and Claude-Henri Lamarque. "A wavelet-Galerkin procedure to investigate time-periodic systems: transient vibration and stability analysis". *Journal of Sound and Vibration* 245.5 (2001), pages 845–875. [hal-00814563].

[18] Gregory Beylkin. "On the representation of operators in bases of compactly supported wavelets". *SIAM Journal of Numerical Analysis* 29.6 (1992), pages 1716–1740. JSTOR: 2158264.

[19] C. Romine and Barry Peyton. "Computing connection coefficients of compactly supported wavelets on bounded intervals". *US Department of Energy* (1997). [hal-01321539].

[20] Gerlind Plonka, Kathi Selig, and Manfred Tasche. "On the construction of wavelets on a bounded interval". *Advances in Computational Mathematics* 4 (1995), pages 357–388. [hal-01322918].

[21] Helle Hein and Lidia Feklistova. "Free vibrations of non-uniform and axially functionally graded beams using Haar wavelets". *Engineering Structures* 33 (2011), pages 3696–3701. [hal-01321538].

[22] John Boyd. *Chebyshev and Fourier spectral methods*. 2nd edition. Dover Publications, 2001. ISBN: 9780486411835.

[23] Stephen Welstead. *Fractal and wavelet image compression techniques*. SPIE - The International Society for Optics and Photonics, 1999. ISBN: 9780819435033.

[24] Andrew Latto, Howard Resnikoff, and Eric Tenenbaum. "The evaluation of connection coefficients of compactly supported wavelets". *Proceedings of the French-USA Workshop on Wavelets and Turbulence, Princeton University*. New York: Springer-Verlag, 1992.

[25] James Doyle. *Wave Propagation in Structures*. New York: Springer, 1997. ISBN: 9781461273042. DOI: 10.1007/978-1-4612-1832-6.

[26] Tonghua Zhang, Yu-Chu Tian, Moses Tadé, and Johan Utomo. "Comments on 'The computation of wavelet-Galerkin approximation on a bounded interval'". *International Journal for Numerical Methods in Engineering* 72.2 (2007), pages 244–251. [hal-01330597].

[27] Jagmohan Humar. *Dynamics of structures*. 3rd edition. CRC Press/Balkema, 2012. ISBN: 9780415620864.