1 Introduction

Numerous physical phenomena, e.g., sound, heat, electrostatics, electrodynamics, elasticity, and quantum mechanics, can be quantitatively described by nonlinear partial differential equations(PDEs), together with a set of additional constraints of initial and boundary conditions, which form initial boundary value problems(IBVPs). The IBVPs play a very important role in science and engineering, and their solution methods have been extensively studied.

Since the introduction by Grossmann and Morlet in 1984, the wavelet theory has been broadly applied in many areas of science and engineering^[1]. One of these successful applications is the numerical solution of PDEs^[2]. These studies have been devoted to developing the space discretization methods for nonlinear boundary value problems(BVPs)^[3]. For example, a modified wavelet Galerkin method has recently been proposed to solve nonlinear BVPs, including the large deflection bending of circular plates^[4], rectangular plates^[5], and Bratu equations^[3], showing much better accuracy and efficiency^[3-5] comparing with conventional ones. The accuracy level is almost independent of the nonlinearity of the equations^{[3, 6]}. When IBVPs are solved with the wavelet-based spatial discretization method, these equations are usually transformed into a set of ordinary differential equations(ODEs) with unknown time-dependent wavelet coefficients or nodal approximations^[7]. These ODEs form the initial value problems(IVPs). To efficiently solve the IVPs, several studies have been devoted to developing the wavelet-based time discretization schemes^[8]. These attempts have included the solution of the vibrations of multi-degrees of freedom(MDOF) systems^[7-10], Burgers equations^[11-12], nonlinear circuit simulation^[13], and nonlinear stiff systems^[14-15], etc. However, all these studies have just treated the IVPs as pseudo-BVPs^[11] by specifying a certain time interval. These pseudo-BVPs can then be readily solved by the aforementioned spatial discretization methods, e.g., the wavelet collocation^{[7, 10-12]} and the wavelet Galerkin method^[16]. The disadvantage of this treatment on time is that the dimension of the resulting algebraic equations can be significantly increased^[7] comparing with most methods in terms of step-by-step integration. A larger time interval usually results in the requirement of much more additional unknown coefficients for wavelet expansion. Therefore, the numerical solutions have been investigated only within a significantly limited time range with most of these wavelet methods^[9]. One main reason for using such time discretization strategies is hypothesized as to control the global error in time domain^{[7, 11]}. However, the simulation of the time response of a complex dynamic system usually requires the balance among the accuracy, stability, and computational cost^[17-18].

In spite of the above progress, the development of a simultaneous space-time wavelet method with adjustable accuracy and high efficiency for general nonlinear IBVPs is still far from complete. One aspect relevant to the present paper is that a wavelet-based, step-by-step, and precisely time-integrating scheme is still lacking. Therefore, in this study, after briefly introducing the properties and construction procedure of a typical wavelet, i.e., the Coiflet, we will first develop a novel scheme of wavelet approximation of square-integrable functions with an adjustable order of accuracy, and then use this approximation scheme to develop a wavelet-based method for the solution of general nonlinear IBVPs. In this method, we will only use one type of Coiflet bases to make the spatial and temporal discretization so that the conflict arising with different approximations in the spatial and temporal dimensions can be effectively avoided. Numerical examples, including the nonlinear vibration of a Duffing oscillator, the Burgers equation^[19-23] in fluid mechanics, and the Klein Gordon equation^[24-27] in quantum mechanics, are considered to demonstrate the accuracy and efficiency of the wavelet-based methods.

2 Brief introduction to Coiflets

Coiflets are a class of compactly supported orthogonal wavelets, whose wavelet function ψ(x) and scaling function ϕ(x) satisfy the vanishing and shifted vanishing moment properties^[28-29], respectively. More specifically, for the Coiflet with a compact support of [0, 3N-1], where N is a positive even integer, we have^[29]

(1)

(2)

(3)

(4)

(5)

where p_k is a set of low-pass filter coefficients of the wavelet system, and is the nth-order moment of the scaling function. As shown by Wang^[29], these properties can be derived to associate with the filter coefficients p_k as follows:

(6)

(7)

(8)

(9)

Equation(6) provides one equation for the determination of the filter coefficients, Eq.(7) provides 3N/2 equations, Eq.(8) provides N independent equations, and Eq.(9) provides N/2 equations. Thus, Eqs.(6)-(9) give 3N independent equations. Once the 3N filter coefficients are obtained, the dyadic values of the scaling wavelet functions and their derivatives and integrals can be easily calculated, following a standard procedure^[29]. Tables 1 and 2 give the values of the scaling function and its derivatives and integrals at the integer points, respectively.

