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

HAL
open science

# The wavelet-Galerkin method for solving PDE's with spatially dependent variables

Simon Jones, Mathias Legrand

► To cite this version:

Simon Jones, Mathias Legrand. The wavelet-Galerkin method for solving PDE's with spatially dependent variables. 19th International Congress on Sound and Vibration (ICSV19), Jul 2012, Vilnius, Lithuania. hal-00719744

HAL Id: hal-00719744

https://hal.science/hal-00719744

Submitted on 20 Jul 2012

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.
# THE WAVELET-GALERKIN METHOD FOR SOLVING PDE'S WITH SPATIALLY DEPENDENT VARIABLES

Simon Jones and Mathias Legrand

Department of Mechanical Engineering, McGill University
817 Sherbrooke Street West, Montreal, Quebec, Canada
e-mail: simon.jones2@mcgill.ca

Discrete orthogonal wavelets are a family of functions with compact support which form a basis on a bounded domain. Use of these wavelet families as Galerkin trial functions for solving partial differential equations (PDE's) has been a topic of interest for the last decade, though research has primarily focused on equations with constant parameters. In the current paper the wavelet-Galerkin method is extended to allow spatial variation of equation parameters. A representative example from the field of vibration illustrates the method: compression waves in a bar with varying elastic modulus. The computed natural frequencies and modeshapes are compared to finite element solutions and show excellent correspondence. The wavelet-Galerkin method is also shown to be an efficient and convenient solution method as the majority of the calculations are performed *a priori* and can be stored for use in solving future PDE's. This efficiency is displayed by performing a stochastic analysis of elastic modulus variation to determine the effect on the frequency response function.

## 1. Introduction

The wavelet-Galerkin method has emerged as an accurate and efficient means of approximating the solution of partial differential equations (PDE's) [1]. Wavelets are well localized, oscillatory functions which provide a basis of $L^2(\mathbb{R})$ and can be modified to a basis of $L^2[a,b]$ where $[a,b]$ is a bounded domain [2]. These localized characteristics of discrete, orthogonal wavelets allow sparse representation of piecewise signals, including transients and singularities, making them useful functions for use in the Galerkin approach when non-smooth or non-periodic solutions are predicted [3].

Discrete, orthogonal wavelets have been used by a number of investigators in a Galerkin approach to solving differential equations. Williams and Amaratunga provide a review of orthogonal wavelet use in engineering [4] and specifically to solutions of linear boundary value problems [5]. Beylkin and Keiser [6] attempt to efficiently capture shock-like responses in nonlinear equations described by the semigroup approach. Restrepo and Leaf [7] look specifically at periodic solutions using orthogonal wavelets, and Pernot and Lamarque [8] investigate transient vibrations and stability analysis. Beylkin [9], Chen et al. [10], and Romine and Peyton [11] all investigate the computation of inner products and other operators of orthogonal wavelets on bounded domains; the exact solution to these operators was paramount in the development of the discrete, orthogonal wavelet-Galerkin method.

In the current investigation, a method for solving one-dimensional partial differential equations with spatially dependent variables is introduced using discrete, orthogonal wavelets, known as Daubechies wavelets [1]. A specific example of axial vibrations of a bar with varying elastic modulus is used to illustrate the approach. Details on how to account for various boundary and loading conditions are given, and the efficiency of the method is shown by performing a stochastic analysis of the variation in elastic modulus.
The paper is broken into four sections. In Section 2 the notation and relevant theory for orthogonal wavelets is reviewed. In Section 3 the orthogonal wavelet-Galerkin method for PDE's with spatially dependent variable is developed using axial vibrations of a bar. Section 4 provides details of the numerical model, and Section 5 presents results validating the method for constant properties and displays efficiency using a stochastic analysis to vary the elastic modulus.

## 2. Notation of orthogonal wavelets

The 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. [10]. The wavelet decomposition is analogous to a filter bank [3]; the scaling functions act as low-pass filters while the wavelets act as high-pass filters.

The wavelet filter coefficients $p_ℓ$ were derived by Daubechies to produce scaling and wavelet functions with specific properties [1, 10], 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) $$

Chen et al. [10] show that a subspace of $L^2[a, b]$ can be formed from the linear spans of the scaling functions at resolution level $J$

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

The span of scaling functions at level $J$ is commonly denoted $V_J$. The characteristics of the scaling function levels allows for a multiresolution analysis of a finite domain on $L^2[a, b]$ by decomposition of the space into a chain of closed subspaces [10]

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

such that

$$ \bigcap_{J=0,1,...} V_J = \{0\} \quad \overline{\bigcup}_{J=0,1,...} V_J = L^2[a,b]. \quad (8) $$
### 3. Orthogonal wavelet-Galerkin method

Discrete, orthogonal scaling functions defined by Eq. (1) are generally highly non-smooth and fractal in nature: as one increases the resolution the shape does not converge but rather continues to increase in complexity. This makes accurately estimating inner products which involve scaling functions difficult when using numerical integration or quadrature. Interestingly, the exact solution to many orthogonal wavelet operators have been derived using the two-scale properties [9, 10, 11] allowing relatively simple implementation of the wavelet-Galerkin scheme. The equation of motion for axial vibrations of a bar is used below as an instructional example (see Fig. 1).

<figure><img src="image_2.png" /><figcaption>Figure 1. Bar used in investigation</figcaption></figure>

#### 3.1 Axial vibrations of a bar

The equation of motion for the elastodynamics of a bar assuming no body forces is [12]

$$ \frac{\partial}{\partial x} \left( EA \frac{\partial y}{\partial x} \right) = \rho A \frac{\partial^2 y}{\partial t^2} + \eta A \frac{\partial y}{\partial t} \qquad (9) $$

where $y(x,t)$ is the axial displacement of the bar, $E$ is the elastic modulus, $\rho$ is the mass per unit volume, $A$ is the cross-sectional area, and $\eta$ is the damping factor. It is assumed the elastic modulus of the bar varies with $x$ (i.e. $E = E(x)$); all other parameters are assumed to remain constant. Furthermore, all loading in this investigation is harmonic, thus the displacement is assumed to be of the form

$$ y(x,t) = Y(x)e^{i\omega t} \qquad (10) $$

where $\omega$ is the frequency of vibration and $i = \sqrt{-1}$. Accounting for the spatial dependency and the harmonic response results in the following equation of motion

$$ E \frac{d^2 Y}{dx^2} + \frac{dE}{dx} \frac{dY}{dx} + \beta Y(x) = 0 \qquad (11) $$

where

$$ \beta = \rho \omega^2 - i \omega \eta. \qquad (12) $$

Without loss of generality, assume the total length of the bar is 1m. Let the domain $x \in [0, 1]$ be discretized by $2^J$ points, hence the distance between discrete points $\Delta x = 1/2^J$.

#### 3.2 Galerkin approximation

In accordance with the Galerkin method [13], a trial solution $u(x) \simeq Y(x)$ is introduced using scaling functions at level $J$ as the test functions

$$ u(x) = \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 J > 0. \qquad (13) $$

At this point in the derivation the variation in the elastic modulus of the bar is known but remains general. The inner product of the scaling function with $E(x)$ is eventually required in the Galerkin
approach. The solution to such an inner product is not generally solvable algebraically, and the fractal nature of the scaling functions make numerical integration prone to numerical error. To allow spatially varying parameters the function $E(x)$ is written in the scaling function domain

$$E(x) \approx \sum_{j=2-L}^{2^J-1} E_j \phi_{j,j}(x) = 2^{J/2} \sum_{j=2-L}^{2^J-1} E_j \phi(2^J x - j) \quad (14)$$

where the coefficients $E_j$ are calculated using the inner product [2]

$$E_j = \int_0^1 E(x) \phi_{j,j}(x) dx \quad (15)$$

The Galerkin residual $R$ is found by substituting Eqs. (13) and (14) into Eq. (11)

$$R = \sum_{j=2-L}^{2^J-1} \sum_{k=2-L}^{2^J-1} E_j u_k \left( \phi_{j,j}^{(0)}(x) \phi_{j,k}^{(2)}(x) + \phi_{j,j}^{(1)}(x) \phi_{j,k}^{(1)}(x) \right) + \beta \sum_{k=2-L}^{2^J-1} u_k \phi_{j,k}^{(0)}(x) \quad (16)$$

where the superscript (n) refers to the order of differentiation. The nth derivative of the scaling function $\phi_{j,j}^{(n)}(x)$ is

$$\phi_{j,j}^{(n)}(x) = 2^{nJ+J/2} \phi^{(n)}(2^J x - j), \quad n = 0, 1, \dots, L/2 - 1. \quad (17)$$

It is important to note the limit of $n = L/2 - 1$ due to the vanishing moment condition [1]; in the current investigation where a second-order derivative is required this condition implies a minimum wavelet order of $L = 6$. 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. [10].

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

$$\int_0^1 R \phi_{J,l}(x) dx = 0 \quad \text{for } l = 2-L, 3-L, \dots, 2^J-1 \quad (18)$$

which using Eq. (16), results in

$$\sum_{j=2-L}^{2^J-1} \sum_{k=2-L}^{2^J-1} E_j u_k (a_{l,j,k} + b_{l,j,k}) + \beta \sum_{k=2-L}^{2^J-1} u_k c_{l,k} = 0 \quad (19)$$

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

$$a_{l,j,k} = \int_0^1 \phi_{J,l}(x) \phi_{j,j}^{(0)}(x) \phi_{j,k}^{(2)}(x) dx \quad (20)$$

$$b_{l,j,k} = \int_0^1 \phi_{J,l}(x) \phi_{j,j}^{(1)}(x) \phi_{j,k}^{(1)}(x) dx \quad (21)$$

$$c_{l,k} = \int_0^1 \phi_{J,l}(x) \phi_{j,j}^{(0)}(x) dx. \quad (22)$$

Chen et al. [10] refer to integrals of form (20) and (21) as three-term connection coefficients, and integrals of form (22) as two-term connection coefficients. Reference [14] suggests improved algorithms for exact calculation of these connection coefficients. The standard notation for the two-term and three-term connection coefficients are

$$\Gamma_k^n(x) = \int_0^x \phi(y) \phi^{(n)}(y-k) dy \quad (23)$$

$$\Omega_{j,k}^{m,n}(x) = \int_0^x \phi(y) \phi^{(m)}(y-j) \phi^{(n)}(y-k) dy. \quad (24)$$
Using this notation allows Eq. (19) to be written in matrix form as

$$[\mathbf{G} + \beta \mathbf{H}] \mathbf{u}_k = \mathbf{0} \qquad (25)$$

where

$$G_{\ell,k} = 2^{5J/2} \sum_{j=2-L}^{2^J-1} E_j \left[ \left( \Omega_{j-\ell,k-\ell}^{0,2} (2^J - \ell) + \Omega_{j-\ell,k-\ell}^{0,2} \ell \right) + \left( \Omega_{j-\ell,k-\ell}^{1,1} (2^J - \ell) + \Omega_{j-\ell,k-\ell}^{1,1} \ell \right) \right] \\ H_{\ell,k} = 2^{J/2} (\Gamma_{k-\ell}^0 (2^J - \ell) + \Gamma_{k-\ell}^0 \ell)$$

and $\Gamma_k^n(x)$ and $\Omega_{j,k}^{m,n}(x)$ can be computed as detailed in Refs. [10, 14]. It should be noted that these connection coefficient matrices are parameter independent, thus they can be stored and reused efficiently for parametric analyses.

### 3.3 Imposing boundary conditions

As shown in Fig. 1, the bar under investigation is clamped at $x = 0$ and subject to a harmonic force at $x = 1$, which is equivalent to

$$u(0) = 0 \qquad (26)$$

$$\frac{du}{dx}(1) = \frac{F}{AE(1)}. \qquad (27)$$

Using Eq. (13) and making use of the compact support of the scaling functions, the boundary conditions can be written in the scaling function domain as

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

$$\sum_{k=2^J-L+2}^{2^J-1} u_k \phi^{(1)}(2^J - k) = \frac{F}{2^{3J/2}AE(1)} \qquad (29)$$

respectively. Replacement of the first and last rows of Eqn. (25) with Eqns. (28) and (29) fully constrains the system, allowing solution via standard matrix inversion techniques.

## 4. Definition of the model

The bar in Fig. 1 has the following properties: length $L_{bar} = 1$ m; density $\rho = 2700$ kg/m³; cross-sectional area $A = 4 \times 10^{-4}$ m²; damping factor $\eta = 1 \times 10^5$ N s/m⁴. The harmonic end load is varied within the range 0 – 10 kHz to determine the point-response function at the end of the bar.

To highlight the efficiency of the wavelet-Galerkin approach, a Monte Carlo simulation of 100 realizations is performed to investigate the model’s sensitivity to stochastic variation of elastic modulus. It is assumed the elastic modulus has the following statistical properties: mean value $\bar{E} = 70$ GPa, standard deviation $\tilde{E} = 7$ GPa, and scale of fluctuation $\delta_E = 5$ m. These statistical properties are used in conjunction with a K-L expansion to produce the realizations of elastic modulus variation along the length of the bar using a modified exponential covariance function; the probability density function (PDF) is described using a log-normal distribution with zero mean, unit variance, and cut-off at negative three standard deviations from mean to ensure positive elastic modulus. Two example realizations of the elastic modulus variation are shown in Fig. 2. Full details on the implementation of the K-L expansion for non-Gaussian PDF’s can be found in Ref. [15].

To act as a comparison, a finite-element (FE) model is constructed using truss elements. This model is used to quantify the accuracy and efficiency of the wavelet-Galerkin approach relative to a known numerical approach.
<figure><img src="image_3.png" /><figcaption>Figure 2. Two example realizations of the variation of elastic modulus with respect to position along the bar</figcaption></figure>

## 5. Results and discussion

The analytic solution for the natural frequencies of an undamped uniform bar undergoing axial vibration is

$$ \omega = \frac{2n + 1}{4L_{\text{bar}}} \sqrt{\frac{E}{\rho}} \quad \text{for } n = 0, 1, 2, \dots \qquad (30) $$

where the frequency $\omega$ is in Hz. For the bar under investigation the first four natural frequencies are: 1272.94 Hz; 3818.81 Hz; 6364.69 Hz; 8910.56 Hz. Both the wavelet-Galerkin (WG) and finite-element (FE) approaches are employed using a harmonic analysis from 0-10 kHz to determine the predicted natural frequencies for an undamped uniform bar (i.e. average elastic modulus used, no stochastic variation). The results are presented in Fig. 3 as percent errors relative to the analytic solutions.

<figure><img src="image_4.png" /><figcaption>Figure 3. Convergence of finite-element (FE) and wavelet-Galerkin (WG) approaches relative to the analytic solution for the first four natural frequencies of a uniform bar under axial vibration</figcaption></figure>

A single FE curve is presented where the percent error relative to the analytic solution is approximately 0.1% or less for the first four modes; a total of 64 degrees of freedom (dof) are required for this convergence tolerance. Three curves representing the WG solutions for different scales are also included: $J = 4$ (16 dof); $J = 5$ (32 dof); $J = 6$ (64 dof). The relative error for the WG approach with $J = 5$ is lower than the FE model with 64 dof; when 64 unknowns are employed for the WG approach the relative error is an order of magnitude smaller than the equivalent FE model. This suggests the wavelet-Galerkin approach requires fewer unknowns to produce equivalent prediction
accuracy compared to the linear finite element approach for this model. As both approaches involve solving a system of linear equations with order equal to the number of unknowns, this equates to faster computational times when performing parametric or Monte Carlo simulations.

To quantify this claim, 100 realizations of the bar with stochastically varying elastic modulus are computed using both the FE and WG methods. The results from the wavelet-Galerkin method ($J = 5$) are presented in Fig. 4. The results from the FE model (64 dof) are virtually equivalent; for example the results for the two realizations shown in Fig. 2 using the FE and WG methods have relative frequency errors of less than 0.15% for the first four modes and absolute amplitude differences of less than 0.1 dB. The results from the stochastic analysis themselves are not the goal of this example, but rather to show the wavelet-Galerkin method can simulate varying realizations of spatially dependent variables accurately and efficiently.

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Monte Carlo results for 100 realizations of stochastically varied elastic modulus: grey lines are individual realizations; black line is mean response</figcaption></figure>

The reduced basis size of the WG method (i.e. 32 unknowns for WG vs. 64 unknowns for FE) is reflected in the computational times for the 100 realization Monte Carlo simulation: 59.8 s for WG; 109.2 s for FE. This relatively low run time for the WG method is not only due to the reduced number of unknowns, but also reflects the problem independence of the connection-coefficients. Since the connection coefficients can be stored and reused for all the simulations, deriving the equations of motion is as computationally inexpensive as using predetermined finite element blocks.

# 6. Conclusions

The wavelet-Galerkin method is introduced as an accurate and efficient means of estimating the solution to partial differential equations with spatially varying parameters. Due to the fractal nature of the scaling functions used as the Galerkin test functions, the spatially varying parameters are transformed into the scaling function domain to allow exact calculations of the required inner products using the connection-coefficient method. Results for an illustrative example of axial vibrations of a bar with varying elastic modulus show the wavelet-Galerkin method requires fewer unknowns than the linear finite element method to predict natural frequencies of equivalent or better numerical accuracy. This reduced size coupled with the problem independence of the connection-coefficients result in efficient calculation of results when performing parametric analyses or Monte Carlo simulations.

# Bibliography

¹ Daubechies, I. *Ten lectures on wavelets*, SIAM, Philadelphia, PA (1992).
² Mallat, S. *A wavelet tour of signal processing*, Academic Press, Burlington, MA (2009).

³ Strang, G., Nguyen, T. *Wavelets and Filter Banks*, Wellesley-Cambridge Press (1996).

⁴ Williams, J. and Amaratunga, K. Introduction to wavelets in engineering, *International Journal for Numerical Methods in Engineering*, **37**, 2365-2388, (1994).

⁵ Amaratunga, K. and Williams, J., Wavelet-Galerkin solution of boundary value problems, *Archives of Computational Methods in Engineering*, **4**(3), 243-285, (1997).

⁶ Beylkin, G. and Keiser, J., On the adaptive numerical solution of nonlinear partial differential equations in wavelet bases, *University of Colorado Technical Report*, (1996).

⁷ Restrepo, J. and Leaf, G., Inner product computations using periodized Daubechies wavelets, *University of California, Los Angeles, Technical Report*, (1993).

⁸ Pernot, S. and Lamarque, C.-H., A wavelet-Galerkin procedure to investigate time-periodic systems: transient vibration and stability analysis, *Journal of Sound and Vibration*, **245**, 845-875, (2001).

⁹ Beylkin, G., On the representation of operators in bases of compactly supported wavelets, *SIAM Journal of Numerical Analysis*, **6**(6), 1716-1740, (1992).

¹⁰ Chen, M., Hwang, C. and Shih, Y., The computation of wavelet-Galerkin approximation on a bounded interval, *International Journal for Numerical Methods in Engineering*, **39**, 2921-2944, (1996).

¹¹ Romine, C. and Peyton, B., Computing connection coefficients of compactly supported wavelets on bounded intervals, *US Department of Energy*, (1997).

¹² Doyle, J., *Wave Propagation in Structures*, Springer-Verlag, New York Inc, (1997).

¹³ Fletcher, C., *Computational Galerkin Methods*, Springer-Verlag, New York, (1984).

¹⁴ Zhang, T., Tian, Y., Tadé, M. and Utomo, J., Comments on ‘The computation of wavelet-Galerkin approximation on a bounded interval’, *International Journal for Numerical Methods in Engineering*, **72**, 244-251, (2007).

¹⁵ Jones, S. and Hunt, H., Predicting surface vibration from underground railways through inhomogeneous soil, *Journal of Sound and Vibration*, **331**(9), 2055-2069, (2012).