## Abstract

For , the eigenfunctions of the nonlinear eigenvalue problem associated with the one-dimensional *q*-Laplacian are known to form a Riesz basis of *L*^{2}(0,1). We examine in this paper the approximation properties of this family of functions and its dual, and establish a non-orthogonal spectral method for the *p*-poisson boundary value problem and its corresponding parabolic time-evolution initial value problem with stochastic forcing. The principal objective of our analysis is the determination of optimal values of *q* for which the best approximation is achieved for a given *p* problem.

## 1. Introduction

The *p*-*Laplace operator* or *p*-*Laplacian*, a generalization of the ordinary Laplace operator, arises naturally in applications from physics and engineering including slow–fast diffusion related to particles (Aronsson *et al.* 1996), superconductivity (Barrett & Prigozhin 2000), wavelet inpainting (Zhang *et al.* 2007), image processing (Kuijper 2007) and game theory (Rossi in press). A typical application in the large *p* limit is a model for slow–fast diffusion for sandpiles (Evans *et al.* 1997). Recently there has been a significant amount of research activity encompassing methods of approximation for solution of nonlinear partial differential equations involving this operator (Barrett & Liu 1993*a*, 1994; Barrett & Prigozhin 2000). The aim of the present paper is to further contribute to this activity by considering the particular case of the one-dimensional *p*-Laplacian and examine in detail the approximation properties of a generalized spectral method described as follows.

Let *p*>1. For let [[*z*]]^{p−1}=*z*|*z*|^{p−2}. By extension from the linear case corresponding to *p*=2, we define the one-dimensional *p*-Laplacian to be the differential operator Δ_{p}*u*=([[*u*′]]^{p−1})′. Here is such that [[*u*′]]^{p−1}∈*H*^{1}(0,1). The corresponding *p*-Poisson boundary value problem is given by
1.1where *g*∈*L*^{2}(0,1). We also consider the related evolution equation
1.2that includes a stochastic forcing term where the noise intensity *ν*≥0 and *W* is a space–time Wiener process. Here
1.3where *β*_{n} are independent scalar Brownian motions, *ψ*_{n} are the eigenfunctions of the covariance operator *Q* and *α*_{n} the corresponding eigenvalues (see Prévôt & Röckner 2007; Liu 2009). For the case of space–time white noise, we have *Q*=*I*. In the case of a deterministic system when *ν*=0, equation (1.1) is the steady-state solution of equation (1.2).

Let *q*>1. The so-called *q*-*sine functions* are defined as the eigenfunctions of the *q*-Laplacian eigenvalue equation (Elbert 1981; Ôtani 1984; Lindqvist 1995):
1.4This family of functions and the corresponding problem (1.4) were studied over 30 years ago by Elbert (1981) and later by Ôtani (1984) and Bennewitz & Saitō (2004) in connexion with the computation of optimal constants in Sobolev-type inequalities. They generalize in a natural fashion the 2-sine basis (corresponding to the linear case), and they have very similar periodicity and interlacing structures.

In Binding *et al.* (2006), analogues of the classical completeness and expansion theorems for the *q*-sine functions were established for *q*≥12/11. Specifically it was shown that they form a Riesz basis of *L*^{2}(0,1). This leads to the following question: what are the approximation properties of this basis, as well as its dual basis, and how they relate to the approximation properties of the standard 2-sine basis?

Below we address this question by examining approximation of the solutions of equations (1.1) and (1.2) via projection methods with a *q*-sine and a dual *q*-sine basis, regarding *q*>1 and *p*>1 as free parameters. A main focus of attention is the determination of optimal values of *q* for which the highest order of convergence is achieved in a given *p*-problem. We demonstrate that standard properties of the 2-sine basis applied to the *p*=2 problem (such as super-polynomial convergence when *g*(*x*) is smooth) are lost, when a *q*-sine basis for *q*≠2 is considered for a *p*≠2 problem. As it turns out, the property of being a basis for the *q*-sine functions conceals a remarkably rich structure that is far from evident given the apparent simplicity of problem (1.4).

Background material on the *p*-sine basis and its dual is considered in §2. There we examine a matrix representation of the Schauder transform introduced in Binding *et al.* (2006). This will be crucial in our subsequent analysis as it gives rise to a stable procedure for constructing numerically both bases.