For a function , we denote

(10)

where m is an integer called the resolution level, and ϕ_{m, k}(x)=2^m/2ϕ(2^mx-k). The sequence of subspaces, , satisfies

(11)

(12)

(13)

(14)

(15)

(16)

(17)

Let f(x)≈ P^mf(x). From Eq.(10), we have

(18)

The expansion coefficients in Eq.(18) can be calculated as follows:

(19)

which have an accuracy with the order N in terms of Eq.(17) and Ref. [29].

Thus, the function f(x) can be represented as follows^[6]:

(20)

where

Then, we have f(x)=P^mf(x), which is any polynomial with the order upto N-1 according to Eqs.(17) and(19).

3 Approximation of interval-bounded functions

If a function is defined on the interval [a, b] with 2^ma and 2^mb as integers, Eq.(20) can be rewritten as follows:

(21)

where ||f(x)-P^mf(x)||₂=O(2^-mN)^[6]. It is noted that x_k=k/2^m may be located outside the interval [a, b], and the function f(x) has no definition in these locations. Therefore, the technique of boundary extension is usually adopted to resolve this problem. In this study, similar to Ref. [29], we consider the power series expansion of the function at each boundary as follows:

(22)

where

(23)

Similarly,

(24)

Equations(23) and(24) mean that, if the boundary derivatives are known, we utilize them in the boundary extension treatment in Eq.(22); otherwise, we will use linear combination of function values at the discrete points in the interval [a, b] to approximate the boundary derivatives. The latter needs to determine the coefficients ζ_{b, i, k}, which have been derived based on the technique of numerical differences^[29]. However, such a treatment is not accurate enough. For example, if f(x)=xⁿ(x∈[a, b], n=1, 2, ..., N-1), we expect that Eq.(22) recovers f(x)=xⁿ(x∈(-∞, +∞)). Of course, this is true if all the boundary derivatives are known and are applied in Eq.(22). However, even if only one of the boundary derivatives is unknown and is given by the traditional numerical difference as shown in Ref. [29], it is not true. Our question is how we can derive the coefficients ζ_{b, i, k} and ζ_{a, i, k} so that Eq.(22) can recover f(x)=xⁿ(x∈(-∞, +∞)) no matter the boundary derivatives are known or not. In the following, we will solve this problem by considering the wavelet approximation of function derivatives. To do this, we substitute Eq.(22) into Eq.(21), take derivatives on both sides of Eq.(21) with respect to x, and consider x=b. Eventually, we have

(25)

where m is assumed to be sufficiently large so that 2^mb-3N+2+M₁>2^ma. This restriction on m is to ensure that there exists a sufficient number of nodal points to approximate the boundary derivatives. We let the right-side of Eq.(25) equal F_b^(j). Then, using F_b^(j) to replace f_b^(j) in the left-side of Eq.(25), we can obtain

(26)

This replacement ensures the consistency between the difference-like approximation of the boundary derivatives as shown in the second line of the right-side of Eq.(23) and the wavelet approximation of the boundary derivatives in Eq.(25). In addition, we can easily verify that, Eq.(26) is exact if f(x) is any polynomial with the order upto N-1, and F_b^(j)=f_b^(j) (j=1, 2, ..., N-1), owing to the property of the wavelet shown in Eq.(17). In order to determine the coefficients ζ_{b, i, k}, we consider Eq.(23), insert into Eq.(26), and take α₂=3N-2-M. Then, we have

(27)

Equation(27) can be rewritten in a matrix form as follows:

(28)

where I is an N × N unit matrix, and

In the above equations, i, j=0, 1, 2, ..., N-1, and k=0, 1, 2, ..., α₂.

Since Eq.(28) is satisfied for any F, from Eq.(28), we have

Theorem 1  If P₁={ζ_{b, i, k}}=(I-B₁)^-1A₁ (i=0, 1, 2, ..., N-1, k=0, 1, 2, ..., α₂) for f(x)=xⁿ(n=0, 1, 2, ..., N-1),

The proof of this theorem is presented in Appendix A. Similarly, we have

Then, we have the following theorem.

Theorem 2  If P₀={ζ_{a, i, k}}=(I-B₀)^-1A₀ (i=0, 1, 2, ..., N-1, k=0, 1, 2, ..., α₁) for f(x)=xⁿ(n=0, 1, 2, ..., N-1),

The proof of this theorem is similar to that in Appendix A. The calculation process of the coefficients ζ_{a, i, k} is shown in Appendix B.

Substituting Eqs.(23) and(24) into Eq.(22) and taking x=k/2^m, respectively, we have

(29)

