A numerical comparison of Chebyshev methods for solving fourth-order semilinear initial boundary value problems
B.K. Muite
Mathematical Institute, University of Oxford, 24-29 St. Giles’, Oxford OX1 3LB, UK

Abstract In solving semilinear initial boundary value problems with prescribed non-periodic boundary conditions using implicit-explicit and implicit time-stepping schemes, both the function and derivatives of the function may need to be computed accurately at each time step. To determine the best Chebyshev collocation method to do this, the accuracy of the Chebyshev differentiation, preconditioned Chebyshev tau and Chebyshev integration methods are compared in the L2 and W 2,2 norms when solving linear fourth order boundary value problems. We find that the best Chebyshev method to use for high resolution computations of solutions to initial boundary value problems is the Chebyshev integration method which uses sparse matrix operations and has a comparable computational cost to a Fourier spectral discretization. Key words: Spectral collocation, fast methods, time-stepping

1

Introduction

The motivation for the comparison of these spectral methods is to compute solutions to high-order semilinear initial boundary value problems found in elastodynamic models for microstructure formation during phase transitions in which a small Ginzburg or capillarity term is added. Typical examples of these semilinear initial boundary value problems can be found in studies by Ahluwalia et al. [1], Hoffmann and Rybka [19] and Truskinovsky [37]. Many
∗ tel +44 (0)1865 273525, fax +44 (0)1865 273583 Email address: muite@maths.ox.ac.uk (B.K. Muite).

Preprint submitted to Elsevier

21 August 2008

previous studies of microstructure formation with non-periodic boundary conditions have used low order finite difference or finite element methods, see for example, Ahluwalia et al. [1], Vainchtein [39], Dondl and Zimmer [10] and the review by Luskin [29]. An objective of this study is to show that Chebyshev collocation methods efficiently simulate multiscale phenomena in regular but non-periodic domains and should be considered as a viable alternative to low order methods. In a typical implicit-explicit (IMEX) or fully implicit time-stepping scheme for an initial boundary value problem, not only is the function required at each time step, but derivatives of the function may also be required to calculate the explicit part of the time step or to perform a fixed point iteration in a fully implicit time stepping scheme. To determine the best method to do this, the accuracy of approximate solutions and derivatives of approximate solutions obtained using different Chebyshev collocation methods in solving linear and initial boundary value problems are compared. The accuracy of solutions to linear boundary value problems are studied because in typical time-stepping schemes, a linear boundary value problem is solved at each time step or fixed point iteration. Chebyshev collocation methods have also been used to obtain high resolution numerical solutions to the KdV, Allen-Cahn and Cahn-Hilliard equations – see, for example, Xu and Tang [41] and Kassam and Trefethen [23]. They may also be useful in examining viscosity solutions to conservation laws, such as Burgers equation regularized by viscosity and dispersion, see, for example, Chen et al. [3] or Hesthaven et al. [17]. Here numerical simulations can indicate the existence of possible vanishing viscosity or vanishing dispersion limits in bounded domains. Kaneda and Ishihara [22] have also used spectral methods to examine scaling laws for turbulent flow with periodic boundary conditions. As explained by Yeung [42] and by Jim´nez and Moser [20], it is also of interest e to examine wall bounded flows. This can be done using Chebyshev spectral methods. In particular, as shown by De la Cruz et al. [9], it is possible to parallelize multidimensional spectral collocation methods efficiently, so that it is feasible to perform large simulations quickly. Following this introduction is a review of previous studies of Chebyshev collocation methods. The next section contains a description and comparison of the accuracy of the Chebyshev spectral integral, preconditioned Chebyshev differentiation and real space Chebyshev differentiation methods in solving linear boundary value problems. An IMEX and a fully implicit time-stepping scheme that can use these Chebyshev spatial discretizations for solving initial boundary value problems are then described. Results of a numerical examination of the spatial and temporal convergence of the IMEX time-stepping scheme for a model problem from the dynamics of phase transformations are summarized. In the final section we show that the fully implicit scheme can be 2

used to simulate problems with stiff nonlinearities for which the IMEX scheme does not converge.

2

Previous Work

There have been many studies of preconditioned Chebyshev and Chebyshev integration methods for boundary value problems, but none of these studies has numerically examined the accuracy of these methods in norms that include derivatives. The preconditioned Chebyshev tau method (hereafter referred to as the preconditioned Chebyshev differentiation to contrast it to the Chebyshev integration method) is described by Gottlieb and Orszag [14, p. 119] and by Canuto et al. [2, p. 173], and has been extended to general orthogonal polynomial expansions by Coutsias et al. [5]. Funaro and Heinrichs [12] and Tuckerman [38] have also discussed similar preconditioning methods for other orthogonal polynomial systems. Both the Chebyshev integration and preconditioned Chebyshev differentiation methods allow for the solution of linear boundary value problems in Chebyshev spectral space in O(N ) operations, and thus by using a Fast Fourier Transform, in O(N log(N )) operations in real space. Unfortunately, the preconditioned Chebyshev differentiation method does not immediately give derivatives of the function, and these must be obtained by differentiation. As has been noted in many studies, numerical differentiation of Chebyshev interpolants is sensitive to errors introduced by finite precision arithmetic in both real space (see for example Trefethen and Trummer [35] or Weideman and Trefethen [40]) and in spectral space if the spectral coefficients are not carefully computed (see for example Coutsias et al. [5] or Hesthaven et al. [17, p. 217]). Greengard [15] showed that, by using the Chebyshev spectral integration method to obtain numerical solutions to two-point linear boundary value problems, it is possible to calculate derivatives without the numerical instability associated with differentiation. Coutsias et al. [5] generalized Greengards results to expansions in other orthogonal polynomials and suggested methods to solve the resulting linear systems efficiently. It should be noted that Coutsias et al. [5] refer to the Chebyshev integration method as the postconditioned Chebyshev method to contrast it to the preconditioned Chebyshev method. However, the term integration better captures the notion that a smoothing operation which does not amplify errors occurs, and so this will be used here. The Chebyshev integration method has so far primarily been used in IMEX schemes for initial boundary value problems with at most two spatial derivatives, for example by Cox and Matthews [6] and by Lundbladh et al. [28]. Clenshaw [4] and Elliot [11] have also used versions of the Chebyshev integration method for linear two-point boundary value problems and for the time-dependent heat equation respectively; however, because their papers were 3

published in 1957 and 1960, the advantages of their method in comparison to other numerical methods in use today are not indicated. Zebib [43] has also shown that, by solving for the highest derivative in a fourthorder nonlinear boundary value problem and then integrating to obtain the lower order derivatives, the accuracy of Galerkin solutions to nonlinear boundary value problems can be improved. When using a full Galerkin method, a nonlinear iterative solution of the resulting equations is required and Zebib [43] found that it was computationally expensive to use a large number of modes when high spatial resolution simulations were required. Mai-Duy [30] and Mai-Duy and Tanner [31] have also compared the accuracy in the L2 norm of Chebyshev collocation differentiation and Chebyshev collocation integration methods to obtain solutions to linear fourth-order boundary value problems. They find that the Chebyshev integration method gives more accurate results than the Chebyshev collocation differentiation method. They do not use a large number of modes nor do they examine the accuracy of the method in norms which include derivatives, so they do not show when the extra effort of using the integration method instead of the collocation differentiation method is justified.

3

Boundary Value problems

In this section, the Chebyshev integration method is described and the differences between the Chebyshev integration, Chebyshev differentiation and preconditioned Chebyshev differentiation methods are explained. It is implicitly assumed that space is discretized using a collocation scheme with Chebyshev polynomials of the first kind discretized at Chebyshev Gauss-Lobatto points. The reason for using this discretization is that it allows the use of the Fast Fourier Transform to calculate integrals and derivatives when computing solutions. 3.1 Spatial Discretization

We explain the differences between the differentiation matrix method, the preconditioned Chebyshev differentiation method, and the Chebyshev integration method by considering the linear boundary value problem uxxxx + uxx + u = f (x) u(−1) = 0, u(1) = 0, (3.1) ux (−1) = 0, ux (1) = 0,

where u is the displacement, x the position and f (x) is a smooth but otherwise unspecified function. The description given here does not include all the details 4

because these can be found in appendices A and B for the pre-conditioned Chebyshev differentiation and Chebyshev integration methods respectively. A full description of the Chebyshev differentiation method can be found in Trefethen [33]. The Chebyshev differentiation matrix method as described by Trefethen [33, p. 145] amounts to solving the following linear matrix equation ¯ (D 4 + D 2 + I)u = f . (3.2)