In §3, we find estimates for the approximation of square integrable functions in terms of their regularity. The dual *q*-sine basis turns out to have very similar approximation properties as the *q*=2 basis (lemma 3.1). On the other hand, however, it is fairly simple to construct smooth functions such that their *q*-sine Fourier coefficients do not decay faster than a power −5/2 for *q*>2 (lemma 3.2). The latter is in stark contrast with the most elementary results in the numerical approximation of solutions of differential equations by orthogonal spectral methods.

Section 4 is devoted to the *p*-Poisson boundary value problem. In theorem 4.2, we find explicit uniform bounds on the distance between any two solutions of equation (1.1), given the distance between the corresponding right-hand sides. We then examine in detail the numerical computation of solutions of equation (1.1) for source terms that are subject to various different regularity constraints. As it turns out, the estimates established in theorem 4.2 appear to be sub-optimal. A more thorough investigation in this respect will be reported elsewhere. See Ebmeyer & Liu (2005) for related results in the context of finite-element approximation of the solutions of equation (1.1), including the higher dimensional case.

In the final §5, we study the numerical approximation of solutions to equation (1.2) both in the deterministic and stochastic systems. We describe our discretization strategy and solve this problem for different values of *p* and *ν*. Our results provide evidence on the performance of the *q*-sine basis for the solution of equation (1.2), by showing the dependence on the parameter *q* of numerically computed *L*^{2} residuals.

## 2. The *q*-sine basis and its dual

The *q*-Laplacian eigenvalue problem (1.4), although nonlinear, has a fairly simple structure. The eigenvalues are found to be *λ*=(*nπ*_{q})^{q}, where . The first eigenfunction *f*_{1}(*x*) associated with the first eigenvalue (*π*_{q})^{q} is strictly increasing in [0,1/2], decreasing in [1/2,1] and it is even with respect to *x*=1/2. It can be extended to an odd function (with respect to *x*=0) in the interval [−1,1] and then to a 2-periodic *C*^{1} function of . If *q*>2, then *f*_{1}′′(*x*) is singular at *x*=1/2. The eigenfunctions *f*_{n}(*x*) associated with the eigenvalues (*nπ*_{q})^{q} satisfy *f*_{n}(*x*)=*f*_{1}(*nx*) for all *n*≥2.

In figure 1, we have plotted the first three eigenfunctions (i) and their corresponding derivatives (ii) for (*a*) *q*=1.4 and (*b*) *q*=10. They typify the case *q*<2 (*a*) and *q*>2 (*b*), respectively. For large *q*, the basis functions *f*_{n} approach zig-zag functions, which are the eigenfunctions of the -Laplace eigenvalue problem.

Below we always assume that the family is normalized by the condition *f*′_{n}(0)=*nπ*_{q} and leave implicit the dependence of *f*_{n} on *q*. In the special case *q*=2, we write , so that is an orthonormal basis.

The Pythagorean identity generalizes to the *q*-sine functions Elbert (1981) as
2.1Integrating this differential expression for small enough *x* leads to the following explicit representation for the inverse function
2.2As we will see below, this representation plays a crucial role in the numerical estimation of *f*_{n}(*x*).

Let the *Schauder transform*, *T*_{q}, be the linear extension of the mapping *e*_{n}↦*f*_{n}. Then *T*_{q}:*L*^{2}(0,1)→*L*^{2}(0,1) is an invertible bounded operator for all *q*≥12/11 (Binding *et al.* 2006; theorem 1). Thus is a Riesz basis of *L*^{2}(0,1) for such a range of the parameter *q*. Further evidence presented in Binding *et al.* (2006) suggests that in fact this is also the case for all *q*>1, but at present this has not been proved rigorously. Unless otherwise specified, we will assume from now on that *q*≥12/11.

The property of a Riesz basis ensures that every *g*∈*L*^{2}(0,1) is represented by a unique series expansion which is convergent in norm. The *q*-*sine Fourier coefficients*, , are given explicitly by *a*_{n}=〈*g*,*f**_{n}〉, where is the basis dual to . Since for all , then . It turns out that for *q*≠2.

The following matrix representation of *T*_{q} is fundamental to our analysis. Let
be the *j*th 2-sine Fourier coefficient of *f*_{1}(*x*). Then the *k*th 2-sine Fourier coefficient of *f*_{n}(*x*) is given by
2.3Hence and, therefore, *T*_{q} has a lower triangular matrix representation in the orthonormal basis (figure 2*a*).