in which^[30]

(30)

In the above equation,

The introduction of β_L and β_R merely aims to assign the boundary conditions to the function f(x). For example, when dⁱf(a) /dxⁱ=0 and other boundary derivatives are unknown, we can simply set β_{L, i}=0 and all the other elements of β_L and β_R equal 1.

Using Eq.(21) to approximate Eq.(29) and making further rearrangement as shown in Ref. [30] yield

(31)

where

(32)

Then, we have

The local accuracy at the boundary regions is expected to be significantly improved comparing with the previous boundary extension treatment based on the numerical difference in Ref. [29].

As an example, we consider the function f(x)=e^x(x∈[0, 1]) to be approximated by different series expansions such as the Fourier series, the Chebyshev orthogonal polynomials, the previous treatment in Ref. [29], and the present wavelet method as shown in Eq.(31) with N=6 and M₁=7. The absolute errors of these approximations are plotted in Fig. 1 when the number of terms in each series is taken as 16. We use to characterize the approximation error, where is the approximation of f(x). Figure 1 shows that the approximation with Eq.(31) is much better than those with the other two methods. The absolute errors obtained by the wavelet method, as a function of the effective "grid size", 2^-m, with different terms of approximation are plotted in Fig. 2. We can see from Fig. 2 that, the errors are scaled as 2^-6m, i.e., the wavelet approximation in terms of Eq.(31) has an order N, which is 6 in this example.

4 Wavelet time-integrating method(WTIM) for the solution of IVPs

In solving IBVPs, one usually needs to spatially discretize the original equations into a set of initial-valued ODEs. These equations can be generally written as follows^[31]:

(33)

where y=(y₁, y₂, ..., y_s)^T is an unknown vector function of the state variables, and f=(f₁, f₂, ..., f_s)^T is a given nonlinear vector function.

To derive the solution method of Eq.(33), we integrate Eq.(33). Then, we have

(34)

By setting a=-∞ and b=j/2^m based on Eqs.(31) and(32), we directly approximate every entry of the function vector f on the time interval as follows:

where j is a positive integer, t_k=k/2^m, and y_k =y(t_k). Based on the above equation, we further consider the approximation of f_i(t, y) for t∈[(j-1)/2^m, j/2^m], and obtain

(35)

Substituting Eq.(35) into Eq.(34) and taking t=j/2^m, we have

(36)

Rearranging Eq.(36) yields

The above equation can be rewritten as follows:

(37)

where

(see Appendix C for the reason why this integral depends only on j-k). Φ_{m, k}(τ) is given in Eq.(32) with a=-∞ and b=j/2^m, and all the elements of β_L and β_R are equal to 1. In solving Eq.(33) with Eq.(37) at the first α₂ steps, we can either set all y⁽ⁱ⁾(0)(i =1, 2, ..., N-1) as unknowns or obtain them from Eqs.(33) or(24) in advance. Then, the values of y_k (k=-α₂, -α₂+1, ..., -1) in Eq.(37) can be obtained through the power series expansion

For the linear ODEs, Eq.(33) can be simplified to

(38)

Applying Eq.(37) to Eq.(38) yields

(39)

where I is a unit matrix, and

(40)

The error estimation of Eq.(40) can be easily obtained based on those of Eqs.(31) and(32) and the relevant results in Ref. [33]. Considering Eq.(40) and defining the approximate solution y_j as follows:

we have

In the case that, the elements of A(t_j) have a finite upper bound, so that the elements of hΓ₀A(t_j) are on the order of O(h) and further ||I-hΓ₀A(t_j)||₂=O(1)^[33]. Then, we have the following error estimation:

which shows that the proposed solution method has an order of convergence N, which can be any positive even integer, as long as the Coiflet with N-1 vanishing moments is adopted.

5 Simultaneous space-time wavelet method(SSTWM) for IBVPs

We consider a class of nonlinear IBVPs as follows:

(41)

where u(x, t) is the unknown function, L⁰ and L¹ denote the differential operators, N(u(x, t)) represents any nonlinear composite function of u(x, t), and i_s, j_s, and s are non-negative integers. For the equations with inhomogeneous boundary conditions, a simple transformation can reduce them into those with homogeneous ones as shown in Ref. [6].

For the solution of Eq.(41), we use the proposed wavelet expansion of Eqs.(31) and(32) to approximate the various functions in this equation, and then adopt the Galerkin method to perform the spatial discretization. To satisfy the boundary conditions in Eq.(41), we set β_{L, is}=β_{L, js}=0 and all the other elements of β_L and β_R to be equal to 1. The resulting modified scaling function basis of Φ_{n, k}(x) in Eq.(32) is specified accordingly, which is re-denoted as Φ_{n, k}(x). The functions and nonlinear terms with x∈ [0, 1] in Eq.(41) can then be approximated by

