# Direct solution of differential equations using the Wavelet-Galerkin method

Rodrigo Burgos, Raul E Silva, Marco Santos

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

## To cite this version:

Rodrigo Burgos, Raul E Silva, Marco Santos. Direct solution of differential equations using the Wavelet-Galerkin method. Mecánica Computacional, 2010, 29, pp.4573-4584. hal-01333839

HAL Id: hal-01333839

https://hal.science/hal-01333839

Submitted on 19 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
DIRECT SOLUTION OF DIFFERENTIAL EQUATIONS
USING THE WAVELET-GALERKIN METHOD

Rodrigo B. Burgos<sup>a</sup>, Raul R. e Silva<sup>a</sup>, Marco A. C. Santos<sup>b</sup>

**Keywords:** Wavelets, Daubechies, Interpolets, Wavelet-Galerkin Method.

**Abstract.** The use of compactly supported wavelet functions has become increasingly popular in the development of numerical solutions for differential equations, especially for problems with local high gradient. Daubechies wavelets have been successfully used as base functions in several schemes like the Wavelet-Galerkin Method, due to their compact support, orthogonality, and multi-resolution properties. Another advantage of wavelet-based methods is the fact that the calculation of the inner products of wavelet basis functions and their derivatives can be made by solving a linear system of equations, thus avoiding the problem of approximating the integral by some numerical method. These inner products were defined as connection coefficients and they are employed in the calculation of stiffness, mass and geometry matrices. In this work, the Wavelet-Galerkin Method has been adapted for the direct solution of differential equations in a meshless formulation. This approach enables the use of a multiresolution analysis. Several examples based on differential equations for beams and plates were studied successfully.

<footnote><sup>a</sup>DEC, Pontifícia Universidade Católica do Rio de Janeiro, Rua Marques de São Vicente, 225, Rio de Janeiro, RJ, Brasil, rburgos@esp.puc-rio.br, http://www.civ.puc-rio.br</footnote>

<footnote><sup>b</sup>LAGEMAR, Universidade Federal Fluminense, Av. Gen. Milton Tavares de Souza, s/n, Niterói, RJ, Brasil, cetale@igeo.uff.br, http://www.igeo.uff.br</footnote>
# 1 INTRODUCTION