### Remark 2.1

The basis of eigenvectors of the -Laplace eigenvalue problem are zig-zag functions (Binding *et al.* 2006; §5). In this case, we can write explicitly. As it turns out,

Let *s*>0. Below we denote by , the Sobolev space of 1-periodic functions *g*∈*L*^{2}(0,1) such that .

### Lemma 2.2

*Let f*_{1} *be the first q*-*sine function as defined above. If* 1<*q*<2, *then* . *If q*>2, *then* . *If q*>4, *then* .

### Proof.

A straightforward argument involving integration by parts yields
Then
2.4Hence *τ*_{q}(*j*)=*O*(*j*^{−2}) as and . From equation (2.1) it follows that *f*′′_{1}(*x*)=*h*(*f*_{1}(*x*)) for .

Let 1<*q*<2. Then
so .

Let *q*>2. Since , then . Moreover, *f*_{1}(*x*)≥2*x* for . Then
If *q*>4, the integral on the right-hand side diverges and so . ■

Evidently the *H*^{s} regularity for *f*_{1}(*x*) found in lemma 2.2 is not optimal. Figure 2*b* shows a numerical estimation of the precise value of *s*(*q*), such that for *r*<*s*(*q*) and for *r*>*s*(*q*). The data for this graph were obtained by computing the decay rate of *τ*_{q}(*j*) for a large truncation of *T*_{q}. A thorough investigation closely related to this lemma in the higher dimensional context can be found in Barrett & Liu (1993*b*), Ebmeyer *et al.* (2005) and references therein.

In the large *q*, we have a limit of *s*=3/2 and this is confirmed by remark 2.1. At *q*=2, we simply have and so *f*_{1} is in for all *s*. According to lemma 2.2, the curve should remain below *s*=3 as *q*→2^{+}. For *q*>3, the graph suggests . However, if *q*<3 it suggests *s*(*q*)>2. As *q*→1, the regularity drops and the limit seems to approach *s*=2. There is an interesting ‘peak’ of regularity around *q*=1.6, which we cannot presently explain.

The 2-sine Fourier coefficients of *f**_{n}(*x*) are given by the *n*-th columns of . This operator has an upper triangular representation in the 2-sine basis. Then *f**_{n}(*x*) are trigonometric polynomials of order *n*. In fact, are parallel for all *q*>1. Unlike the *q*-sine functions, not all dual *q*-sine functions have the same periodicity structure.

We now describe a stable numerical procedure for computing these two bases. The *q*-sine functions can be approximated by first estimating *f*_{1} using equation (2.2). Numerical integration yields on [0,1/2]. Although the integral is singular at *y*=1, the value is known, so we do not need to consider quadrature points too close to this singularity. In our numerical procedure, we chose a fine uniform grid and apply a cumulative Simpson's rule. This gives an approximation of *f*_{1} for *x*∈[0,0.5] with a controlled tolerance and by symmetry we obtain for *x*∈[0,1]. Note that is given on a non-uniform grid. The Pythagorean identity (2.1) immediately yields the derivative *f*′_{1} of *f*_{1} and hence we can use it to approximate the former with a , defined also on the non-uniform grid.

Once we have constructed and , we use periodicity and symmetry to find corresponding approximations of *f*_{n} and *f*′_{n} for *n*>1. This involves considering scaled copies and , to form and . Here, the number of non-uniform grid points on [0,1] grows with each *n*=2:*N*.

To obtain *f*_{n} and *f*_{n}′ on a uniform grid *x*_{j}=*jh* for *j*=0:*J*, rather than the non-uniform grid that arises from the numerical integration, we have considered numerical interpolation by piecewise cubic polynomials. This gives and defined at *x*=*x*_{j}. The former is the approximated basis and the latter the corresponding derivatives that we use for further computation.

Figure 1 was generated with an implementation of the numerical scheme just described on an uniform grid with *J*=4000 points (*h*=2.5×10^{−4}). The integral (2.2) was approximated with 2×10^{5} points.

The dual basis is found from an *N*×*N* truncation, , of the Schauder transform. In practice, we first compute *τ*_{q}(*j*) and assemble . Then, we define approximations (*f**_{n})^{h} for *n*=1:*N*, as the trigonometric polynomials whose *k*th 2-sine Fourier coefficients are the (*n*,*k*) entry of the matrix .