Here D k is the real space Chebyshev differentiation matrix of order k, I is the identity matrix, u is a vector with the approximate solution values for u at the nodal points, and f is a vector with the forcing function values at the nodal points. As explained by Trefethen [33, p. 58] and by Canuto et al. [2, p. 88], the formulas for the entries in the matrix D can be derived by differentiating interpolating polynomials and evaluating the derivatives at the nodal points. In all computational experiments considered here, D was obtained using the function cheb.m which can be found in Trefethen [33, p. 58]. Note that the solution of the linear system in eq. (3.2), only gives u and not its ¯ derivatives. Note also that since D 4 = (D 1 )4 , the highest order differentiation matrix is modified to ensure that the solution satisfies clamped boundary conditions. To be precise, as explained by Trefethen [33, p. 146] we restrict our search for polynomial interpolants that satisfy the boundary conditions by obtaining solutions in polynomials that have (1 − x2 ) as a factor. Thus, if xi ∈ [0, 1] is the ith Chebyshev Gauss-Lobatto nodes and x is a vector with the node locations then ¯ D 4 = diag(1 − x2 )D 4 − 8diag(x)D 3 − 12D 2 × diag 1/(1 − x2 ) , where diag(x) is a diagonal matrix whose entries are the quadrature points. As shown in Fig. 1(a), the resulting differentiation matrices are full, so solving the large linear systems using these matrices can take some time. To obtain the preconditioned differentiation method as described by Gottlieb and Orszag [14, p. 119], eq. (3.2) is transformed into an equation for the coeffcients of the truncated Chebyshev expansions for u and f (x), denoted by ˆ u and f respectively. The resulting infinite system of equations is truncated, ˆ to obtain ˆ ˆ ˆ D4 + D2 + I u = f ˆ (3.3) ˆ where D k is the k th order Chebyshev differentiation matrix in Chebyshev space (the entries of these spectral differentiation matrices can be found in Hesthaven et al. [17, p. 258] or Coutsias et al. [5]). The resulting equations are then multiplied by the highest order integration matrix to obtain the sparse 5

(a) Discretized real space bound- (b) Discretized Chebyshev space ary value problem differentiation preconditioned boundary value ˜ operator D 4 + D 2 + I. problem differentiation operator 0 2 4 ˆ Pˆ + Pˆ + Pˆ + LBC1. D D D

(c) Discretized Chebyshev space (d) Discretized Chebyshev space ˜ differentiation operator D. boundary value problem integra˜ ˜ ˜ tion operator S 0 + S 2 + S 4 + ˆ LBC2. Fig. 1. Sparsity patterns with 256 Chebyshev modes

system of equations
4 2 0 4 ˆ ˆ Pˆ + Pˆ + Pˆ + LBC1 u = Pˆ f + RBC1. D D D ˆ Dˆ

(3.4)

Here Pˆ is the preconditioned Chebyshev differentiation matrix of order k D in which the four rows for the coeffcients of the four highest order Chebyshev polynomials have been set to zero, u is the vector with coefficients for each ˆ ˆ Chebyshev polynomial for the approximation of u and f is the vector with the coefficients for each Chebyshev polynomial for the approximation of f (x). The reason for setting these rows to zero is that the equations for the highest modes are used to enforce the boundary conditions exactly, instead of satis6

k

ˆ fying the differential equation. Thus LBC1 contains the coeffcients for the ˆ boundary conditions in Chebyshev space and RBC1 contains the values of 4 2 0 ˆ these boundary conditions. Full expressions for Pˆ , Pˆ , Pˆ , RBC1 D D D ˆ and LBC1 can be found in appendix A. Note again that, this linear system is solved to find an approximate solution to u, and thus this solution needs to be differentiated numerically to obtain approximations of the derivatives of the exact solution. This differentiation can either be done in spectral space by using the spectral differentiation matrices, which, if there are N modes, requires O(N 2 ) operations, or by transforming to real space and using the Fast Fourier Transform to differentiate the resulting series, (see for example, Trefethen [33, p. 78]) which requires O(N log N ) operations. The sparsity pattern for the preconditioned matrix operator is shown in Fig. 1(b), which shows that the solution of this system of equations requires O(N ) operations. To obtain derivatives using spectral differentiation matrices requires O(N 2 ) operations, because in this case, as shown in Fig. 1(c), the matrix has a full upper triangular structure. Greengard [15] explains that, the Chebyshev integration method amounts to solving for the highest order derivative once eq. (3.2) is transformed into an equation for the truncated Chebyshev expansions of u and f (x) ˆ ˆ ˆ ˆ ˆ ˆ S 0 + S 2 + S 4 + LBC2 uxxxx = f + RBC2. ˆ (3.5)

ˆ In this equation, LBC2 is a matrix with the equations that the coeffcients ˆ of the Chebyshev expansion should satisfy, RBC2 is a vector with the values of these boundary conditions and uxxxx is a vector with the coeffcients of the ˆ ˆ ˆ truncated series expansion for uxxxx . The matrix LBC2 and vector RBC2 fix the four coeffcients obtained from the indefinite integral of uxxxx by using the boundary conditions. This linear system is solved to find uxxxx , which is then ˆ integrated to find u and its first three derivatives. Full details on how to do this ˆ are in appendix B. An important observation in this derivation is that, because the Chebyshev basis is a polynomial basis, the four integration constants, c1 , c2 x, c3 x2 and c4 x3 only involve combinations of the low order Chebyshev polynomials, and hence the integration matrices remain sparse. This is not true for a non-periodic function whose highest derivative is expanded in a Fourier series with lower order derivatives being obtained by integration. Finally, in the preconditioned Chebyshev method, if the truncated expansion for u has N + 1 modes then a (N + 1) × (N + 1) linear system is solved. In the ˆ Chebyshev integration method however, four further equations are obtained because of the integration constants, thus if the truncated expansion for uxxxx ˆ has N + 1 modes, then a (N + 5) × (N + 5) matrix system is solved. The typical sparsity pattern for the matrix obtained when using the Chebyshev integration method matrix is shown in Fig. 1(d) and is similar to the sparsity pattern obtained using the preconditioned Chebyshev differentiation method 7

shown in Fig. 1(b) — in both cases the top four rows are full and the diagonal contains five dense bands. 3.2 Numerical Results

This section contains a numerical comparison of the accuracy in the L2 and W 2,2 norms of the Chebyshev differentiation, preconditioned Chebyshev differentiation and Chebyshev integration methods for solving two linear boundary value problems. When solving the linear system for the Chebyshev integration method in this numerical comparison, a solution which satisfies all the boundary conditions is obtained. Variations of the Chebyshev integration method in which solutions of the homogeneous boundary value problem are added to a particular solution of the inhomogeneous boundary value problem, as described by Coutsias et al. [5] are not reported here. This is because for the linear boundary value problems considered here, these modifications do not affect the accuracy of the numerical solutions that are obtained. Numerical comparisons of the accuracy of approximate derivatives of the solution obtained by the preconditioned Chebyshev differentiation method are obtained using the Fast Fourier Transform and obtained using spectral differentiation matrices are included, because in this case the results differ. Calculations in this section were done in MATLAB 7.3 on a laptop with a 2GHz Intel Core Duo processor and 2Gb of RAM. MATLABs backslash was used to obtain solutions to the resulting linear systems of equations. We examine numerical approximations to the following two boundary value problems, with the associated boundary conditions and given exact solutions, uxxxx + 2uxx + u = cos(x), u(−1) = 0, u(1) = 0, ux (−1) = 0, ux (1) = 0, (3.6)

u(x) = and

4x cos2 (1) sin(x) − cos(x) {sin(2) − 2 + x2 [2 + sin(2)]} , 8 [2 + sin(2)]

50−4 uxxxx − u = 10, u(−1) = 0, u(1) = 0, ux (−1) = 0, ux (1) = 0,

(3.7)

u(x) =

10 sinh(50) cos(50x) + 10 sin(50) cosh(50x) − 10. sin(50) cosh(50) + sinh(50) cos(50)

8

0.05

10

0

0.04

10

−5

0.03 Error 0.02 10
−15

10

u

−10

0.01

0 −1

−0.5

0 x

0.5

1

10

−20

10

0

10 10 10 Number of Chebyshev Modes

1

2

3

10

4

(a) Exact solution to eq. (3.6).

(b) Numerical convergence for eq. (3.6).
10
5

5 0

10 −5 Error −10 −15 10 −20 −25 −1 10 10 u

0

−5

L2 error Integration W2,2 error Integration
−10

L2 error Differentiation W W
2,2

error Differentiation error Preconditioned 1
1 2 3 4

