# The computation of Wavelet-Galerkin approximation on a bounded interval

Ming-Quayer Chen, Chyi Hwang, Yen-Ping Shih

► To cite this version:

Ming-Quayer Chen, Chyi Hwang, Yen-Ping Shih. The computation of Wavelet-Galerkin approximation on a bounded interval. International Journal for Numerical Methods in Engineering, 1996, 39 (17), pp.2921-2944. 10.1002/(SICI)1097-0207(19960915)39:173.0.CO;2-D . hal-01330583

HAL Id: hal-01330583

https://hal.science/hal-01330583

Submitted on 11 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
THE COMPUTATION OF WAVELET-GALERKIN
APPROXIMATION ON A BOUNDED INTERVAL

MING-QUAYER CHEN,

*Department of Chemical Engineering, National Cheng Kung University, Tainan 701, Taiwan*

CHYI HWANG,

*Department of Chemical Engineering, National Chung Cheng University, Chia-Yi 621, Taiwan*

YEN-PING SHIH

*Department of Electrical Engineering, National Taiwan Ocean University, Keelung 200, Taiwan*

This paper describes exact evaluations of various finite integrals whose integrands involve products of Daubechies' compactly supported wavelets and their derivatives and/or integrals. These finite integrals play an essential role in the wavelet-Galerkin approximation of differential or integral equations on a bounded interval.

KEY WORDS: wavelet orthogonal bases; wavelet-Galerkin method; Burgers' equation; numerical method; approximation

# 1. INTRODUCTION

In a recent paper, Daubechies¹ constructed a family of orthonormal bases of compactly supported wavelets for the space of square-integrable functions, $L^2(\mathbf{R})$. Due to the fact that they possess several useful properties, such as orthogonality, compact support, exact representation of polynomials to a certain degree, and ability to represent functions at different levels of resolution, Daubechies' wavelets have gained great interest in the numerical solutions of ordinary and partial differential equations.²⁻¹¹ In each of these papers the wavelet discretization of differential equations is based on the Galerkin approach. As a result, the wavelet-Galerkin scheme involves the evaluation of connection coefficients¹², ¹³ to approximate derivatives as well as non-linear terms. The connection coefficients are integrals with integrands being products of wavelet bases and their derivatives. Due to the derivatives of compactly supported wavelets being highly oscillatory, it is difficult and unstable to compute the connection coefficients by the numerical evaluation of integral. In order to overcome this problem, dedicated algorithms have been devised¹²⁻¹⁶ for the exact evaluation of connection coefficients.

We note that the connection coefficients and the associated computation algorithms developed in References 12–16 are essentially based on an unbounded domain. It is therefore not surprising that the above-mentioned applications of the wavelet-Galerkin method are limited to the cases where the problem domain is unbounded or the boundary condition is periodic. In order to apply the wavelet-Galerkin method to the solution of finite-domain problems, the authors¹⁷, ¹⁸ have derived algorithms for computing some finite integrals of products of wavelets and their derivatives or
integrals. As motivated by the successful construction of orthogonal interval wavelets,¹⁹,²⁰ we present in this paper algorithms for computing useful finite integrals of wavelets for the wavelet-Galerkin method on a bounded interval. More precisely, let $\phi(x)$ denote the scaling function of Daubechies' wavelets, we consider the computation of the following functionals:

$$ \phi^{(n)}(x) = \frac{d^n \phi(x)}{dx^n} \qquad (1) $$

$$ \theta_n(x) = \underbrace{\int_0^x \int_0^{y_n} \dots \int_0^{y_2}}_{n\text{-tuple}} \phi(y_1) dy_1 \dots dy_{n-1} dy_n \qquad (2) $$

$$ M_k^m(x) = \int_0^x y^m \phi(y-k) dy \qquad (3) $$

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

$$ \Lambda_k^{m,n}(x) = \int_0^x y^m \phi^{(n)}(y-k) \phi(y) dy \qquad (5) $$

$$ \Upsilon_k^{m,n}(x) = \int_0^x y^m \theta_n(y-k) \phi(y) dy \qquad (6) $$

$$ \Omega_{jk}^{m,n}(x) = \int_0^x \phi(y) \phi^{(m)}(y-j) \phi^{(n)}(y-k) dy \qquad (7) $$

In the above equations, $j, k \in \mathbb{Z}$, $m, n, x \in \mathbb{Z}^+ + \{0\}$, and $\mathbb{Z}$ and $\mathbb{Z}^+$ denote the sets of integers and positive integers, respectively.

Before attempting to present the computational algorithms for the above integrals, we first give in the following section a brief description of the construction of Daubechies' class of orthonormal wavelets and the basic properties of these wavelets.

## 2. DAUBECHIES' ORTHONORMAL WAVELETS

The family of compactly supported orthonormal wavelets constructed by Daubechies¹ includes members from highly localized to highly smooth. Each wavelet member is governed by a set of $L$ (an even integer) coefficients $\{p_k: k = 0, 1, \dots, L-1\}$ through the two-scale relation

$$ \phi(x) = \sum_{j=0}^{L-1} p_j \phi(2x-j) \qquad (8) $$

and the equation

$$ \psi(x) = \sum_{j=2-L}^{1} (-1)^j p_{1-j} \phi(2x-j) \qquad (9) $$

where $\phi(x)$ and $\psi(x)$ are called scaling function and mother wavelet, respectively. The fundamental support of the scaling function $\phi(x)$ is in the interval $[0, L-1]$ while that of the corresponding wavelet $\psi(x)$ is in the interval $[1 - L/2, L/2]$.
The coefficients $p_k$ in the two-scale relation (8) are called the wavelet filter coefficients. Daubechies¹ established these wavelet filter coefficients to satisfy the following conditions:

$$ \sum_{j=0}^{L-1} p_j = 2 \qquad (10) $$

$$ \sum_{j=0}^{L-1} p_j p_{j-m} = \delta_{0,m} \qquad (11) $$

$$ \sum_{j=2-L}^{L-1} (-1)^j p_{1-j} p_{j-2m} = 0 \quad \text{for integer } m \qquad (12) $$

$$ \sum_{j=0}^{L-1} (-1)^j j^m p_j = 0, \quad m = 0, 1, \dots, L/2 - 1 \qquad (13) $$

where $\delta_{0,m}$ is the Kronecker delta function. Correspondingly, the constructed scaling function $\phi(x)$ and wavelet $\psi(x)$ have the following properties:

$$ \int_{-\infty}^{\infty} \phi(x) dx = 1 \qquad (14) $$

$$ \int_{-\infty}^{\infty} \phi(x-j) \phi(x-m) dx = \delta_{j,m} \qquad (15) $$

$$ \int_{-\infty}^{\infty} \phi(x) \psi(x-m) dx = 0 \quad \text{for integer } m \qquad (16) $$

$$ \int_{-\infty}^{\infty} x^k \psi(x) dx = 0, \quad k = 0, 1, \dots, L/2 - 1 \qquad (17) $$

It is noted that the property (17) is equivalent to that the elements of the set $\{1, x, \dots, x^{L/2-1}\}$ are linear combinations of $\phi(x-k)$, integer translates of $\phi(x)$. The exact expression for such linear combinations is given by Reference 12:

$$ \sum_{l=-\infty}^{\infty} l^n \phi(x-l) = x^n + \sum_{j=1}^{n} (-1)^j \binom{n}{j} M_j^{\phi} x^{n-j}, \quad n = 0, 1, \dots, L/2 - 1 \qquad (18) $$

In the above equation, $M_j^\phi$ denotes the jth moment of $\phi(x)$, which can be computed by the recursive relation:¹³

$$ \begin{aligned} M_k^\phi &= \int_{-\infty}^{\infty} x^k \phi(x) dx \\ &= \frac{1}{2^{k+1}-2} \sum_{i=0}^{L-1} \sum_{l=1}^{k} \binom{k}{l} p_i i^l M_{k-l}^\phi \end{aligned} \qquad (19) $$

with the initial condition $M_0^\phi = 1$.

Denote by $L^2(\mathbb{R})$ the space of square-integrable functions on the real line. Let $V_j$ and $W_j$ be the subspaces generated, respectively, as the $L^2$-closure of the linear spans of $\phi_{j,k}(x) = 2^{j/2}\phi(2^j x - k)$ and $\psi_{j,k}(x) = 2^{j/2}\psi(2^j x - k)$, $k \in \mathbb{Z}$. Then the condition (16) implies that

$$ V_{j+1} = V_j \oplus W_j \qquad (20) $$

The above relation further implies²¹

$$ V_0 \subset V_1 \subset \dots \subset V_j \subset V_{j+1} \qquad (21) $$
and

$$V_{j+1} = V_0 \oplus W_0 \oplus W_1 \oplus \cdots \oplus W_j \quad (22)$$

where $\oplus$ denotes the orthogonal direct sum. For the Daubechies wavelets, we have the following orthogonal properties:

$$\int_{-\infty}^{\infty} \phi_{j,k}(x) \phi_{j,l}(x) \, dx = \delta_{k,l} \quad (23)$$

$$\int_{-\infty}^{\infty} \psi_{j,k}(x) \psi_{l,m}(x) \, dx = \delta_{j,l} \delta_{k,m} \quad (24)$$

$$\int_{-\infty}^{\infty} \phi_{j,k}(x) \psi_{j,m}(x) \, dx = 0 \quad (25)$$

It is noted that there are no explicit expressions for calculating the values of the scaling function $\phi(x)$ and the corresponding mother wavelet $\psi(x)$ at an arbitrary point of $x$. However, the function values of $\phi(x)$ and $\psi(x)$ at the dyadic points $k/2^j$ for integers $j$ and $k$ can be recursively computed from the two-scale relations (8) and (9) provided that $\phi(1), \phi(2), \dots, \phi(L-2)$ have been obtained. The algorithm for obtaining $\phi(k)$ for $k=1, 2, \dots, L-2$ will be given in the next section.

### 3. COMPUTING THE DERIVATIVES AND INTEGRALS OF WAVELETS

In this section, we present computational algorithms for evaluating the derivatives and integrals of the scaling function $\phi(x)$. Once the values of scaling function $\phi(x)$ are obtained, the mother wavelet $\psi(x)$ can be computed from the relation (9).

#### 3.1. Evaluation of the nth derivative of $\phi(x)$

Denote by $\phi^{(n)}(x)$ the nth derivative of the scaling function $\phi(x)$:

$$\phi^{(n)}(x) = \frac{d^n \phi(x)}{dx^n} = \frac{d}{dx} \phi^{(n-1)}(x), \quad \phi^{(0)}(x) = \phi(x) \quad (26)$$

It follows from (18) that the derivatives $\phi^{(n)}(x)$, $n = 0, 1, \dots, L/2 - 1$, exist. Hence, applying the two-scale relation (8) to the above equation, we have

$$\phi^{(n)}(x) = 2^n \sum_{k=0}^{L-1} p_k \phi^{(n)}(2x-k) \, dx \quad (27)$$

This is the two-scale relation for $\phi^{(n)}(x)$, which is analogous to the two-scale relation (8) for $\phi(x)$. Obviously, this two-scale relation can be used to compute the values of $\phi^{(n)}(x)$ at all dyadic points $x = k/2^j$, $j = 1, 2, \dots$, provided that the values $\phi^{(n)}(k)$, $k = 1, 2, \dots, L-2$, are given.

To obtain the values $\phi^{(n)}(x)$ at integer points, we substitute $x = 1, 2, \dots, L-2$ into the two-scale relation (27) to give the homogeneous linear system of equations

$$2^{-n} \Phi = P \Phi \quad (28)$$

where

