On discontinuous solutions of hyperbolic equations

Computer methods in applied mechanics and engineering

ELSEVIER

Comput. Methods Appl. Mech. Engrg. 116 (1994) 103-112

On discontinuous solutions of hyperbolic equations Knut S. Eckhoff Department of Mathematics, University of Bergen, All~gt. 55, N-5007 Bergen, Norway

Abstract

A modified Fourier method which is applicable to numerical simulationsof discontinuoussolutions of hyperbolicequations is presented. Truncated Fourier series representations are de-truncated at every timestep by utilising step-functions in the reconstruction of the discontinuoussolutions. The de-truncated expansions are used to eliminate the numerical dispersion and diffusion which normally ruin the discontinuoussolutions in a few timesteps when traditional Fourier methods are applied.

I. Introduction

Spectral methods are in many cases known to work extremely well for problems involving smooth functions [1, 2]. For problems involving functions of finite regularity, however, and in particular for problems involving discontinuous functions, there are still a number of issues to be settled. It is for instance well known that the traditional spectral methods have to be modified if satisfactory results are to be obtainable for linear hyperbolic equations with discontinuous initial data [3], as well as for non-linear hyperbolic equations [4]. The approximations obtained by truncated Fourier series expansions for discontinuous functions are highly oscillatory, and the oscillations are particularly bad close to the discontinuities. This oscillatory behaviour near a discontinuity is known as the Gibbs phenomenon, and is regarded as the major cause of the problems connected with application of spectral methods to problems involving discontinuous functions. Most of the recipes suggested in the literature for curing these problems for spectral methods are intended at eliminating, or at least considerably reducing, the oscillatory behaviour of the truncated expansions by various ways of filtering. Even though they often produce acceptable results, most of the suggested filtering techniques do not seem to take advantage of the full power of the spectral methods [5]. At the same cost, they so far only seem at best to be able to produce results of accuracy comparable to that obtained by traditional finite difference methods. The exceptions to this rule are the promising results obtained by utilising step-functions in the reconstruction of discontinuous functions, an approach which was initiated by Gottlieb et al. [6], and which has been further developed in [7-10]. That approach is in accordance with the idea introduced by Lax [11], namely that the oscillatory behaviour obtained by a high-order method such as a spectral method should contain enough information for us to be able to reconstruct the proper non-oscillatory discontinuous function by a post-processing filter. The idea of introducing step-functions was originally only utilised as a cosmetic post-processing operation, intended at displaying the relevant physical results obtained by the calculations at a few preselected values of time. In a recent work [9], however, an approach is attempted where the discontinuous function is reconstructed at every timestep. T h e purpose of the present paper is to advocate a new approach to the numerical solution of initial value problems with discontinuous solutions which is based on some of the same underlying ideas applied in [9], but with radically different 0045-7825/94/$07.00 © 1994 Elsevier Science B.V. All rights reserved SSDI: 0045-7825(93)E0188-E

104

Knut S. Eckhoff / Comput. Methods Appl. Mech. Engrg. 116 (1994) 103-112

ways of thinking at some crucial points. The objective is to take further advantage of the inherent potential of spectral methods, and thus improve upon the results that can be obtained in applications. We illustrate our approach by considering initial value problems with discontinuous initial data for linear hyperbolic equations with smooth coefficients, i.e. problems of the same type as those considered in [3]. We shall furthermore show how the ideas presented can be applied to non-linear problems such as those considered in [9].

2. Linear hyperbolic equations We first consider the application of Fourier methods to the problem of obtaining accurate numerical solutions of well-posed initial value problems with discontinuous initial data for linear hyperbolic systems of partial differential equations of the form

c)u Ou O---f+ A(x, t) ~ + B(x, t)u + c(x, t) = 0 ,

(1)

where u = {u~ . . . . . urn} are the dependent variables, while A, B are given m × m matrices with smooth 2"rr-periodic coefficients and c is a given m-dimensional vector with smooth 2xr-periodic components. For this purpose, we look closely at the following model equation:

Ou Ou O---t-+ a(x, t) ~ + b(x, t)u + c(x, t;v) = 0 ,