L2 error Preconditioned
2,2

−15

W2,2 error Preconditioned 2
0

−0.5

0 x

0.5

1

10

10 10 10 Number of Chebyshev Modes

10

(c) Exact solution to eq. (3.7).

(d) Numerical convergence for eq. (3.7).

Fig. 2. Convergence for linear boundary value problems. In the legend, Integration — Chebyshev integration method where the bordered banded matrix solved in one step; Differentiation — real space collocation differentiation method; Preconditioned 1 — preconditioned differentiation method with derivatives obtained using spectral differentiation matrices; and Preconditioned 2 — preconditioned differentiation method with derivatives obtained using the Fast Fourier Transform.

The exact solutions are plotted in Figs. 2(a) and 2(c). Figures 2(b) and 2(d) show the convergence of numerical solutions to these linear boundary value problems obtained using Chebyshev differentiation, preconditioned Chebyshev differentiation and Chebyshev integration methods respectively. Trefethen and Trummer [35] have shown that using large real space Chebyshev differentiation matrices in numerical calculations gives results that are sensitive to errors introduced by finite precision arithmetic. The figures show results using as many modes as could be used with the available amout of memory. Trefethen and Trummer [35] also note that differentiation of the Chebyshev interpolant of a function using the Fast Fourier Transform is sensitive to errors introduced by finite precision arithmetic, thus the numerical instability observed here cannot be alleviated by using this method to cal9

culate derivatives. Greengard [15] showed that, provided the resulting linear systems are solved accurately, by using an integration matrix formulation the numerical instabilities that occur when obtaining the solution and derivatives of the solution to two-point boundary value problems are avoided. The figures confirm these previous results when extended to fourth order boundary value problems and show that the real space Chebyshev differentiation matrix method is a poor method when a large number of modes are used in both L2 and W 2,2 norms. The preconditioned Chebyshev differentiation method is a good method when accurate results in the L2 norm are required. If derivatives are obtained using the Fast Fourier Transform, the preconditioned Chebyshev method is a poor method to use when accurate approximations of the derivative are required. Finally, the figures show that accurate results in the W 2,2 norm can be obtained using the preconditioned Chebyshev differentation method when derivatives are found using spectral differentiation matrices in spectral space. As explained by Hesthaven et al. [17, p. 220], the reason for the difference in the accuracy of real space differentiation and spectral space differentiation is that, even though both differentiation matrices are ill conditioned, the accurately computed spectral expansion coefficients decay sufficiently rapidly to ensure that when they are multiplied by the spectral differentiation matrices, the results remain accurate. If the Fast Fourier Transform is used to obtain the spectral expansion coefficients, a computational error is introduced into the calculated spectral expansion coefficients which is then amplified upon differentiation (see Higham [18, p. 451] for an analysis of errors introduced by a simple Fast Fourier Transform algorithm).

Since the Chebyshev integration method uses more modes, N + 5, as opposed to N +1 for the preconditioned Chebyshev differentiation method, it is surprising that the Chebyshev integration method is slightly less accurate than the preconditioned Chebyshev differentiation method when more than 100 modes are used. In particular, Canuto et al. [2, p. 177] suggest that the Chebyshev integration method should be more accurate because it has more degrees of freedom. Figure 2(d) shows that the preconditioned Chebyshev differentiation method is less accurate when there are less than 100 modes, because here the extra degrees of freedom give better resolution. When more than 100 modes are used, the preconditioned method seems to give a matrix which can be solved more accurately. As we have not done a full error analysis of the solution procedure used by MATLABs “backslash” for these matrix systems, the reason for this difference in accuracy is unclear. It is likely that since the preconditioned and Chebyshev integration linear systems have similar matrix 0 structures, the faster decay of the terms on the right hand side, Pˆ f in eq. Dˆ ˆ (3.4) as opposed to f in eq. (3.5), makes it easier to solve the preconditioned matrix system more accurately than the Chebyshev integration matrix system (see, for example, Higham [18, p. 120]). 10

4

Initial boundary value problems

In this section, the performance of an IMEX time-stepping scheme that uses the Chebyshev integration method and pre-conditioned Chebyshev differentiation methods to solve semilinear initial boundary value problems is examined. The accuracy and computational cost of two different implementations of the Chebyshev integration method and of the preconditioned Chebyshev differentiation method for IMEX time-stepping schemes are compared. The modification made to the Chebyshev integration method is to combine a particular solution of the inhomogeneous boundary value problem solved at each time step to solutions of the homogeneous boundary value problem presolved before timestepping began so that as suggested by Coutsias et al. [5], only a banded and not a bordered banded linear system is solved at each timestep. The comparison also includes a fully implicit time-stepping scheme using the Chebyshev integration spatial discretization. This section ends with an example of a semilinear equation with a stiff nonlinear term for which only the fully implicit scheme can be used. 4.1 An implicit-explicit temporal discretization

A simple second order time-stepping scheme is constructed for the equation ρutt − βuxxt = γ 2 u3 − ux x
x

−

2

uxxxx ,

(4.1)

where x ∈ [−1, 1] and t ∈ [0, 1] are the spatial and variables and where ρ, β, γ and are real positive constants. This equation is used as a simplified dynamic model for phase transformations and has an exact traveling wave solution obtained by Truskinovskii [36] which we rewrite as u(x, t) = √ 2 log cosh   γ
  

κx − ωt κ2

ω2 2
2

ρ−

β2 6
2

+

γ 2 κ2 2

  ωβ(κx − ωt) − √ . 2  3κ2 2

(4.2) In this solution, κ is the wavenumber and ω the wave frequency. The boundary conditions, u(−1, t), u(1, t), ux (−1, t) and ux (1, t), and initial conditions u(x, 0) and ut (x, 0) are obtained from this exact solution. Since single domain spectral collocation methods cannot approximate discontinuous functions well, only smooth solutions are used for the numerical comparison, for which the inequality, ω 2 (ρ − β 2 /6 2 ) + γ 2 κ2 ≥ 0, must hold. To construct a second-order IMEX time-stepping scheme, time dependent terms are discretized using second-order finite difference approximations centered at the next time step. This is done using second-order backward differentiation approximations for the terms utt and uxxt , and an Adams-Bashforth or 11

forward extrapolation method for the nonlinear term (u3 − ux )x . The resulting x time-stepping scheme is 2uj+1 − 5uj + 4uj−1 − uj−2 3uj+1 − 4uj + uj−1 xx xx − β xx 2 δt 2δt = 2γ 2 (uj )3 − uj − γ 2 (uj−1 )3 − uj−1 − 2 uj+1 . x x x x xxxx ρ
x x

(4.3)

In this equation, the superscript on u denotes the time step at which the function is evaluated. This time discretized eq. (4.3) is rearranged to obtain the new iterate uj+1 in terms of the previous two iterates uj , uj−1 and uj−2 uj+1 uj+1 − 3β xx + 2 uj+1 xxxx δt2 2δt j j−1 j−2 5u − 4u + u 4uj − uj−1 xx =ρ − β xx 2 δt 2δt + 2γ 2 (uj )3 − uj − γ 2 (uj−1 )3 − uj−1 x x x x 2ρ
x

(4.4)
x

.

This shows that at each time step, a linear boundary value problem is solved and so the Chebyshev integration method can be used. The descriptions of the implementation of Chebyshev integration methods in IMEX timestepping schemes in the literature differ (see, for example, Coutsias et al. [5], Cox and Matthew [6] and Lundbladh et al. [28]). For completeness, a description of the implementation of the numerical scheme in eq. (4.4) using the Chebyshev integration method is now provided. A similar approach is used for the preconditioned Chebyshev differentiation method and so this is not included. The highest spatial derivative term in eq. (4.4) is approximated N j by a truncated Chebyshev series, uj n=0 an Tn (x). Lower derivative xxxx ≈ approximations are found by integration as explained in appendix B. The vector g j is then defined as the sum of the terms in the scheme that were ˆ computed in the previous three time steps in Chebyshev space g j := ρ ˆ 5ˆj − 4ˆj−1 + uj−2 u u ˆ 4ˆj − uj−1 u ˆxx − β xx 2 δt 2δt

+ 2γ 2 3(uj )2 uj − uj − γ 2 3(uj−1 )2 uj−1 − uj−1 , x xx xx x xx xx where the vectors uj , and uj are the vectors of coefficients of the truncated ˆ ˆxx Chebyshev polynomials which approximate the solution and its second derivative at the j th time step. The vectors uj , and uj are the real space approximax xx tions to ux and uxx at the j th time step at the N + 1 grid points.The nonlinear 2 term, 3 (uj ) uj − uj , is projected onto the space of N + 1 Chebyshev polyx xx xx nomials by using a Fast Fourier Transform, see, for example, Trefethen [33, p. 75]. In doing so, aliasing errors are neglected since these are expected to be small if the solution is sufficiently well resolved, as explained by Canuto et 12