$$\Phi = [\phi^{(n)}(1) \ \phi^{(n)}(2) \ \dots \ \phi^{(n)}(L-2)]^T \quad (29)$$
(the superscript T denotes transpose) and P is the $(L-2) \times (L-2)$ matrix

$$ P = [p_{2j-k}]_{1 \le j,k \le L-2} \qquad (30) $$

with $j$ being the row index and $k$ the column index. Equation (28) indicates that the unknown vector $\Phi$ is the eigenvector of the matrix $P$ corresponding to the eigenvalue $2^{-n}$. The values of $\phi^{(n)}(1), \phi^{(n)}(2), \dots, \phi^{(n)}(L-2)$ can be determined uniquely by first finding the eigenvector of the matrix $P$ which corresponds to the eigenvalue $2^{-n}$, and then normalizing with the condition

$$ \sum_{k=1}^{L-2} (-k)^n \phi^{(n)}(k) = n! \qquad (31) $$

which is obtained by differentiating (18) $n$ times and then letting $x = k$.

Once the values of $\phi^{(n)}(x)$ for $x = 1, 2, \dots, L-2$ are obtained, the relation

$$ \phi^{(n)}\left(\frac{k}{2^j}\right) = 2^n \sum_{l=0}^{L-1} p_l \phi^{(n)}\left(\frac{k}{2^{j-1}} - l\right) \qquad (32) $$

and the facts that $\phi^{(n)}(x) = 0$ for $x \le 0$ and $x \ge L-1$ allow one to determine the values of $\phi^{(n)}(x)$ at $x = k/2^j$ for $k = 1, 3, 5, \dots, 2^j(L-1)-1$ and $j = 1, 2, \dots$. Here, it should be pointed out that the above algorithm for generating $\phi^{(n)}(k)$, $k = 1, 2, \dots, L-2$, also holds for $n=0$, which corresponds to the scaling function $\phi(x)$ itself. For $n=0$, (27) is exactly the same as the two-scale relation (8) for $\phi(x)$. Moreover, the relation (18) with $n=0$ is called the resolution of identity of the scaling function $\phi(x)$.

### 3.2. Evaluation of multiple integrals of $\phi(x)$

Denote by $\theta_n(x)$ the $n$-tuple integrals of $\phi(x)$, i.e.

$$ \begin{aligned} \theta_n(x) &= \underbrace{\int_0^x \int_0^{y_n} \dots \int_0^{y_2}}_{n\text{-tuple}} \phi(y_1) dy_1 \dots dy_{n-1} dy_n \\ &= \int_0^x \theta_{n-1}(y) dy \end{aligned} \qquad (33) $$

The algorithm for computing $\theta_n(x)$ at dyadic points $k/2^j$ has been derived by authors in a recent paper.¹⁷ For completeness, we repeat the derivation here.

Applying the two-scale relation (8) to (33), we have

$$ \begin{aligned} \theta_n(x) &= \sum_{k=0}^{L-1} p_k \int_0^x \int_0^{x_n} \dots \int_0^{x_2} \phi(2x_1 - k) \, dx_1 \dots \, dx_n \\ &= 2^{-n} \sum_{k=0}^{L-1} p_k \theta_n(2x - k) \end{aligned} \qquad (34) $$

This is the two-scale relation for $\theta_n(x)$, which is analogous to the two-scale relation (8) for $\phi(x)$. Hence, the values of $\theta_n(x)$ at dyadic points $k/2^j$ can be computed recursively for $j = 1, 2, \dots$ provided that the values of $\theta_n(x)$ at integer points $x = 1, 2, \dots$, are pre-computed. Since the values of $\theta_n(x)$ for $x \ge L-1$ do not vanish, the computation of $\theta_n(x)$ for $x = 1, 2, \dots$ is quite different from that of $\phi(x)$.
Here, we first examine some properties of $\theta_n(x)$ for $x \ge L-1$. Let us recall the facts that $\phi(x)$ vanishes for $x \ge L-1$, and that $\int_0^{L-1} \phi(x) dx = 1$. Hence, for $x \ge L-1$, we can obtain from the last equation of (33) that

$$ \theta_1(x) = \int_0^{L-1} \phi(y) dy + \int_{L-1}^{x} \phi(y) dy = 1 \quad (35) $$

$$ \begin{aligned} \theta_2(x) &= \int_0^{L-1} \theta_1(y) dy + \int_{L-1}^{x} \theta_1(y) dy \\ &= \theta_2(L-1) + \int_{L-1}^{x} 1 dy \\ &= \theta_2(L-1) + (x - L + 1) \end{aligned} \quad (36) $$

$$ \begin{aligned} \theta_3(x) &= \int_0^{L-1} \theta_2(y) dy + \int_{L-1}^{x} [\theta_2(L-1) + (y-L+1)] dy \\ &= \theta_3(L-1) + \theta_2(L-1)(x-L+1) + \frac{1}{2!} \theta_1(L-1)(x-L+1)^2 \end{aligned} \quad (37) $$

By induction, we can write $\theta_n(x)$ for $x \ge L-1$ as follows:

$$ \theta_n(x) = \sum_{j=0}^{n-1} \frac{(x-L+1)^j}{j!} \theta_{n-j}(L-1) \quad (38) $$

where $\theta_1(L-1) = 1$, but the values of $\theta_{n-j}(L-1)$ for $j = 0, 1, \dots, n-2$ are still to be determined. For determining $\theta_n(L-1)$ for $n = 2, 3, \dots$, we go back to the two-scale relation (34) for $\theta_n(x)$. Substituting $x = L-1$ into the equation yields

$$ \theta_n(L-1) = 2^{-n} \sum_{k=0}^{L-1} p_k \theta_n(2L-2-k) \quad (39) $$

Since $2L-2-k \ge L-1$ for $k = 0, 1, \dots, L-1$, from (38) it follows that

$$ \theta_n(2L-2-k) = \sum_{j=0}^{n-1} \theta_{n-j}(L-1) \frac{(L-1-k)^j}{j!} \quad (40) $$

As a result, $\theta_n(L-1)$ can be expressed as

$$ \begin{aligned} \theta_n(L-1) &= 2^{-n} \sum_{k=0}^{L-1} p_k \sum_{j=0}^{n-1} \theta_{n-j}(L-1) \frac{(L-1-k)^j}{j!} \\ &= 2^{-n} \left[ \theta_n(L-1) \sum_{k=0}^{L-1} p_k + \sum_{j=1}^{n-1} \sum_{k=0}^{L-1} p_k \theta_{n-j}(L-1) \frac{(L-1-k)^j}{j!} \right] \end{aligned} \quad (41) $$

Applying the relation (10), we then have the recursive formula

$$ \theta_n(L-1) = \frac{1}{2^n-2} \sum_{j=1}^{n-1} \left( \sum_{k=0}^{L-1} p_k \frac{(L-1-k)^j}{j!} \right) \theta_{n-j}(L-1) \quad (42) $$

for computing $\theta_n(L-1)$ for $n = 2, 3, \dots$, starting with $\theta_1(L-1) = 1$. Having obtained $\theta_n(x)$ for an integer $x \ge L-1$, the values of $\theta_n(x)$ for $x = 1, 2, \dots, L-2$ can be determined from the following linear system of equations:

$$ (I - 2^{-n} P) \Theta_n = c \quad (43) $$
where

$$
\Theta_n = [\theta_n(1) \ \theta_n(2) \ \dots \ \theta_n(L-2)]^T \quad (44)
$$

$$
\mathbf{c} = [c_1 \ c_2 \ \dots \ c_{L-2}]^T \qquad (45)
$$

$$
c_i = \sum_{\substack{2-i \ge L-1 \\ k=0,1,\dots,L-1}} p_k \theta_n(2i-k) \quad (46)
$$

and the matrix P is defined by (30). The equations of the linear system (43) come from (34) by letting $x = 1, 2, \dots, L-2$.

4. COMPUTING THE CONNECTION COEFFICIENTS OF WAVELETS OVER A BOUNDED INTERVAL

In solving a differential-integral equation of the form

$$
f \left( x, \frac{dy}{dx}, \frac{d^2 y}{dx^2}, \dots, \int^x y dx_1, \int \int^x y dx_2 dx_1, \dots \right) = 0 \quad (47)
$$

defined on a bounded interval $x \in [a, b]$ by the wavelet-Galerkin method, the dependent variable $y(x)$ is approximated by the wavelet series

$$
\hat{y}(x) = \sum_{k=N_a}^{N_b} \hat{y}_k \phi_{J,k}(x) \quad (48)
$$

By the Galerkin approach, the coefficients $\hat{y}_k$ are determined from the following equations:

$$
\int_a^b \phi_{J,k}(x) f \left( x, \frac{d\hat{y}}{dx}, \frac{d^2 \hat{y}}{dx^2}, \dots, \int^x \hat{y} dx_1, \int \int^x \hat{y} dx_2 dx_1, \dots \right) = 0 \quad (49)
$$

for $k = N_a, N_a + 1, \dots, N_b$. The construction of these equations involves the computation of finite integrals in (3)-(7). As $x \to \infty$, these integrals are called the connection coefficients of the wavelets.¹³ In this section, algorithms will be derived for computing the connection coefficients of wavelets over a bounded interval.

4.1. Evaluation of $M_k^m(x) = \int_0^x y^m \phi(y-k) dy$

Let the integral of the product of $y^m$ and $\phi(y-k)$ be denoted by

$$
M_k^m(x) = \int_0^x y^m \phi(y-k) dy \quad (50)
$$

Performing integration by parts successively $m$ times on the above integral, we have

$$
\begin{align*}
M_k^m(x) &= \int_0^x y^m \phi(y-k) \, dy \\
&= x^m \theta_1(x-k) - m \int_{-k}^{x-k} (y+k)^{m-1} \phi(y) \, dy \\
&\vdots \\
&= \sum_{i=0}^m (-1)^i \frac{m!}{(m-i)!} x^{m-i} \theta_{i+1}(x-k) + (-1)^{m+1} m! \theta_{m+1}(-k) \tag{51}
\end{align*}
$$
Hence, the integral $M_k^n(x)$ can be computed exactly in terms of $\theta_m(x)$, which is the multiple integral of the scaling function $\phi(x)$ defined in Section 3.1.

### 4.2. Evaluation of integrals $\Gamma_k^n(x) = \int_0^x \phi^{(n)}(y-k)\phi(y)dy$

Denote the integral of the product of the scaling function $\phi(x)$ and its nth-order derivative $\phi^{(n)}(x-k)$ by

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

The evaluation of $\Gamma_k^n(x)$ plays an important role in applying the wavelet-Galerkin method to solve differential equations.¹⁸ Since there are no explicit expressions for representing the scaling function $\phi(x)$ and its derivatives, the values of the integral $\Gamma_k^n(x)$ cannot be computed directly from (52). An algorithm has been developed by Beylkin¹² to compute $\Gamma_k^n(\infty)$. The authors¹⁷ developed an algorithm for computing the values of the integrals $\Gamma_k^n(x)$ at dyadic points $k/2^j$ for positive integers $j$ and $k$. To make this paper self-contained, we derive the algorithm in this subsection.

As we begin, we note that $\Gamma_k^n(x)$ has the following properties:

$$ \Gamma_k^n(x) = \Gamma_k^n(L-1) \quad \text{for } x \ge L-1 \quad (53) $$

$$ \Gamma_k^n(x) = 0 \quad \text{for } |k| \ge L-1 \text{ or } x \le 0 \text{ or } x \le k \quad (54) $$

$$ \Gamma_{-k}^{n}(L-1) = (-1)^{n} \Gamma_{k}^{n}(L-1) \quad \text{for } k \ge 0 \quad (55) $$