The use of wavelet-based numerical schemes has become increasingly popular in the last two decades. Wavelets have several properties that are especially useful for representing solutions of differential equations (DE's), such as orthogonality, compact support and exact representation of polynomials of a certain degree. Their capability of representing data at different levels of resolution allow the efficient and stable calculation of functions with high gradients or singularities, which would require a dense mesh or higher order elements in a Finite Element analysis (Qian and Weiss, 1992).

A complete basis of wavelets can be generated through dilation and translation of a mother scaling function. Although many applications use only the wavelet filter coefficients of the multiresolution analysis, there are some which explicitly require the values of the basis functions and their derivatives, such as the Wavelet Finite Element Method (WFEM) (Ma et al., 2003).

Compactly supported wavelets have a finite number of derivatives which can be highly oscillatory. This makes the numerical evaluation of integrals of their inner products difficult and unstable. Those integrals are called connection coefficients and they appear naturally when applying a numerical method for the solution of a DE. Due to some properties of wavelet functions, these coefficients can be obtained by solving an eigenvalue problem using filter coefficients.

The most commonly used wavelet family is the one developed by Ingrid Daubechies (1988). All the mathematical foundation for the wavelet theory was formulated for Daubechies wavelets and then extended to other families.

Working with dyadically refined grids, Deslauriers and Dubuc (1989) obtained a new family of wavelets with interpolating properties, later called Interpolets. Unlike Daubechies wavelets, Interpolets are symmetric, which is especially interesting in numerical analysis. The use of Interpolets instead of Daubechies wavelets considerably improves the method's accuracy (Burgos et al., 2008).

The use of wavelets as interpolating functions in numerical schemes such as the Galerkin Method holds some promise due to their multiresolution properties. Accuracy can be improved by increasing either the level of resolution or the order of the wavelet used. For detection of singularities, the increase in the level of resolution seems to work better.

Two examples were used for validating the proposed method in a one-dimensional scheme. First, a beam with a concentrated load was used to test the method's ability to capture singularities. In a second example, the critical loads and buckling modes of a doubly clamped beam were calculated at different levels of resolution. In order to evaluate the method's capability to solve two-dimensional problems, it was then applied for a thin plate with excellent results.

# 2 WAVELET THEORY

Multiresolution analysis using orthogonal, compactly supported wavelets has been successfully applied in numerical simulation. Wavelets are localized in space, which allows local variations of the problem to be analyzed at various levels of resolution. In the following expression, known as the two-scale relation, $a_k$ are the filter coefficients of the wavelet scale function $\varphi$ and $N$ is the wavelet order.

$$ \varphi(x) = \sum_{k=0}^{N-1} a_k \varphi(2x-k) = \sum_{k=0}^{N-1} a_k \varphi_k(2x). \quad (1) $$
In general, there are no analytical expressions for wavelet functions, which can be obtained using iterative procedures like Eq. (1). In order to comply with the requirements of orthogonality and compact support, wavelets present, in general, an irregular fractal-like shape. Figure 1 shows Daubechies wavelets at four different orders. It's evident that the higher the order of the wavelet, the greater its smoothness.

<figure><img src="image_3.png" /><figcaption>Figure 1: Daubechies Wavelets: (a) N = 4, (b) N = 6, (c) N = 8, (d) N = 10</figcaption></figure>

## 2.1 Wavelet properties

The set of properties summarized in Eq. (2) is valid for Daubechies wavelets but can be adapted to other wavelet families, such as Delaurier-Dubuc Interpolets. Some of these properties, like compact support and unit integral, are required for the use of the wavelet family in numerical methods. Others, like orthogonality, are desirable but not extremely necessary.

$$ 
\begin{align} 
\text{supp}(\varphi) &= [0 \ N-1], \notag \\ 
\int_{-\infty}^{+\infty} \varphi(x)dx &= 1, \notag \\ 
\int_{-\infty}^{+\infty} \varphi(x-i)\varphi(x-j)dx &= \delta_{i,j}, \tag{2} \\ 
x^m &= \sum_{k=-\infty}^{+\infty} c_k^m \varphi(x-k), \quad m \le N/2-1. \notag 
\end{align} 
$$
The last expression in Eq. (2) derives from the vanishing moments property, which states that a set of shifted and scaled Daubechies wavelets of order N is capable of representing exactly an $N/2 - 1$ degree polynomial.

## 2.2 Wavelet Derivatives

In the process of solving a DE using numerical methods, derivatives of the basis functions tend to appear. As there is no analytical expression for wavelets, derivatives are obtained in dyadic grid points and the refinement of the solution depends on the level of resolution needed (Lin et. al, 2005). The scale relation can be differentiated d times, generating the following expression:

$$ \varphi^{(d)}(x) = 2^d \sum_{i=0}^{N-1} a_i \varphi^{(d)}(2x-i). \quad (3) $$

Applying Eq. (3) to integer points results in the following system of equations shown in matrix form. In Eq. (4), A represents the filter coefficients matrix, I is the identity matrix and $\Gamma^{(d)}$ is the vector containing derivative values at integer points of the grid.

$$ \begin{gathered} (2^d \mathbf{A} - \mathbf{I})\Gamma^{(d)} = 0, \\ \mathbf{A} = [a_{2i-k}]_{0 \le i, k \le N-1}. \end{gathered} \quad (4) $$

Eq. (4) is an eigenvalue problem which, for unique solution, has to be normalized using the so-called moment equation, derived from the wavelet property of exact polynomial representation. This equation is given by Latto et al. (1992) and provides a relation between derivative values at integer points.

$$ \begin{gathered} d! = \sum_{i=0}^{N-1} M_i^d \varphi^{(d)}(x-i), \\ M_i^j = \frac{1}{2^{j+1}-2} \sum_{k=0}^{j} \binom{j}{k} i^{j-k} \sum_{l=0}^{k-1} \binom{k}{l} M_0^l \left( \sum_{i=0}^{N-1} a_i i^{k-l} \right). \end{gathered} \quad (5) $$

Once the derivative is obtained at integer values, the scale relation can be applied for any $x = k/2^j$.

$$ \varphi^{(d)}\left(\frac{k}{2^j}\right) = 2^d \sum_{i=0}^{N-1} a_i \varphi^{(d)}\left(\frac{k}{2^{j-1}} - i\right), \quad k = 1, 3, 5, \dots, N-2. \quad (6) $$

## 2.3 Connection Coefficients

Assuming that a function f is approximated by a series of interpolating scale functions, the following may be written:

$$ f(\xi) = \sum_k \alpha_k \varphi_k(\xi). \quad (7) $$

The process of solving a differential equation by some numerical method requires the calculation of the inner products of the basis functions and their derivatives. These inner products are defined as connection coefficients ($\Gamma$):

$$ \Gamma_{i,j}^{d_1, d_2} = \int \varphi^{(d_1)}(\xi - i) \varphi^{(d_2)}(\xi - j) d\xi. \quad (8) $$
The values for the limits of the integral in Eq. (8) depend on which method is used to impose boundary conditions. In this work, the limits are given by $[0 \ 2^m]$, where $m$ is the wavelet level of resolution. This method allows the use of Lagrange multipliers to deal with boundary conditions, similarly to what is usually done in a meshless scheme (Nguyen et al., 2008). Connection coefficients at level $m$ can be obtained through the calculation at level 0, thus avoiding its recalculation while increasing the level of resolution. Wavelet dilation and translation properties allow the calculation of connection coefficients within the interval $[0 \ 1]$ to be summarized by the solution of an eigenvalue problem based only on filter coefficients (Zhou & Zhang, 1998).

$$ \left( \mathbf{P} - \frac{1}{2^{d_1+d_2-1}} \mathbf{I} \right) \Gamma^{d_1,d_2} = 0, \qquad (9) \\ P_{i,j:k,l} = a_{k-2i}a_{l-2j} + a_{k-2i+1}a_{l-2j+1}, $$

Since Eq. (9) leads to an infinite number of solutions, there is the need for a normalization rule that provides a unique eigenvector. This unique solution comes with the inclusion of an adapted version of the moment equation mentioned before (Latto et al., 1992).

$$ \sum_{i} \sum_{j} M_{i}^{k} M_{j}^{k} \Gamma^{d_{1},d_{2}}_{i,j} = \frac{(k!)^{2}}{(k-d_{1})!(k-d_{2})!(2k-d_{1}-d_{2}+1)}. \qquad (10) $$

## 2.4 Delauriers-Dubuc Interpolets

The basic characteristics of interpolating wavelets require that the mother scaling function satisfies the following condition (Shi et al., 1999):

$$ \varphi(k) = \delta_{0,k} = \begin{cases} 1, & k = 0 \\ 0, & k \neq 0 \end{cases}, \quad k \in \mathbb{Z}. \qquad (11) $$

The filter coefficients for Delauriers-Dubuc Interpolets can be obtained by an autocorrelation of the Daubechies filter coefficients. Interpolets satisfy the same requirements as other wavelets, specially the two-scale relation, which is fundamental for their use as interpolating functions in numerical methods. Figure 2 shows the Interpolet IN6 (autocorrelation of DB6, Daubechies wavelet of order 6). Its symmetry and interpolating properties are evident. There is only one integer abscissa which evaluates to a non-zero value.

<figure><img src="image_4.png" /><figcaption>Figure 2: Interpolet IN6 scaling function with its full support</figcaption></figure>
All expressions used for the calculation of derivatives, connection coefficients and moments of Daubechies wavelets can be applied to Interpolets. Of course, due to the correlation, the support [0 N-1] in the expressions for Daubechies becomes [1-N N-1] for Interpolets.

## 3 WAVELET-GALERKIN METHOD

The numerical solution of differential equations is one of the possible applications of the wavelet theory. The Wavelet-Galerkin Method (WGM) results from the use of wavelets as interpolating functions in a traditional Galerkin scheme. In the following sections, the WGM will be applied to solve typical DE's for structures like beams and plates.

### 3.1 Beam subjected to axial load

The equation of a beam subjected to an axial load is given by:

$$ \frac{\partial^4 w}{\partial x^4} + \frac{P}{EI} \frac{\partial^2 w}{\partial x^2} = 0. \qquad (12) $$

Stiffness and geometry matrices can be obtained by substituting the displacement w by a series of interpolating functions. Adimensional coordinates ξ within the interval [0 1] are used in wavelet space, leading to the subsequent expressions:

$$ \begin{aligned} \bar{k}_{i,j} &= \int_{0}^{1} \varphi_{i}''(\xi) \varphi_{j}''(\xi) d\xi = \Gamma_{i,j}^{2,2}, \\ \bar{g}_{i,j} &= \int_{0}^{1} \varphi_{i}'(\xi) \varphi_{j}'(\xi) d\xi = \Gamma_{i,j}^{1,1}. \end{aligned} \qquad (13) $$

As done in the Finite Elements Method (FEM), the critical loads and buckling modes can be obtained by solving an eigenvalue problem of the form:

$$ \left\{ \begin{bmatrix} \mathbf{\Gamma}^{2,2} & -\mathbf{G}^T \\ -\mathbf{G}^T & \mathbf{0} \end{bmatrix} - \frac{P}{EI} \begin{bmatrix} \mathbf{\Gamma}^{1,1} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} \end{bmatrix} \right\} \left\{ \begin{matrix} \mathbf{\alpha} \\ \mathbf{\lambda} \end{matrix} \right\} = \left\{ \begin{matrix} \mathbf{0} \\ \mathbf{0} \end{matrix} \right\}, \qquad (14) $$