al. [2, p. 163]. A sparse, bordered banded system of equations, ˆ uxxxx ˆ ˆ j+1 Lˆj+1 = g j + RBC2 , where 3β ˆ2 2ρ ˆ4 ˆ ˆ S0 − S + 2 S + LBC2, 2δt δt is then obtained. In these equations: uj+1 is the vector of coefficients of the ˆxxxx Chebyshev polynomials which approximate the solution for the fourth derivaj ˆˆ tive at the j th time step; g j is a vector with entries gn ; the matrices S 0 , S 2 ˆ 4 ˆ ˆ and S are the zeroth, second and fourth order integration matrices; LBC2 j+1 ˆ is a matrix with the equations for the boundary conditions; and RBC2 is a vector with the values of the boundary conditions. Explicit expressions for ˆˆˆ ˆ ˆ j+1 are given in appendix B. S 0 , S 2 , S 4 , LBC2 and RBC2 ˆ L :=
2

(4.5)

In typical implementations of time-stepping schemes, eq. (4.4) is multiplied ˆ by δt2 . Here by not doing so, the smallest non-zero matrix entries in S 4 do 2 ˆ0 not become too small relative to the entries in S . Finally, even though 2 ˆ ˆ multiplies the term S 0 , the non-zero entries of S 0 are of order one, and so small values of can be used and the linear system still solved using finite precision floating point arithmetic. At each time step, the sparse linear system given by eq. (4.5) is solved and used to advance in time in the space of Chebyshev polynomials. It is important to ˆ note that, for the linear operator L in eq. (4.5) used in this implicit time step to be invertible, it must include the integration matrix for the highest derivaˆ ˆ ˆ tive term, S 0 . This is because only the combination S 0 + LBC2 has full rank. Thus the integration matrix formulation as described here cannot be used in an explicit time-stepping scheme. Coutsias et al. [5] explain that, it is possible to restrict the linear operator to a subspace where it is invertible, which would allow an explicit time-stepping scheme to be used. However, stability constraints make explicit time-stepping schemes unattractive for the numerical solution of stiff problems that occur in models for phase transformations with a fourth order spatial derivative. When solving a system arising from the Chebyshev integration method with different right hand sides repeatedly as is done in IMEX time-stepping, Coutsias et al. [5] suggest that it may be faster to pre-solve a homogeneous boundary value problem first to obtain a set of solutions which can then be added to satisfy a variety of boundary conditions. This is because at each time step one can then solve a well conditioned banded diagonal matrix system. The particular solution of the well conditioned inhomogeneous boundary value problem does not satisfy the boundary conditions because the top four rows of the matrix are chosen to simplify the structure of the matrix. In the implementation considered here, the first four diagonal entries are one, and all other 13

entries in the first four rows are zero. However, solutions of the homogeneous boundary value problem can then be added to a particular solution of the inhomogeneous boundary value problem to obtain the solution that satisfies the boundary conditions. Cox and Matthews [6] have implemented a variant of this method using analytic solutions for the homogeneous boundary value problem when solving the Navier-Stokes equations but, possibly unaware of Coutsias et al.’s [5] results, did not choose the matrices that are inverted to solve the inhomogeneous boundary value problem to have a banded structure. Lundbladh et al. [28] use two different particular solutions to the inhomogeneous boundary value problem and take a linear combination of these so as to satisfy the boundary values. As explained by Cox and Matthews [6], Lundbladh et al.’s method requires more computations, because two and not one linear system are solved at each time step. When using the Chebyshev integration method to do time-stepping, as indicated in eq. (4.4), a linear boundary value problem is solved at each timestep. This can be done in several ways, which change the structure of the linear system that is solved. If the full matrix system is solved at each time step, a bordered banded system is obtained. As explained by Coutsias et al. [5], to solve a banded matrix system at each time step, one first presolves the linear boundary value problem in eq. (4.4) with a right hand side that is zero to obtain four independent solutions to the homogeneous problem,
2 hi uxxxx

−

2ρ 3β hi uxx + 2 uhi = 0, 2δt δt

(4.6)

where i can be equal to 1, 2, 3 or 4. The four corresponding sets of boundary conditions are uh1 (−1) = 1, uh2 (−1) = 0, uh3 (−1) = 0, uh4 (−1) = 0, uh1 (1) = 0, uh2 (1) = 1, uh3 (1) = 0, uh4 (1) = 0, uh1 (−1) = 0, xx uh2 (−1) = 0, xx uh3 (−1) = 1, xx uh4 (−1) = 0, xx uh1 (1) = 0, xx uh2 (1) = 0, xx uh3 (1) = 0, xx uh4 (1) = 1. xx (4.7) (4.8) (4.9) (4.10)

These solutions can either be computed numerically, as suggested by Coutsias et al. [5], or analytically, as suggested by Cox and Matthews [6] before timestepping begins. To obtain the particular solution to eq. (4.4), the first four diagonal entries in the matrix that are used to fix the boundary conditions, ˆ LBC2 in eq. (4.5), are set to be one. All other entries in these rows are zero. Similarly, the four entries for the boundary conditions in the right hand side ˆ j+1 in eq. (4.5) are set to zero. vector RBC2 In the numerical examples computed using the bordered banded matrices, the solution of the matrix system was done by the sparse LU factorization implemented in MATLAB. The banded linear systems were also solved by a sparse LU factorization using MATLABs backslash, but the default settings 14

were changed to ensure that a band solver was used. Other solvers such as KLU and UMFPACK 5.0.2 (this is a slightly newer version of UMFPACK than that which is installed in MATLAB 7.3) which are part of Suitesparse [7] were also tested and found to give similar results. Since the best sparse solution method is architecture and solver setting dependent, we have not done an extensive numerical comparisons of sparse solvers. An introduction to current sparse solver technology can be found in Davis [8]. In starting the multi-step time-stepping scheme, the functions u0 , u−1 and u−2 need to be obtained. The initial conditions give u0 and u0 . To obtain t u−1 and u−2 , we follow Kress [24] and use a Taylor expansion, the initial velocity and the partial differential equation to extrapolate backwards in time and still have O(δt2 ) accuracy. Thus, u−1 := u0 − δtu0 + t and u
−2

δt2 0 u 2 tt

(4.11)

:= u −

0

2δtu0 t

(2δt)2 0 + utt . 2

(4.12)

The initial acceleration, u0 , is not given explicitly by the initial conditions, tt but by differentiating the initial conditions u0 and u0 with respect to x, u0 , t x u0 and u0 can be calculated and then the partial differential equation xxxx xxt used to find that, u0 = ρ−1 βu0 + γ 2 (u0 )3 − u0 tt xxt x x
x

−

2

u0 xxxx .

(4.13)

This can then be used to calculate u−1 and u−2 so that a second order accurate numerical approximation of the true solution can be calculated.

4.2

A fully implicit temporal discretization

Time stepping schemes for variational problems that nearly conserve momentum or nearly conserve energy are often used in computational solid mechanics where it is expected that dissipation is small. Reviews of these methods can be found in Hairer, Lubich and Wanner [16, p. 204], Kane et al. [21], Leok [25], Lew [26], Lew et al. [27] and Marsden and West [32]. In most solid mechanics applications, explicit time integration schemes are used because it is expensive to solve the resulting systems. Here, a fully implicit scheme is used because stiff nonlinear and higher order capillarity terms can make explicit and IMEX schemes unfavorable. Note that, because a Chebyshev spatial discretization is used, a full analysis of the scheme would show that it is variational in a 15

weighted energy space. However, because well resolved spectral spatial discretizations give very accurate approximations of analytic solutions, it is expected that the scheme will retain the good energy dissipation properties that variational integration schemes have.

To introduce a stiff nonlinear term, we consider the equation,

ρutt − β uxt 1 + 2u2 + u4 x x =γ
2

x

u3 x

− ux

x

−

2

uxxxx + f (x, t),

(4.14)

for x ∈ [−1, 1] and t ∈ [0, 1]. The forcing function, f (x, t), the initial and boundary conditions are chosen so that eq. (4.14) has the exact solution

u = sin [πxt] + cos [4π (x − 3t)] .

(4.15)

The dissipation term in eq. (4.14),

uxt 1 + 2u2 + u4 x x

x

,