$$ \Gamma_{-k}^{n}(x) = (-1)^{n} \Gamma_{k}^{n}(L-1) \quad \text{for } x - k \ge L-1 \quad (56) $$

The properties (53) and (54) come from the fact that the support of $\phi(x)$ is in the interval $[0, L-1]$, which does not overlap with that of $\phi^{(n)}(x-k)$ for $|k| \ge L-1$. The last two properties are derived from (52) by performing integration by parts for $n$ times.

Now, applying the two-scale relation (27) and (8) to (52), we have

$$ \begin{aligned} \Gamma_k^n(x) &= 2^n \sum_{i=0}^{L-1} \sum_{j=0}^{L-1} p_i p_j \int_0^x \phi^{(n)}(2y - 2k - i) \phi(2y - j) dy \\ &= 2^{n-1} \sum_{i=0}^{L-1} \sum_{j=0}^{L-1} p_i p_j \int_0^{2x-j} \phi^{(n)}(y - 2k - i + j) \phi(y) dy \\ &= 2^{n-1} \sum_{i=0}^{L-1} \sum_{j=0}^{L-1} p_i p_j \Gamma_{2k+i-j}^n (2x-j) \end{aligned} \quad (57) $$

This relation allows one to compute the exact values of $\Gamma_k^n(x)$ at dyadic points $x = k/2^j$ provided that the values of $\Gamma_k^n(x)$ at $x = 1, 2, \dots, L-1$ are given. To determine the values of $\Gamma_k^n(x)$ for integers $k$ and $x$, we first compute the values of $\Gamma_k^n(L-1)$ for $k = 0, 1, \dots, L-2$. Letting $x = L-1$, we obtain from (53) and (57) that

$$ \begin{aligned} \Gamma_k^n(L-1) &= 2^{n-1} \sum_{i=0}^{L-1} \sum_{j=0}^{L-1} p_i p_j \Gamma_{2k+i-j}^n (2L-2-j) \\ &= 2^{n-1} \sum_{i=0}^{L-1} \sum_{j=0}^{L-1} p_i p_j \Gamma_{2k+i-j}^n (L-1) \end{aligned} \quad (58) $$
For $n > 0$, the equations obtained by substituting $k = 0, 1, \dots, L-2$ into (58) constitute the following homogeneous system for unknowns $\Gamma_k^n(L-1)$, $l = 0, 1, \dots, L-2$:

$$ \Gamma^n(L-1) = \mathbf{D} \Gamma^n(L-1) \qquad (59) $$

$$ \Gamma^n(L-1) = [\Gamma_0^n(L-1) \ \Gamma_1^n(L-1) \ \dots \ \Gamma_{L-2}^n(L-1)]^T \qquad (60) $$

$$ \mathbf{D} = [d_{l,m}]_{0 \le l,m \le L-2} \qquad (61) $$

$$ d_{l,m} = 2^{n-1} \left( \sum_{\substack{0 \le i \le L-1 \\ 0 \le j \le L-1 \\ 2l+i-j=m}} p_i p_j + (-1)^n \sum_{\substack{0 \le i \le L-1 \\ 0 \le j \le L-1 \\ 2l+i-j=-m}} p_i p_j \right) \qquad (62) $$

It is noted that in the formulation of (59) the property (56) has been used.

As indicated by (59), the vector $\Gamma^n(L-1)$ is the eigenvector of the matrix $\mathbf{D}$ that corresponds to the unity eigenvalue. However, the homogeneous linear system of equations in (59) does not admit the vector $\Gamma^n(L-1)$ to be determined uniquely. The required additional constraint or normalization condition can be derived from (18). The derivation of such a normalization condition is given below.

The nth-order derivative of (18) is

$$ \sum_{k=-\infty}^{\infty} k^n \phi^{(n)}(x-k) = n! \qquad (63) $$

Multiplying both sides of this equation by $\phi(x)$ and then taking integration from $x = -\infty$ to $\infty$, we obtain

$$ \begin{aligned} \sum_{k=-\infty}^{\infty} k^n \int_{-\infty}^{\infty} \phi^{(n)}(x-k) \phi(x) dx &= \sum_{k=-\infty}^{\infty} k^n \Gamma_k^n(L-1) \\ &= n! \int_{-\infty}^{\infty} \phi(x) dx \\ &= n! \end{aligned} \qquad (64) $$

Substituting the relations $\Gamma_k^n(L-1) = 0$ for $|k| \ge L-1$ and $\Gamma_{-k}^n(L-1) = (-1)^n \Gamma_k^n(L-1)$ into the above equation, we have

$$ \sum_{k=1}^{L-2} k^n \Gamma_k^n(L-1) = \frac{n!}{2} \quad \text{for } n > 0 \qquad (65) $$

which is the desired normalization condition for the eigenvector of the matrix $\mathbf{D}$ associated with the unity eigenvalue.

Once the values of $\Gamma_k^n(L-1)$ are obtained, we can determine the values of $\Gamma_k^n(x)$ for $x = 0, 1, \dots, L-2$ and $|k| \le L-2$. To this end, we note from (53)-(56) that there are only $(L-2)^2$ independent members in the set $\{\Gamma_k^n(x) : x = 1, 2, \dots, L-2; x - L + 2 \le k \le x - 1\}$. Let the independent members be packed in the vector

$$ \Gamma^n = [\Gamma^n(1) \ \Gamma^n(2) \ \dots \ \Gamma^n(L-2)]^T \qquad (66) $$

where

$$ \Gamma^n(x) = [\Gamma_{x-L+2}^n(x) \ \Gamma_{x-L+3}^n(x) \ \dots \ \Gamma_{x-1}^n(x)], \quad x = 1, 2, \dots, L-2 \qquad (67) $$

Then the equations obtained from (57) by substituting $k$ by $x - L + 2$, $x - L + 3$, ..., $x - 1$ and $x$ by $1, 2, \dots, L-2$ can be put into the following matrix-vector form:

$$ (2^{1-n}\mathbf{I} - \mathbf{Q})\Gamma^n = \mathbf{d} \qquad (68) $$
where

$$
Q = \begin{bmatrix}
Q_{1,1} & Q_{1,2} & \cdots & Q_{1,L-2} \\
Q_{2,1} & Q_{2,2} & \cdots & Q_{2,L-2} \\
\vdots & \vdots & \ddots & \vdots \\
Q_{L-2,1} & Q_{L-2,2} & \cdots & Q_{L-2,L-2}
\end{bmatrix} \tag{69}
$$

$$
Q_{i,j} = [q_{i,j,k,m}]_{1 \le k,m \le L-2} \quad (70)
$$

$$
q_{i,j,k,m} = p_{2i-j} p_{L-1-2k+m} \tag{71}
$$

$$
\mathbf{d} = [d(1) \ d(2) \ \dots \ d((L-2)^2)] \quad (72)
$$

$$
d(mL - 2m + k) = \sum_{i,j \in \mu(k,m,L)} p_i p_j \Gamma_{2k+i-j}^{n} (L-1) \quad (73)
$$

and the index set $\mu(k,m,L)$ is given by

$$
\mu(k,m,L) := \{(i,j) : (L-1 \le 2m-2k-i \text{ or } L-1 \le 2m-j), 0 \le i,j \le L-1\} \quad (74)
$$

It is important to note that the eigenvalue set of the matrix Q includes 2⁻ⁿ, m = 0, 1, ..., L−2 with multiplicities L/2 − |((L−2)/2) − m|, but does not include 2. Hence, the vector Γ⁰ can be obtained by solving the linear system (68). For n > 0, however, the matrix (2¹⁻ⁿI − Q) is singular and the rank deficiency is n. To have a linear independent system of equations for the unknown vectors Γⁿ, n = 1, 2, ..., we seek in the following for additional relations among the members of Γⁿ.

Multiplying both sides of (63) by $\phi(y)$ and then taking integration from $y=0$ to $y=x$, we obtain

$$
\sum_{k=-\infty}^{\infty} k^n \int_{0}^{x} \phi(y) \phi^{(n)}(y-k) dy = n! \int_{0}^{x} \phi(y) dy \quad (75)
$$

or

$$
\sum_{k=-\infty}^{\infty} k^n \Gamma_k^n(x) = n! \theta_1(x) \quad (76)
$$

For an integer $x = m$, the above equation can be written as

$$
\sum_{k=m-L+2}^{m-1} k^n \Gamma_k^n(m) = n! \theta_1(m) \quad (77)
$$

where $\theta_1(x)$ is defined in Section 3.2. The above equation can be written in the form

$$
[(m - L + 2)^n (m - L + 3)^n \cdots (m - 1)^n] \Gamma^n(m) = n! \theta_1(m) \quad (78)
$$

Hence, the linear independent system of equations for $\Gamma_n$ can be finally obtained as follows: (i)
For $i = 1,2,...,n$, replace the $m$th row ($m \in \{1,2,...,L-2\}$) of the submatrix $Q_{i,i}$ in (69) by
$[(i-L+2)^n (i-L+3)^n \cdots (i-1)^n]$ and the $m$th row of $Q_{i,j}, j \neq i$ by a zero row vector,
(ii) Replace the $m$th element of the subvector of $d$ that corresponds to the submatrix $Q_{i,i}$ by $n! \theta_1(i)$.

4.3. Evaluation of integrals $\Lambda_k^{m,n}(x) = \int_0^x y^m \phi^{(n)}(y-k) \phi(y) dy$

In solving a differential equation with non-constant coefficients by the wavelet-Galerkin method,
it is required to compute the following integral:

$$
\Lambda_k^{m,n}(x) = \int_0^x y^m \phi^{(n)}(y-k) \phi(y) \, dy, \quad m, n \ge 0 \quad (79)
$$

In this subsection, we devote an effort to the evaluation of this integral.
Since $\phi^{(n)}(x) = 0$ for $x \ge L - 1$ or $x \le 0$, it is easy to verify that the integral $\Lambda_k^{m,n}(x)$ has the following properties:

$$ \Lambda_k^{m,n}(x) = \Lambda_k^{m,n}(L-1) \quad \text{for } x \ge L-1 \text{ or } x-k \ge L-1 \quad (80) $$

$$ \Lambda_k^{m,n}(x) = 0 \quad \text{for } x \le 0 \text{ or } x \le k \text{ or } |k| \ge L-1 \quad (81) $$

Applying the two-scale relations (8) and (27) for $\phi(x)$ and $\phi^{(n)}(x)$, respectively, to the integral (79), we obtain the following two-scale relation for $\Lambda_k^{m,n}(x)$:

$$ \Lambda_k^{m,n}(x) = 2^{n-m-1} \sum_{i=0}^{L-1} \sum_{j=0}^{L-1} \sum_{l=0}^{m} p_i p_j \binom{m}{l} j^l \Lambda_{2k-j+i}^{m-l,n} (2x-j) \quad (82) $$

Since $\Lambda_k^{0,n}(x) = \Gamma_k^n(x)$, the above relation allows one to recursively calculate $\Lambda_k^{m,n}(x)$ for $m = 1, 2, \dots$. However, for fixed integers $m$ and $n$, the values of $\Lambda_k^{m,n}(x)$ at various integers $k$ and $x$ have to be determined simultaneously. The formulation of a linear system of equations for solving $\Lambda_k^{m,n}(x)$ at integer points of $k$ and $x$ is given below.

Substituting $x = L - 1$ into (82), we obtain:

$$ \begin{aligned} 2^{1+m-n} \Lambda_k^{m,n}(L-1) = & \sum_{i=0}^{L-1} \sum_{j=0}^{L-1} p_i p_j \Lambda_{2k-j+i}^{m,n}(L-1) \\ & + \sum_{i=0}^{L-1} \sum_{j=0}^{L-1} \sum_{l=1}^{m} x p_i p_j \binom{m}{l} j^l \Lambda_{2k-j+i}^{m-l,n}(L-1) \end{aligned} \quad (83) $$