where the matrix G is associated with boundary conditions and λ is a vector of Lagrange multipliers. The main difference in relation to the FEM is that the unknowns in vector α are the interpolating coefficients of the basis functions instead of nodal displacements. In fact, there is no need to establish nodal coordinates.

### 3.2 Thin Plate

The bending of a thin plate with thickness t is modeled by the following DE:

$$ D \left( \frac{\partial^4 w}{\partial x^4} + 2 \frac{\partial^4 w}{\partial x^2 \partial y^2} + \frac{\partial^4 w}{\partial y^4} \right) = q(x, y), \quad D = \frac{Et^3}{12(1 - v^2)}. \qquad (15) $$

Displacement w(x,y) is modeled using bi-dimensional wavelets, which are products between one-dimensional wavelets:

$$ w(x, y) = \sum_{i} \sum_{j} d_{ij} \phi(x - i) \phi(y - j). \qquad (16) $$
As done in a tradicional Galerkin approach, Eq. (16) is substituted in the DE and integrated, leading to a system of equations which contains the wavelets' connection coefficients (Chen et al., 2004).

$$kd = f,$$

$$k = D \left[ \Gamma_{[0,2^m]}^{00} \otimes \Gamma_{[0,2^m]}^{22} + \nu (\Gamma_{[0,2^m]}^{20} \otimes \Gamma_{[0,2^m]}^{02} + \Gamma_{[0,2^m]}^{02} \otimes \Gamma_{[0,2^m]}^{20}) + \Gamma_{[0,2^m]}^{22} \otimes \Gamma_{[0,2^m]}^{00} + 2(1-\nu^2)\Gamma_{[0,2^m]}^{11} \otimes \Gamma_{[0,2^m]}^{11} \right], \quad (17)$$

