1 Introduction

Nonlinear differential equations (DEs) exist widely in science and engineering areas ^[1-9]. The so-called p-Laplacian equation is such a nonlinear DE, and has been successfully used to model the problems of thermology ^[7], fluid dynamics ^[8], and structural mechanics ^[9]. A general form of the two-dimensional p-Laplacian equation can be expressed as follows:

(1)

where is a known function, and is defined by

(2)

It can be seen from Eqs. (1) and (2) that, this is a class of second order elliptic boundary value problems in the divergence form with gradient nonlinearity. When p=2 and g(x, y)=0, Eq. (1) degenerates into the normal Laplacian equation. The p-Laplacian equations of Eq. (1) not only has broad applications, but also represents a benchmark DE for verifying the numerical methods. Huang et al. ^[10] emphasized that the p-Laplacian equations hold most of the very basic difficulties for the numerical solutions of general nonlinear DEs, and many existing numerical techniques actually did not work well in their solutions. Under this situation, there have been many attempts to solve the equations. Chow ^[11] considered the application of the finite element method (FEM), and derived the error estimates in the energy norm. Barrett and Liu ^[12] studied the continuous piecewise linear finite element approximation (FEA), and numerically solved several radial symmetric problems. Ainsworth and Kay ^[13] further considered the FEA of the p-Laplacian equations defined on a polygonal domain, and derived the error estimates in terms of the mesh size and the parameter p. Huang et al. ^[10] proposed the preconditioned descent algorithms (PDA) with the FEA. Zhou et al. ^[14] brought up the preconditioned hybrid conjugate gradient algorithm (PHCGA), and adopted the FEA. In addition to the FEA, Andreianov et al. ^[1] used the finite volume approximation (FVA) on the rectangular meshes to discretize the p-Laplacian equations. Oberman ^[15] built a scheme with the finite difference method (FDM). Lefton and Wei ^[16] adopted the penalty method associated with finite elements. Pezzo et al.^[17] constructed the interior penalty discontinuous Galerkin method to approximate the p-Laplacian equations. These studies have achieved significant advances in the numerical solution of the p-Laplacian equations. However, as has been pointed out by Huang et al. ^[10] and emphasized by Zhou et al. ^[14], there still exist difficulties in the conventional methods to successfully solve the p-Laplacian equations.

As one of the developing potent mathematical tools, the wavelet theory has been successfully used to solve the partial DEs (PDEs) in many fields of science and engineering ^[18-24] by establishing wavelet-based numerical methods. He and Han ^[18] proposed a wavelet FDM to solve the elastic wave equation. Ding et al. ^[19] introduced a wavelet multiscale method for the Maxwell equation inversion. Xiang et al. ^[20] provided a wavelet-based FEM to solve the Poisson equation. Recently, our group has developed a modified wavelet Galerkin method (MWGM) ^[25], which has been successfully used for solving different types of nonlinear integral equations ^[26] and DEs ^[27-29] with applications to large deflection bending and vibration for plates ^[30-32]. The associated numerical results have shown that such an MWGM has much better accuracy and efficiency than most existing conventional numerical methods. As inspired by the progresses in wavelet-based numerical techniques, we propose a new wavelet approximation scheme to improve the MWGM, and then use it to discretize the p-Laplacian equations, so that a convenient and high accurate wavelet method can be established to numerically solve the p-Laplacian equations.

2 Wavelet approximation of interval-bounded functions

Coiflets are a family of compactly supported orthogonal wavelets. For this type of wavelets, the scaling function φ (x) has a compact support [0, 3 N-1], and satisfies

(3)

(4)

(5)

where N is a positive even integer, , and the filter coefficients p _k can be obtained through the following relations ^[25]:

(6)

(7)

(8)

(9)

Once the filter coefficients are obtained by solving Eqs. (6) -(9), the values of φ (x) at the dyadic points can be obtained from Eq. (3). During the construction of φ (x), N and M₁ can be artificially chosen ^[25]. The filter coefficients p₀, p₁, …, p₁₇ for N=6 and M₁=7 are listed in Table 1 ^[25].

According to the multi-resolution analysis of the wavelet theory ^[25], the function space L²(R) can be divided into a series of nested subspaces