It follows from (81) that $\Lambda_k^{m,n}(L-1)$ vanishes for $|k| \ge L-1$. Then, substituting $k = 2-L$, $3-L, \dots, L-2$ into (83), we have the following linear system of equations for the unknowns $\Lambda_k^{m,n}(L-1)$:

$$ (2^{1+m-n}\mathbf{I} - \mathbf{A}) \Lambda_k^{m,n}(L-1) = \mathbf{b} \quad (84) $$

Here $\mathbf{I}$ is the $(2L-3) \times (2L-3)$ identity matrix, and

$$ \Lambda_k^{m,n}(L-1) = [\Lambda_{2-L}^{m,n}(L-1) \ \Lambda_{3-L}^{m,n}(L-1) \ \dots \ \Lambda_{L-2}^{m,n}(L-1)]^T \quad (85) $$

$$ \mathbf{A} = \begin{bmatrix} a_{2-L,2-L} & a_{2-L,3-L} & \cdots & a_{2-L,L-2} \\ a_{3-L,2-L} & a_{3-L,3-L} & \cdots & a_{3-L,L-2} \\ \vdots & \vdots & \ddots & \vdots \\ a_{L-2,2-L} & a_{L-2,3-L} & \cdots & a_{L-2,L-2} \end{bmatrix} \quad (86) $$

$$ a_{r,s} = \sum_{\substack{i=j=s-2r \\ 0 \le i,j \le L-1}} p_i p_j \quad (87) $$

$$ \mathbf{b} = [b(2-L) \ b(3-L) \ \dots \ b(L-2)]^T \quad (88) $$

$$ b(k) = \sum_{i=0}^{L-1} \sum_{j=0}^{L-1} \sum_{l=1}^{m} p_i p_j \binom{m}{l} j^l \Lambda_{2k-j+i}^{m-l,n}(L-1) \quad (89) $$

Since the matrix $\mathbf{A}$ has the eigenvalues $2^{1-i}$ for $i = 0, 1, \dots, L-1$, the matrix $(2^{1+m-n}\mathbf{I} - \mathbf{A})$ is singular for $0 \le n-m \le L-1$. Hence, for the case of $n \le m$, we can directly solve the linear system (84) to find $\Lambda_k^{m,n}(L-1)$ for $k = 2-L, 3-L, \dots, L-2$. As for the case of $n > m$, we have to look for an additional equation which allows, along with (84), one to determine the vector $\Lambda_k^{m,n}(L-1)$. To this end, we apply the relation (18) to obtain:

$$ \sum_{k=-\infty}^{\infty} k^{n-m} \int_{-\infty}^{\infty} y^m \phi^{(n)}(y) \phi(y-k) dy = \sum_{j=0}^{n-m} (-1)^j \binom{n-m}{j} M_j \phi \int_{-\infty}^{\infty} y^{n-j} \phi^{(n)}(y) dy \quad (90) $$
By performing integration by parts on the last integral of (90) for $n-j$ times, we have

$$
\begin{align*}
\int_{-\infty}^{\infty} x^{n-j} \phi^{(n)}(x) \, dx &= (-1)^{n-j} (n-j)! \int_{-\infty}^{\infty} \phi^{(j)}(x) \, dx \\
&= \begin{cases} (-1)^{n-j} (n-j)! [\phi^{(j-1)}(\infty) - \phi^{(j-1)}(-\infty)], & j \ge 1 \\ (-1)^n n! \int_{-\infty}^{\infty} \phi(x) \, dx, & j = 0 \end{cases} \\
&= \begin{cases} 0, & j \ge 1 \\ (-1)^n n!, & j = 0 \end{cases} \tag{91}
\end{align*}
$$

Substituting (80), (81), and (91) into (90), we have

$$
\sum_{k=2-L}^{L-2} \sum_{j=0}^{m} k^{n-m+j} \Lambda_{-k}^{m-j,n} (L-1) = (-1)^n n! \quad \text{for } n \ge m \quad (92)
$$

Making a rearrangement yields

$$
\sum_{k=2-L}^{L-2} k^{n-m} \Lambda_{-k}^{m,n}(L-1) = (-1)^n n! - \sum_{j=1}^{m} \sum_{k=2-L}^{L-2} k^{n-m+j} \Lambda_{-k}^{m-j,n}(L-1) \quad (93)
$$

This is the desired additional relation for the vector $\Lambda^{m,n}(L-1)$. Hence, the non-singular system of equations for $\Lambda^{m,n}(L-1)$ can be formed by replacing any equation in (84) by (93).

Having obtained the values of $\Lambda_k^{m,n}(x)$ for $x = L-1$, we are to determine the values of $\Lambda_k^{m,n}(x)$ for $x = 1, 2, \dots, L-2$ and $|k| \le L-2$. If $m$ and $n$ are specified, there are only $(L-2)^2$ independent unknowns to be determined simultaneously. These independent unknowns are packed into the vector

$$
\Lambda^{m,n} = [\Lambda^{m,n}(1) \ \Lambda^{m,n}(2) \ \dots \ \Lambda^{m,n}(L-2)]^T \quad (94)
$$

where

$$
\Lambda^{m,n}(x) = [\Lambda_{x-L+2}^{m,n}(x) \ \Lambda_{x-L+3}^{m,n}(x) \ \dots \ \Lambda_{x-1}^{m,n}(x)], \quad x = 1, 2, \dots, L-2 \quad (95)
$$

By letting $k = x - L + 2, x - L + 3, \dots, x - 1$ and $x = 1, 2, \dots, L - 2$ in (82), we obtain the following system of equations for the unknown vector $\Lambda^{m,n}$:

$$
(2^{1+m-n}I - Q)\Lambda^{m,n} = e \quad (96)
$$

where

$$
e = [e(1) \ e(2) \ \dots \ e((L-2)^2)]^T \quad (97)
$$

$$
e(xL - 2x + k) = \sum_{i,j \in v_{m,n}(k,x)} p_i p_j \Lambda_{2k+i-j}^{m,n}(L-1) \\
+ \sum_{i=0}^{L-1} \sum_{j=0}^{L-1} \sum_{l=1}^{m} p_i p_j \binom{m}{l} j^l \Lambda_{2k-j+i}^{m-l,n}(2x-j) \quad (98)
$$

The matrix **Q** is defined by (69) and the index set $v_{x,n}(k,x)$ is given by

$$
v_{m,n}(k,x) := \{(i,j): (L-1 \le 2x - 2k - i \text{ or } L-1 \le 2x - j), 0 \le i,j \le L-1\} \quad (99)
$$

As mentioned before, the eigenvalue set of the matrix **Q** contains $2^{-i}$, $i = 0, 1, \dots, L-2$, with multiplicities $L/2 - |((L-2)/2) - i|$. Hence, for $n > m$ the matrix $(2^{1+m-n}\mathbf{I} - \mathbf{Q})$ is singular and its rank deficiency is $n-m$. To have a linear system of independent equations for $\Lambda^{m,n}$, we seek in the following for additional relations among the elements of $\Lambda^{m,n}$.
Multiplying both sides of (63) by $y^m\phi(y)$ and then taking integration from $y=0$ to $y=x$, we obtain

$$ \sum_{k=-\infty}^{\infty} k^n \int_{0}^{x} y^m \phi^{(n)}(y-k) \phi(y) dy = n! \int_{0}^{x} y^m \phi(y) dy \quad (100) $$

For integer $x=i$, the above equation can be written as

$$ \sum_{k=i-L+2}^{i-1} k^n \Lambda_k^{m,n}(i) = n! M_0^m(i) \quad (101) $$

The equations in (101) with $i=1,2,...,n-m$ can be used together with (96) to construct a linear system of independent equations for the unknown vector $\Lambda_k^{m,n}$. The procedure of constructing such a system has been described in Section 4.2.

At this point, we have derived algorithms for computing the values of $\Lambda_k^{m,n}(x)$ at integer points $x$. Once these values are obtained, the values of $\Lambda_k^{m,n}(x)$ at dyadic points $x = k/2^j$ can be recursively obtained from the two-scale relation (82).

### 4.4. Evaluation of integrals $\Upsilon_k^{m,n}(x) = \int_0^x y^m \theta_n(y-k) \phi(y) dy$

In solving an integro-differential equation with non-constant coefficients by the wavelet-Galerkin method, it is required to compute the integral

$$ \Upsilon_k^{m,n}(x) = \int_0^x y^m \theta_n(y-k) \phi(y) dy, \quad n > 0 \quad (102) $$

Since $\phi(x) = 0$ for $x \ge L-1$ or $x \le 0$, it is easy to verify that the integral $\Upsilon_k^{m,n}(x)$ has the following properties:

$$ \Upsilon_k^{m,n}(x) = \Upsilon_k^{m,n}(L-1) \quad \text{for } x \ge L-1 \quad (103) $$

$$ \Upsilon_k^{m,n}(x) = 0 \quad \text{for } x \le 0 \quad \text{or} \quad x \le k \quad \text{or} \quad k \ge L-1 \quad (104) $$

With the substitutions of the two-scale relations (34) and (8) for $\theta_n(x)$ and $\phi(x)$, respectively, into (102), we obtain

$$ \Upsilon_k^{m,n}(x) = 2^{-(m+n+1)} \sum_{i=0}^{L-1} \sum_{j=0}^{L-1} \sum_{l=0}^{m} p_i p_j \binom{m}{l} j^l \Upsilon_{2k-j+i}^{m-k,n}(2x-j) \quad (105) $$

It is noted that for $k \le 1-L$ and $x \ge 0$, we have that $x-k = x+L-1 \ge L-1$ and the integral $\theta_n(x-k)$ is explicitly given by (42). As a result, the integral $\Upsilon_k^{m,n}(x)$ for $k \le 1-L$ can be explicitly evaluated by

$$ \Upsilon_k^{m,n}(x) = \sum_{j=0}^{n-1} \theta_{n-j}(L-1) \left[ \sum_{l=0}^{j} \binom{j}{l} \frac{(1-L)^l}{j!} M_0^{m+j-l}(x) \right] \quad (106) $$

Having derived formulas for computing $\Upsilon_k^{m,n}(x)$ for $x > L-1$, we can solve the independent unknown $\Upsilon_k^{m,n}(x)$, $x = 1,2,...,L-1$, from (105). For fixed integers $m$ and $n$, the independent unknowns of $\Upsilon_k^{m,n}(x)$ are packed in the vector

$$ \Upsilon^{m,n} = [\Upsilon^{m,n}(1) \ \Upsilon^{m,n}(2) \ \dots \ \Upsilon^{m,n}(L-1)]^T \quad (107) $$

where

$$ \Upsilon^{m,n}(k) = [\Upsilon_{k-L+1}^{m,n}(k) \ \Upsilon_{k-L+2}^{m,n}(k) \ \dots \ \Upsilon_{k-1}^{m,n}(k)], \quad k = 1, 2, \dots, L-1 \quad (108) $$
The total number of elements of the vector $\Upsilon^{m,n}$ is $(3L-4)(L-1)/2$. It is noted that the required linear system of equations for solving this unknown vector is formed from (105), and the system matrix is full rank. In constructing such a system of equations, the values of the integral $\Upsilon_k^{m-i,n}(x)$ in (105) for $i < m$ are known. Hence, the vectors $\Upsilon^{m,n}$ can be computed recursively for $m = 0, 1, \dots$ and $n = 1, 2, \dots$ starting from $\Upsilon^{0,1}$.

### 4.5. Evaluation of three-term connection coefficients

The integrals

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