(42)

(43)

(44)

where Φ_{n, k}(x) implies that the boundary derivatives of the corresponding functions are not specified, or all the elements of β_L and β_R are equal to 1.

Substituting Eqs.(42)-(44) into Eq.(41) yields

(45)

Applying the Galerkin method to Eq.(45), we have

(46)

where

The entries in these matrices, i.e., a_lk, b_lk, c_lk, and e_lk, can be obtained exactly through the procedure suggested by Wang^[29]. Then, Eq.(46), as an IVP, can be readily solved with the proposed WTIM. However, most interestingly, we note that, in the semi-discretization system(46) of the nonlinear IBVP problem, A, B, C, E, and G are all constant matrices, which are completely independent of the unknown vector U(t) and time t. Thus, in the subsequent time integration for solving Eq.(46), the matrix generated in the spatial discretization does not need updating, which implies that a fully decoupling between the spatial and temporal discretizations is achieved in the present wavelet formulation.

6 Numerical examples

Four different dynamic problems are considered to demonstrate the accuracy and efficiency of the proposed wavelet-based methods, including the free vibration equation of a single oscillator, the nonlinear vibration equation of a Duffing oscillator, the Burgers equation, and the Klein-Gordon equation in quantum mechanics.

Problem 1 Free vibration equation of a single oscillator

We consider the free vibration of a single oscillator with the dimensionless governing equation as follows:

(47)

We take ξ=4π² and the initial condition as

The analytical solution of Eq.(47) can be immediately obtained as x(t)=cos (2πt). We define

Then, Eq.(47) can be changed into

(48)

Figure 3 shows the displacement response of this oscillator, obtained with the proposed WTIM with the Coilfet of N=6 and M₁=7, the central difference method, the Newmark-β method, and the Wilson-θ method, respectively. It can be seen from Fig. 3 that, the proposed wavelet method has the smallest period elongation among all the classical methods. Figure 4 plots the absolute error of the WTIM as a function of the step size h, which shows that the WTIM has an accuracy order of 6.5, i.e., N+1/2. It also shows that, the accuracy order of the WTIM is unchanged, whether we set all y⁽ⁱ⁾(0) (i=1, 2, ..., N-1) as unknowns or obtain them from Eqs.(33) or (24) in advance.

Problem 2  Free vibration equation of the Duffing oscillator

The free vibration of the Duffing oscillator represents the vibration of an unforced pendulum with nonlinear restoring force. We choose this problem to test the proposed WTIM because the analytical solution of this nonlinear equation can be obtained by the following Jacobi-elliptic functions^[18]:

(49)

If ω_n=1, the analytical solution is where k²=η/(2η+2), and cn(·) is the Jacobi-elliptic function^{[18, 32]}. Denote

Then, Eq.(49) can be rewritten as follows:

(50)

We use the proposed WTIM to solve Eq.(50), where we choose the Coiflets with N=4, M₁=7 and N=6, M₁=7, respectively. Figures 5 and 6 show the absolute errors of x(t) as a function of the step size h at different moments and under different η, respectively. It can be seen from Figs. 5 and 6 that, the errors decay as h⁶ for N=6 and as h⁴ for N=4. This example again demonstrates that, the convergence order of this proposed wavelet method is N. This accuracy remains unchanged when the intensity of nonlinearity of Eq.(50), characterized by the parameter η, increases.

Problem 3  Burgers' equation

Burgers' equation is a fundamental and typical nonlinear partial differential equation in fluid mechanics, and has been widely used to describe the wave propagation phenomena in acoustics and hydrodynamics. In this study, we consider the following one-dimensional case:

(51)