differs from that in eq. (4.1) because it is nonlinear. This unusual dissipation term is chosen because in frame indifferent viscoelastic models, the dissipation term is also nonlinear (see, for example, Lew [26] and Lew et al. [27]). This equation will allow us to show that for these types of equations it is favorable to use a fully implicit scheme. 16

Equation (4.14) is approximated using the following temporal discretization,

uj+1 − 2uj + uj−1 δt2  uj+1 + uj β uj+1 − uj  xx x xx 1+2 x − 2 δt 2 ρ uj + uj−1 β uj − uj−1  xx xx x 1+2 x − 2 δt 2 − 2β uj+1 x − δt uj x
   

2

uj+1 + uj x x + 2 uj + uj−1 x x + 2 + 2 uj  x
 

4

 

2

4

 

uj+1 x

+ 2

uj x

3

+
3

uj+1 x

uj+1 + uj xx xx 2

uj − uj−1  uj + uj−1 x x x − 2β x δt 2 γ 2  uj+1 + uj x x = 2 2 + − γ 2
2 2

uj + uj−1  uj + uj−1 x xx xx +x 2 2

x



3

uj+1 + uj  x x − 2
3



uj x

+ 2

uj−1 x

−

uj x

+ 2

uj−1 x


x

4

j−1 uj+1 + 2uj xxxx xxxx + uxxxx +

f x, tj+1/2 + f x, tj−1/2 2

.

(4.16)

This time discretization can be obtained by approximating the LagrangeD’Alembert variational principle using the generalized midpoint rule, see Kane et al. [21], Lew [26] and Marsden and West [32] for further details. An error analysis using finite difference approximations in time at uj can be used to show that this scheme is second order accurate.

To implement this scheme numerically, the nonlinear terms are obtained by using fixed point iterations. We shall let uj+1,k+1 denote iterate k + 1 for u at timestep j + 1. Then, in these iterations a linear boundary value problem is solved to obtain uj+1,k+1 , uj+1,k+1 and uj+1,k+1 until uj+1,k+1 −uj+1,k converges xx xxxx 17

to a specified tolerance in the W 2,2 norm. The fixed point iteration scheme is
2 uj+1,k+1 uj+1,k+1 + uj+1,k+1 − β xx δt2 2δt 4 xxxx j j−1 2 j−1 2u − u u j−1 =ρ 2uj − β xx − xxxx + uxxxx 2 δt 2δt 4  2 j+1,k j j+1,k + uj uj+1,k + uj β uxx − uxx  ux x x x 2 + + 2 δt 2 2

ρ

4

 

β uj − uj−1  uj + uj−1 xx xx x + 2x 2 δt 2 uj+1,k − uj  uj+1,k + uj x x x + 2β x δt 2 uj − uj−1  uj + uj−1 x x x + 2β x δt 2
j+1,k + uj γ 2  ux x + 2 2



2

uj + uj−1 x x + 2
3

4

 



uj+1,k + uj  uj+1,k + uj x xx xx +x 2 2






3

uj + uj−1  uj + uj−1 xx xx x +x 2 2

x



3

uj+1,k + uj  x x − 2 uj x
j−1/2

+

γ 2

2



uj x

+ 2

j−1 ux

3

−

+ 2 .

uj−1 x

 
x

+

f x, t

j+1/2

+ f x, t 2

(4.17)

In starting the fixed point iterations, the first iterates for u, ux and uxx are obtained by extrapolation from the previous two solutions, for example, uj+1,1 = 2uj − uj−1 . For this fully implicit scheme we will only use the Chebyshev integration formulation with a bordered banded matrix system. In the numerical examples, convergence of the iterative scheme was deemed to have been achieved once uj,k+1 − uj,k 2,2 < 10−14 .
W

Since the solution procedure for the linear systems found at each iterate is the same as for the IMEX scheme detailed in the previous subsection, that discussion is not repeated here. Similarly, the second order accurate starting values are obtained using extrapolation as shown in eqs. (4.11), (4.12) and (4.13).

4.3

Numerical convergence results for initial boundary value problems.

In this section, the numerical convergence of the Chebyshev spectral integration and preconditioned Chebyshev differentiation methods to the smooth exact one dimensional traveling wave solution of eq. (4.1) with the exact so18

4 3 2 u 1 0 −1 1 1 0 x −1 0 t 0.5

Fig. 3. Exact traveling wave solution to eq. (4.1). A close-up view of the traveling kink is shown, away from this region the solution is almost planar.

lution eq. (4.2) is examined. The convergence of the fully implicit scheme to the solution of eq. (4.14), which has nonlinear dissipation term, is also demonstrated. A numerical solution to eq. (4.1) was simulated in MATLAB for x ∈ [−1, 1] and from t = 0 to t = 1 using boundary conditions from the exact solution. The exact solution was also used to provide the initial conditions u0 and u0 t and from these, the startup values for the numerical approximation scheme were obtained by extrapolation. The parameter values in the simulations were ρ = 2, = 0.1, β = 0.001, ω = 1 and κ = 1, and a plot of the exact solution is in Fig. 3. Smaller values of and larger values of ρ were also tried. However, the numerical results showed that the Chebyshev integration scheme which combines homogenous and particular solutions performs even more poorly in these cases, while the other schemes retained the same trends. We have therefore chosen to present a comparison using these values, even though, in typical applications ρ will be larger and smaller. In these simulations, the maximum difference in the L2 and W 2,2 norms between the exact solution and the numerical solution during the time interval of simulation were calculated and are plotted in Fig. 4. Figures 4(b) and 4(c) show the results for simulations done with 1024 Chebyshev modes to assess how the accuracy and computation time changed as the timestep was reduced. Figures 4(c) and 4(d) show the results of simulations done with a fixed timestep of 10−3 to assess how the accuracy and computation time changed as the number of modes increased. In measuring the computational time in these simulations, only the time required for time-stepping was measured, because the main interest is in estimating computation times for simulations which require many timesteps and for which the setup time is a negligible proportion of the total time required. Figure 4(a) shows that the preconditioned Chebyshev method is as accurate as the Chebyshev integration method when the full linear system is solved, but is computationally more expensive. This is expected because the Chebyshev 19

10 10 10 Error 10 10 10

4

10

3

2

0

−2

Time for computation
−4 −3 −2

10

2

10

1

−4

10

0

−6

10 −5 10

−8

10

10 Time step

10

10 −5 10

−1

10

−4

10 Time step

−3

10

−2

(a) Error vs. timestep with 1024 (b) Computation time vs. timestep Chebyshev modes. with 1024 Chebyshev modes.
10
60

L!([0,1],L2),bbi L ([0,1],W 10
40

10

3

! !

2,2

),bbi Time for computation 10
2

L ([0,1],L ),bia L!([0,1],W2,2),bia L ([0,1],L ),sd
! 2

2

bbi bia sd bbiImp

Error

10

20

L!([0,1],W2,2),sd L!([0,1],L2),bbiImp L ([0,1],W
! 2,2

10

1

),bbiImp

10

0

10

0

10

−20

10

1

10 10 Number of modes

2

3

10

4

10

−1

10

1

10 10 Number of modes

2

3

10

4

(c) Error vs. number of Chebyshev (d) Computation time vs. number of modes with a time step of 10−3 . Chebyshev modes with a time step of 10−3 . Fig. 4. Temporal and spatial convergence results for an initial boundary value problem. In the legend, bbi – Chebyshev integration method using the IMEX time stepping scheme where the linear systems are solved by LU factorization of the bordered banded matrix; bia – Chebyshev integration method using the IMEX time stepping scheme where the homogeneous equation is pre-solved analytically and then the inhomogeneous equations are solved by sparse LU factorization at each time step; sd – preconditioned differentiation method using the IMEX time stepping scheme where derivatives are obtained by spectral differentiation matrices at each time step; bbiImp - Chebyshev integration method using the fully implicit time stepping scheme and the linear systems are solved by LU factorization of the bordered banded matrix.

differentiation matrices in spectral space are dense upper triangular matrices and the Chebyshev integration matrices are sparse. Surprisingly, the actual computational cost of solving the bordered banded matrix system that immediately satisfies the boundary conditions was less than first pre-solving to obtain homogeneous solutions and then repeatedly solving banded systems for particular solutions to which linear combinations of the solutions to the homogeneous equation were added to satisfy the boundary conditions. The reason for the difference is that the computational cost of adding the homogeneous 20