for $0 \le m \le n \le L/2 - 1$ and $j, k, m, n, x \in \mathbb{Z}$, play an important role in solving non-linear differential equations with the wavelet-Galerkin method.²² In Reference 13, $\Omega_{j,k}^{m,n}(\infty)$ is called the three-term connection coefficient of wavelets. In order to extend the wavelet-Galerkin method to finite-domain problems, we devote an effort to compute the values of $\Omega_{j,k}^{m,n}(x)$ at dyadic points $x = k/2^j$.

As we begin, we note that $\Omega_{j,k}^{m,n}(x)$ has the following properties:

$$ \Omega_{j,k}^{m,n}(x) = 0 \quad \text{for } |j|, |k|, \text{ or } |j-k| \ge L-1 \quad (110) $$

$$ \Omega_{j,k}^{m,n}(x) = 0 \quad \text{for } x-j, x-k, \text{ or } x \le 0 \quad (111) $$

$$ \Omega_{j,k}^{m,n}(x) = \Omega_{j,k}^{m,n}(L-1) \quad \text{for } x-j, x-k, \text{ or } x \ge L-1 \quad (112) $$

The formal properties (110) and (111) come from the fact that the supports of $\phi^{(i)}(x)$ for $i = 0, 1, 2, \dots$ are all in the interval $[0, L-1]$.

Now, substituting two-scale relations (8) and (27) into (109), we have

$$ \Omega_{j,k}^{m,n}(x) = 2^{m+n-1} \sum_{i_u=0}^{L-1} \sum_{i_b=0}^{L-1} \sum_{i_c=0}^{L-1} p_{i_u} p_{i_b} p_{i_c} \Omega_{2j+i_b-i_u, 2k+i_c-i_c}^{m,n} (2x - i_a) \quad (113) $$

Since $\Omega_{j,k}^{m,n}(x) = \Omega_{j,k}^{m,n}(L-1)$ for $x \ge L-1$, the values of $\Omega_{j,k}^{m,n}(L-1)$ for integers $j, k, m, n$ can be first computed. For this purpose, we note from (110)-(112) that, for fixed integers $m$ and $n$, there are $3L^2 - 9L + 7$ unknown $\Omega_{j,k}^{m,n}(L-1)$, where $k \in [j+2-L, j-2+L] \cap [2-L, L-2]$ and $j = 2-L, \dots, 0, \dots, 2-L$. The equations in (113) associated with these unknowns constitute a linear homogeneous system as follows:

$$ v = 2^{1-m-n} S v \quad (114) $$

where

$$ v = [v_{2-L} \ v_{3-L} \ \dots \ v_{L-2}]^T \quad (115) $$

$$ v_j = [\Omega_{j,v}^{m,n}(L-1) \ \Omega_{j,v+1}^{m,n}(L-1) \ \dots \ \Omega_{j,\mu}^{m,n}(L-1)] \quad (116) $$

$v = \max(j+2-L, 2-L)$, $\mu = \min(j+L-2, L-2)$, and the elements of matrix $S$ involve the triple products of the form $p_{i_u} p_{i_b} p_{i_c}$.

It is found by computation that the matrix $S$ has eigenvalues $2^{1-k}$, $k = 0, 1, \dots, L-2$, and that the multiplicity of the eigenvalue $2^{1-k}$ is $k+1$. Obviously, (114) indicates that the vector $v$ is the eigenvector of the matrix $S$ corresponding to the eigenvalue of $2^{(1-m-n)}$. However, since the rank lost of $S$ at this eigenvalue is $m+n+1$, it is not sufficient to determine the vector $v$ uniquely from (114). To have a system of linearly independent equations for solving the unknown elements of
the vector v, we need extra equations. As before, the required extra equations can be derived from the moment equation (63). By first multiplying both sides of (63) by $\phi^{(m)}(y-j)$ and $\phi^{(n)}(y-k)$, respectively, and then taking integration from $y=0$ to $y=x$, we have

$$ \sum_{k} k^n \Omega_{j,k}^{m,n}(x) = n! \Gamma_{j}^{m}(x) \quad (117) $$

$$ \sum_{j} j^m \Omega_{j,k}^{m,n}(x) = m! \Gamma_{k}^{m}(x) \quad (118) $$

For fixed integers $m$ and $n$, we can form an independent system of equations by deleting the equations of (114) which correspond to the unknowns $\Omega_{j,0}^{m,n}$, $j = 2, 3, \dots, 2-L+n$, and to the unknowns $\Omega_{k,k}^{m,n}$, $k = 1, 2, \dots, m$, and then adding the equations in (117) for $j = 0, 1, \dots, n$, and the equations in (118) for $k = 1, 2, \dots, m$, with $x = L-1$. It is noted that such a replacement of equations is by no means unique, and a rigorous validation of this rank remedy is difficult. In practice, however, the rank of the resultant system of equations can be checked by a numerical computation, and the obtained values of $\Omega_{j,k}^{m,n}(L-1)$ can be verified by using equations in (117) and (118) for $2-L \le j, k \le L-2$.

Having obtained the values of $\Omega_{j,k}^{m,n}(L-1)$ for $j, k \in \mathbb{Z}$, we can proceed to determine the values of $\Omega_{j,k}^{m,n}(x)$, for $x = 1, 2, \dots, L-2$ and $j, k \in \mathbb{Z}$. It can be shown from the properties of $\Omega_{j,k}^{m,m}$ in (110)-(112) that, for fixed integers $m$ and $n$, there are only $(L-2)^3$ independent unknowns among $\Omega_{j,k}^{m,n}(x)$, $x = 1, 2, \dots, L-2, j, k \in \mathbb{Z}$. These $(L-2)^3$ independent unknowns are given by $\Omega_{j,k}^{m,n}(x)$ with $x-L+2 \le j, k \le x-1$, and $x = 1, 2, \dots, L-2$. The linear system of $(L-2)^3$ equations constructed from (113) by letting $x-L+2 \le j, k \le x-1$, and $1 \le x \le L-2$ has the form

$$ w = Rw + f \quad (119) $$

where $w$ is the vector formed by the independent unknowns $\Omega_{j,k}^{m,n}(x)$, $x-L+2 \le x-1, x = 1, 2, \dots, L-2$, the vector $f$ contains elements $\Omega_{j,k}^{m,n}(L-1)$. We have found by numerical computation that the matrix $R$ has a unity eigenvalue with multiplicity $K$ given by

$$ K = \begin{cases} \sum_{i=1}^{m+n} i & \text{if } m+n \le L/2 \\ \sum_{i=L/2+1}^{m+n} \left( \frac{3L}{2} - 2i \right) + \frac{(L+2)L}{8} & \text{if } L/2 < m+n \le L-2 \end{cases} \quad (120) $$

In other words, (119) is singular and it cannot be used to determine the unknown vector $w$ without adding extra equations. To remedy the rank deficiency of the linear system (119), we again need the moment equations in (117) and (118). Instead of using the approach of equation replacement in Sections 4.2 and 4.3 or the approach of equation argument in Reference 13, we propose here to combine moment equations to some equations in (119). The moment equations in (117) for $j = 0, 1, \dots, x-1$ and $x = 1, 2, \dots, L-2$ are respectively combined with the equations in (119) that correspond to unknowns $\Omega_{j,j}^{m,n}(x)$, and the moment equations in (118) for $k = x-L+3, \dots, x-2$ and $x = 1, 2, \dots, L-2$ are respectively combined with the equations in (119) corresponding to the unknowns $\Omega_{x-1,k}^{m,n}(x)$. Although the rank lost of $R$ is $K$, we combine $2K$ moment equations with (119) to ensure the resultant system is of full rank. As illustrated by numerical computation, this manner of equation combination indeed allows us to form a system of linear independent equations for the unknown vector $w$. In addition, the moment equations in (117) and (118) can be used to verify the computational results.
<table><caption>Table I. Daubechies' wavelet filter coefficients for $L = 6$</caption><thead><tr><th>$i$</th><th>$p_i$</th></tr></thead><tbody><tr><td>0</td><td>0.470467207784E+00</td></tr><tr><td>1</td><td>0.114111691583E+01</td></tr><tr><td>2</td><td>0.650365000526E+00</td></tr><tr><td>3</td><td>-0.190934415568E+00</td></tr><tr><td>4</td><td>-0.120832208310E+00</td></tr><tr><td>5</td><td>0.498174997316E-01</td></tr></tbody></table>

<table><caption>Table II. The values of $\phi(x)$ $\phi^{(1)}(x)$ and $\theta_1(x)$</caption><thead><tr><th rowspan="2">$x$</th><th colspan="3">$L=6$</th></tr><tr><th>$\phi(x)$</th><th>$\phi^{(1)}(x)$</th><th>$\theta_1(x)$</th></tr></thead><tbody><tr><td>0.0</td><td>0.00000000E+00</td><td>0.00000000E+00</td><td>0.00000000E+00</td></tr><tr><td>0.5</td><td>0.60517847E+00</td><td>0.15416762E+01</td><td>0.14131460E+00</td></tr><tr><td>1.0</td><td>0.12863351E+01</td><td>0.16384523E+01</td><td>0.60074157E+00</td></tr><tr><td>1.5</td><td>0.44112248E+00</td><td>-0.24468283E+01</td><td>0.10529082E+01</td></tr><tr><td>2.0</td><td>-0.38583696E+00</td><td>-0.22327582E+01</td><td>0.10967114E+01</td></tr><tr><td>2.5</td><td>-0.14970591E-01</td><td>0.12730265E+01</td><td>0.98506614E+00</td></tr><tr><td>3.0</td><td>0.95267546E-01</td><td>0.55015936E+00</td><td>0.98548673E+00</td></tr><tr><td>3.5</td><td>-0.31541303E-01</td><td>-0.37227297E+00</td><td>0.10033183E+01</td></tr><tr><td>4.0</td><td>0.42343456E-02</td><td>0.44146491E-01</td><td>0.99965909E+00</td></tr><tr><td>4.5</td><td>0.21094451E-03</td><td>0.43985356E-02</td><td>0.99999151E+00</td></tr><tr><td>5.0</td><td>0.00000000E+00</td><td>0.00000000E+00</td><td>0.10000000E+01</td></tr></tbody></table>

## 5. NUMERICAL RESULTS

In this section we present computed values of functionals defined in (1)-(7) for the Six-coefficient Daubechies scaling function. For $L=6$, the coefficients $p_k$, $k=0,1,\dots,5$, are listed in Table I.

Table II lists the values of the scaling function $\phi(x)$, the derivative of $\phi(x)$, $\phi^{(1)}(x)$, and the integral of $\phi(x)$, $\theta_1(x)$, for $x = 0, 0.5, 1.0, 1.5, \dots, 5.0$. Table III lists the values of $M_k^m(x)$ for $k=0, m=1,2,3$ and $x=1,2,\dots,5$. The values of non-vanishing functionals of $\Gamma_k^1(x)$, $\Lambda_k^{1,1}(x)$, $\Upsilon_k^{1,1}(x)$ and $\Omega_{j,k}^{0,1}(x)$ at $x=1,2,3,4,5$, are listed in Tables IV-VII, respectively. It is noted that in computing the values in Tables II-VII, the 12 decimal digits of Daubechies' wavelet coefficients are used. We have verified the obtained values by using the moment equations that are not employed to remedy the rank deficiency, and the equation residues are all in the order of $10^{-9}$. Moreover, the computed three-term connection coefficients $\Omega_{j,k}^{0,1}(5)$ agree with those listed in Reference 13.

## 6. AN APPLICATION EXAMPLE

The non-linear parabolic equation

$$ \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} = \frac{1}{Re} \frac{\partial^2 u}{\partial x^2} \qquad (121) $$