for which we have

so that a function f(x)∈ L²(R) can be approximated by the projection of the function from L²(R) to V _J as follows ^[25]:

(10)

where J is the so-called resolution level ^[25], and

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

(11)

and

(12)

in which N≥ n, and C is a constant depending on f(x) and φ (x). We note that the property of the Coiflet allows any polynomials with the order up to N-1 to be exactly represented by Eq. (11).

Since the scaling function basis is defined on the whole real line, Eq. (11) cannot be directly applied to the functions bounded on a finite interval. When approximating f(x)∈ L²[0, 1], the summation index k in Eq. (11) satisfies

at the same time, such that the function values of f((M₁ + k)/2 ^J) are undefined. To resolve this conflict, one usually needs to use the zero, symmetric, or periodic extension of the function near the boundaries of the interval. However, these treatments may sometimes introduce artificial "jump" to the values of the function and/or its derivatives close to the boundaries, and eventually result in numerical instability or large approximation error. Wang ^[25] has considered a natural extension treatment by using the Taylor series expansion, in which the boundary derivatives are approximated by the numerical difference methods of the third-order. Although such a treatment can effectively reduce the boundary error when approximating the interval-bounded functions, the relevant approximation formula is not able to exactly represent all the interval bounded polynomials with the order up to N-1, since only low order numerical difference approximations to the boundary derivatives are adopted. In order to overcome this drawback, in this study, we consider the technique of the Lagrange polynomial interpolation.

For the function f(x)∈ L²[0, 1], we choose m+1 nodal points, i.e., x₀, x₁, …, x _m ∈ [0, 1], and the corresponding function values f(x₀), f(x₁), …, f(x _m). Then, we have the Lagrange polynomial interpolation to the function at x close to x₀ as follows:

(13)

Particularly, if the nodal points are spaced equally with the distance 1/2 ^J, i.e.,

we have

(14)

Based on Eq. (14), we can construct a new function with the boundary extension to the original f(x) as follows:

(15)

It is easy to verify that, for any polynomials with the order up to m, we have

Similar to Ref. [25], using Eq. (11) to approximate and letting m= N-1, we have

(16)

where

(17)

(18)

(19)

(20)

The approximation formula of Eq. (16) can exactly represent any polynomials on [0, 1] and with the order up to N-1.

Equation (16) is also valid for the approximation of multi-dimensional functions. Considering the two-dimensional function u(x, y), where (x, y)∈ [0, 1]× [0, 1], treating u(x, y) as the function of x, and then applying Eq. (16), we have

(21)

or

(22)

Then, treating each term of as the function of y, and using Eq. (16) once again, we have

(23)

Combining Eqs. (22) and (23), we have

(24)

or

(25)

where

(26)

(27)

We note that the errors of the approximations in Eqs. (16) and (24) are still on the same order as that in Eq. (12), i.e., O (2^{- JN}). To verify the accuracy of this approximation, we take

as examples. Figure 1 shows the relative L² error, defined by , of f(x)=tanh x, where is the approximation of f, by using the Fourier series, various orthogonal polynomials, and Eq. (16) with N=6 and M₁=7, respectively. Under the condition of 17 expansion terms, Fig. 1 shows that the proposed Coiflet-based method has a much higher accuracy than those of the conventional methods. Figure 2 illustrates the relative L² error at x=0.5 by using Eq. (16) for the resolution levels J= 4, 5, 6, and 7, respectively. We can see from Fig. 2 that the error is exactly on the order of O(2^{- NJ}), where N=6.

Figure 3 shows the relative L² error in approximating u(x, y) by using Eq. (24) with J=4. Figure 4 illustrates the cross-section of the error profile at x=0.5. It can be seen from Figs. 3 and 4 that the maximum error can be on the order of 10^-8, and the boundary error can be on the order of 10^-10. This indicates that the proposed approximating scheme is very accurate.

3 Solution method of two-dimensional general p-Laplacian equations

Using Eq. (24) to approximate the two-dimensional function u(x, y) ((x, y)∈ [0, 1]× [0, 1]) and then performing derivative calculations yield

(28)

Considering α =1, β =0 and α =0, β =1 and rewriting Eq. (28) in the matrix form yield

(29)