solutions to the particular solution outweighs the savings in solving a banded matrix system. Figures 4(a) and 4(b) show that although for a given timestep the fully implicit scheme is computationally more costly than the IMEX schemes, it has better L2 accuracy. In particular, for a time step of 10−4 , the implicit scheme is five times more accurate in the L2 norm than the IMEX scheme but it only takes twice as long to compute. The results in the W 2,2 norm are not so clear since for large time steps, the IMEX method is more accurate but for small timesteps the implicit scheme is more accurate. Again, examining Figure 4(a) we find that combining solutions of the homogeneous boundary value problem with the solution of the inhomogeneous boundary value problem solved at each timestep gives poor accuracy when the time step is less than 10−3 . The reason for the poor accuracy in the time dependent case is that rounding errors obtained when satisfying the boundary conditions accumulate at each time step. The reason for the poor accuracy is not that the diagonal banded portion of the matrix system is not strictly diagonally dominant. This is because all three implementations of the Chebyshev integration method discussed here were equally accurate for the linear boundary value problem eq. (3.7), which has a similar structure to the boundary value problem solved at each time step when the time step is small (see Golub and Van Loan [13] for an analysis of solving banded matrix systems). We conclude that for boundary value problems, there is no time integration, so the effect of rounding errors created by enforcing the boundary conditions does not grow. As explained by Higham [18, p. 79], it is possible to reduce these rounding errors by compensated summation when obtaining multiplicative coefficients for the homogeneous solutions and adding them to the particular solution. Unfortunately, this is not a portable method because it depends on how floating point operations are implemented. Figures 4(c) and 4(d) show how the computational accuracy and computational time vary as the number of modes is increased. Again, it is fastest to use the Chebyshev integration method and solve the bordered banded matrix system rather than to obtain a particular solution to which solutions of the homogeneous boundary value problem are added. When more than 256 Chebyshev modes are used, it is more computationally costly to use the preconditioned Chebyshev method with spectral differentiation in the IMEX time-steppingschemes. A particularly interesting feature in Fig. 4(c) is that, as the number of modes is increased, a smaller timestep is not required to ensure the stability of the IMEX schemes. Greengard [15] has hinted that the Chebyshev integration method may overcome the time step restrictions for explicit time stepping schemes using the real space Chebyshev differentiation method pointed out 21

10

4

2 1

10

2

10 Error
0 −1 −2 1 1 0 x −1 0 t 0.5 u

0

L!([0,1],L2) 10
−2

L!([0,1],W2,2)

10

−4

10 −5 10

−6

10

−4

10 Time step

−3

10

−2

(a) Exact solution to eq. (4.14) given (b) Error vs. timestep with 1024 by eq. (4.15) . Chebyshev modes when approximating eq. (4.14). Fig. 5. Temporal convergence results for the fully implicit scheme eq. (4.16) to the initial boundary value problem with nonlinear dissipation given in eq. (4.14).

by Trefethen and Trummer [35] and Weideman and Trefethen [40] (see Trefethen and Embree [34, p. 287] for a recent review of stability of spectral spatial discretizations of initial boundary value problems). The Chebyshev spectral integration matrices used here are not normal, but their eigenvalues and singular values are clustered near the origin in such a manner that, as the size of the discretization grows, the largest eigenvalues and singular values are bounded independent of the size of the discretization. This would suggest that, provided the resulting linear systems are solved accurately, IMEX time-stepping methods using Chebyshev integration matrices will have stability properties that are independent of the spatial resolution. Finally, Fig. 5(a) shows the exact solution to eq. (4.14) given by eq. (4.15) plotted over the simulation time. Figure 5(b) demonstrates second order convergence of the implicit scheme in eq. (4.16). In these simulations 1024 Chebyshev √ modes were used and we set the physical parameters to ρ = 10000, = 0.001 and β = 0.00001. Although the L∞ ([0, 1], W 2,2 ) error appears to be large, the exact solution has an L∞ ([0, 1], W 2,2 ) norm of 2.5x104 , and so the relative error in the numerical approximation is small. Numerical experiments showed that when the damping parameter, β, is large, very small time steps are needed to ensure that the iterative scheme converges. In most solid mechanics applications however, damping is small and so this limitation is not very restrictive. Numerical tests showed that the IMEX scheme in eq. (4.4) was unstable for reasonable timesteps when applied to eq. (4.14), even when the forcing term was set to zero. In applying the IMEX scheme, the nonlinear part of the dissipation was approximated using 22

the second order extrapolation formula, so that uxxt 2u2 + u4 + 4uxt ux + u3 x x x ≈ 2 uxxt − uxxt 2u2 x 2u2 x + + u4 x u4 x + 4uxt ux + + 4uxt ux + u3 x u3 x
t=jδt t=(j−1)δt t=(j−2)δt

j−2 3uj−1 − 4uxx + uj−3 xx xx ≈2 2u2 + u4 x x 2δt 3uj−1 − 4uj−2 + uj−3 x x x +4 ux + u3 x 2δt

t=(j−1)δt

−

j−2 3uxx

−

4uj−3 xx

+

uj−4 xx

+4

3uj−2 x

2δt − 4uj−3 + uj−4 x x 2δt

2u2 + u4 x x ux + u3 x
t=(j−2)δt

.

5

Conclusion

The Chebyshev integration method allows one to obtain high resolution simulations to semilinear initial boundary value problems on finite nonperiodic intervals with time dependent boundary conditions. If less than 100 modes are required, the standard Chebyshev differentiation method is easiest to use. If high resolution simulations are required, and the function only needs to be approximated correctly in the L2 norm, then the preconditioned Chebyshev differentiation method is the best one to apply. When a wide range of scales are present in smooth solutions to high-order initial boundary value problems and spatial derivatives of the function need to be estimated, then the Chebyshev integration method is the appropriate one to use for numerical investigation of these solutions.

Acknowledgements

The author thanks J.M. Ball, L. Castell, N. Condette, M. Leok, P. Plech´˘, ac C. Ortner, U. Salman, E. S¨li, W.T. Tee, L.N. Trefethen and two anonymous u readers for very helpful comments during the course of this work. In particular, one of the readers is thanked for suggesting examination of the accuracy of solutions obtained using the preconditioned Chebyshev differentiation method with derivatives obtained using spectral differentiation matrices. The support of the MULTIMAT Marie Curie Research Training Network (MRTN-CT-2004505226) is gratefully acknowledged. 23

Appendices

A

Preconditioned Chebyshev differentiation matrices

To obtain the fourth order preconditioned Chebyshev differentiation matrices, we extend Gottlieb and Orszag [14, p. 119] second order preconditioned Chebyshev differentiation matrix approach. The formulas for the fourth order preconditioned Chebyshev differentiation matrices are similar to those for the second order preconditioned Chebyshev differentiation matrices, but, as they are not available elsewhere, we give them here. The presentation is similar to that in Canuto et al. [2, p. 173] and Gottlieb and Orszag [14, p. 119] for second order linear boundary value problems. We consider the linear boundary value problem Au + Buxx + Cuxxxx = g, (A.1)

where A, B and C are constants and g is a smooth but otherwise unspecified function. The associated boundary conditions are u(±1) = α and ux (±1) = β. We expand g, u, uxx and uxxxx in terms of Chebyshev polynomials,
∞

g=
n=0 ∞

gn Tn (x), an Tn (x),
n=0 ∞

(A.2) (A.3) (A.4)

u= uxx =

an Tn (x)
n=0

and
∞

uxxxx =
n=0

an Tn (x).

(A.5)

We equate modal coefficients in eq. (A.1), Aan + Ban + Can = gn . (A.6)

We now consider approximations of u, uxx and uxxxx obtained by finite Chebyshev series with N + 1 terms. Rather than use the spectral differentiation matrices to obtain a relationship between an , an , and an , the preconditioned matrices can be found from a relationship between an , an , and an by using the recursion relation 2nan = cn−1 an−1 − an+1 , (A.7)

24

where cn = 0 for n < 0 or n > N, cn = 2 for n = 0, cn = 1 otherwise. This recursion relation can be found in Canuto et al. [2, p. 87] or Gottlieb and Orszag [14, p. 161]. Using eq. (A.7) repeatedly we find that an = cn−1 an−1 an+1 − 2n 2n (A.8)

an =

cn−1 cn−2 cn−1 cn 1 an−2 − + an + a (A.9) 4n(n − 1) 4n(n − 1) 4n(n + 1) 4n(n + 1) n+2 an = cn−1 cn−2 cn−3 a 8n(n − 1)(n − 2) n3 c2 cn−1 cn−2 cn cn−1 − + 2 n−1 + a 8n(n − 1)(n − 2) 8n (n − 1) 8(n + 1)n2 n−1 cn cn+1 cn−1 + + a + 2 (n − 1) 2 8n 8(n + 1)n 8(n + 2)(n + 1)n n+1 1 a − 8(n + 2)(n + 1)n n+3