is known as the Burgers equation,²³ which is one of the simplest partial differential equations describing both non-linear propagation and diffusive effects. It represents a first step in the hierarchy of approximations of the Navier-Stokes equation. Depending on the magnitude of the viscous term,
<table><caption>Table III. The values of $M_0^m(x)$ for $m = 1, 2, 3$</caption><thead><tr><th rowspan="2">x</th><th colspan="3">L=6</th></tr><tr><th>M<sub>0</sub><sup>1</sup>(x)</th><th>M<sub>0</sub><sup>2</sup>(x)</th><th>M<sub>0</sub><sup>3</sup>(x)</th></tr></thead><tbody><tr><td>1</td><td>0.40662888E+00</td><td>0.30755931E+00</td><td>0.24746348E+00</td></tr><tr><td>2</td><td>0.10138604E+01</td><td>0.10542175E+01</td><td>0.11618193E+01</td></tr><tr><td>3</td><td>0.77363652E+00</td><td>0.53674234E+00</td><td>0.52903635E-01</td></tr><tr><td>4</td><td>0.81603481E+00</td><td>0.66267182E+00</td><td>0.42355381E+00</td></tr><tr><td>5</td><td>0.81740117E+00</td><td>0.66814467E+00</td><td>0.44546004E+00</td></tr></tbody></table>

<table><caption>Table IV. The values of $\Gamma_k^1(x)$ for $L = 6$</caption><thead><tr><th>x</th><th>k</th><th>Γ<sub>k</sub><sup>1</sup>(x)</th><th>x</th><th>k</th><th>Γ<sub>k</sub><sup>1</sup>(x)</th><th>x</th><th>k</th><th>Γ<sub>k</sub><sup>1</sup>(x)</th></tr></thead><tbody><tr><td rowspan="4">1</td><td>-3</td><td>-0.96071829E-02</td><td rowspan="4">3</td><td>-1</td><td>-0.74488223E+00</td><td rowspan="4">5</td><td>-4</td><td>-0.34246575E-03</td></tr><tr><td>-2</td><td>0.24682915E+00</td><td>0</td><td>0.45379527E-02</td><td>-3</td><td>-0.14611872E-01</td></tr><tr><td>-1</td><td>-0.10642085E+01</td><td>1</td><td>0.73437630E+00</td><td>-2</td><td>0.14520548E+00</td></tr><tr><td>0</td><td>0.82732896E+00</td><td>2</td><td>-0.12428316E+00</td><td>-2</td><td>-0.74520548E+00</td></tr><tr><td rowspan="4">2</td><td>-2</td><td>0.14376046E+00</td><td rowspan="4">3</td><td>-1</td><td>0.89648439E-05</td><td rowspan="4">0</td><td>0</td><td>-0.11162444E-16</td></tr><tr><td>-1</td><td>-0.77113404E+00</td><td>0</td><td>0.74528563E+00</td><td>1</td><td>0.74520548E+00</td></tr><tr><td>0</td><td>0.74435080E-01</td><td>1</td><td>-0.14539422E+00</td><td>2</td><td>-0.14520548E+00</td></tr><tr><td>1</td><td>0.56789284E+00</td><td>2</td><td>0.15053970E-01</td><td>3</td><td>0.14611872E-01</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>4</td><td>0.34246575E-03</td></tr></tbody></table>

<table><caption>Table V. The values of $\Lambda_k^{1,1}(x)$ for $L = 6$</caption><thead><tr><th>x</th><th>k</th><th>Λ<sub>k</sub><sup>1,1</sup>(x)</th><th>x</th><th>k</th><th>Λ<sub>k</sub><sup>1,1</sup>(x)</th><th>x</th><th>k</th><th>Λ<sub>k</sub><sup>1,1</sup>(x)</th></tr></thead><tbody><tr><td rowspan="4">1</td><td>-3</td><td>0.61238997E-02</td><td rowspan="4">3</td><td>-1</td><td>-0.33884086E+00</td><td rowspan="4">5</td><td>-4</td><td>-0.39622225E-05</td></tr><tr><td>-2</td><td>0.16013363E+00</td><td>0</td><td>-0.48571617E+00</td><td>-3</td><td>-0.67621931E-03</td></tr><tr><td>-1</td><td>-0.74526083E+00</td><td>1</td><td>0.10513436E+01</td><td>-2</td><td>0.19212883E-01</td></tr><tr><td>0</td><td>0.57900506E+00</td><td>2</td><td>-0.26736083E+00</td><td>-2</td><td>-0.12104326E+00</td></tr><tr><td rowspan="4">2</td><td>-2</td><td>0.36873826E-01</td><td rowspan="4">4</td><td>0</td><td>-0.49996340E+00</td><td rowspan="4">0</td><td>0</td><td>0.10224223E+01</td></tr><tr><td>-1</td><td>-0.39613509E+00</td><td>1</td><td>0.10853692E+01</td><td>1</td><td>-0.12104326E+00</td></tr><tr><td>0</td><td>-0.33492810E+00</td><td>2</td><td>-0.33106420E+00</td><td>2</td><td>0.19212883E-01</td></tr><tr><td>1</td><td>0.69350982E+00</td><td>3</td><td>0.44926126E-01</td><td>3</td><td>-0.67621931E-03</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>4</td><td>-0.39622225E-05</td></tr></tbody></table>

<table><caption>Table VI. The values of $\Upsilon_k^{1,1}(x)$ for $L = 6$</caption><thead><tr><th>x</th><th>k</th><th>Υ<sub>k</sub><sup>1,1</sup>(x)</th><th>x</th><th>k</th><th>Υ<sub>k</sub><sup>1,1</sup>(x)</th><th>x</th><th>k</th><th>Υ<sub>k</sub><sup>1,1</sup>(x)</th></tr></thead><tbody><tr><td rowspan="6">1</td><td>-4</td><td>0.40660602E+00</td><td rowspan="6">3</td><td>-4</td><td>0.77363681E+00</td><td rowspan="6">4</td><td>1</td><td>-0.76349559E-01</td></tr><tr><td>-3</td><td>0.40660225E+00</td><td>-3</td><td>0.77388859E+00</td><td>2</td><td>0.58249586E-02</td></tr><tr><td>-2</td><td>0.40271563E+00</td><td>-2</td><td>0.76364511E+00</td><td>3</td><td>-0.14233365E-03</td></tr><tr><td>-1</td><td>0.44173660E+00</td><td>-1</td><td>0.83517453E+00</td><td>5</td><td>0.81740146E+00</td></tr><tr><td>0</td><td>0.10494506E+00</td><td>0</td><td>0.40735866E+00</td><td>-3</td><td>0.81749716E+00</td></tr><tr><td>2</td><td>0.10146510E+01</td><td>1</td><td>-0.48100476E-01</td><td>-2</td><td>0.80964464E+00</td></tr><tr><td rowspan="6">2</td><td>-4</td><td>0.10080259E+01</td><td rowspan="6">4</td><td>2</td><td>-0.16157050E-02</td><td rowspan="6">0</td><td>-1</td><td>0.87377189E+00</td></tr><tr><td>-3</td><td>0.10080259E+01</td><td>-4</td><td>0.81603510E+00</td><td>1</td><td>-0.46400092E+00</td></tr><tr><td>-2</td><td>0.10401149E+01</td><td>-3</td><td>0.81613080E+00</td><td>1</td><td>-0.13185340E+00</td></tr><tr><td>-1</td><td>0.10491431E+01</td><td>-2</td><td>0.80824884E+00</td><td>2</td><td>0.28977257E-01</td></tr><tr><td>0</td><td>0.56252533E+00</td><td>-1</td><td>0.87069022E+00</td><td>3</td><td>-0.11265687E-02</td></tr><tr><td>1</td><td>0.17523787E-01</td><td>0</td><td>0.50032473E+00</td><td>4</td><td>-0.81106612E-05</td></tr></tbody></table>
<table><caption>Table VII. The values of $\Omega_{j,k}^{0,1}(x)$ for $L = 6$</caption><thead><tr><th>x</th><th>j</th><th>k</th><th>Ω<sub>j,k</sub><sup>0,1</sup>(x)</th><th>x</th><th>j</th><th>k</th><th>Ω<sub>j,k</sub><sup>0,1</sup>(x)</th><th>x</th><th>j</th><th>k</th><th>Ω<sub>j,k</sub><sup>0,1</sup>(x)</th></tr></thead><tbody><tr><td rowspan="10">v 1</td><td>-3</td><td>-3</td><td>-0.71559141E-03</td><td rowspan="10">3</td><td>1</td><td>-4</td><td>0.31710703E+00</td><td rowspan="10">5</td><td>-1</td><td>-2</td><td>0.87462277E-01</td></tr><tr><td>-2</td><td>-2</td><td>0.13550268E-02</td><td>-3</td><td>-3</td><td>-0.13345958E+00</td><td>-1</td><td>-1</td><td>-0.24844590E+00</td></tr><tr><td>-1</td><td>-1</td><td>0.78711878E-05</td><td>2</td><td>-2</td><td>0.52989782E-02</td><td>0</td><td>0</td><td>0.31358775E+00</td></tr><tr><td>0</td><td>0</td><td>-0.61485101E-03</td><td>-1</td><td>-1</td><td>-0.11566110E-01</td><td>1</td><td>1</td><td>-0.12527255E+00</td></tr><tr><td>-2</td><td>-3</td><td>0.52989782E-02</td><td>0</td><td>0</td><td>0.25492658E-01</td><td>2</td><td>2</td><td>-0.61773222E-02</td></tr><tr><td>-2</td><td>-2</td><td>-0.11566110E-01</td><td>1</td><td>1</td><td>-0.19359467E-01</td><td>3</td><td>3</td><td>0.32901969E-04</td></tr><tr><td>-1</td><td>-3</td><td>0.25492658E-01</td><td>4</td><td>0</td><td>2</td><td>0.25336150E-07</td><td>0</td><td>-4</td><td>0.24592627E-05</td></tr><tr><td>0</td><td>0</td><td>-0.19359467E-01</td><td>-4</td><td>-4</td><td>0.49689213E+00</td><td>-3</td><td>-3</td><td>0.23763642E-02</td></tr><tr><td>-1</td><td>-3</td><td>-0.22426219E-01</td><td>-3</td><td>-3</td><td>0.32901910E-01</td><td>-2</td><td>-2</td><td>0.93548371E-01</td></tr><tr><td>-2</td><td>-2</td><td>0.60684150E-01</td><td>-2</td><td>-2</td><td>0.14542464E-02</td><td>-1</td><td>-1</td><td>-0.62717551E+00</td></tr><tr><td>-1</td><td>-1</td><td>-0.17563529E+00</td><td>1</td><td>-1</td><td>-0.24844521E+00</td><td>0</td><td>0</td><td>0.29322842E-10</td></tr><tr><td>0</td><td>0</td><td>0.13782300E+00</td><td>0</td><td>0</td><td>0.31359528E+00</td><td>1</td><td>1</td><td>0.49689180E+00</td></tr><tr><td rowspan="5">0</td><td>-3</td><td>-2</td><td>0.82360955E-02</td><td>-4</td><td>-4</td><td>-0.12529008E+00</td><td>2</td><td>2</td><td>0.32902689E-01</td></tr><tr><td>-2</td><td>-2</td><td>-0.19635413E+00</td><td>-3</td><td>-3</td><td>-0.61351029E-02</td><td>3</td><td>3</td><td>0.14522972E-02</td></tr><tr><td>-1</td><td>-1</td><td>-0.91407419E+00</td><td>2</td><td>-2</td><td>-0.16454414E-01</td><td>4</td><td>4</td><td>0.15241273E-05</td></tr><tr><td>0</td><td>0</td><td>0.70948150E+00</td><td>-1</td><td>-1</td><td>0.37778887E-01</td><td>1</td><td>-3</td><td>-0.44564810E-03</td></tr><tr><td>1</td><td>1</td><td>-0.13345958E+00</td><td>0</td><td>0</td><td>-0.46700711E-01</td><td>-2</td><td>-2</td><td>-0.20741515E-01</td></tr><tr><td rowspan="10">2</td><td>-2</td><td>-2</td><td>-0.16454414E-01</td><td>1</td><td>1</td><td>0.20316776E-01</td><td>-1</td><td>-1</td><td>0.87462280E-01</td></tr><tr><td>-1</td><td>-1</td><td>-0.37778887E-01</td><td>3</td><td>2</td><td>-0.71559141E-03</td><td>0</td><td>0</td><td>-0.24844590E+00</td></tr><tr><td>0</td><td>0</td><td>-0.46700711E-01</td><td>-4</td><td>-4</td><td>0.13550268E-02</td><td>1</td><td>1</td><td>0.31358775E+00</td></tr><tr><td>1</td><td>1</td><td>0.20316776E-01</td><td>-3</td><td>-3</td><td>0.78711888E-05</td><td>2</td><td>2</td><td>-0.12527255E+00</td></tr><tr><td>-1</td><td>-2</td><td>0.87355550E-01</td><td>-2</td><td>-2</td><td>-0.61485101E-03</td><td>3</td><td>3</td><td>0.61773222E-02</td></tr><tr><td>-1</td><td>-1</td><td>-0.24981584E+00</td><td>0</td><td>0</td><td>-0.76206365E-06</td><td>4</td><td>4</td><td>0.32901969E-04</td></tr><tr><td>0</td><td>0</td><td>0.31710703E+00</td><td>5</td><td>-4</td><td>-0.76206365E-06</td><td>0</td><td>0</td><td>0.13394071E-03</td></tr><tr><td>1</td><td>1</td><td>-0.13345958E+00</td><td>0</td><td>0</td><td>-0.44639911E-06</td><td>2</td><td>-2</td><td>0.49255216E-02</td></tr><tr><td>0</td><td>-2</td><td>0.94031790E-01</td><td>-4</td><td>-4</td><td>0.19493567E-05</td><td>-1</td><td>-1</td><td>-0.16451344E-01</td></tr><tr><td>-1</td><td>-1</td><td>-0.62003156E+00</td><td>-3</td><td>-3</td><td>0.48873854E-06</td><td>0</td><td>0</td><td>0.37810269E-01</td></tr><tr><td>0</td><td>0</td><td>-0.19146537E-01</td><td>-2</td><td>-2</td><td>-0.12296325E-05</td><td>1</td><td>1</td><td>0.37810269E-01</td></tr><tr><td>1</td><td>1</td><td>0.54276748E+00</td><td>-3</td><td>-1</td><td>-0.32455570E-04</td><td>2</td><td>2</td><td>-0.46774186E-01</td></tr><tr><td rowspan="5">3</td><td>-2</td><td>-2</td><td>-0.22426219E-01</td><td>0</td><td>0</td><td>-0.72614859E-03</td><td>3</td><td>3</td><td>0.20491688E-01</td></tr><tr><td>-1</td><td>-1</td><td>0.60684150E-01</td><td>1</td><td>1</td><td>0.12518005E-02</td><td>4</td><td>4</td><td>-0.13589007E-03</td></tr><tr><td>0</td><td>0</td><td>-0.17563529E+00</td><td>2</td><td>2</td><td>0.24982636E-03</td><td>3</td><td>-1</td><td>-0.32455570E-04</td></tr><tr><td>1</td><td>1</td><td>0.13782300E+00</td><td>-4</td><td>-4</td><td>-0.11881821E-02</td><td>0</td><td>0</td><td>-0.72614859E-03</td></tr><tr><td>0</td><td>0</td><td>-0.24844521E+00</td><td>-3</td><td>-3</td><td>0.44515936E-03</td><td>1</td><td>1</td><td>0.12518005E-02</td></tr><tr><td rowspan="5">0</td><td>-1</td><td>-1</td><td>0.31359528E+00</td><td>-2</td><td>-2</td><td>0.13394071E-03</td><td>2</td><td>2</td><td>0.24982636E-03</td></tr><tr><td>0</td><td>0</td><td>-0.12529008E+00</td><td>-1</td><td>-1</td><td>0.49255216E-02</td><td>3</td><td>3</td><td>-0.11881821E-02</td></tr><tr><td>2</td><td>2</td><td>-0.61351029E-02</td><td>0</td><td>0</td><td>-0.16451344E-01</td><td>4</td><td>4</td><td>0.44515936E-03</td></tr><tr><td>-1</td><td>-1</td><td>-0.62715213E+00</td><td>2</td><td>2</td><td>0.37810269E-01</td><td>0</td><td>0</td><td>-0.76206365E-06</td></tr><tr><td>0</td><td>0</td><td>0.28821311E-03</td><td>-4</td><td>-4</td><td>-0.46774186E-01</td><td>1</td><td>1</td><td>-0.44639911E-06</td></tr><tr><td rowspan="3">1</td><td>-1</td><td>-1</td><td>0.49612985E+00</td><td>-3</td><td>-3</td><td>0.20491688E-01</td><td>2</td><td>2</td><td>0.19493567E-05</td></tr><tr><td>2</td><td>2</td><td>0.34806875E-01</td><td>-2</td><td>-2</td><td>-0.13589007E-03</td><td>3</td><td>3</td><td>0.48873857E-06</td></tr><tr><td>0</td><td>0</td><td>-0.24981584E+00</td><td>-1</td><td>-1</td><td>-0.44564810E-03</td><td>4</td><td>4</td><td>-0.12296325E-05</td></tr></tbody></table>