(2)

where u = u(x, t) may denote any one of the unknown functions u, . . . . . u,, in (1). The coefficients a(x, t), b(x, t), c(x, t; v) in (2) are regarded in the following as known 2rr-periodic functions. Even though (2) in our context is a single linear hyperbolic equation, the considerations in the following can also be related to cases where the inhomogeneity c(x, t; v) in (2) depends on the other unknown functions determined by the other differential equations in (1) as well. In such cases it is assumed that the proper solution v(x, t) of those differential equations is substituted into c. Thus, (2) may be regarded as one of the equations in the given linear system of partial differential equations (1). The basic underlying assumption in such cases is clearly that the associated initial value problem is well posed; i.e. that a uniquely determined solution is known to exist and therefore in principle can be substituted into c. If we consider initial value problems with discontinuous initial data, the function c will then normally be a linear combination of discontinuous functions and the first derivative of discontinuous functions, while the coefficients, a, b in the following will be assumed to be smooth functions. In this work, we consider discontinuous 2~r-periodic solutions of the differential equations (1), (2) which are piecewise smooth on [0, 2"rr]. According to the classical theory [12], differentiation of Fourier series associated with discontinuous functions cannot be accomplished by termwise differentiation. Discontinuous solutions of the differential equations (1), (2) are not solutions in the classical sense, however, they are generalized solutions and should therefore be considered as generalized functions. In the context of generalized functions [13], differentiation of Fourier series can be accomplished by termwise differentiation in formally the same way as for smooth functions. Therefore in the following we consider .the occurring functions as generalized functions whenever necessary. In our handling of the partial differential equation (2), we represent the solution as well as the coefficients in (2) at each instant t by truncated Fourier series expansions. The emphasis is on the maintenance of accurate approximations for the occurring Fourier coefficients, as the time integration goes on. In order to achieve the stated goal, we start by considering the complete Fourier series for the functions involved. Thus, we assume that

u(x, t) ~ ~

•k(t) e 'k" .

Substitution of (3) into (2) formally leads, for k = 0, - 1 , +-2 . . . . .

(3) to

Knut S. Eckhoff / Comput. Methods Appl. Mech. Engrg. 116 (1994) 103-112

105

dtl k

dt (t) -I-~'k(t ) -t- ~k(t) -t- C.k(t)

= O,

(4)

where we have set Ou

+~

a(x, t)-~x (x, t) ~ ~, fk(t) e i*x ,

(5 /

+oe

b(x, tlu(x, t ) ~ E c(x, t;v) ~ ~

~,(t) e i*x ,

(61

e,(t) e ikx

(7)

k~-~

With the emphasis on the maintenance of accurate approximations for the Fourier coefficients that occur, we see from ( 4 ) - ( 7 ) that we have to look carefully at the way the Fourier coefficients associated with the product of two functions are calculated in the spectral method we are applying. In fact, the two middle terms (5) and (6) in (4) are such terms, and it is of vital importance that the calculation of those terms is in accordance with the stated goal. In the cases where (2) is one component of the system of Eq. (1), the function c is also normally obtained as a linear combination of terms of the types (5), (6). Let us therefore consider two arbitrarily given functions +~ t/(X)

E

+oo

a^ ke ikx ,

v(x)-

k= - ~

~

6ke ikx ,

(8)

k= - ~

and form the product W(X) = a(x)o(x)

~

E

E

a m l ) n e i("l+'O~ ~

171= - o0 n ~ - ~

E

I b k e ikx

(9)

k=_~

Formally, we then have, for k = 0, - 1 , +-2 . . . . . -t-oo

if. = E

a,,,v. = E

/;.i + / i = k

+~

a.,Vk-,,, = d~o6, + E [a.,Vk-., + a-.,6k+,,,l.

/I'1 = - o ¢

(10)

~1=1

which constitute the traditional convolution sums [1]. On the assumption that one of the functions involved, a(x) or v(x), is a smooth function, the traditional convolution sums (10) are known to be valid even if the other function is a generalized function [13]. Let us proceed with our reasoning by considering the case where the time t is fixed and where a(x) is one of the coefficients a, b in the linear hyperbolic equation (2) and v(x) is either the unknown function u or its derivative Ou/ax. Since the coefficients a, b in (2) are assumed to be smooth functions, the solution of an initial value problem for (2) will depend continuously on those coefficients a, b [14], and truncated Fourier series expansions are known to provide good approximations for a, b. Provided that the even integer M is chosen large enough, the solutions of (2) will therefore change only slightly if we approximate a(x) by its truncated Fourier series expansion PMa(X) defined by M/2-1

PMa(x) =

~

?~keikX.

(11)

k=-M/2+l

In our context it is important to note that this conclusion is valid both for smooth solutions and for solutions of finite regularity of Eq. (2). On this background, we conclude that for smooth functions a(x) it is advisable in (10) to introduce the approximation M/2-1

ffk = ~06k + ~

[a,,6k-,. + ~-,,6k+?,]"

(12)

m=1

When we apply a numerical procedure, we are clearly forced to truncate any expansion which enters into the constructed algorithm. Thus, in some way we have to truncate the expansion for the unknown

106

K n u t S. E c k h o f f

/ C o m p u t . M e t h o d s A p p l . M e c h . Engrg. 116 (1994) 1 0 3 - 1 1 2

function (3) also, and the most readily apparent way of doing this clearly is to assume that for some finite even integer N, we have N/2-1

u(x, t) ~

t) =

PNll(X,

~

hk(t ) e ik~

(13)

k=-N/2+[

With this assumption, we see from (4) that we have to calculate accurate approximations for fk(t) and +-(N/2- 1) if we are to be able to advance the Fourier coefficients occurring in (13) accurately in time. Those calculations must clearly be consistent with the convolution sums (12) in order to be relevant in our context. In the exceptional case where M = 2 , which corresponds to the case where a(x) is no longer dependent on x, i.e. a(x) =- c~0, we have that rbk = hov k for every k involved. Hence, we can conclude that for the cases where the coefficients a(x, t) and b(x, t) in (2) are independent of x, it seems perfectly possible to maintain accurate approximations for the occurring Fourier coefficients by just adopting an accurate time integration scheme for (4). Numerical experiments confirm this expectation beyond doubt. On the other hand, when (2) has coefficients a(x, t) and b(x, t) which do depend on x, we are faced with the disturbing fact that although all occurring coefficients dj in (12) may be known, we cannot a priori assume that the coefficients 6; are known for 1/I ~> N/2 since those coefficients are directly related to the corresponding coefficients for the unknown function (13). Thus, at least some of the higher order coefficients fk(t), ~,k(t) in (4) will depend on Fourier coefficients which are not readily available. In particular, for the case M = N, we may write (12) in the following way for k = - 1 , ---2. . . . . +-(N/2 1):

~k(t) for k = 0 , -+1, +2 . . . . .

N/2-lk[- I

~k =

~

N/2-1

a,,,6k-,,, +

m= - N I 2 + l k l + l

~

[a,,,Ok-,,, + a-,,,Ok .... ] ,

(14)

m=N/2-lk]

while we have N/2- 1

~,, =

Z

h,,,6_,,,.

(15)

m = -N/2+ 1

By the traditional Galerkin method, the expression (15) is left unchanged, while in the last sum in (14), only one of the two terms appearing in the brackets is retained since the other is unavailable. The error made by utilising the Galerkin approximation for the convolution sums (14) is substantially much larger for functions of finite regularity than for smooth functions. In fact. in the extreme case where v(x) is the derivative of a discontinuous function, we have

6 k = O(1)

as k---~-+~.

(16)

The two terms appearing in the brackets in the last sum in (14) will therefore be of the same order in that case, and the error made by the traditional Galerkin method can be expected to be intolerable at least for the highest order modes. This expectation is confirmed beyond doubt by numerical experiments when v(x) is the derivative of a discontinuous function, and also when v(x) is a discontinuous function, both for the Galerkin approximation and for the straight forward collocation approximation [16]. If our goal of maintaining accurate approximations for the occurring Fourier coefficients shall be reached, there is on this background, at least for discontinuous functions, a need to do better than the Galerkin method and the straightforward collocation method in approximating (14). For solving (2) in the case with smooth coefficients, the most apparent alternative seems to be to try to calculate as accurate approximations for the lacking Fourier coefficients in the convolution sums (14) as possible. For the solution of (2) this means that at each time step, we would have to try to determine an accurate approximation for P2N_2tt(X, t) solely from knowledge of Puu(x, t). Such a de-truncation of the truncated Fourier series (13) is not as impossible as it may appear at first glance. In fact, as we shall see in the following section, a suitable de-truncation can be obtained as a direct result of the reconstruction of discontinuous functions described in [15].

Knut S. Eckhoff / Comput. Methods Appl. Mech. Engrg. 116 (1994) 103-112

107

3. The method of de-truncation

As in the preceding section, we assume that a(x) is a smooth function for which the truncated Fourier series expansion PNa(x) is known. We assume furthermore that u(x)'is a discontinuous function for which PNu(x) is known, and that we are interested in obtaining an accurate approximation for the truncated Fourier series expansion PN{a(x)u(x)} for the product a(x)u(x) of those two functions. In [15], it is shown how knowledge of PNu(x) can be used to approximately reconstruct the function u(x) in the following way: M

u(x) = v(x) + ~ , A y ( x ; 3't),

(17)

j=l

where v(x) is a continuous 2"rr-periodic function which is approximately given by an Nth order truncated Fourier series expansion, A t, 3'j are some constants for j = 1 . . . . . M, and V ( x ; / 3 ) is the 2-rr-periodic extension of the following function:

I

-~7 ( / 3 - ~ " - x ) ,

V(x;/3) =

if - , r < x < / 3 , (18)

t.-~-~-w(/3+~r

x),

if/3
where/3 is a parameter at our disposal such that - r r ~
1 (/'k(/3)==2-~ll-i~e -ik~ ,

k = ±1, ±2 . . . . .

(19)

Hence it follows from (17) that M

def

A t e -ik*, = 2-trik(t~, - 6,) = Ck ,

k = -+1, ±2 . . . . .

(20)

j=l

Since it follows from our assumptions that 6 k = O(Ik1-2) and tJk = O ( [ k l - ' ) as k---->---~, approximate equations for the constants A t, 3'j in (17) follows from (20). In particular, if the locations of the discontinuities 3,t, j = 1 . . . . . M are known with an accuracy of the order O(N-2), Eq. (20) with t3k neglected constitutes for k = N / 2 - M , N / 2 - M + 1 . . . . . N / 2 - 1 a linear algebraic system of equations which determine the associated jumps A t with an accuracy at least of the order O(N -1) as N---* oo. For linear systems of hyperbolic equations, this completes the reconstruction, since the discontinuity locations can be determined from the initial values by the method of characteristics [14]. When the discontinuity locations are not known, however, as will be the case for instance for non-linear problems, those locations can be determined with sufficient accuracy from (20) by the method provided in [15]. In fact, if we seek an algebraic equation of degree M Z

M

-~-xIzM-I~-x2zM-2~-''''~-XM_IZJI-XM=O.

with roots z~ = e -~j, j = 1 . . . . . M, the coefficients X t, j = 1 . . . . . inhomogeneous linear algebraic equation: C k + C k _ t X 1 + Ck_2X 2 + ." "Ck_MX M = 0 ,

(21) M in (21) will satisfy the following

(22)

for every integer k ~ 0, 1,2 . . . . . M. This is used in [15] to determine accurate approximations for the discontinuity locations Yr" The reconstructed function (17) gives us an opportunity to de-truncate the known expansion PNu(x) to any desired order, since by (17), (19), we obtain the approximate Fourier coefficients

Knut S. Eckhoff / Comput. Methods Appl. Mech. Engrg. 116 (1994) I03-I12

108

Uk -- 2"rrik

j=l

A t e -ik~q ,

k =

+-N/2, +-(N/2 + 1) . . . . .

(23)

where the error is of the order O([k[-2) as k---~ -+oo. In particular, an approximation for P21v_2u(x) can thus be easily determined, and the dispersion-free convolution sums (14) can therefore be utilised both for (5) and (6). As we know from the discussion in [15], the accuracy of the reconstruction (17), and hence of the de-truncation (23), is determined by the smoothness of the exact function v(x) in (17). If that function v(x) is infinitely smooth, the de-truncation will by the above procedure be determined with spectral accuracy, while we will have finite order accuracy if v(x) is continuously differentiable a finite number of times only. In the latter case the de-truncation can be made more accurate by the method described in the following section. From a computational point of view, it is important to note that the de-truncated convolution sums (14) can be calculated efficiently by the use of FFT [16]. If the procedure is applied at each timestep in the advancing of the solution of (2), it should with this background be possible to obtain discontinuous solutions with substantially improved accuracy in the cases where the coefficients in (2) are smooth functions. Numerical experiments fully confirm this expectation, and have also shown that the described approach is fairly robust in practical applications. Further details on this approach will be given elsewhere.

4. Improved methods of de-truncation For cases where the exact continuous function v(x) in (17) is of finite regularity, the method of de-truncation described in the preceding section can be improved. For that purpose we now let V0(x;/3) = V(x;/3) denote the 2~r-periodic step-function defined by (18). Based on that function we may then successively define a sequence of functions V,,(x;/3) by

V.(x;/3)=fV,,_,(x;/3)dx,

n=l.2

.....

(24)

where the constants of integration are successively determined by

f r Vn(x;/3)dx=O,

n=l,2

.....

(25)

From this definition, it follows that for each n = 1,2 . . . . , V~(x; fl) is a 2~r-periodic function of finite regularity with "'(P)" v . t x ; f l ) continuous everywhere for p = 0 , . . ,n-l, but with v"(")',tx;/3) only piecewise continuous with jump-discontinuities of magnitude + 1 at x =/3 + 2m~r, m = 0, - 1, -+2 . . . . . For the higher order derivatives V.tP) (x;/3), -p ~>n + 1 there are no jumps at those singularity locations. More specifically, it follows from (18) and (24) that V.(x;/3) is a piecewise polynomial function for each n = 0 , 1,2 . . . . . i.e. Vn(x; ~) is a polynomial of degree n + 1 on each interval (/3 +2m~r, /3 + 2(m + 1)~r), m = 0, -+1, -+2 . . . . . In particular, it follows from (18), (24) and (25) that V.(x;/3) = U.(x-/3), where U.(~) is the 2~r-periodic extension of the following function for n = 0, 1 . . . . : U.(~:) =

(n + 1)! Bn÷l

'

when 0 < s¢ < 2"rr '

(26)

where Bj(x), j = 1, 2 . . . . are the Bernoulli polynomials [17]. We note that U,(~:) is an odd function when n is an even number and U,(~:) is even when n is odd. If we now consider a 2~r-periodic function u(x) which is known to be piecewise smooth on [-'rr, ~r], the interval [-~r, ~] can be divided into a finite number of subintervals on which u(x) is smooth. At the endpoints of those subintervals, however, the function u(x) a n d / o r some (or all) of its derivatives may

Knut S. Eckhoff / Comput. Methods Appl. Mech. Engrg. 116 (1994) 103-112

109

have jump-discontinuities. Thus, the assumption that the 2~r-periodic function of finite regularity u(x) is piecewise smooth on [-'rr, "rr], is clearly equivalent with the assumption that u(x) for any given integer Q can be written in the following way for some finite integer M: Q

M

u(x) = v(x) + ~ ~ ATV.(x; ?t),

(27)

n=0j=l

where o(x) is some Q times continuously differentiable 2~r-periodic function which is piecewise smooth on [-~r, xr], the functions V.(x; Yt) are given by (18), (24), (25) with /3 = 3't, and A~, 3't are some constants for j = 1 . . . . . M, n = 0, 1 , . . . , Q. From the properties of the functions V.(x;/3) discussed above, it is clear that the points x = Yt, J = 1 . . . . . M, are the locations for the singularities of the function u(x) given by (27), and A ° are the associated jumps function u(x) itself, A~ are the associated • 2 jumps for u'(x), A t are the associated jumps for u"(x), and so on until, finally, A are the associated jumps for u(°)(x) at those singularity locations. Clearly, there will be no loss of generality by assuming that - ~ ~ ~'1 < 3'2 < " " " < ~'M < ~. The sum on the right hand side in (27) is then seen to constitute a 2~r-periodic piecewise polynomial function of degree at most Q + 1, i.e. it is a polynomial of degree at most Q + 1 on each of the intervals ('/M-2~r, Yl), (Y1,'/2), (~%,Y3). . . . , (~M-~,'/M), (~'M, Y1+2~) • The problem of determining good approximations for the function v(x) and the constants A~', ~j in (27) solely from knowledge of a finite n u m b e r of Fourier coefficients associated with the function u(x) can be solved in an analogous way as for the case Q = 0 considered in the preceding section. In fact, from (18), (19), (24) and (25) we see that, for n = 0, 1, 2 . . . . , we have A

(

A

)0(/3)= 0,

(

v. )k(/3)

e

-ik~

2rr(ik)-+l ,

k = +-1, -+2 . . . . .

(28)

Hence, it follows from (27) that o 1 M .=0 2"rr(ik) "+1 jl= At e-ikvj = /~k

k = 0, -+1, +-2 . . . . .

(29)

Since it follows from our assumptions that t3k = O(]kl - ° - 2 ) and t~k = O((Ik l-1) as k--~ +_oo, approximate equations for the constants A~,3,~ in (27) follows from (29). In particular, if the locations of the discontinuities Yr J = 1 . . . . . M are known with an accuracy of the order O ( N - ° - 2 ) , Eq. (29) with Vk neglected constitutes for k = N/2 - M(Q + 1), N / 2 - M(Q + 1) + 1 . . . . . N/2 - 1 a linear algebraic system of equations which determine the associated jumps A~ with an accuracy at least of the order O(N -0-1+') as N---> oo for each j = 1 , . . . , M and n = 0, 1 . . . . . Q. For linear systems of hyperbolic equations this completes the reconstruction, since the discontinuity locations can be determined from the initial values by the method of characteristics [14]. When the discontinuity locations are not known, those locations can be determined with sufficient accuracy from (29) by an improvement of the method provided in [15]• More details on this improved reconstruction will be given elsewhere [18]• An approximated reconstruction of the function u(x) on the form (27) can on this background be calculated from PNu(x). The reconstructed function u(x) given by (27) gives us the opportunity to de-truncate the known expansion PNu(x) to any desired order, since by (28), (27), we obtain the approximate Fourier coefficients Q

M

t~k = n=0 ~ 2"tr(ik) 1 "+1 j=l ~ Aj,, e -ikvj ,

k =

+-NI2, -+(NI2 + 1) . . . . .

(30)

With this improved de-truncated expansion, we may calculate the convolutions considered in the preceding sections with improved accuracy. The resulting solutions of finite regularity of the considered linear hyperbolic equations with variable coefficients will consequently appear with improved accuracy•

Knut S. Eckhoff / Comput. Methods Appl. Mech. Engrg. 116 (1994) 103-112

110

5. Non-linear hyperbolic equations We conclude this paper by giving some remarks on how it may be possible to make the application of spectral methods more feasible for discontinuous solutions of non-linear problems in the form of conservation laws au a a~- + ~x f(U,X, t) = O , (31) where u = {u~ . . . . . u,,,) are the dependent variables and where f ( u , x , t ) is a given m-dimensional vector function which is assumed to be smooth with respect to all its arguments and 2rr-periodic with respect to the variable x. We restrict ourselves to discussing 2rr-periodic solutions u(x, t) of (31) which for each t > 0 are piecewise smooth on the interval [-rr,-rr]; we therefore need to design a procedure which enables us to calculate accurate truncated Fourier series approximations for the occurring non-linear terms for such functions u(x, t). The procedure for the calculation of such terms will depend on the type of term we are dealing with, and will in general be more complicated than the procedure we are going to describe in the following for a special type of non-linear term. The procedure we are going to describe can be fairly easily modified, however, to cover most of the types of non-linear terms occurring in applications in view of the discussion in the preceding sections. Further discussion on that problem will be given elsewhere. We limit our discussion here to the case where we have two scalar functions v(x), w(x) which are known to be of finite regularity and for which only the truncated Fourier series PNV(X) and PNw(x) are known. The non-linear term we consider in the following is simply the product of those two functions

v(x)w(x). In order to make the application of spectral methods feasible in this connection, we clearly have to design a procedure which enables us to calculate an accurate truncated Fourier series approximation PN {O(X)W(X)} solely from knowledge of the truncated Fourier series PNO(X) and PNW(X). The procedure of de-truncation described in the preceding sections may also play a decisive role here. In fact, the truncated Fourier series expansion for v(x) can be de-truncated, and an accurate approximation

v(x) = V(x) + PP(x)

(32)

can be established, where V(x) is given by an Nth order truncated Fourier series expansion and PP(x) is a piecewise polynomial function. Similarly, the truncated Fourier series expansion for w(x) can be de-truncated, and an accurate approximation

w(x) = W(x) + PQ(x) ,

(33)

can be established, where W(x) is given by an Nth order truncated Fourier series expansion and PQ(x) is a piecewise polynomial function. Both the piecewise polynomial functions PP(x) and PQ(x) are obtained by the method described in the preceding sections as linear combinations of the special family of piecewise polynomial functions V,,(x; ~). If we now consider the product of the two functions given by (32) and (33), we obtain an approximation on the form

v(x)w(x) : V(x)W(x) + V(x)PQ(x) + W(x)PP(x) + PP(x)PQ(x) .

(34)

In (34) we can regard V(x) and W(x) as smooth functions which can be identified with their Nth order truncated Fourier series expansions. Hence the Nth order truncated Fourier series expansion for the first term on the right hand side in (34) can be calculated by the usual Fourier Galerkin or Fourier collocation method. The second and third term on the right hand side in (34) are terms of the same type as those considered in the preceding sections in connection with solutions of finite regularity of linear equations with variable coefficients. Hence, the N-th order truncated Fourier series expansions for the second and third term on the right hand side in (34) can be calculated by the procedure described earlier in this work. In this connection, it should be added that in addition to eliminating aliasing, it will also be more economical if we combine the calculation of the first term on the right hand side in (34) with the calculation of the second or the third term.

Knut S. Eckhoff / Comput. Methods Appl. Mech. Engrg. 116 (1994) 103-112

111

It remains to design a p r o c e d u r e for the calculation of the N t h o r d e r truncated Fourier series expansion for the last term on the right hand side in (34). In this connection, it suffices to r e m e m b e r that since both functions PP(x) and PQ(x) are piecewise polynomial, so also is their p r o d u c t PP(x)PQ(x), which therefore can be easily expressed as a linear c o m b i n a t i o n of the special family of piecewise polynomial functions V,(x; [3). F r o m the knowledge of such a linear c o m b i n a t i o n , and in view of (28), it is a straightforward calculation to obtain any desired truncated Fourier series expansion for the function which constitutes the last term on the right hand side in (34). We have thus established a p r o c e d u r e by which we are able at each timestep to handle at least special types of non-linear terms in (31) with spectral accuracy in principle, in the advancing of solutions of finite regularity. In [16] that procedure has been shown to work with impressive accuracy for the following non-linear scalar conservation law: -~-/t(x, t) + ~

= 0,

(35)

w h e r e u = u(x, t) d e n o t e s the u n k n o w n function to be determined. Eq. (35) is often referred to as the inviscid B u r g e r s ' equation, and seems to be accepted as the universal model p r o b l e m for scalar c o n s e r v a t i o n laws [4].

Acknowledgment This p a p e r is partly based on work d o n e while the author was engaged at the S I N T E F Multiphase Flow L a b o r a t o r y , T r o n d h e i m , N o r w a y . The research was s u p p o r t e d by N T N F C o n t r a c t No. OT45.23626 u n d e r the P r o g r a m for multiphase flow technology ( P R O F F ) 1989-91. T h e support from T h e R o y a l N o r w e g i a n Council for Scientific and Industrial Research ( N T N F ) , A m o c o , Elf, Esso, N o r s k H y d r o , Saga, Shell, Statoil, Total, A k e r Engineering, F r a m o Engineering, N o r w e g i a n A p p l i e d T e c h n o l o g y , K v a e r n e r Engineering, and S I N T E F is thankfully acknowledged.

References [1] C. Canuto, M.Y. Hussaini, A. Quarteroni and T.A. Zang, Spectral Methods in Fluid Dynamics (Springer, New York, 1988). [2[ D. Gottlieb and S.A. Orszag, Numerical Analysis of Spectral Methods: Theory and Applications (SIAM, Philadelphia, PA, 1977). [3] A. Majda, J. McDonough and S. Osher, The Fourier method for nonsmooth initial data, Math. Comp. 32 (1978) 1041-1081. [4] E. Tadmor, Convergence of spectral methods for nonlinear conservation laws, SIAM J. Numer. Anal. 26 (1989) 30-44. [5] L.S. Mulholland and D.M. Sloan, The effect of filtering on the pseudospectral solution of evolutionary partial differential equations, J. Comput. Phys. 96 (1991) 369-390. [6] D. Gottlieb, L. Lustman and S.A. Orszag, Spectral calculations of one-dimensional inviscid compressible flows, SlAM J. Sci. Statist. Comput. 2 (1981) 296-310. [7] S. Abarbanel and D. Gottlieb, Information content in spectral calculations, in: E.M. Murman and S.S. Abarbanel, eds., Progress in Scientific Computing, Vol. 6, Proc. U.S.-Israel Workshop, 1984 (Birkhfiuser Boston Inc., 1985) 345-356. [8] S. Abarbanel, D. Gottlieb and E. Tadmor, Spectral methods for discontinuous problems, NASA-CR-177974, ICASE Report no. 85-38 (1985); also in: K.W. Morton and M.J. Baines, eds., Numerical Methods for Fluid dynamics II, Proc. Conf., Reading, 1985 (Clarendon Press, Oxford, 1986) 129-153. [9] W. Cai, D. Gottlieb and C.-W. Shu, Essentially nonoscillatory spectral Fourier methods for shock wave calculations, Math. Comp, 52 (1989) 389-410. [10] D. Gottlieb, Spectral methods for compressible flow problems, in: Soubbaramayer and J.P. Boujot, eds., Proc. 9. Internat. Conf. Numer. Methods Fluid Dynamics, Saclay, France, 1984, Lecture Notes in Physics 218 (Springer, Berlin, 1985) 48-61. [111 P.D. Lax, Accuracy and resolution in the computation of solutions of linear and nonlinear equations, in: C. de Boor and G.H. Golub, eds., Recent Advances in Numerical Analysis, Proc. Symposium Univ. of Wisconsin-Madison (Academic Press, New York, 1978) 1(}7-117. [12] A. Zygmund, Trigonometric Series, Vol. 1 (Cambridge Univ. Press, Cambridge, 1968). [13] D.C. Champeney, A Handbook of FOURIER THEOREMS (Cambridge Univ. Press, Cambridge, 1987). [14] R. Courant and D. Hilbert, Methods of Mathematical Physics, Vol. 11: Partial Differential Equations (Wiley, New York, 1962).

112

Knut S. Eckhoff / Comput. Methods Appl. Mech. Engrg. 116 (1994) 103-112

[15] K.S. Eckhoff, Accurate and efficient reconstruction of discontinuous functions from truncated series expansions. Math. Comp. 61 (1993) 745-763. [16] K.S. Eckhoff, On the handling of linear and nonlinear terms in spectral methods, SINTEF report no. STFll F91053 (1991). [17] A. Erd61yi, W. Magnus, F. Oberhettinger and F.C. Tricomi, Higher Transcendental Functions (McGraw-Hill, New York, 1953). [18] K.S. Eckhoff, Accurate reconstructions of functions of finite regularity from truncated Fourier series expansions, Math. Comp. (1994) submitted.