$$f = \int_{0}^{1} q(x) \Phi^{T} dx.$$

The symbol $\otimes$ indicates Kronecker product. The system is solved using the stiffness matrix provided by Eq. (17) and imposing essential boundary conditions with Lagrange multipliers as done before for one-dimensional DE's.

## 4 EXAMPLES

Figure 3 shows a simple example of a beam subjected to a concentrated load at its midpoint. This example was formulated in order to verify the ability of the wavelet method to deal with singularities, since the load generates a discontinuity in the shear force diagram.

<figure><img src="image_5.png" /><figcaption>Figure 3: Beam with concentrated load</figcaption></figure>

This example is easily solved by dividing the beam in two elements and applying the load as a nodal force. In this work, since degrees of freedom don't have a fixed position, the load is transformed into the wavelet space:

$$q(\xi) = P\delta\left(\xi - \frac{1}{2}\right) \rightarrow \int_{0}^{1} q(\xi)\varphi(\xi - i)d\xi = P\varphi\left(\frac{1}{2} - i\right). \quad (18)$$

The example was solved using the IN8 Interpolet at different levels of resolution and the results for bending moment and shear force diagrams are shown in Figures 4 and 5.

It is clear that higher levels of resolution are necessary in order to capture the singularity that occurs where the load is applied. Nevertheless, results are considerably good, since the solution is obtained in wavelet space and no discretization was performed. The discontinuity in the slope of the bending moment is captured even for a low level of resolution.