the Burgers' equation behaves as an elliptic, parabolic or hyperbolic partial differential equation. Therefore, Burgers' equation has been widely used as a model equation for testing and comparing computational techniques.

When the viscous coefficient $1/Re$ in (121) is small, Burgers' equation involves phenomena that change rapidly from point to point. Problems of this type require mathematical representations that can respond to a local variation. The compactly supported wavelets derived by Daubechies¹ are good bases for 'local' problems. They have been used successfully to solve Burgers' equation with periodic boundary condition.²²,²⁴ Due to the lack of accurate finite-domain connect coefficients,
the wavelet bases have not been used to solve finite-domain problems with non-periodic boundary conditions. In the following, we apply the proposed interval wavelet-Galerkin approximation scheme to find the solution of the finite-domain Burgers' equation.

Consider the Burgers' equation (121) with the following boundary conditions:

$$u(-1, t) = 1 \qquad (122)$$

$$u(1, t) = 0 \qquad (123)$$

Let $u(x,t)$ be approximated by the $J$th-level wavelet series

$$u(x,t) \triangleq \sum_{k=2-L-2^J}^{2^J-1} u_{J,k}(t) \phi_{J,k}(x) = 2^{J/2} \sum_{k=2-L-2^J}^{2^J-1} u_{J,k}(t) \phi(2^J x - k) \qquad (124)$$

where $J > 0$. Substituting (124) into (121) and then applying the Galerkin discretization scheme, we have

$$\sum_{k=2-L-2^J}^{2^J-1} a_{l,k} \frac{du_{J,k}(t)}{dt} + \sum_{k=2-L-2^J}^{2^J-1} \sum_{l=2-L-2^J}^{2^J-1} b_{l,k,l} u_{J,k}(t) u_{J,l}(t) \\ = \frac{1}{Re} \sum_{k=2-L-2^J}^{2^J-1} c_{l,k} u_{J,k}(t), \quad l = 2-L-2^J, 3-L-2^J, \dots, 2^J-1 \qquad (125)$$

In the above differential equations, the coefficients $a_{l,k}$, $b_{l,k,l}$ and $c_{l,k}$ are given by

$$a_{l,k} = \int_{-1}^{1} \phi_{J,l}(x) \phi_{J,k}(x) dx = \Gamma_{k-l}^{0}(2^J - l) - \Gamma_{k-l}^{0}(-l - 2^J) \qquad (126)$$

$$b_{l,k,l} = \int_{-1}^{1} \phi_{J,l}(x) \phi_{J,k}(x) \frac{d\phi_{J,l}(x)}{dx} dx = 2^{3J/2} [\Omega_{k-l,i-l}^{0,1}(2^J - l) - \Omega_{k-l,i-l}^{0,1}(-l - 2^J)] \qquad (127)$$

$$c_{l,k} = \int_{-1}^{1} \phi_{J,l}(x) \frac{d^2 \phi_{J,k}(x)}{dx^2} dx = 2^{2J} [\Gamma_{k-l}^{2}(2^J - l) - \Gamma_{k-l}^{2}(-l - 2^J)] \qquad (128)$$

where $\Gamma_k^n(x)$ and $\Omega_{J,k}^{m,n}(x)$ are defined by (4) and (5), respectively. For $u(x,t)$ to satisfy the boundary conditions (122) and (123), the expansion coefficients $u_{J,2-L-2^J}(t)$ and $u_{J,2^J-1}(t)$ must satisfy the following relations:

$$u_{J,2-L-2^J}(t) = \frac{1}{\phi(L-2)} [2^{-J/2} - \sum_{i=1}^{L-3} \phi(i) u_{J,-i-2^J}(t)] \qquad (129)$$

$$u_{J,2^J-1}(t) = \frac{-1}{\phi(1)} \sum_{i=2}^{L-2} \phi(i) u_{J,2^J-i}(t) \qquad (130)$$

As a result, there are $2^{J+1} + L - 4$ independent expansion coefficients $u_{J,k}(t)$, $k = 3 - L - 2^J$, $4 - L - 2^J, \dots, 2^J - 2$, that are described by the differential equations in (125) with $u_{J,2-L-2^J}(t)$ and $u_{J,2^J-1}(t)$ being replaced by (129) and (130), respectively.

The initial conditions for the differential equations in (125) are derived from the initial condition $u(x,0)$ of the problem. Here, we consider the discontinuous initial condition

$$u(x, 0) = \begin{cases} 1 & \text{for } -1 \le x \le 0 \\ 0 & \text{for } 0 < x \le 1 \end{cases} \qquad (131)$$

and the continuous initial condition

$$u(x, 0) = \frac{1}{2}(1-x) \qquad (132)$$
For the case of the discontinuous initial condition specified by (131), the initial conditions $u_{j,k}(0)$ satisfy

$$ 
\begin{aligned}
\sum_{k=2-L-2^J}^{2^J-1} a_{l,k} u_{j,k}(0) &= \int_{-1}^{0} \phi_{j,l}(x) \, dx \\
&= \frac{1}{2^{J/2}} [\theta_1(-l) - \theta_1(-l - 2^J)] \\
&l = 2 - L - 2^J, 3 - L - 2^J, \dots, 2^J - 1
\end{aligned} 
\quad (133) $$

From these equations and relations in (129) and (130), we can obtain the initial conditions $u_{j,k}(0)$, $k = 3-L-2^J, 4-L-2^J, \dots, 2^J-2$. Similarly, for the case of the continuous initial condition given by (132), the initial conditions $u_{j,k}(0)$ satisfy

$$ 
\begin{aligned}
\sum_{k=2-L-2^J}^{2^J-1} a_{l,k} u_{j,k}(0) &= \frac{1}{2} \int_{-1}^{0} (1-x) \phi_{j,l}(x) \, dx \\
&= \frac{1}{2^{J/2+1}} [\theta_1(2^J-l) - \theta_1(-l-2^J) - 2^{-J} M_l^1(2^J) + 2^{-J} M_l^1(-2^J)] \\
&l = 3 - L - 2^J, 4 - L - 2^J, \dots, 2^J - 2
\end{aligned} 
\quad (134) $$