(A.10)

an =

cn−1 cn−2 cn−3 cn−4 a 16n(n − 1)(n − 2)(n − 3) n−4 cn−1 c2 cn−1 cn−2 cn−3 n−2 + − 16n(n − 1)(n − 2)(n − 3) 16n(n − 1)2 (n − 2) c2 cn−2 cn cn−1 cn−2 n−1 + + a 2 (n − 1)2 16n 16(n + 1)n2 (n − 1) n−2 c2 cn−1 cn n−1 + + 16n(n − 1)2 (n − 2) 16n2 (n − 1)2 cn cn−1 cn cn−1 + + 2 (n − 1) 16(n + 1)n 16(n + 1)n2 (n − 1) c2 cn+1 cn n + + a 2 n2 16(n + 1) 16(n + 1)2 n2 n cn−1 cn − + 2 (n − 1) 16(n + 1)n 16(n + 1)2 n2 cn+1 cn+2 + + a 16(n + 2)(n + 1)2 n 16(n + 3)(n + 2)(n + 1)n n+2 1 + a (A.11) 16(n + 3)(n + 2)(n + 1)n n+4

25

Equations (A.8)–(A.11) are now used to obtain a sparse matrix system for eq. (A.1). Equation (A.11) is of the form, an = µ1 an−4 + µ2 an−2 + µ3 an + µ4 an+2 + µ5 an+4 (A.12)

where after some simplification and using the fact that we are interested in coefficients for 4 ≤ n ≤ N to eliminate cn , we find that µ1n := cn−4 , 16n(n − 1)(n − 2)(n − 3) 1 , 4(n + 1)n(n − 1)(n − 3) (A.13)

µ2n := −

(A.14)

µ3n :=

3 , 8(n + 2)(n + 1)(n − 1)(n − 2) 1 4(n + 3)(n + 1)n(n − 1)

(A.15)

µ4n := − and µ5n =

(A.16)

1 . 16(n + 3)(n + 2)(n + 1)n

(A.17)

Eliminating an from eq. (A.6) by using eq. (A.11) gives, A (µ1n an−4 + µ2n an−2 + µ3n an + dn+2 µ4n an+2 + dn+4 µ5n an+4 ) + B µ1n an−4 + µ2n an−2 + µ3n an + dn+2 µ4n an+2 + dn+4 µ5n an+4 + Can = (µ1n gn−4 + µ2n gn−2 + µ3n gn + dn+2 µ4n gn+2 + dn+4 µ5n gn+4 ) , which holds for 4 ≤ n ≤ N and for which dn = 0 for n < 0 or n > N, dn = 1 otherwise. We now eliminate an−4 , ... ,an+4 from eq. (A.18). To do so eq. (A.9) is used along with the observation that µ1n an−4 + µ2n an−2 + µ3n an + µ4n an+2 + µ5n an+4 = ν1n an−2 + ν2n an + ν3n an+2 , (A.18)

(A.19)

26

where after some simplification and again using the fact that we are interested in the case 4 ≤ n ≤ N to eliminate cn , we find that ν1n = 1 , 4n(n − 1) (A.20)

ν2n = − and

1 2(n − 1)(n + 1)

(A.21)

ν3n = The final matrix system is

1 . 4n(n + 1)

(A.22)

A (µ1n an−4 + µ2n an−2 + µ3n an + dn+2 µ4n an+2 + dn+4 µ5n an+4 ) + B (ν1n an−2 + ν2n an + dn+2 ν3n an+2 ) (A.23) + Can = (µ1n gn−4 + µ2n gn−2 + µ3n gn + dn+2 µn gn+2 + dn+4 µ5n gn+4 ) , (A.24) which holds for 4 ≤ n ≤ N . This gives a system of N + 1 unknowns and N + 3 equations. The last four equations are obtained from the boundary conditions, u(±1) = α and ux (±1) = β. We use the values of the Chebyshev polynomials and their first derivatives at x = ±1, which are Tn (±1) = (±1)n , dTn (±1) = (±1)n+1 n2 dx and can be found in Hesthaven et al. [17, p. 258]. The equations for the boundary conditions are then
N

α= and
N

(±1)n an
n=0

β=

(±1)n+1 n2 an .
n=0

(A.25)

This set of equations corresponds to the matrix system ˆ ˆ ˆ APˆ + B Pˆ + C Pˆ + LBC1 u = Pˆ g + RBC1, D D D Dˆ
0 2 4 0

(A.26)

in which u is a vector with the coefficients for the Chebyshev series of u, the ˆ ˆ ˆ equations for the boundary conditions are in LBC1 and RBC1 contains the values of these boundary conditions. 27

B

Fourth order Chebyshev integration matrices.

To obtain the Chebyshev integration matrices, we use the following indefinite integral identities (see for example Hesthaven et al. [17, p. 257]):

T0 (x) = T1 (x),

T1 (x) =

T2 (x) 4

and

Tn (x) =

Tn+1 (x) Tn−1 (x) − , 2(n + 1) 2(n − 1)

where Tn (x) is the nth Chebyshev polynomial. Suppose
∞

uxxxx =
n=0

bn Tn (x).

(B.1)

Then by using the indefinite integral identities we find that
∞ b2 bn−1 − bn+1 T1 (x) + Tn (x), 2 2n n=2

uxxx = e3 + b0 −

(B.2)

uxx =e2 + e3 −
∞

+
n=3

b1 b3 b4 b0 b2 + T1 (x) + −+ T2 (x) 8 8 4 6 24 bn−2 bn bn+2 − + Tn (x), 4n(n − 1) 2(n − 1)(n + 1) 4n(n + 1)

(B.3)

ux =e1 + e2 −

b0 b2 b4 e3 b1 3b3 b5 + − T1 (x) + − + − T2 (x) 8 12 48 4 24 64 192 b0 b2 b4 b6 + − + − T3 (x) 24 32 80 480 ∞ bn−3 3bn−1 + − 8(n + 1)n(n − 2) n=4 8n(n − 1)(n − 2) 3bn+1 bn+3 + − Tn (x), (B.4) 8(n + 2)n(n − 1) 8(n + 2)(n + 1)n

28

and e3 b1 3b3 b5 + − + T1 (x) 8 48 128 384 e2 b0 11b2 b4 b6 + − + − + T2 (x) 4 24 384 120 1920 b1 3b3 b5 b7 e3 − + − + T3 (x) + 24 128 320 576 5760 b0 b2 b4 b6 b8 + − + − + T4 (x) 192 240 480 1680 13440 ∞ bn−4 bn−2 + − 4(n + 1)n(n − 1)(n − 3) n=5 16n(n − 1)(n − 2)(n − 3) bn+2 3bn − + 8(n + 2)(n + 1)(n − 1)(n − 2) 4(n + 3)(n + 1)n(n − 1) bn+4 + Tn (x). (B.5) 16(n + 3)(n + 2)(n + 1)n

u =e0 + e1 −

In these equations, ei are constants of integration and bi are the coefficients of the Chebyshev series for uxxxx . To implement the method numerically, this infinite system of equations must be truncated. If there are N + 1 modes for uxxxx , then the linear system above has N + 5 coefficients, one for each Chebyshev mode and 4 coefficients to be satisfied by the boundary conditions. To show how this method is used, and to define the integration matrices, consider the linear boundary value problem

Au + Buxx + Cuxxxx = g,

where A, B and C are constants and g is a smooth but otherwise unspecified function. The boundary conditions are u(±1) = α and ux (±1) = β. After transforming the functions into the space of Chebyshev polynomials, the resulting equation for the integration matrices is ˆ ˆ ˆ ˆ ˆ ˆ AS 4 + B S 2 + C S 0 + LBC2 uxxxx = g + RBC2, ˆ

ˆ in which the spectral integration matrix of order k is denoted by S k . This system of equations is used to define each spectral integration matrix. This 29

system of equations is explicitly given by

Ae0 + Be2 + Cb0 = g0 , b1 3b3 b5 b1 b3 e3 − + + B e3 − + + Cb1 = g1 , A e1 − + 8 48 128 384 8 8 b0 11b2 b4 b6 b4 e2 b0 b2 A − + − + +B −+ + Cb2 = g2 , 4 24 384 120 1920 4 6 24 b1 3b3 b5 b7 b3 b5 b1 e3 − + − + +B − + + Cb3 = g3 , A 24 128 320 576 5760 24 16 48 b2 b4 b6 b8 b4 b6 b0 b2 A − + − + +B − + + Cb2 = g4 , 192 128 480 1680 13440 48 30 80 and for 4 < n ≤ N ,