Denote

Then, from Eq. (24), we have

(30)

Equation (30) can be rewritten in the matrix form as follows:

(31)

where

(32)

Equating Eqs. (29) and (31), multiplying Φ _{J, m}(x)Φ _{J, n}(y), and performing integration from 0 to 1 for both x and y, we have

(33)

where

From Eq. (33), we have

(34)

where P is symmetric.

Then, we have

(35)

where symbols and are the so-called Hadamard products ^[34] of matrices, which are defined as the element-wise multiplication and the power operation, respectively, e.g.,

Define

Then, we have

where

Finally, the general p-Laplacian equation shown in Eq. (1) can be discretized as follows:

(36)

where H and S are functions of U, G={ g(k/2 ^J, l/2 ^J)}, and k, l=0, 1, …, 2 ^J. Denote F(U)= A^T RB+ B^T SA+ G. Then, we can rewrite Eq. (46) as follows:

(37)

which is a matrix equation depending on the unknown u(k/2 ^J, l/2 ^J). With the classic Newton iteration method, we can obtain the solution of Eq. (37).

4 Numerical examples

Example 1 Huang et al. ^[10] considered an axial symmetric p-Laplacian equation by assigning

to Eq. (1). For such a p-Laplacian equation, the exact solution can be given by

(38)

Huang et al. ^[10] solved this equation by combining the PDA iteration method and the FEA. By using the meshes with 1 601 unknowns, and after 9 iteration steps, the obtained result has the L² norm error defined by , where is the approximate solution, i.e., 5.097 8× 10^-4.

To apply the proposed wavelet Galerkin method, we use Eq. (16) to approximate the unknown function u(r). The boundary conditions can be embedded into the approximation series by adjusting the summation index from 0~2 ^J to 1~2 ^J-1. Define U={ u _k} as the exact solution and as the approximate solution. We consider both the l² and l^∞ errors, i.e., and , which are listed in Table 2 with the resolution levels 4, 5, 6, 7, respectively. It can be seen that the wavelet method can reach the similar error level 5.8151× 10^-4 by using only 15 unknowns.

Example 2 Bermejo and Infante ^[33] provided the benchmark test problem of Eq. (1) associated with

Obviously, u(x, y)=10 x(1- x) y(1- y) is the exact solution. By proposing the full approximation storage multigrid algorithm (FASMA), Bermejo and Infante ^[33] obtained the finite element solution for this problem. The best result in Ref. ^[33] is achieved by using 9 409 nodes. The relative L² error defined by is estimated to be 5.382 0× 10^-4. In addition, three pre-and post-smoothing steps have been adopted during the solution to improve the computing efficiency and accuracy ^[33].

When using the proposed wavelet Galerkin method, Eq. (24) is adopted to approximate the unknown function u(x, y), where the range of the summation index changes to 1~2 ^J-1 so that the zero boundary condition can be satisfied. Table 3 shows the relative l² error defined by . It can be seen from Table 3 that the present wavelet results are much more accurate than that obtained by the FASMA with the additional pre-and post-smoothing treatments ^[33].

Example 3 Andreianov et al. ^[1] considered Eq. (1) associated with

The exact solution is u(x, y)= e ^x sin (3 π x) sin (3 π y). Figure 5 shows the relative l^p errors defined by and obtained by using the proposed wavelet Galerkin method, the admissible finite volume scheme ^[1], and the non-admissible finite volume scheme ^[1]. From the figure, we can see that the mesh size h for the wavelet method can be regarded as 1/2 ^J, where J is the resolution level, and the results of the proposed wavelet method are much better than those of the admissible finite volume scheme and the non-admissible finite volume scheme.

5 Conclusions

In this paper, we have considered a new boundary extension technique in terms of the Lagrange interpolating polynomial to improve the accuracy of the Coiflet-based approximation of functions defined on an interval. The obtained approximation formula can exactly represent any interval bounded polynomials with the order up to N-1, where the Coiflet with the compact support [0, 3 N-1] is adopted. Combining such an approximation algorithm with the conventional Gelerkin method, we establish the solution procedure of the challenging p-Laplacian equations. The numerical examples indicate that the proposed wavelet method is much more accurate and efficient than several other numerical methods.