where $a_{l,k}$, $\theta_n(x)$ and $M_k^n(x)$ are defined in (126), (2) and (3), respectively. Hence, the initial conditions $u_{j,k}(0)$, $k = 2-L-2^J, 4-L-2^J, \dots, 2^J-1$, can be obtained from solving the linear equations of (134) with the additional relations of (129) and (130).

Having obtained the differential equations and initial conditions for $u_{j,k}(t)$, we can utilize a numerical integration scheme to find the expansion coefficients $u_{j,k}(t)$. In actual computation, we have invoked the subroutine DIVPAG in IMSL mathematical library²⁵ to solve the initial-value problems. In calling this subroutine, the tolerance parameter and the initial time step are respectively set to $10^{-6}$ and $0.1$, while all the other parameters are set to their default values. By choosing $J = 5$ and $L = 6$, the time evolutions of the waveforms of the wavelet-series solution of Burgers' equation are plotted in Figures 1-4. Figures 1 and 2 show the solutions of Burgers' equation with a linear initial condition for $Re = 1$ and $Re = 100$, respectively. The results shown in Figures 3 and 4 are the solutions evolved from a discontinuous initial condition. From Figures 1 and 3, it can be observed that for a low Reynolds number or high viscous coefficient, the solution to Burgers' equation evolves to a smooth steady state, irrespective of continuous or discontinuous initial conditions. As can be seen from Figures 2 and 4, the solution evolutions of Burgers' equation with low viscous coefficients are quite different from those with high viscous coefficients. For the discontinuous initial condition, the solution of Burgers' equation with $Re = 100$ is a shock wave. Before the shock wave touches the right boundary, its penetration is just like in an infinite domain. However, when the shock wave touches the right boundary, its action like the billow beats the reef and generates the spray. In this situation, we cannot find the steady-state solution since the problem becomes singular.

It is noted that there is no exact solution for the finite-domain Burgers' equation. In order to verify the numerical solution, it is often compared with that of the infinite domain. For the purpose of demonstrating that the solutions shown in Figures 1-4 are true approximate solutions of the finite-domain Burgers' equation, we compare in Table VIII the wavelet-series solutions $u(x, 0.92)$ with the solution $u_p(x, 0.92)$ obtained by Galerkin approximation using a ninth-degree polynomial as a trial solution²³ and the solution $u_\infty(x, 0.92)$ of the infinite-domain Burgers' equation. Obviously, the local bases used in this paper can respond well to the sharp change of the shock wave.
<figure><img src="image_2.png" /><figcaption>Figure 1. Numerical solutions of Burgers' equation with $Re = 1$ and linear initial condition (132)</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 2. Numerical solutions of Burgers' equation with $Re = 100$ and linear initial condition (132)</figcaption></figure>

## 7. CONCLUSIONS

The wavelet-Galerkin method has been shown in the literature to be a powerful tool for the numerical solution of partial differential equations. The computation of wavelet-Galerkin approximation relies heavily on the evaluation of connection coefficients, which are integrals with their integrands
<figure><img src="image_4.png" /><figcaption>Figure 3. Numerical solutions of Burgers' equation with $Re = 1$ and discontinuous initial condition (131)</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4. Numerical solutions of Burgers' equation with $Re = 100$ and discontinuous initial condition (131)</figcaption></figure>

being the product of wavelet bases and their derivatives or integrals. We have described the algorithms for the exact evaluations of connection coefficients for Daubechies' compactly supported orthonormal wavelets on a bounded interval. These algorithms extend the application scope of the wavelet-Galerkin method to finite-domain problems. They also play an essential role in applying the newly developed interval wavelets$^{19,20}$ to the numerical solution of partial differential equations. This role has been demonstrated by a numerical example.
<table><caption>Table VIII. Numerical solutions of finite-domain Burgers' equation</caption><caption>Re = 100, t = 0.92</caption><thead><tr><th>x</th><th>W - G, J = 4</th><th>W - G, J = 5</th><th>PA<sup>23</sup></th><th>ES<sup>23</sup></th></tr></thead><tbody><tr><td>-1.0</td><td>1.0000</td><td>1.0000</td><td>1.0000</td><td>1.0000</td></tr><tr><td>-0.9</td><td>0.8508</td><td>0.8722</td><td>0.9956</td><td>1.0000</td></tr><tr><td>-0.8</td><td>0.8809</td><td>0.8727</td><td>1.0456</td><td>1.0000</td></tr><tr><td>-0.7</td><td>0.8810</td><td>0.8727</td><td>1.0672</td><td>1.0000</td></tr><tr><td>-0.6</td><td>0.8811</td><td>0.8728</td><td>1.0402</td><td>1.0000</td></tr><tr><td>-0.5</td><td>0.8817</td><td>0.8737</td><td>0.9831</td><td>1.0000</td></tr><tr><td>-0.4</td><td>0.8860</td><td>0.8781</td><td>0.9303</td><td>1.0000</td></tr><tr><td>-0.3</td><td>0.8981</td><td>0.8913</td><td>0.9128</td><td>1.0000</td></tr><tr><td>-0.2</td><td>0.9204</td><td>0.9157</td><td>0.9444</td><td>1.0000</td></tr><tr><td>-0.1</td><td>0.9487</td><td>0.9465</td><td>1.0159</td><td>1.0000</td></tr><tr><td>0.0</td><td>0.9752</td><td>0.9742</td><td>1.0963</td><td>1.0000</td></tr><tr><td>0.1</td><td>0.9920</td><td>0.9914</td><td>1.1411</td><td>1.0000</td></tr><tr><td>0.2</td><td>0.9989</td><td>0.9983</td><td>1.1057</td><td>1.0000</td></tr><tr><td>0.3</td><td>0.9948</td><td>0.9980</td><td>0.9613</td><td>0.9998</td></tr><tr><td>0.4</td><td>0.9260</td><td>0.9610</td><td>0.7099</td><td>0.9764</td></tr><tr><td>0.5</td><td>0.1090</td><td>0.1156</td><td>0.3933</td><td>0.1861</td></tr><tr><td>0.6</td><td>-0.0270</td><td>0.0029</td><td>0.0905</td><td>0.0015</td></tr><tr><td>0.7</td><td>0.0063</td><td>0.0000</td><td>-0.1017</td><td>0.0000</td></tr><tr><td>0.8</td><td>0.0000</td><td>0.0000</td><td>-0.1154</td><td>0.0000</td></tr><tr><td>0.9</td><td>0.0000</td><td>0.0000</td><td>0.0091</td><td>0.0000</td></tr><tr><td>1.0</td><td>0.0000</td><td>0.0000</td><td>0.0000</td><td>0.0000</td></tr></tbody></table>

W-G: wavelet-Galerkin; PA: polynomial approximation; ES: exact solution

ACKNOWLEDGEMENT

This work was supported by the National Science Council of the Republic of China under Grant NSC-82-0402-E-006-247.

REFERENCES

1. I. Daubechies, 'Orthonormal bases of compactly supported wavelets', Commun. Pure Appl. Math., 41, 909-996 (1988).

2. R. Glowinski, W. M. Lawton, M. Ravachol and E. Tenenbaum, 'Wavelet solution of linear and nonlinear elliptic, parabolic and hyperbolic problems in one space dimension', in *Computing Methods in Applied Sciences and Engineering*, Chapter 4, 1990, pp. 55-120.

3. C. Brislawn and I. G. Rosen, 'Wavelet based approximation in the optimal control of distributed parameter systems', Numer. Functional Anal. Optim., 12, 33-77 (1991).

4. E. Bacry, S. Mallat and G. Papanicolaou, 'A wavelet based space-time numerical method for partial differential equations', Math. Modeling Numer. Anal., 26, 793-834 (1992).

5. C. Zhiqian and E. Weinan, 'Hierarchical method for elliptic problems using wavelet', Commun. Appl. Numer. Methods, 8, 819-825 (1992).

6. K. Amaratunga, J. R. Williamm, S. Qian and S. Yokoyama, 'Wavelet based hierarchical solutions of partial differential equations', Proc. Complas III, 3rd Int. Conf. on Computational Plasticity, Fundamentals and Applications, Barcelona, Spain, April 1992.

7. S. Jaffard and Ph. Laurecot, 'Orthonormal wavelets, analysis of operators, and applications to numerical analysis', in C. K. Chui (ed.), *Wavelets—A Tutorial in Theory and Applications*, 1992, pp. 543-601.

8. J. C. Xu and W. C. Shann, 'Galerkin-wavelet methods for two-point boundary value problems', Numer. Math., 63, 123-144 (1992).

9. S. Qian and J. Weiss, 'Wavelets and the numerical solution of boundary value problems', Appl. Math. Lett., 6, 47-52 (1993).

10. S. Qian and J. Weiss, 'Wavelets and the numerical solution of partial differential equations', J. Comput. Phys., 106, 155-175 (1993).

11. K. Amaratunga, J. R. Williamm, S. Qian and J. Weiss, 'Wavelet-Galerkin solutions for one dimensional partial differential equations', Int. j. numer. methods eng., 37, 2703-2716 (1994).
12. G. Beylkin, 'On the representation of operators in bases of compactly supported wavelets', *SIAM J. Numer. Anal.*, **29**, 507-537 (1993).

13. A. Latto, H. L. Resnikoff and E. Tenenbaum, 'The evaluation of connection coefficients of compactly supported wavelets', in Y. Maday (ed.), *Proc. French-USA Workshop on Wavelets and Turbulence*, Princeton University, June 1991, New York, Springer, 1994.

14. W. Dahmen and C. A. Micchelli, 'Using the refinement equation for evaluating integrals of wavelets', *SIAM J. Numer. Anal.*, **30**, 507-537 (1993).

15. L. Jameson, 'On the wavelet based differentiation matrix', *J. Scientific Comput.*, **8**, 267-305 (1993).

16. K. McCormick and R. O. Wells, Jr., 'Wavelet calculus and finite difference operators', *Math. Comput.*, **63**, 155-173 (1994).

17. M. Q. Chen, C. Hwang and Y. P. Shih, 'A wavelet-Galerkin method for solving population balance equations', *Comput. Chem. Eng.*, in revision.

18. M. Q. Chen, C. Hwang and Y. P. Shih, 'A wavelet-Galerkin method for solving Stefan problems', *J. Chinese Inst. Chem. Engrs.*, in revision.

19. A. Cohen, I. Daubechies and P. Vial, 'Wavelets on the interval and fast wavelet transform', *Appl. Comput. Harmonic Anal.*, **1**, 355-367 (1994).

20. C. K. Chui and E. Quak, 'Wavelets on a bounded interval', in D. Braess and L. L. Schumaker (eds.), *Numerical Methods in Approximation Theory*, Birkhauser, Basel, 1992, pp. 53-75.

21. S. Mallat, 'Multiresolution approximation and wavelets', *Trans. Amer. Math. Soc.*, **315**, 69-88 (1989).

22. A. Latto and E. Tenenbaum, 'Les ondelettes a support compact et la solution numerique de l'équation de Burgers', *C. R. Acad. Sci. Paris*, **311**, 903-909 (1990).

23. C. A. J. Fletcher, 'Burgers' equation: a model for all reason', in J. Noye (ed.), *Numerical Solutions of Partial Differential Equations*, North-Holland, New York, 1982.

24. J. Liandrat and Ph. Tchamitchian, 'Resolution of the 1D regularized Burgers equation using a spatial wavelet approximation', *NASA Contractor Report 187480: ICASE Report No. 90-83*, December 1990, pp. 36.

25. *IMSL: User's Manual for IMSL MathLibrary: Fortran Subroutines for Mathematical Applications*, Version 1.1, 1989.