### Remark 2.3

In our numerical approximation of the set of basis functions, we are careful to fully resolve oscillations on the basis function *f*_{N}, taking at least 20 mesh points per wavelength and for most computations 100 per wavelength. For *N*≤50, we use *h*=(100*N*)^{−1} and so resolve each oscillation in the basis function with 100 spatial points. For *N*>50, we use *h*=(20*N*)^{−1} and so resolve with 20 points. We also examined convergence of orthogonality of the basis and dual in the spatial discretization and noted *O*(*h*^{2}) for *q*≈10 through to *O*(*h*^{4}) for *q*≈1.4.

## 3. Approximation of source terms

Any given *g*∈*L*^{2}(0,1) can be approximated by either
for large *N*. These two expansions converge as in the norm of *L*^{2}(0,1) and also pointwise for almost all *x*∈[0,1]. Unlike in the linear case corresponding to *q*=2, the rates of decrease of ∥*g*−*g**_{N}∥ and ∥*g*−*g*_{N}∥ can be very different when *q*≠2.

Since the dual basis comprises trigonometric polynomials, on the one hand, we can formulate the following natural statement.

### Lemma 3.1

*Let g*∈*L*^{2}(0,1). *For all*

### Proof.

By definition, . Since the matrix associated with is an upper triangular (see §2), According to equation (2.4), we have . Hence Thus, ■

Therefore, the *q*-sine dual expansion of any converges super-polynomially fast. On the other hand, however, it is not difficult to construct examples of smooth functions *g* with a subsequence of *q*-sine Fourier coefficients decaying slowly.

### Lemma 3.2

*Let g*(*x*)=*e*_{1}(*x*). *If j is prime, then* .

### Proof.

Since *T*_{q} has an (infinite) lower triangular matrix representation in the orthonormal basis , we can find the entries of the corresponding matrix representation of by pivoting and forward substitution (Gaussian elimination). It is readily seen that is necessarily the lower triangular and its diagonal should be constant and equal to *τ*_{q}(1)^{−1}.

