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

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.


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 .R/ and bounded domains L 2 OEa; 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 semiorthogonal 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 shocklike 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 "wraparound" 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-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.

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 fp`W`D 0; 1; : : : ; L 1g, where L is an even integer. The fundamental two-scale equations in wavelet theory are defined as where .x/ is the scaling function and .x/ is the wavelet function, with fundamental support over the finite intervals OE0; L 1 and OE1 L=2; L=2, respectively. These equations can be used to determine the value of the scaling and wavelet function at dyadic points x D n=2 J , n D 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 the coefficients sum to two the scaling function and its translates are orthogonal the scaling and wavelet functions are orthogonal the wavelet function has L=2 vanishing moments A basis of L 2 .R/ can be formed from the linear spans of the scaling and wavelet functions at resolution level J , respectively [6] 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 OEa; b by decomposition of the space into a chain of closed subspaces with care taken to properly account for functions at the boundaries [16] Noting that W J is the orthogonal complement to Equation (10) implies

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.

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.
It is assumed the loading and displacements are harmonic Here y is the transverse beam displacement, E is the elastic modulus, I is the second moment of area, is the mass per unit volume, A is the cross-sectional area, L b is the length of the beam, Á is the damping factor, ! is the frequency of vibration and i D p 1. The beam in this investigation is assumed rectangular in crosssection 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 D 1=12bh 3 and A D 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 The total length of the beam is L b D 2m. Let the domain x 2 OE0; L b be discretized by 2 J points, hence the distance between discrete points x is given by is convenient in the following sections to consider a domain of OE0; 1, thus a change of variable is introduced

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 where J > 0. The limits of summation are selected such that an approximation of u.x/ can be made for the whole domain x 2 OE0; 1. When x D 0, the summation in Equation (19) involves scaling functions evaluated at integer values k. The scaling function .z/ has only L 2 non-zero integer values corresponding to z D 1; 2; : : : ; L 2 [13]. Therefore, the only non-zero J;k .x/ correspond to k D 2 L; 3 L; : : : ; 1.
Similarly when x D 1 the only non-zero values correspond to k D 2 J L C 2; 2 J L C 3; : : : ; 2 J 1. Therefore the required summation spans the values k D 2 L; 3 L; : : : ; 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 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 OEa; 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 C 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 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 D 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 D 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/ D 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) where the superscript (n) refers to the order of differentiation. The nth derivative of the scaling function .n/ J;k .x/ is defined as By amalgamating Equations (8), (18) and (26), and accounting for the vanishing moment condition in Equation (7), it is possible to write It is important to note the limit of n D L=2 1 due to the vanishing moment condition; in the current investigation where a fourthorder derivative is required this condition implies a minimum scaling function order of L D 10. A method for determining values of the scaling function and its derivatives at dyadic points x D 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] Z 1 0 R J;`. x/ dx D 0;`D 2 L; 3 L; : : : ; 2 J 1 (28) which, using Equation (24), results in for`D 2 L; 3 L; : : : Chen et al. refer to triple product integrals of this form as threeterm 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 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.

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 D 0 and pinned at x D 1, which is equivalent to u.0/ D 0, u ;x .0/ D 0, u.1/ D 0, and E O I u ;xx .1/ D 0 respectively, where O I is the moment of inertia at x D 1. Consider the clamped boundary conditions at x D 0; using Equation (19) these can be written in the scaling function domain as Noting that .x/ and .1/ .x/ have non-zero values for integers x D 1; 2; : : : ; L 2 only [13], these equations simplify to Now consider the pinned boundary condition at x D 1. If the end is free of any applied moments, the BC's can be written Accounting again for the non-zero values of .x/ and .2/ .x/, these can be simplified to If instead there is a moment M.t / D M e i !t applied at x D 1, as depicted in Figure 1, the respective BC becomes

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 D 0 and the damping factor Á D 0. Equation (39) is a generalized eigenvalue problem inˇ.!/: the eigenvalues can be used to determine the natural frequencies of the system ! 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 where u m k are the master scaling function participation factors and u s k are the slave factors which will be constrained. O where O C > D OEI > ; C > is the constraint matrix. I is an identity matrix of order .2 J C L 6/ and C is a 4 .2 J C L 6/ matrix derived from Equations (43), (44), (47) and (48). The slave participation factors should be selected such that the condition number of 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 C is given in the Appendix.
Substituting Equation (51) into (50) and pre-multiplying by the transpose of O C results in the constrained system of equations 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).

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), This system can be solved using standard matrix inversion for participation factors u m k . The solution y.x; t / for the beam can be found by solving for u.x/ using Equation (19), finding y.x; t / D u.x/e i !t , and finally transforming back to the x domain using x D L b x.

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 [27] is also constructed using the beam parameters listed in Table 1 to act as a comparison model. 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 D 5; 6; 7; 8.  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.

Free Vibration
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 preand 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 D 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 purposed wavelet-Galerkin method.

Forced Vibration
The frequency response function for the displacement of the midspan of the beam (x D 1 m) due to a 10 N m harmonic moment at x D 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 D 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.

Conclusion
The discrete orthogonal wavelet-Galerkin method is introduced as an efficient method of solving partial differential equation  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 pinnedclamped 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. Here a sample constraint matrix C from Section 3.4 is constructed using the boundary conditions given in Equations (43) or u s k D Cu m k .