In a second example, critical loads for a doubly clamped beam were obtained by solving an eigenvalue problem using stiffness and geometry matrices, as in Eq. (14). Results at different levels of resolution are shown in Table 1. All values are normalized by EI/L². Figure 6 shows the shape of the first three buckling modes obtained for the mentioned example. Results were normalized.
<figure><img src="image_6.png" /><figcaption>Figure 4: Bending moment using IN8</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 5: Shear force using IN8</figcaption></figure>

<table><caption>Table 1. Critical loads obtained at different levels of resolution for a doubly clamped beam</caption><thead><tr><th>MODE N°.</th><th>EXACT</th><th>LEVEL 0</th><th>LEVEL 2</th><th>LEVEL 4</th><th>LEVEL 6</th></tr></thead><tbody><tr><td>1</td><td>39.4784</td><td>39.4788</td><td>39.4788</td><td>39.4784</td><td>39.4784</td></tr><tr><td>2</td><td>80.7629</td><td>81.7817</td><td>80.7779</td><td>80.7629</td><td>80.7629</td></tr><tr><td>3</td><td>157.9137</td><td>159.1386</td><td>158.6607</td><td>157.9145</td><td>157.9137</td></tr></tbody></table>

The same example was analyzed using the IN6 Interpolet in order to verify the sensitivity to wavelet order. Figure 7 shows the relative error obtained for the third critical load at different levels of resolution for each Interpolet.
<figure><img src="image_8.png" /><figcaption>Figure 6: Buckling modes for the doubly clamped beam</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 7: Relative error in third critical load at different levels of resolution</figcaption></figure>

Finally, to test the possibility of extending the method to two-dimensional problems, a thin plate was modeled using the equations developed in previous sections. Figure 8 shows a square plate with all edges clamped subjected to a concentrated load applied at its center.

The plate was modeled using the IN6 Interpolet at level 3, leading to a total number of 289 degrees of freedom. The result for the central displacement was $w = -0.00557 PL^2/D$ which represents an error of 0.5% when compared to the exact solution $w = -0.00560 PL^2/D$. Results were extremely good, considering that a FE mesh using 32x32 plate elements with 12 degrees of freedom each gives an error of 0.7% in the central displacement.

Figure 9 shows the results for the bending moments $M_x$, $M_y$ and the twisting moment $M_{xy}$. Displacements and moments distribution were obtained using the wavelet's second derivatives. The errors in the bending moments $M_x$ and $M_y$ at the center point were 4%.
<figure><img src="image_10.png" /><figcaption>Figure 8: Clamped plate subjected to a concentrated load at the center</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 9: Results for moments $M_x$, $M_y$, $M_{xy}$ and displacement $w$</figcaption></figure>