Assume that *j* is prime. According to equation (2.3), the only non-zero entries in the *j*th row of *T*_{q} are *τ*_{q}(*j*) in the first position and *τ*_{q}(1) in the *j*th position. Therefore, the only non-zero entries in the *j*th row of (*T*_{q})^{−1} are −*τ*_{q}(*j*)/*τ*_{q}(1) in the first position and *τ*_{q}(1)^{−1} in the *j*th position. As the Fourier sine coefficients of are obtained from the *j*th column of , the desired conclusion follows. ■

By virtue of lemma 2.2, the prime *q*-sine Fourier coefficients of for *q*>2 cannot decrease faster than *j*^{−5/2} in the large *j* limit. Observe that this is in stark contrast with the most elementary results in the numerical approximation of solutions of differential equations by orthogonal spectral methods.

### Remark 3.3

The finite set of basis functions *f*_{n} and dual for *n*=1:*N* generate corresponding *N*-dimensional subspaces *V* _{N} and of *L*^{2}(0,1). Instead of computing directly with these non-orthogonal bases, one can apply the Gram–Schmidt algorithm in order to obtain orthonormal bases of these subspaces. This has a numerical advantage of not needing to store both the basis and dual. We also considered this approach, however, we found little advantage in terms of accuracy.

Let us now consider various numerical tests on the approximation of regular functions by and . Once the bases have been obtained on uniform mesh, we can examine their approximation properties numerically by looking at the decay of suitable residual for benchmark sources *g*∈*L*^{2}(0,1).

Figure 3 shows the typical outcomes of an experiment to determine the dependence in *q* of the *L*^{2}(0,1) residual. We have fixed here *N*=40, varied *q*=1:100, and computed residuals in the approximation of the following four functions:
3.1in figure 3*a*–*d*, respectively. In the figure, we include both the basis and its dual, as well analogous calculations with the orthogonalized bases of *V* _{40} and . In figure 3*a*, we see an optimal *q*_{opt}=10 as expected, in figure 3*b*–*d*, we increase the regularity of *g* and see *q*_{opt} decreases with values of 4.25, 2.9, 2.55 for figure 3*b*–*d*, respectively.

This experiment gives a general insight about the *q*-behaviour of *L*^{2}-residuals in the approximation of functions with different degrees of regularity by *q*-sine bases and their duals. For the simple functions considered, our calculations indicate the following general behaviour of the *q*-sine basis:

— as

*q*→1 the residual deteriorates,— there is always a single minimum corresponding to an optimal

*q*=*q*_{opt},— as , the residual curve becomes asymptotically constant with no local maximum for

*q*>2.

In contrast, the approximation error is almost constant in *q* for the *q*-sine dual basis. This is indeed a consequence of lemma 3.1 and the fact that are trigonometric polynomials. Our tests indicate that there does not seem to be a clear advantage in orthogonalizing the basis or the dual basis for errors with an order of magnitude above 10^{−5}.

In table 1, we have estimated *α*>0 such that ∥*g*−*g*_{N}∥<*βN*^{−α}, where *β*>0 is independent of *N*. The data indicate that for , *q*_{opt}≈2. Moreover, as *q* increases, we should expect *α* to decrease and stabilize always below 1.5 for . This is indeed suggested by lemma 2.2, figure 2*b* and lemma 3.2, and it is confirmed by the last row of the table.

In figure 4*a*, we have computed the residual ∥*g*−*g*_{10}∥ (red) and ∥*g*−*g*_{20}∥ (blue) for randomly generated constrained to *g*(0)=*g*(1)=0. To construct a random function *g* such that the norm is finite, we find , where *a*_{j}=(1+*j*^{2})^{−s/2}|*j*|^{−1/2−δ} for some small *δ*>0 and *β*_{j} are independent identically distributed *N*(0,1). Two realizations of functions *g* obtained in this way are shown in the inserts in figure 4. Note that *q*_{opt} is achieved at different places but close to 2. In figure 4*b*, we examine *q*_{opt} over 200 realizations of functions in . For any fixed realization and value of *N*, *q*_{opt}≠2 although in the limit, this optimal parameter appears to be close to 2. For *N*=100, we also show the results of 1000 realizations. In this case, the mean value is closer to 2 although the variation is still large. In figure 4*c*, we show the mean values where functions are taken in for *s*=−0.5, 0, 0.5, 1 and 2 with 200 realizations. For *s*<2, the variability in *q*_{opt} is far larger. For fixed *N*, *q*_{opt}>2.

The observed outcome of this experiment strongly support the conjecture that *q*_{opt} typically approaches two as the regularity of *g* is increased.

## 4. Numerical solution of the *p*-Poisson equation

We now address the question of approximating the solutions of equation (1.1) by means of a *q*-sine and dual basis. In view of lemmas 3.1 and 3.2, we begin by determining uniform estimates on how sensitive this solution is under perturbations of the right-hand side. Analogous questions have certainly been considered in more general frameworks, however, here we focus on the explicit calculation of the constants involved.

A key ingredient in the estimates presented below is the fact that equation (1.1) can be integrated explicitly. Let the Volterra operator
Note that *V g*∈AC[0,1] for all *g*∈*L*^{1}(0,1). Furthermore *V* :*L*^{s}(0,1)→*L*^{t}(0,1) is a contraction operator for all and its norm can be explicitly determined (Bennewitz & Saitō 2004, theorem 1.1). Let
Then *h*_{g}(*γ*) is a continuous function, decreasing in *γ*, for all fixed *g*∈*L*^{1}(0,1). Let
be the unique root such that *h*_{g}(*γ*_{0}(*g*))=0. Then
4.1is the unique solution of equation (1.1).

We firstly establish concrete Hölder estimates on *h*_{g}(*γ*). Without further mention in this section, we will fix *r*=1/(*p*−1) and denote by ∥⋅∥_{s} the norm of *L*^{s}(0,1). In the case *s*=2, we will continue suppressing the sub-index.

### Lemma 4.1

*Let g*∈*L*^{1}(0,1) *and* −∥*g*∥_{1}≤*γ*≤*μ*≤∥*g*∥_{1}. *Then*

### Proof.

Suppose first that 0<*r*≤1. From the graph of [[*z*]]^{r} for |*z*|≤*M* it is readily seen that [[*z*]]^{r}−[[*w*]]^{r}≥*rM*^{r−1}(*z*−*w*) for all −*M*≤*w*≤*z*≤*M*. Then
for *M*=2∥*g*∥_{1}, and *γ* and *μ* as in the hypothesis.

In a similar fashion, let *r*>1. Then (*z*−*w*)^{r}≤2^{r−1}([[*z*]]^{r}−[[*w*]]^{r}) for all −*M*≤*w*≤*z*≤*M*. Indeed, if 0≤*w*<*z*, a straightforward argument shows that
if *w*<*z*≤0,
and if *w*<0<*z*,
achieved when *w*=−*z*. Thus
■

### Theorem 4.2

*Let u and* *be solutions of equation (1.1) with corresponding sources g and* *. Let* *. Then
*

### Proof.

Let 0<*r*≤1 and *s*=2/*r*≥2. By virtue of equation (4.1)
Note that |[[*z*]]^{r}−[[*w*]]^{r}|≤2^{1−r}|*z*−*w*|^{r} for all 0≤|*w*|≤|*z*|. Thus
According to lemma 4.1,
This ensures the first statement.

Let *r*>1. In a similar fashion as before, we see that
Lemma 4.1 and similar arguments as for the previous case yield
This completes the proof. ■

The right-hand side bound in the above theorem approaches 2 as *r*→0 independently of the value of . In fact for . Since *l*_{r}(*τ*)→0 for almost all *τ*∈[0,1] and , the Dominated Convergence Theorem yields as *r*→0. This is a well-known property of the -Laplacian (see Rossi in press). As mentioned in §1, the approach considered in Ebmeyer & Liu (2005) in the context of finite-element approximation of the solutions of equation (1.1) may provide an insight on whether the constants found in the above theorem are optimal.

We now describe how the *p*-Poisson problem is discretized. The strong formulation (1.1) leads to the following weak formulation using integration by parts and the boundary conditions:
4.2for any absolutely continuous test function *v*. We expand both *u* and *v* using a basis , where *ϕ*_{n} is either *f*_{n} or . After truncation this leads to the following nonlinear system of equations for unknown coefficients :
4.3The case *p*=2 evidently reduces to the standard linear system to solve for the Poisson equation.

Numerically the right-hand side of equation (4.3) is approximated via quadrature rules. The nonlinear system may then be solved, for example, by Newton's method. Since the Jacobian is a full matrix in this case, a banded approximation appears to provide sufficient accuracy. However, the results presented below were found using the trust-region dogleg method with a full Jacobian matrix as implemented in Matlab.

In figure 5, we have used this scheme to solve the *p*-Laplacian problem with *p*=5 and examine the *L*^{2} error taking *N*=40. In figure 5*a*, we have fixed a solution where we clearly expect and observe that *q*=2 is the optimal basis. In figure 5*b*, we have fixed *u*(*x*)=*f*_{1}(*x*) for *p*=5 and so we observe the *q*=5 as the optimal basis.

In figure 6, we take *g*=*g*_{b} from equation (3.1) that turns out to be a typical form of forcing for sandpile problems. We solve the *p*-Laplacian problem for (*a*) *p*=1.8, (*b*) *p*=3, (*c*) *p*=5 and *p*=10 with *N*=40 and examine the *L*^{2} error in the solution as we vary the *q* basis. To estimate this error, we take as exact the solution with 2-sines and 2*N* modes. We observe that the optimal basis for representing the solutions is no longer the standard *q*=2 for *p*=3,5,10. For problems with moderate *p* (for example, *p*=3), we see two distinct minima. For larger *p* problems, however, the *q*=2 basis becomes less competitive.

In table 2, we give estimates of *q*_{opt} for *g*=*g*_{b} and *g*=1 for *N*=20,40,60 and *N*=80 modes. The changes in *q*_{opt} with *N* are explained by the interchange of the two minima in figure 6. For *N*=100, this interchange occurs for *p*=5 and *p*=10 as illustrated in figure 7*b*.

## 5. The time-dependent *p*-Laplacian

In this final section, we consider the evolution equation (1.2) both in the deterministic (*ν*=0) and the stochastically forced regime (*ν*≠0). Slow–fast diffusion is often taken in the large *p* limit, as a model for sandpile growth. For further details including arguments about the validity of the modelling (see, for example, Aronsson *et al.* 1996; Evans *et al.* 1997; Falcone & Finzi Vita 2006; Andreu *et al.* 2009; Caboussat & Glowinski 2009; Igbida 2010). Our purpose here is to consider the time-dependent problem for different values of *p* and examine the choice of basis in the formation of the sandpile.

We discretize the weak form of the equation and truncate to solve the time-dependent version of equation (4.3) given by
5.1This expression is then discretized in time to get a nonlinear system of equations to solve for at each step. Below we consider an Euler's method that reduces to
where Δ*β*_{m} are independent identically distributed random variables with mean zero and variance Δ*t* (recall (1.3)), and we have assumed the eigenfunctions of *Q* are now given by *e*_{n}. For the case of stochastic forcing, we take *M*≫*N*.

In figure 8, we consider *ν*=0 with *g*=*g*_{b} and we plot the time evolution for (*a*) *p*=1.8 (*q*_{opt}=2), (*b*) *p*=3 (*q*_{opt}=2.95), (*c*) *p*=5 (*q*_{opt}=3.9) and (*d*) *p*=10 (*q*_{opt}=5.75). In each case, the evolution found numerically quickly converges towards the steady state.

Figure 8 was obtained by choosing in each case the predicted *q*=*q*_{opt} basis for the steady state. Intuitively this choice of basis should be near optimal for large *t*. On the other hand, however, there is no reason to presume that it is so for small values of *t*. In figure 9, we examine the spatial error at each step for four different bases with *N*=40. For comparison, we have chosen *p*=10 and *q*∈{1.8,2,*q*_{opt},10}. As the solution approaches the steady state, in each case, the basis *q*_{opt} has the smallest error as is expected (in fact a gain of almost an order of magnitude in accuracy compared with *q*=2 is observed). It is remarkable however that, in the transient regime, these bases appear to be far less accurate than other choices of *q*, such as *q*=2 and *q*=1.8.

Let us now consider the stochastically forced case, *ν*>0. Figure 10 shows plotted solution with time-dependent noises that are *H*^{1} in space and white in time for (*a*–*c*), and white both in space and time for (*d*–*f*). The forcing corresponds to *ν*=0.2. We expect that on average, the noisy solution is simply that of the deterministic system (which has a unique stable solution). We see in figure 10*a*,*b*,*d*,*e*, the time evolution of the solution for one particular realization of the noise. This should be compared with the deterministic case in figure 8*d*.

In figure 10*c*,*f*, we plot the evolution of the error for fixed time step with *q*∈{1,2,*q*_{opt},10}. Similar to the deterministic case, in the transient regime, the error for *q*_{opt} and *q*=10 is far higher than that for *q*=2 or *q*=1.8. Where the noise effects are large (such as for the spatially correlated noise in *H*^{1} in figure 10*c*), the optimal basis is no longer clear. However, where the effect is small (such as the white noise in figure 10*f*), we clearly see in the time-dependent evolution about the steady state that the *q*=4.85 outperforms the other bases (and in particular *q*=2). For sandpile-type problems, the introduction of a spatially white noise is a natural choice. It is interesting to note the small effect on the dynamics in this realization.

At present, it is unclear whether the use of a *q*-sine basis with *q*≠2 provides any real computational advantage over the natural choice *q*=2 for the solution of equation (1.2). There is clearly a computational overhead in obtaining the former for *q*≠2 that impacts significantly on efficiency. However, we have observed that the nonlinear problem (4.3) is solved faster in the optimal basis. This is certainly worth further investigation. For example, for *N*=40, we observe approximately a 20 per cent speed up in the nonlinear solve over the standard 2-sine basis. Thus, if solving many fixed *p*-problems, the corresponding optimal *q*-sine basis may not only be more accurate but also more efficient. In our current implementation, an average approximation to the optimal *q*-sine basis may be precomputed for a given regularity of *g*. Lemma 2.2 gives some preliminary indication on how to solve the problem of *a priori* determining the optimal basis. Our numerical results suggest that for large *p* we expect an optimal basis with *q*>2 for a right-hand side with a discontinuous derivative.

## Acknowledgements

We kindly thank A. Vignes and B. Rynne for discussions related to this paper. The first author acknowledges support from the Université Paris Dauphine, where part of this research was carried out.

- Received September 17, 2010.
- Accepted March 18, 2011.

- This journal is © 2011 The Royal Society