where the initial boundary condition is u(x, 0)=2π sin(π x)/(100Re+Recos(πx)) for Case a^{[19, 23]} and u(x, 0)=sin (πx)for Case b^[20-22]. The exact solutions for these two cases can be found in Refs. [19]-[23]. Equation(51) is solved with the proposed simultaneous space-time wavelet method, where the spatial resolution level is taken as n=4. The time step size h and the Reynolds number are taken as 1/64 and 200 for Case a and 1/256 and 10 for Case b. By contrast, many researchers have solved Eq.(51) with different numerical methods, such as the differential quadrature method^[19], the methods based on the multiquadratic quasi-interpolation^[20] and spline interpolation^[21], the quadratic B-spline finite element method^[22], and the cubic B-splines collocation method^[23]. The time derivative is usually discretized with the forward difference methods^[19-21] or integration method based on the Runge-Kutta method^[23]. Table 3 lists the maximum absolute errors of the numerical solution at t=1, the proposed wavelet method with N=6 and M₁=7, and several other methods^[19-23]. It can be observed from Table 3 that, the wavelet method is much more accurate than those proposed in Refs. [19]-[23]. The solutions of Burgers' equation can develop a sharp transition region with significantly small width when the initial condition is smooth^[23]. This situation poses a challenge for conventional numerical methods. For the case in which u(x, 0)=sin(2πx) and Re=100, Fig. 7 shows the numerical results on u(x, t) as a function of x at different moments t, where the spatial resolution level n=6, and the time step size h=1/2⁸. We can see from Fig. 7 that, the numerical solution of Burgers' equation clearly presents a sharp shock wave front after a certain propagation time, implying the competency of the proposed wavelet method in solving nonlinear problems with extremely large gradients.

Problem 4  Klein-Gordon equation

The nonlinear Klein-Gordon equation plays an important role in mathematical physics, particularly in investigating the behavior and interaction of the solutions in condensed matter physics, nonlinear wave equations, and dispersive phenomena in relativistic physics. A one-dimensional case can be given by

(52)

When f(x, t)=6xt(x²-t²)+x⁶t⁶, the analytical solution^[24-27] is u(x, t)=x³t³. Eq.(52) has been solved by the mixed finite difference and collocation methods^[24-27], where the time discretization is implemented with the finite difference^{[24-25, 27]} or Runge-Kutta method^[26], and the spatial discretization is implemented with the collocation method associated with the spline basis^[24-25] and the radial basis functions^[26-27], respectively. To use the proposed wavelet method, we define

(53)

Equation(53) can then be changed into

(54)

Taking the resolution level as n=4 in the space domain, Eq.(54) can be discretized into a system of first-order nonlinear ODEs with 34 degrees of freedom. The maximum absolute errors of the numerical solution obtained by the proposed wavelet method and the three other methods^[24-27] have been given in Table 4, which shows that the proposed SSTWM has much better accuracy even under a larger time step than that of the finite difference collocation, finite element collocation, and radial basis function meshless methods^[24-27].

7 Conclusions

Based on the unique property of the Coiflets in the approximation of functions, a space-time fully decoupled discretization method with high precision is proposed to solve nonlinear IBVPs. With the wavelet Galerkin method, the PDEs are usually transformed into a set of ODEs with unknown nodal parameters depending on time. Then, a novel wavelet-based step-by-step time-integrating scheme is constructed for the ODEs, which has the Nth-order convergence. Most promisingly, this order of convergence can be adjusted to any even number when the Coiflet with the order of the vanishing moments of this number minus one is adopted. Different numerical examples show that the proposed wavelet-based solution method has much better accuracy and efficiency than most existing numerical methods.

Appendix A

Proof  We first consider the following relation:

(A1)

Then, we consider

(A2)

For f(x)=xⁿ(n=0, 1, 2, ..., N-1), Eq.(A2) becomes

(A3)

or

(A4)

For fⁱ(x)=xⁿ(i, n=0, 1, 2, ..., N-1), considering Eq.(21) and noting Eq.(17), we have

(A5)

Inserting x=b yields

(A6)

or

(A7)

From Eqs.(A1), (A3) or(A4), and(A7), we have

(A8)

where F_d={2^-imf⁽ⁱ⁾(b)}, f(x)=xⁿ, and i, n=0, 1, 2, ..., N-1. From Eq.(A8), we can further obtain

(A9)

or

(A10)

Then, for f(x)=xⁿ(n=0, 1, 2, ..., N-1), if we define

Eq. (A9) can be exactly satisfied.

Appendix B

To derive the coefficients ζ_{a, i, k}, substituting Eq.(23) into Eq.(22), taking the derivatives on both sides of Eq.(10) with respect to x, and then taking x=a, we have

(B1)

where we assume that m is sufficiently large so that 2^ma+M₁-1 < 2^mb. Using F_a^(j) to replace f_a^(j) in Eq.(B1), we have

(B2)

Inserting into Eq.(B2) and taking α₁=M₁-1, we have

(B3)

Equation(B3) can be rewritten as

where

(B4)

Appendix C

When a=-∞ and b=j/2^m, Eq.(32) becomes

(C1)

Define

Then, we have

(C2)

We define

(C3)

Then, Eq.(C2) can be rewritten as follows:

(C4)

Equation(C4) shows that, for any j and k, the integral, i.e., , depends only on j-k. Therefore, we have

(C5)