Different types of boundary conditions and loadings were tested for a square plate and the values obtained for central displacement are summarized in Table 2. Results were compared with exact solutions given by Timoshenko and Woinowsky-Krieger (1959).
<table><caption>Table 2: Results for different types of boundary conditions and loadings</caption><thead><tr><th>Boundary Conditions and Loading Type</th><th>Exact</th><th>WGM</th><th>Error</th></tr></thead><tbody><tr><td>Clamped / Uniform</td><td>0.00126 qL<sup>4</sup>/D</td><td>0.00126 qL<sup>4</sup>/D</td><td>0.4 %</td></tr><tr><td>Clamped / Concentrated</td><td>0.00560 PL<sup>2</sup>/D</td><td>0.00557 PL<sup>2</sup>/D</td><td>0.5 %</td></tr><tr><td>Simply Supported / Uniform</td><td>0.00406 qL<sup>4</sup>/D</td><td>0.00406 qL<sup>4</sup>/D</td><td>0.1 %</td></tr><tr><td>Simply Supported / Concentrated</td><td>0.01160 PL<sup>2</sup>/D</td><td>0.01156 PL<sup>2</sup>/D</td><td>0.3 %</td></tr></tbody></table>

## 5 CONCLUSIONS

This work presented the formulation and validation of the Wavelet-Galerkin Method using Deslauriers-Dubuc Interpolets. It was also shown that wavelets have the ability of capturing discontinuities without the need to place nodes where they occur.

As in the traditional FEM and other numerical methods, the accuracy of the solution can be improved either by increasing the level of resolution or the wavelet order. Sometimes, lower order wavelets at higher resolution can give better results than higher order wavelets at lower resolutions.

For two-dimensional problems, results for displacements and bending moments were extremely good, although only regular geometry problems were studied. The extension of the method to irregular geometries is still a challenge.

Since the unknowns of the method are interpolation coefficients instead of nodal displacements, it is possible to obtain a smooth representation of bending moments even with a reduced number of degrees of freedom.

All matrices involved can be stored and operated in a sparse form, since most of their components are null, thus saving computer resources. Due to the compact support of wavelets, the sparseness of matrices increases along with the level of resolution.

## REFERENCES

* Burgos, R. B., Cetale Santos, M. A. and Silva, R. R., The use of wavelets and interpolets as interpolating functions in a Finite Element analysis, *Proceedings of the 29th CILAMCE - Iberian Latin American Congress on Computational Methods in Engineering*, 2008.

* Chen, X., Yang, S., Ma, J. and He, Z., The construction of wavelet finite element and its application, *Finite Elements in Analysis and Design*, 40: 541–554, 2004.

* Daubechies, I., Orthonormal bases of compactly supported wavelets, *Communications in Pure and Applied Mathematics*, 41: 909-996, 1988.

* Deslauriers, G. and Dubuc, S., Symmetric iterative interpolation processes, *Constructive Approximation*, 5: 49-68, 1989.

* Latto, A., Resnikoff, L. and Tenenbaum, E., The evaluation of connection coefficients of compactly supported wavelets, *Proceedings of the French-USA Workshop on Wavelets and Turbulence*, 76-89, 1992.

* Lin, W., Kovvali, N. and Carin, L., Direct algorithm for computation of derivatives of the Daubechies basis functions, *Applied Mathematics and Computation*, 170: 1006-1013, 2005.

* Ma, J., Xue, J., Yang, S. and He, Z., A study of the construction and application of a Daubechies wavelet-based beam element, *Finite Elements in Analysis and Design*, 39: 965–975, 2003.
Nguyen, V. P., Rabczuk, T., Bordas, S. and Duflot, M., Meshless methods: a review and computer implementation aspects, Mathematics and Computers in Simulation, 79: 763-813, 2008.

Qian, S and Weiss, J., Wavelets and the numerical solution of partial differential equations, *Journal of Computational Physics*, 106: 155-175, 1992.

Shi, Z., Kouri, D. J., Wei, G. W. and Hoffman, D. K., Generalized symmetric interpolating wavelets, *Computer Physics Communications*, 119: 194-218, 1999.

Timoshenko, S. and Woinowsky-Krieger, S., *Theory of Plates and Shells*, McGraw-Hill, 1959.

Zhou, X. and Zhang, W., The evaluation of connection coefficients on an interval, *Communications in Nonlinear Science & Numerical Simulation*, 3: 252-255, 1998.