bn−2 bn−4 − 16n(n − 1)(n − 2)(n − 3) 4(n + 1)n(n − 1)(n − 3) 3bn bn+2 + − 8(n + 2)(n + 1)(n − 1)(n − 2) 4(n + 3)(n + 1)n(n − 1) bn+4 + 16(n + 3)(n + 2)(n + 1)n + Cbn = gn . A

Here bn = 0 for n > N and gn are the Chebyshev series expansion coefficients for g. To obtain the equations that fix the last four coefficients, the boundary conditions, u(±1) = α and ux (±1) = β, are imposed. This gives the following equations e2 e3 7b0 5b1 47b2 9b3 − ± + 4 12 192 384 1920 640 b4 b5 b6 b7 b8 − ± − ± + 160 1152 13440 5760 13440 N bn−2 bn−4 + (−1)n − 16n(n − 1)(n − 2)(n − 3) 4(n + 1)n(n − 1)(n − 3) n=5 3bn bn+2 + − 8(n + 2)(n + 1)(n − 1)(n − 2) 4(n + 3)(n + 1)n(n − 1) bn+4 + 16(n + 3)(n + 2)(n + 1)n

α =e0 ± e1 +

30

and β =e1 ± e2 + b0 b1 5b2 3b3 b4 b5 b6 − ± + − 12 24 96 64 120 192 480 N 3bn−1 bn−3 − + (±1)n 8n(n − 1)(n − 2) 8(n + 1)n(n − 2) n=4 bn+3 3bn − , + 8(n + 2)n(n − 1) 8(n + 2)(n + 1)n e3 4

ˆ where again, bn = 0 for n > N . These equations are in the matrix LBC2 and ˆ vector RBC2. Once the matrix system is solved, we can use the coefficients to find the functions and their derivatives using the truncated versions of eqs. (B.1) to (B.5).

References [1] R. Ahluwalia, T. Lookman, and A. Saxena. Dynamic strain loading of cubic to tetragonal martensites. Acta Materialia, 54:2109–2120, 2006. [2] C. Canuto, M.Y. Hussaini, A. Quarteroni, and T. Zang. Spectral methods fundamental in single domains. Springer, 2006. [3] G.-Q. Chen, Q. Du, and E. Tadmor. Spectral viscosity approximations to multidimensional scalar conservation laws. Math. Comp., 61(204):629– 643, 1993. [4] C.W. Clenshaw. The numerical solution of linear differential equations in Chebyshev series. Proc. Camb. Phil. Soc., 53:134–149, 1957. [5] E.A. Coutsias, T. Hagstrom, and D. Torres. An efficient spectral method for ordinary differential equations with rational function coefficients. Mathematics of Computation, 65(214):611–635, 1996. [6] S.M. Cox and P.C. Matthews. A pseudospectral code for convection with an analytical/numerical implementation of horizontal boundary conditions. Int. J. Numer. Meth. Fluids, 25:151–166, 1997. [7] T. Davis. Suitesparse, www.cise.ufl.edu/research/sparse/. [8] T. Davis. Direct methods for sparse linear systems. SIAM, 2006. [9] R. de la Cruz, B.K. Muite, and H. Servat. Strong linear scaling for spectral simulations of time dependent semilinear partial differential equations on Marenostrum. In 8th World Congress on Computational Mechanics, 2008. [10] P. Dondl and J. Zimmer. Modeling and simulation of martensitic phase transitions with a triple point. J. of Mech. Phys. Solids, 52:2057–2077, 2004. [11] D. Elliot. A method for the numerical integration of the one-dimensional heat equation using Chebyshev series. Proc. Cambridge Philos. Soc., 57:823–832, 1960. [12] D. Funaro and W. Heinrichs. Some results about the pseudospectral 31

[13] [14] [15] [16]

[17] [18] [19]

[20] [21]

[22] [23] [24] [25] [26] [27] [28]

[29] [30]

[31]

approximation of one-dimensional fourth-order problems. Numer. Math., 58:399–418, 1990. G. Golub and C. Van Loan. Matrix Computations. The John Hopkins University Press, 1989. D. Gottlieb and S.A. Orszag. Numerical Analysis of Spectral Methods: Theory and Applications. SIAM, 1977. L. Greengard. Spectral integration and two-point boundary value problems. SIAM J. Numer. Anal., 28(4):1071–1080, 1991. E. Hairer, C. Lubich, and G. Wanner. Geometric Numerical Integration. Springer Series in Computational Mathematics. Springer, second edition, 2006. J.S. Hesthaven, S. Gottlieb, and D. Gottlieb. Spectral methods for timedependent problems. Cambridge University Press, 2007. N. Higham. Accuracy and stability of numerical algorithms. SIAM, 2006. K.-H. Hoffmann and P. Rybka. On convergence of solutions to the equation of viscoelasticity with capillarity. Comm. Partial Differential Equations, 25:1845–1890, 2000. J. Jim´nez and R. Moser. What are we learning from simulating wall e turbulence? Phil. Trans. Roy. Soc. A., 365:715–732, 2007. C. Kane, J.E. Marsden, M. Ortiz, and M. West. Variational integrators and the Newmark algorithm for conservative and dissipative mechanical systems. Int. J. Numer. Meth. Engng., 49:1295–1325, 2000. Y. Kaneda and T. Ishihara. High resolution direct numerical simulation of turbulence. J. of Turbulence, 7(20):1–17, 2006. A.-K. Kassam and L.N. Trefethen. Fourth-order time stepping for stiff pdes. SIAM J. Sci. Comput., 26(4):1214–1233, 2005. W. Kress. Error estimates for deferred correction methods in time. Appl. Numer. Math., 57:335–353, 2007. M. Leok. Foundations of computational geometric mechanics. PhD thesis, Caltech, 2004. A. Lew. Variational time integrators in computational solid mechanics. PhD thesis, Caltech, 2003. A. Lew, J.E. Marsden, M. Ortiz, and M. West. Variational time integrators. Int. J. Numer. Meth. Engng., 60:153–212, 2004. A. Lundbladh, D.S. Hennigson, and A.V. Johansson. An efficient spectral integration method for the solution of the Navier-Stokes equations. Technical Report FFA-TN 1992-28, Aeronautical Research Institute of Sweden, 1992. M. Luskin. On the computation of crystalline microstructure. Acta Numerica, 5:191–258, 1996. N. Mai-Duy. An effective spectral collocation method for the direct solution of high-order odes. Commun. Numer. Meth. Engng., 22:627–642, 2006. N. Mai-Duy and R. Tanner. A spectral collocation method based on integrarted Chebyshev polynomials for two-dimensional biharmonic 32

[32] [33] [34] [35] [36] [37]

[38] [39]

[40]

[41] [42] [43]

boundary-value problems. J. Comp. and Appl. Math., 201:30–47, 2007. J.E. Marsden and M. West. Discrete mechanics and variational integrators. Acta Numerica, 10:357–514, 2001. L.N. Trefethen. Spectral Methods in MATLAB. SIAM, 2000. L.N. Trefethen and M. Embree. Spectra and Pseudospectra. Princeton, 2005. L.N. Trefethen and M.R. Trummer. An instability phenomenon in spectral methods. SIAM J. Numer. Anal., 24(5):1008–1023, 1987. L. Truskinovskii. Equilibrium phase interfaces. Sov. Phys. Dokl., 27(7):551–552, 1982. L. Truskinovsky. Kinks versus shocks. In J.E. Dunn, R. Fosdick, and M. Slemrod, editors, Shock induced transitions and phase structures in general media, volume 52 of IMA Vol. Math. Appl., pages 185–229, 1993. L.S. Tuckerman. Transformations of matrices into banded form. J. Comput. Phys, 84:360–376, 1989. A. Vainchtein. Stick-slip interface motion as a singular limit of the viscosity-capillarity model. Mathematics and Mechanics of Solids, 6(3):323–341, 2001. J.A.C. Weideman and L.N. Trefethen. The eigenvalues of second-order spectral differentiation matrices. SIAM J. Numer. Anal., 25(6):1279– 1298, 1988. C. Xu and T. Tang. Stability analysis of large time-stepping methods for epitaxial growth models. SIAM J. Numer. Anal., 44:1759–1779, 2006. P. Yeung. Lagrangian investigations of turbulence. Annu. Rev. Fluid Mech., 34:115–142, 2002. A. Zebib. A Chebyshev method for the solution of boundary value problems. J. Comput. Phys., 53:443–455, 1984.

33



Add Comment

To add a comment, you need to login or signup (you can signup without leaving this page)