NASA Conmctor Report 177974 

lCASE REPORT NO. 85-38 

leASE 

NASA-CR-177974 

19860004484 

SPECTRAL METHODS FOR DISCONTINUOUS PROBLEMS 

Saul Abarbanel 

David Gottlieb 

Eitan Tadmor 

NASA Contract No. NAS1-17070 

August 1985 

r I , ..... , '''11 r ~ '11 : )' '\ r ~ rf ,. ~ ,. \" , , I I ~ l. • t I .. I. .. 

• ,-J·r.I~1 t'Ui 

~:O\I }' II 'ORr- f< 'f I,,;~/) 

INSTITUTE FOR CO~WUTER APPLICATIONS IN SCIENCE AND ENGINEERING 

NASA Langley Research Center, Hampton, Virginia 23665 

Operated by the Universities Space Research Association 

NJ\SI\ 

National Aeronautics and 

Space Administration 

Langley Research Center 

Hampton Virginia 23665 

111111111111111111111111111111111111111111111 

NFOI005 
1 Report No NASA CR-177974 I 2 Government Accession No 3 ReciPient's Catalog No 

ICASE Report No. 85-38 

4 Title and Subtitle 5 Report Date 

SPECTRAL METHODS FOR DISCONTINUOUS PROBLEMS 1001:; 

6 -I'irlormlllg Organization Code 

7 Author(sl 8 Performing Organization Report No 

Saul Abarbanel, David Gottlieb, and Eitan Tadmor 85-38 

10 Work Unit No 

9 Performing Organization Name and Address 

Institute for Computer Applications in Science 

and Engineering 11 Contract or Grant No 

Mail Stop 132C, NASA Langley Research Center NAS1-17070 

H::Imnt"nn VA ?~n6'l 13 Type of Report and Period Covered 

12 Sponsoring Agency Name and Address -

National Aeronautics and Space Administration - Washington, D.C. 20546 ln~tc:I~L.~genCY""t'odI ~ 

505-31-83-01 

15 Supplementary Notes Submitted to Conference on 

Langley Technical Monitor: Numerical Methods for Fluid 

J. C. South Jr. Dynamics, Reading 

Final Reoort 

16 Abstract 

We show that spectral methods yield high-order accuracy even when applied 

to problems with discontinuities, though not in the sense of pointwise 

accuracy. Two different procedures are presented which recover pointwise 

accurate approximations from the spectral calculations. 

~ 

--

, 

. 

17 Key Words (Suggested by Author(sll 18 Distribution Statement 

spectral methods, discontinuities, 64 - Numerical Analysis 

shock recovery 

Unclassified - Unlimited 

19 Secunty aa55lf (of thiS reportl 20 Security aaSSlf (of thiS pagel 21 No of Pages 22 Price 

Unclassified Unclassified 31 A03 

N-305 For sale by the National Technical Information Service, Springfield, Virginia 22161 NASA-Langley, 1985 
SPECTRAL METHODS FOR DISCONTINUOUS PROBLEMS 

Saul Abarbanel 

David Gottlieb 

Eitan Tadmor* 

Tel-Aviv University, Tel-Aviv, Israel 

and 

Institute for Computer Applications in Science and Engineering 

Abstract 

We show that spectral methods yield high-order accuracy even when applied 

to problems with discontinuities, though not in the sense of pointwise 

accuracy. Two different procedures are presented which recover pointwise 

accurate approximations from the spectral calculations. 

Research was supported in part by the National Aeronautics and Space 

Administration under NASA Contract No. NASl-17070 while the authors were in 

residence at lCASE, NASA Langley Research Center, Hampton, VA 23665. 

Additional support was also provided in part by the Air Force Office of 

Scientific Research under Contract No. AFOSR 83-0089 for the first and second 

authors, and by the National Science Foundation under Grant No. DMS85-03294 

and Army Research Office under Grant No. DAAG29-85-K-0190 for the third 

author. 

* Bat-Sheva Foundation Fellow 

i 
1. INTRODUCTION 

Consider the evolution partial d1fferentiation equation ut = Lu, on a 

finite interval, where L is a hyperbolic operator. The solution u has a 

projection PN u on a finite subspace (which may for example consist of the 

first N modes in a Galerkin method, or 

interval), and a numerical approximation 

N collocating points in the 

generated by some spectral 

method. For linear operators it is known from the Lax equivalence theorem 

that if the scheme is consistent and stable, then uN ~pproximates PN u in 

some appropriate norm. If u is smooth, then the theorem implies that uN 

approximates the solution u in the same sense. 

In practice, one looks at the point values of uN at the grid points and 

takes it as an approximation to the values of the true solution u at these 

points. We shall call this approach the realization of the computed solution 

via its grid-points value. The aims of the paper are: 1) demonstrate that 

when u is a complicated function, this realization will not produce 

acceptable results; 2) to suggest different ways for the realization of the 

solution in such cases. 

The following examples give a very clear illustrat10n of the misleading 

results that may be obtained by pointwise realization. 

Example 1 

Consider the equation 

o " x " 2'1r 

(1) 
-2-

where u(x) and uO(x) are periodic functions and ~(x) is a discontinuous 

function. If we expand uO(x) in Fourier series we get 

CD 

E k=-CD 

where 

a = ~ f2 n -ikx k 2n 0 uO(x)e • 

The solution u(x) is thus given by 

CD 

u(x) = E ak e

ikt 

k=-CD 

ikx 

e 

(2a) 

(2b) 

Suppose that (1) is solved numerically by the Fourier-Galerkin method, namely 

we seek a trigonometric polynomial of the form 

that satisfies 

(

auN _ aUN eikX) = 

a t ax' 0 , 

N 

E ak k=-N 

ikx 

e 

From (3) it is clear that 

and 

-N ... k ... N 

-N ... k ... N 

(3) 

(4) 
-3-

yielding the solution 

Therefore 

ikt 

e 

ikx 

e (5) 

Equat10n (5) implies that uN(x,t), obtained from the numerical solution (3), 

coincides with PN u(x, t), the Galerkin projection of u, thus yielding the 

best possible convergence of uN to PN u. However, since the Fourier series 

of u(x,t) converges very slowly, the point values will not 

approximate well u(xJ,t). In general, one would witness the Gibbs phenomenon 

of overshoot in the neighborhood of the discontinuity and global oscillations 

allover the domain. In fact, even the init1al approximation, uN(x,O), 

displays the same behavior in relation to uO(x). 

In the second example we show that the same phenomenon occurs even if the 

numerical initial point values do approximate the true initial point values to 

a high degree of accuracy. 

Example 2 

Consider the equation (1) where uO(x) is the saw-tooth function 

uO(x,x) = {Ax 

A(2x - n) 

x < x 

(6) 

x > x 

for some k, 0 < k < 2N-l, x = ~ (k+ liz) •• 
-4-

In the pseudospectral Fourier method we seek a trigonometr1c polynomial 

(7) 

such that 

aVN aVN --=-- at ax 

at the points j 0, ••• ,2N-l (8a) 

(8b) 

Since vN is a polynomial of degree N, (8a) implies that 

aVN aVN --=--

at ax (9) 

for all x Moreover, from (8b) it is clear that vN(x,O) is the 

(unique) trigonometric polynomial of order N that interpolates uO(x) at the 

points Xj' j = 0, ••• ,2N-l, thus 

(10) 

where 2N-l 

= 2N~ L 

R. j=O 

(11 ) 

Performing (11) we get 

AN [k - N + .5] (12) 
where 

A_1f_ 

2NcR, 

-5-

-id (k+l) N 

2 1 - e ... n 1 + ictn -"_I. _ ------1~·1f~R,~---- N 

-N1 - e 

~ = 2, It I * N. 

The numerical solution vN(x,t) of (9), (10) is 

and upon manipulating (12), (13) one gets 

t * a (13) 

(14) 

(15) 

The trigonometric interpolant FN(X,x) collocates uO(x,x) at the grid 

points Xj. However, in between the grid points it oscillates. If we read 

the values of vN(x,t) at the grid points, then by (14) 

and unless t = ~ for some integer N m, we will get solution that looks 

oscillatory. Thus, even though the initial approximation looks smooth at the 

grid points, when it evolves in time the oscillations will present themselves 

at the points Xj. 
-6-

The conclusion one might draw from the above examples is that spectral 

methods (or any h~gher-order methods) are useless when applied to 

discontinuous function. A different approach is to look at a different 

realization of the numerical solution rather than the pointwise one. We will 

argue that high-order accurate information is contained in the numerical 

solution and demonstrate how that information can be extracted in such a way 

that accurate pointwise approximation to the true solution can be obta~ned. 

2. INFORMATION AND HOW TO EXTRACI IT 

Consider the linear equation 

u = Lu t 

u(O) = uo 

(16) 

where L is a linear hyperbolic operator with variable coefficients and uO 

is a discontinuous function. For simplicity, we will restrict ourselves to a 

periodic, one (space) dimensional problem though the results are more general, 

(see Gottlieb and Tadmor [2]). Let v be the solution of the auxiliary 

problem 

v = - L * v 

t 

(17) 

v(a) = vo' 

where va is a coo function. B f ecause 0 th e h b yper 0 Ii c it yo, f L (17) i s a 

well-posed problem. In Lemma 1 we quote the well-known Green's identity. 
-7-

Le1llD.a 1: Let u( t) and v( t) be the solutions of (16) and (17) at some 

level t, then 

(18) 

Assume now that (16) and (17) are discretized by the Fourier-Galerkin 

method. That is, we seek UN and vN that are trigonometric polynomials of 

degree N such that for every k, Ikl ( N 

(auN ikx) ~ - L UN' e = 0 (19a) 

(uN(O) - uO' e

ikX) = 0, (19b) 

(ikX aVN * ) e , (at + L VN) = 0 (19c) 

e

ikx , (vN(O) - Vo) = O. (19d) 

We have also a Green identity for UN and vN. 

L~ 2: 

(20) 

Proof: Since vN(t) and uN(t) are Nth-order trigonometric polynomials we 

use (19a) and (19c) to get 
-8-

and therefore 

which implies (20). 

We will proceed by showing the relation of the RHS of (20) to that of 

(18) • 

~3: 

(21) 

where 

(22) 

for every s. 

Proof: From (19b) it is clear that 

(23) 

Also, 
-9-

and since Vo is a C~ function, 

for every s. (24) 

Now 

and in view of (23) and (24), 

where 

and this proves the Lemma. 

From Lemmas 1 - 3 we can conclude: 

Theorem 1: Let u(t) and v(t) be the solutions of (16) and (17), 

be the solutions of the FourierGalerkin approximations of (16) and (17). Then 

for every s. (25) 

The proof is an immediate consequence of (18), (20), and (21). 
-10-

Assume now that the Fourier-Galerkin method described in (19c) and (19d) 

is stable, then vN(t) approximates vet) within spectral accuracy, that is 

nvo 

UvN(t) - v(t)n = e:2 < K __ s 

Ns- 1 

We can, therefore, replace vN(t) in (25) and get 

(uN(t),v(t») = (u(t),v(t») + e: 

where £ is spectrally small. We use now the fact that every C~ function 

vet) can be obtained from some Vo in (17). This is, in fact, one of the 

definitions of hyperbolicity. We can, therefore, state: 

Theore. 2: Let u(t) be the (nonsmooth) solution of (16) and let uN(t) be 

the solution of the spectral Galerkin approximation to (16). Then for any C~ 

function vet) 

(26) 

where e: is spectrally small. 

Thus, uN(t) approximates weakly u(t) within spectral accuracy. It is 

in this sense that contains a highly accurate information about 

u(t). We will show later how to use this information in order to obtain 

spectral accurate approximation to the grid-point values of u(t). 
-11-

We turn now to the pseudospectral Fourier case. Here we need some 

preprocessing of the initial data in order to prove the same result as in 

Theorem 2. 

Theore1l 3: Let uN(x, t) be a trigonometric polynomial of order N that 

satisfies 

at x = x J' 

_ 1TJ 

Xj - N ' j '"' O, ••• ,2N-l 

(27) 

Ikl .. N, 

is the solution of the pseudospectral Fourier scheme, but 

initially uN(x,O) is obtained by the Galerkin projection). 

Then for every smooth function u(x,t) 

r 

2N-l =I21T W ~ ~(Xj ,t) vex. ,t) u(x,t) v(x,t)dx + e: 

J=O J 0 

(28) 

where e: is spectrally small, provided that the pseudospectral approximation 

is stable. 

Proof: Let vN be the solution of the pseudospectral Fourier approximation 

of (l7a) and let vN(O) be the Galerkin projection of vo' that is 

Ikl .. N. (29) 

From (27) and the analog equation for vN, one gets 
-12-

(30) 

From the exactness of the trapezoidal rule for polynomials of degree 2N, we 

conclude 

(31) 

Note that the initial functions are not the 

interpolants of Uo and vo as in the usual pseudospectral methods but rather 

the Galerkin approximation to and We recall now Lemma 3 and 

equality (18) to establish (28). The proof is thus completed. 

It is interesting to note the way in which the information is contained. 

The interpolant of Uo looks smooth at the grid points, whereas the Galerkin 

approximation of Uo looks oscillatory on the grid points. It means that in 

order to preserve the information one has to require initially osc11latorylooking solution. The information is preserved 1n the structure of the 

oscillations. 

We will show now a way of using (26) and (28) in order to construct a 

better approximation to then the one given by 

uN(x,t) is given by either the Galerkin method or the pseudospectral method). 

From (28) and (26) it is clear that in order to get a good approximation 

to u(y,t) at some point y (0,21T), we need to find a function vy(x,t) 

such that 

21T j[ u(x,t) v (x,t)dx = u(y,t) + £1' 

o y 
-13-

where £1 is spectrally small. By (26) we will have 

2'IT j[ ~(x,t) v (x,t)dx = u(y,t) + £ + £1 

o y 

(32) 

for the Galerkin method and 

2N-1 

; E ~(x. ,t) v (x. ,t) = u(y,t) + £ + £1 

j=O J y J 

(33) 

for the pseudospectral method. 

For conveniency we will shift the interval [O,2'IT] to [-'IT,'IT]. Let 

p(x) be a CD C -function vanishing outside the interval [-'IT ,'IT] satisfying 

p(O) = 1. (34) 

Let Dp(X) be the Dirichlet kernel, namely 

D (y) =_1 

P 2 'IT 

~ eikx 1 sin(p+ liz) y L.J = ZiT sin (y/2) (35) 

Ikl<p 

We set now 

(36) 

One can prove (see [2 ]) that 

J

'IT e u(x)~ ,P(y-x)dx = u(y) + £2 (37) 

where £2 is spectrally small. 

Thus, it is possible to extract accurate pointwise values from uN(x). 
-14-

3. NUMERICAL RESULTS 

In this section we demonstrate the efficacy of the smoothing procedure 

outlined above. As a test function we have chosen the piecewise em-function 

{ 

. x S1n 2 

f(x) = 

-sin i 

(38) 

'II' ( X ( 2'11' • 

.. Denote its spectral approximation by fN(x), and let be the 

pseudospectral approximation to f(x). It is evident from the first column of 

.. Tables I and III that fN(yv) - the spectral approximation sampled at 

y = v'll'/N - V 

do not approximate f(y) within spectral accuracy. v 

In fact, 

.. .. the error committed by f I28 (yv) 1S only half of that committed by f 64 (yv). 

Regarding the pseudospectral approximation, fN(x), it, of course, collocates 

the exact values at the sampling grid points, tN(yv) = f(yv); yet, in between 

approximate with1n 

first-order accuracy only, as shown in the first column of Tables II and IV. 

In order to construct our regularization kernel in (36), we define the 

cut-off function p(~) = Pa(~) to be 

a~2 

~2 _ 1 I ~ I < 1 

o otherwise 

namely, P (~) is a e

m_ function whose support is the interval I~I < 1. 

a 

~ to be used in (36) is of the form 

,I,a ,p(y) = _1_ P (a-1 y) sin (p+ liz )y/a ~ 2'11'a a sin y/2a 

(39) 

(40) 
-15-

,. 

The post-processing procedure of the spectral approximation fN involves 

convoluting fN against ~9,P, namely 

f(x) ~ 1 [21T f

N

(y)p(X9

Y) sin (p+ liz )(x-y)/9 dy 

21T9 0 sin (x-y) /29 (41) 

where x 1S a fixed point of interest. (In practice, we use the trapezoidal 

rule to evaluate the right-hand-side of (41) taking a large number of 

quadrature points.) 

The parameter 9 was chosen as 

9 (42) 

this guarantees that ~ is so locahzed that it does not interact with 

regions of discontinuity. 

It should be noted, in this stage, that 1f 9 was so chosen to be the 

same for each x, (and not as in (42», the formula (41) admits a simpler 

form; that is, if 

then 

co 

~ iky LJ O'k e 

k=-co 

N ,. 

f(x) ~ 1: f(k)O'k 

k=-N 

ikx 

e 

This procedure can be carried out efficiently in the Fourier space. 

(43) 

(44) 

Next, we turn to the post-processing for the pseudospectral approximation 

,. 

fN(x) which is s1mpler than (41). In fact, in this case 
-16-

2 2N-1 e 

f(x) ..... 2N'Ir E f(y )111 ,p(x-y ). 

\1=0 \I \I 

(45) 

Note that carrying out the smoothing procedure defined in (45) does not 

involve any extra evaluation of f(y) in points other than y , in contrast 

\I 

to spectral smoothing procedure in (41). As before, the parameter e was 

chosen according to (42). We have yet to determine the parameters p and 

a. The parameter p must be equal to Na for 0 < a < 1, in order to 

assure infinite accuracy. (In our computations, a ... 8.) Finally, we feel 

that a is problem dependent and we chose a = 10. We have not tuned the 

parameters to get optimal results; further tuning may improve the quality of 

our filtering procedure. 

In Tables I, II, III, and IV we give the results of the smoothing 

procedure at several points in the domain. The pointwise values are now 

recovered with high accuracy. The first column in each table indicates the 

points in which the procedure was performed. We limited ourselves to four 

points in the interval (O,'Ir) because of the symmetry of the function f(x). 

A 

The second column gives either the spectral approximacion fN(x) or the 

pseudospectral approximation 1N(x), N = 128 in Table I and II and N = 64 in 

Tables III and IV. The third column gives the smoothed results, when filtered 

by (41) on (45), at the same points as in column I. 

The results indicate the dramatic improvement obtained by the smoothing 

procedure. Moreover, note that the error committed by 1128 is 

better than the one committed by 164 (or £64) only by a factor of 2 whereas 

after the post-processing the error improves by a factor of 104• 
-17-

Table I. Results of smoothing of the spectral 

approximation of f(x), N = 128 

'lrV 

X =-

v 8 

v equals 

2 

3 

4 

5 

I f (x ) - fN(x ) I 

v v 

3.2 (-3) 

5.2 (-3) 

7.8 (-3) 

1.1 (-2) 

Table II. Same as Table I for the pseudospectral approximation 1N(x). 

x 1 =!. (v+ liz ) v+ /2 8 

v equals 

2 

3 

4 

5 

5 (-3) 

8.1 (-3) 

1.2 (-2) 

1.8 (-2) 

A * 

If - fN 1/11 

at x - x 

v 

5.8 (-10) 

7.9 (-10) 

6.3 (-10) 

1.1 (-10) 

* If - 1N 1/11 

at x = xv+ liz 

7 (-10) 

7.9 (-10) 

6.4 (-10) 

1.2 (-10) 
-18-

Table III. Results of smoothing of the spectral 

approximation of f(x), N = 64 

lrV 

x =- If(x) - fN(x )1 v 8 v v 

v equals 

2 6.4 (-3) 

3 1 (-2) 

4 1.5 (-2) 

5 2.3 (-2) 

Table IV. Same as Table III for the pseudospectral approximation, 1N(x). 

xv+ 1/2 = i (v+ 1/2 ) 

v equals 

2 

3 

4 

5 

1 (-2) 

1.6 (-2) 

2.4 (-2) 

3.6 (-2) 

,. 

If - f * ljIl N 

at x = x 

v 

4.8 (-6) 

5.9 (-6) 

7.7 (-6) 

8.9 (-6) 

4.1 (-6) 

6 (-6) 

7.8 (-6) 

8.9 (-6) 
-19-

4. A DIFFERENT METHOD FOR EXTRACTING INFOBHATION 

In this section we would like to present a different approach for 

extracting the information from an oscillatory solution. The idea is to 

subtract from the solution those oscillations that correspond to the saw-tooth 

function discussed in Example 2. This leads to the following procedure: 

Let iR.x 

e , be the solution of the pseudospectral 

approximation to a hyperbolic problem. We try to find an unknown smooth 

function and a (oscillatory) saw-tooth function FN(x-t,xs ) with an unknown 

jump 2nA at an unknown location Xs such that 

r

N- 1 

H = L 

j=O 

uN(x. ,t) - AFN(X.,x ) - c - J J s 

(46) 

is minimized. Note that we have 2p + 3 unknowns in (46): A, xs ' c and 

2P values of b

t (l * 0). 

The conditions for local minima of H are found from the following 

2p + 3 equations: 

2N-l p iR.x. aH o ==) L 2 - cF. - F. L. bt -= e J = 0 (47) ujFj - AFj aA J J 

j=O t=-p 

vO 

where Fj .. FN(xj ,xs ), uj = ~(Xj ,t). Also, 
-20-

2N-1 P u.x. aH o ==) ~ -AF ~ b! e 

J -= u. - c - 0 

ac J J (48) 

J=O !=-p 

U:O 

2N-1 P ix.! aH o ==) ~ F: - AF: F. - cF: - F: ~ b! e J 

as 

u. 

J J J J J J 

o (49) 

j=O !=-p 

vO 

N L aaR.(s) u.x. 

where F: = aFN(x.,x )/as = • e J. and 

J J s as , 

!=-N 

aH A 

~b = 0 ==) b = u - Aa , a m m m Iml 1,2, ••• ,p (50) 

m 

where 

Substituting (50) into (47), (48), and (49) we get, respectively: 

uo - AaO - c = 0 (51) 

(52) 

(53) 

where a'(s) = aaR. (x)/as. Next, we combine (52) and (53) to get a single 
-21-

nonlinear equation for s: 

o (54) 

where all sums run over p < It I < N. 

Equation (54) is solved iteratively for s. Having found s, one 

immediately obtains from Example 2 all the at(s)'s. Then from (50) we have 

the bm's, and A from (52). Finally, having A we find c from (51). 

The minimum thus obtained may be a local one while we are seeking a global 

minimum. This means that in practice one searches for the global minimum. 

We now give an example that illustrates the efficacy of the procedure. We 

solve the following problem: 

{

sin ~ 

uN(x,O) = 2 

-sin ~ 2 

o < x < 2~, t > 0 (55) 

(56) 

(57) 

We ran the problem on several grids and exhibit here the numer~cal results for 

the case N = 8 (i.e., 16 subintervals in the domain (0,2~». The 

unadulterated results at t = ~/2N on the grid points are shown in Figure 1. 
-22-

.~ 

-.5 

-1.0 

Figure 1 

x 

o 

exact solut~on 

unsmoothed pseudospectral solut~on 

smoothed solut~on 

(N = 8) 
-23-

Table V 

j exact solution error 1 = error 2 = error 1 

lexact-unsmoothedl I exact-smoothed I error 2 

0 9.80 x 10-2 5.86 x 10-5 5.86 x 10-5 1.00 

1 9.80 x 10-2 1.24 x 10-2 5.86 x 10-5 211 

2 2.90 x 10-1 2.57 x 10-2 6.30 x 10-5 408 

3 4.71 x 10-1 4.13 x 10-2 7.33 x 10-5 563 

4 6.34 x 10-1 6.15 x 10-2 9.30 x 10-5 661 

5 7.73 x 10-1 9.11 x 10-2 1.31 x 10-4 695 

6 8.82 x 10-1 1.43 x 10-1 2.16 x 10-4 662 

7 9.57 x 10-1 2.70 x 10-1 4.42 x 10-4 611 

8 -9.95 x 10-1 1.00 x 100 1.10 x 10-2 91 

9 -9.95 x 10-1 2.68 x 10-1 1.34 x 10-3 200 

10 -9.57 x 10-1 1.42 x 10-1 4.42 x 10-4 321 

11 -8.82 x 10-1 9.07 x 10-2 2.16 x 10-4 420 

12 -7.73 x 10-1 6.12 x 10-2 1.32 x 10-4 464 

13 -6.34 x 10-1 4.11 x 10-2 9.30 x 10-5 442 

14 -4.71 x 10-1 2.55 x 10-2 7.32 x 10-5 348 

15 -2.90 x 10-1 1.22 x 10-2 6.30 x 10-5 194 

We then post-processed these ~(xj,n/2N) values according to the procedure 

described above. The filtered values are shown on the same graph, and the 

errors listed in Table V are computed before and after processing. The 

dramatic improvement is evident. 

Next we demonstrate the procedure in the case of the Euler equation for 

gas dynamics. Because the physical problem involves inflow, outflow, and noflow boundary condl.tions, periodicity could not be imposed and we use the 

Tchebyshev, rather than Fourier, pseudospectral method. 
-24-

The physical problem is that of a wedge, inserted as a zero angle of 

attack, into a uniform supersonic flow of an ideal gas with y = 1.4. An 

oblique shock develops in time and the flow reaches, after a while, a steady 

state. The time-dependent Euler equations in two-space dimensions were 

discretized by the pseudospectral Tchebyshev method in space with an 8x8 grid 

and a modified Euler scheme was used for the time discretization. Since we 

are interested in the steady state only, the accuracy for the time integration 

is of little importance. In order to be sure that a steady state is reached, 

the code was run until all physical quantities did not change to 11 

significant figures over a span of 100 time steps. The values of the density 

in the steady state at the grid points together with the grid points 

themselves are given in Table VI. 

Table VI. 

p y 

1.862 1.851 1.869 1. 871 1.837 1.865 1.892 1.885 1.878 1 • 

1.862 1.870 1.867 1.820 1.870 1. 954 1.899 1.803 1.759 .961 

1.862 1.854 1.852 1.904 1.877 1.770 1.782 1.864 1.900 .853 

1.862 1.871 1. 876 1.812 1.838 1. 969 1. 975 1.884 1.841 .691 

1.862 1.848 1.842 1.935 1.899 1.703 1. 710 1.890 1.984 .5 

1.862 1.883 1. 894 1.729 1.832 2.429 2.994 3.255 3.316 .308 

1.862 1.808 1.810 2.387 3.133 3.375 3.224 3.054 3.002 .146 

1.862 2.115 2.868 3.288 3.176 2.965 3.006 3.136 3.187 .038 

1.862 3.083 3.046 2.975 3.087 3.108 3.024 3.013 3.016 0 

x 0 .038 .146 .308 .5 .691 • 853 .961 1 • 
-25-

Note that the raw data in Table VI seems to indicate roughly the same y-shock 

locat10n at Xo 1, xl = .961, and x2 = .853, namely between the grid points 

Y4 = .3086 and Y5 = .500. This means that because of the coarse Tchebyshev 

grid the shock location cannot be resolved to better than 20% of the domain. 

In fact, the correct shock locations at those x-stations are y = .434 for 

xo, y = .417 for xl and y = .370 for x2' 

In the present case it is not necessry to employ a saw-tooth piecewise 

smooth function, as was done in the previous section, because there is no need 

to preserve periodicity. Instead, we subtract from the oscillatory data an 

expansion of the Heaviside function, S(Y'Ys): 

S(Y'Ys) (58) 

where d1 , the state ahead of the shock, and d2, the magnitude of the 

discontinuity, are constant. The description here of 

independent of x, has to do with the fact that the two-dimensional results of 

the pseudospectral algorithm were post-processed at fixed x-stations. The 

expansion of S(Y'Ys) is given by 

where Tt(Y) is the Tchebyshev polynomial of order t, 

Tt(y) = cos[t cos1(y)], and 

(59) 
-26-

1 = (s + Z)/N 

[

'Il't 1) 'Il't At (s) = sin N (s + -r) /N sin '2N 1 ( t ( N-1 

1 ~(s) = sin[(s + -r»)/2N. 

If s is an integer, then on the grid points, y. 

J cos('Il'J/N). 

S(y.,y ). J s (60) 

The L2-norm which we wish to minimize is now, at any given x-station: 

N 1 Q<N 

H = :E - [PN(Y.) - d1 - d2 SN(Y"y ) -:E b n T n (yJ)]2 (61) 

j=O c j J J s t=1 ~ ~ 

1 ( j ( N-1 

(62) 

j = 0, N 

Differentiating (61) with respect to the parameters d1' d2' sand 

b

t (1 ( t ( P < N), using the orthogonality relations for the Tchebyshev 

polynomials and manipulations similar to those used in the previous section, 

we get p + 3 nonlinear algebraic equations which are completely analogous to 

(50) - (53). They are: 

1,2, ••• ,p. (63) 

o (64) 
where 

-27-

N N 

E A2 = 0 ct At Pt - d2 t~1 cR. 

t=p+l R. 

N N 

E cR. A' Pt - d E c t At A' 

R.=p+1 R. 2 R.=p+1 t 

N 

P., = N

2 "c1

. p(yj)T.,(yJ

.) ,., c

t ~ ,., 

j=O J 

A' =.L A (s). as t 

(65) 

= 0 (66) 

(67) 

(68) 

Again, we combine (65) and (66) into a single nonlinear equation for the shock 

location index, s: 

(69) 

where all the sums are from R. = p+1 to R. = N. 

The procedure for extracting the shock location, jump magnitude and smooth 

part of the solution from the raw data p(x,y j) (given in Table VI) is 

exactly the same as described above for the Fourier problem. 

For the wedge-flow problem considered here, this procedure applied in the 

case of a coarse net (N = 8), located the shock with an error only in the 

fourth significant figure. The smooth part was recovered to within 1% at the 

worst field point. 
-28-

Conclusion 

We have demonstrated that the realization of a numerical solut10n via its 

gr1d-point value may be misleading when the true solution has a complicated 

structure which is not resolved by the grid. We have shown, however, that the 

numerical solution does contain highly accurate information about the solution 

and we suggested two ways of extracting this information. 

The analysis outlined in this chapter is a linear one (though the 

procedure was applied also to nonlinear problems). However, in [28] Lax has 

argued that more information about the solution is contained in high 

resolution schemes even in the nonlinear case. In fact, using notions from 

information theory, Lax has shown that the €-capacity of the set of 

approximate solutions is closer to the €-capacity of the set of the 

projections of exact solutions if the numerical scheme is a high-order scheme. 

In the area of d1gital filters one always processes the data in order to 

overcome the Gibbs phenomenon. If we look at the initial conditions as an 

input s1gnal and at the final result as the output signal, the idea of 

filtering is a natural one. 
29 

REFERENCES 

[1] GOTTLIEB, D., M. Y. HUSSAINI, and S. A. ORSZAG, "Introduction to the 

Proc. of the Symposium on Spectral Methods," SIAM CBMS Series, 1983. 

[2] GOTTLIEB, D. and E. TADMOR, "Recovering pointwise values of 

discontinuous data within spectral accuracy," ICASE Report No. 85-3, 

NASA CR-172535, 1985. 

[3] LAX, P. D., "Accuracy and resolution in the computat1on of solutions of 

linear and nonlinear equations," Recent Advances in Numerical Analysis, 

Proc. Symp. Mathematical Research Center, University of Wisconsin, 

Academic Press, 1978, pp. 107-117. 

[4] MAJDA, A. and S. OSHER, "The Fourier method for nonsmooth initia data," 

Math. Comp., Vol. 32, 1978, pp. 1041-1081. 

[5] MOCK, M. S. and LAX, P. D., "The computation of discontinuous solutions 

of linear hyperbolic equations," Comm. Pure Appl. Math., Vol. 31, 1978, 

pp. 423-430. 
End of Document