Convergence to steady state of solutions of Burgers' equation

Applied Numerical North-Holland

Mathematics

161

2 (1986) 161-179

CONVERGENCE TO STEADY STATE OF SOLUTIONS OF BURGERS’ EQUAnON

*

Gunilla KREISS Royal Institute of Technology, Stockholm, Sweden

Heinz-Otto KREISS California Institute of Technology, Pasadena, CA, U.S.A.

Consider the initial-boundary value problem for Burgers’ equation. It is shown that its solutions time, to a unique steady state. The speed of the convergence depends on the boundary conditions exponentially slow. Methods to speed up the rate of convergence are also discussed.

converge, in and can be

1. Introduction In many gasdynamical problems one tries to calculate the steady state solution by solving the corresponding time dependent problem. One hopes that for t -+ 00 the solution converges to a unique steady state. Recently, Salas, Abarbanel and Gottlieb [l] considered the initial-boundary value problem 4 + t(u*),

t>o,

=f(x),

O
4x9 0) = g(x).

(I -1)

They used f(x)

= sin x cos x,

g(x)=b

sin x,

O
and showed that the solution u( x, t) of the above problem converges to a steady state u(x), as t + co, but that u(x) depends on the initial data. In this paper we consider the viscous problem u,,+ $(u2), = %,

+f(x),

t>,o,

O
C>O,

(1.2a)

with initial and boundary conditions 4%

0) = g(x),

and the corresponding

U(0, t) = a,

(1.2b)

steady state problem

:(y2)x=Eyx,+f(X), y(0) = a,

U(1, t) = b,

OGX
f>O,

(1.3)

y(l) = b.

* Research was partially supported by the Office of Naval Research under NOOO14-83-K-0422 and National Science Foundation Grant DMS-8312264. Additional support was provided by the National Aeronautics and Space Administration under NASA Contract No. NASl-17070 while the authors were in residence at the Institute for Computer Applications in Science and Engineering, NASA Langley Research Center, Hampton, VA 236655225. 0168-9274/86/$3.50

0 1986, Elsevier Science Publishers

B.V. (North-Holland)

G. Kreiss. H.-O. Krew

162

/ Convergence to stea& smte

For simplicity we restrict ourselves to two cases: (I) a>O3b, a>, -b, f(x)=O, (2) c( = b = 0, f is such that there exists an (Ywith 0 < cx < 1 such that f(x) > 0 for 0 < x < N, f(x) -c0 for (Y-Cx < 1, f(0) = f(1) = 0, f,(O) afo > 0 and f,(l) >fo. We will show that (1.3) has a unique solution and discuss the properties of _Y(x). We shall also show that in all cases we consider, the limit of v(x) as 6 + 0 exists. Thus. if lim u(x, t) =.v(x) 1+x exists, we obtain a unique steady state solution of the inviscid equation (1 .l) if we first let t + x and then E + 0. This is in contrast to the procedure in [l], where the two limit procedures are taken in the reverse order. We shall prove that the eigenvalues of the eigenvalue problem A+=

-(Y+),++ru~

(1.4)

G(O)=+(l)=O,

are all negative. Therefore, the solution of (1.2) converges to the solution of (1.3) provided u( x, 0) = g(x) is sufficiently close to -v(x). In another paper we shall prove that u( x, f) converges to y(x) as t + cc for arbitrary initial data. The speed of convergence is determined by the eigenvalues, A,, of (1.4).We shall show that the eigenvalue distribution depends on f(x) and on a, h in the following way: There is a constant c > 0 which does not depend on 6 such that (1)

if a>

-b,

f-OthenO>

(2)

if a = -b,

(3)

if u = b = 0,

(4)

if a=b=O,

f= 0

then -A,

'f(x)

J0 s

-c/E>~,

"'f(x)

>h2>

= O(e-‘I’)

dx # 0, th en -c>h, dx = 0, th en -A,

...;

> 0, -C./C > A, > A, > . . . : >A,> =O(ee”‘)>O.

...;

(1.5)

-c>hz>X,>

....

We expect a reasonable speed of convergence in the first and third case, while in the second and fourth case the speed should be extremely slow due to the eigenvalue -A, = O(e”‘). This is confirmed by numerical experiments. We see that at first u(x, t) quite rapidly approaches the same limit as the inviscid equation (1.1) which consists of solutions of the stationary equation G2),

=f(x)

connected by a shock. Once the viscous shock has been formed, the solution of (1.2) becomes quasi-stationary and the shock creeps extremely slowly to the ‘right’ position. We can explain the behavior, because by linearizing around the quasistationary solution we find that the eigenvalues of the corresponding eigenvalue problem have a similar distribution as earlier. If -A, = O(e-‘I’) then the speed of convergence is so slow that the above method to calculate the steady state is impractical, see Figs. 1 and 3. However, we can use the same technique as Hafez, Parlette and Salas in [2] to speed up the convergence. See Figs. 2 and 4. Unfortunately, not only the speed of convergence but also the condition number of the stationary problem deteriorates. We have to calculate with O(e”‘) decimals to obtain correct results. To avoid an excessive number of decimals we have used a quite large E in our numerical calculations. The situation becomes much better in a two-dimensional case, which we discuss in the last section. Now there is a whole sequence of eigenvalues -j_L,,=o(j%),

.j=l,

2 ,....

G. Kreiss, H.-O. Kreiss / Convergence to steady state

163

-. 5

1 Fig. 1. Convergence in time without convergence acceleration. Numerical solutions at different time stages for the case r=0.05, f=O, a=l, b=-1, u(x,0)=1+2(e-2x - l)/(l - em2). Between each curve there are 200 time steps = 40 time units. The calculation is made with time step k = 0.2 and N = 50 grid points

close to zero. However, they are only algebraically and not exponentially close to zero. We indicate how to modify the procedure to accellerate the speed of convergence. We believe that the viscous model (1.2) better explains what happens in actual calculations than the inviscid equation (1.1). Practically all numerical methods have some viscosity built in. Also, from a physical point of view, the solution we are interested in is the limit of solutions of a viscous equation. Finally we want to point out that the appearance of small eigenvalues has also been observed by Brown, Kath, Kreiss and Henshaw, Naughton (private communication).

2. Uniqueness, existence and properties of the steady state solution We start with uniqueness,

which can be proven

by standard

Lemma 2.1. If the steady equation (1.3) has a solution, Proof. Let u, u be two solutions. t(pw)X

= fWXX>

then it is unique.

Then w = u - u is the solution

p = 2.4+ u,

techniques.

w(0) = w(1) = 0.

of (2.1)

If w # 0 then the zeros of w are isolated. Let X with 0 < X G 1 be the first zero to the right of x = 0. Without restriction we can assume that w > 0 for 0 < x < JE, i.e. w,(O) 2 0 and w,(X) < 0.

164

G. Kreiss, H.-O. IQeiss / Convergence to steady state

1.0

.5

0

-. 5

1.0 0

.5

.25

.75

1.0

Fig. 2. Convergence in time with convergence acceleration. Numerical solutions at different time stages for the same case in Fig. 1. Between each curve there are 15 time steps and one extrapolation step. The same time step, k = 0.2, and number of grid points, N = 50, are used.

Integration

of (2.1) gives us -c(~w,(x)~+~w,(o)~)=~[w,]~=~[pw]~=o.

Thus w,(O) = w,(X) = 0. We can consider (2.1) as an initial value problem w(0) = w,(O) = 0 whose solution is w(x) = 0, and the lemma is proved. 0 We shall now discuss the properties of the solution. a > 0 > 13, a > -b. Integrating (1.3) gives us ey,=+y2-c,

O
y(0) = a.

Let us start

with initial

with the case f(x)

data

= 0,

(2.2)

The constant c has to be determined so that y(1) = b. We necessarily have c = 2d2 > ia2, because with c G :a2, y, > 0 for all x, and y(1) = b cannot be satisfied. We can solve equation (2.2) explicitly. This is done by writing (2.2) in the form Y(X)

Ja

2E i.e.

dy Xdl, jj2- d2 = / ,,

G. Kreiss, H.-O. Kreiss / Convergence to steady state

165

1.0

.5

C

-. 5

1. [II 0

I

I

I

.25

.50

.75

J 1.0

Fig. 3. Convergence in time without convergence acceleration. Numerical solutions at different time stages for the case c = 0.04, f = $T sin(lrx) cos(nx), a = b = 0, u( x, 0) = isin( Between each curve there are 100 time steps. The calculation is made with time step k = 0.1 and N = 50 grid points.

Therefore Y(1) = b implies d = a + O(e-‘/‘), l-re Y(X)

=

a

-a(1

and

-x)/c

1 + 7 e-4’-X)/’

’

with 7 = s.

(2.3)

Away from the boundary layer at x = 1 we have Y(x) = a + O(e-a(‘--X)/f). Thus, for c + 0, Y(x) converges to a for 0 Q x < 1. If a = -b we consider (2.2) on the interval 0 < x ,< $, with boundary conditions Y(0) = a, Y(i) = 0 and obtain a solution Y,(x) of the form (2.3). The solution on the whole interval is given by Ylbh

Y(X) = i

-y,(l-x),

ifO
+,

if$
In Figs. 5 and 6 we have plotted y(x) for two different sets of boundary values. Consider case 2, where f only vanishes at x = 0, ar, 1 and a = b = 0. Without restrictions we can assume that /0

f

0x

20.

(2.4)

166

G. Kreiss, H.-O. Kreiss / Conoergence to steu& state

I

1.01 0

.25

I

I

i

.50

.75

1.0

Fig. 4. Convergence in time with convergence acceleration. Numerical solutions at different time stages for the same case as in Fig. 3. Between each curve there are 20 time steps and one extrapolation step. The same time step, k = 0.1. and number of grid points, N = 50, are used.

If this is not true, we transform %=1-x,

f”=

The new problem

by introducing

new variables,

jj= -y.

-f,

satisfies (2.4).

Lemma 2.2. Lety(x) Y,(l)

the problem

be the solution of (1.3)

cYX(O)

Proof. Integration

G K,>

4

= O~;~l . .

F(x) { lb,(x)

= J(f(C;)

d.$ and h(x)

I> + (h,(O)

= {2F(x).

Then

I.

of (1.3) gives

e(y, -r,(O))

= iY’-

where y,(O) is determined U, =y,(O)

F,

by y(l) = 0. If u = y - h, then u is the solution

- h, + c-‘uh

(2.5)

Y(O) = 0, of

+ +c-‘u*,

u(0) = 0.

Assume that y,(O) > K,. It follows that y,(O) - h,(x) is positive and thus u and U, are positive for all x > 0. In particular u(1) > 0 and y(l) = u(l) + h(l) > 0, which contradicts y(l) = 0. Thus y, (0) < K, . Also %(l)

= v,(O)

This proves the lemma.

- fo) q

G 40).

G. Kreiss, H.-O. Kreiss / Convergence to steady state

167

.5-

I .5

0 Fig. 5. The solution

1.0

of (1.2) when f = 0, a = 1, b = 0 and c = 0.05.

Lemma 2.3. Let y(x) be the solution of (1.3) and let E be sufficiently small. If F(1) > 0 then y(x) > 0 for 0 < x < 1 andy( x) has exactly one maximum. If F(1) = 0 then there exists an X with 0 -C X < 1 such that y( x) > 0 for 0 < x -CX, andy( x) < 0 for X < x -C 1. Also y( x) has exactly one minimum and one maximum. In both cases 1y(x) 1 < max ) F(x) 1. Proof. At extrema

y,, =

y, = 0 and

-c-If

=

CO =0 1 >O

Thus y cannot have a minimum three possibilities, namely

forO
(24 Since y(0) = y(l) = 0 there are only

y>O

forO
y has exactly one maximum,

(2.7a)

y
forO
(2.7b)

y>O

forO
y has exactly one minimum, O
y-z0

forF
y has exactly one maximum

and one minimum.

We shall prove that if F(1) > 0 then (2.7b) and (2.7~) are not possible, then (2.7a) and (2.7b) are not possible. Let F(1) > 0. Suppose (2.7b) holds. Then Y,(O) G 0,

Y,(l) a 0.

(2.7~)

and that if F(1) = 0

G. Kreiss, H.-O. Kreiss / Convergence to steady state

1.c

.5

0

-. 5

1.0 0

.5

Fig. 6. The solution

1.0

of (1.2) when f = 0, a =l,

b = - 1 and c = 0.05.

BY (2.5) 0 < r(y,(I)

-y,(O))

This is a contradiction,

= -F(l)

< 0.

(2.8)

so (2.7b) cannot hold. Now suppose (2.7~) is valid. Then y,(O) > 0 and

by P.8) y,(O) >, r-‘P(1). If E is small enough this is impossible by Lemma 2.2. Let F(1) = 0. Assume that (2.7a) or (2.7b) are valid. By (2.8) y,(O) =y,(l), possible if y,(O) = y,(l) = 0. Differentiating (1.3) gives us

which is only

EYXXX =YYX, + (YX)’ -f,. Thus Y(0) =y,(O)

=y,,(O)

= 09

y,,,(O)

CO,

Y(1) =y,(I) =y,,(I) = 09 Y,,,(I) < 0. This implies that y must change sign at least once, which contradicts the assumption, and therefore (2.7~) must hold. It remains to show that 1y(x) 1 is bounded by max 1F(x) I. Since y(0) = y(l) = 0, the maximum absolute value of y is found at a local extrema, where y, = 0. Thus, from (2.5) it follows that

I Y(X) I G cm,?] I JIX) - EYX(O)I G oy:I I JI-4 I. . . . This finishes the proof.

0

G. Kreiss, H.-O. Kreiss / Convergence to steady state We can use the usual singular perturbation detail, see e.g. [3].

169

methods to discuss the behavior of the solution in

Theorem 2.1. Let F(1) > 0, assume that (1.3) has a solution and that c is sufficient& small. Then y( x) has a boundary layer at x = 1. For 1 - 0( e 1log(e) I) Q x < 1, y(x) is close to w(x) which is the solution of ewX=+w2-F(l),

-oo
w(l) = 0.

(2.9)

In any interval 0 < x0 < x < 1 - O(c 1log(c) 1) y(x) = h(x) + q(x,

c),

h(x) = IliF(x)

=:xg(x),

(2.10)

where u1 and its derivatives are bounded independently of C. For 0 < x < x0 < a we have y(x) = h(x) + et.@),

5Z= x/k,

where u and the derivatives d’u/dY converges to h(x) for 0 Q x < 1.

(2.11)

are bounded independently of e. Thus, for e + 0, y(x)

Proof. We indicate only the proof of (2.11). In the proof we shall use Ii,- I2 and I to denote the

intervals 0 G R 4 1, 1 G 2 G x,/6

and 0 d 2 < x,/h,

respectively. We shall also use

where I is an interval. We introduce a new variable in (1.3), y(x) = h(x) f EU(X/&). This gives us ~~~-(Zg(_x)+&u)u~-h~u= u(O)=O,

-h,,,

O&%x,/&, (2.12)

u(x,/J;)=u,,

where u. = ul(xo, E) is bounded independently of c. From x0 < (Yand the assumption f,(O) a fo > 0 it follows that h,(x) 2 ho > 0 for 0 < x
Iul G Ihxx/h.xl ~(l/ho)IIh,,(x)II,=:~. Thus IIu II I Gmaxb,, 4.

(2.13)

Next we want to estimate I]uz I]I. First we consider the interval Ii = [0, 11. By (2.12) and (2.13) there are constants C, and C2 such that

II uL%II I,

G G

II u2 II I, +

c2.

It is well known (see [4]) that one can estimate I(uz ]II, in terms of )Iu ]I1,, and I]ufp ]II,, i.e. for every constant S there is a constant C( 8) such that

II u; II I, G6 II UP2II I, + C(S) II u II I,* Thus for S = i(C,)-’ we obtain a bound for (1uzz ]I1,, which gives us a bound for ]Iu? I]r,. Especially, I u,(l) I is bounded.

G. Krelss, H.-O. Kreiss / Convergence to steucfy state

170

In the remaining

interval

F> F(h)

I, = [l, x,/6],

we have

= Q-,(0)(1 + O(h)).

Thus zg+&u=J2F/&+&>/~+o(&). i.e. for sufficiently

small 6

zg+&Gt/m. At local extrema

of up, upi = 0 and we have, by (2.12)

Thus

II UYII I, G max(Id)

I3 l4GJr)

I3 PI3

and U; is bounded independently of E in the whole interval. By differentiating (2.12) bounds higher derivatives of u can be obtained. It is also clear that as E + 0, y(x) converges to h(x). This finishes the proof. q

for

If F(1) = 0 then the solution switches at X from m + O(C) to - m + O(E). In each subinterval 0 < x < X and X < x < 1 the local behavior of the solution is of the same type as in the first case. As E -+ 0, y(x) converges to h(x) for 0 G x -CX and to -h(x) for X < x < 1. In general, the position of X can only be obtained by detailed calculation. However, if f(x) is antisymmetric around x = : then X = i. This is the only case we consider. We shall now discuss the existence of a solution. For this we need two lemmata. Lemma 2.4. For sufficiently Proof. By integrating

large E the steady state equation

(1.3) twice, we can write the equation

(1.2) has a solution. in the form

f~‘y2(t) dt - ~‘F(5) d5 + co= 0, or after the change of variable

y = nJ

For n = 0 the above equations have a unique sufficiently small n_ This proves the lemma. 0 Lemma 2.5. Let p(x)

be a smooth function.

A+= -(P44X++&

Therefore

the same

is true

for all

Consider the eigenvalue problem

444=+(1)=@

The eigenvalues are real and negative.

solution.

(2.14)

G. Kreiss, H.-O. Kreiss / Convergence to steady state

Proof.

We introduce

a new variable

171

G(x) by

and obtain xlc, = C\c/?..- crc/ =: L$J, #(O) =$(l)

44

= :A(4

(2.15)

+ (1/4+(x))*,

= 0.

(2.15) is selfadjoint and therefore the eigenvalues are real. Let + f 0, X be a solution of (2.14) and let 2 be the first zero of \c/ to the right of x = 0. We can assume that 4 > 0 for 0 -Cx < 5Z. Thus r&(O) 2 0 and &,(a) G 0, and integration of (2.14) gives us X/$(x)

dx=<[&];
0

It follows that h i 0. If X = 0, the only possible X -C 0, which proves the lemma. •I

solution

of (2.14) would

be c$(x) = 0. Thus

Now we can prove Theorem 2.2. The equation (1.3) has a unique solution for all c > 0. Proof. We have already shown that (1.3) has a solution for sufficiently large E. We will now employ continuation in c to prove existence for all E > 0. Assume we have shown existence for c > Z. We want to show that there is a solution for c = E. By Lemma 2.3 the solutions of (1.3) are uniformly bounded for 2 < E G E + 1. Therefore the same is true for the first three derivatives. Thus we can select a sequence of solutions Y(X, %)9

v = 1, 2,. . .‘)

lim E, = E, Y’oo

such that d’ lim yv(x, YAW dxJ

d’ 6”) = PY(X, dx’

and Y( x, E) is the desired solution.

r),

Linearizing

(Y(X, i)SY)X = @Y),,

j=

0, 1,2,

the equation

+ (C - E)Y(X, E),

around

Y( x, C) gives us

6Y(O) = 6Y(l) = 0.

By the previous lemma X = 0 is not an eigenvalue of the above equation solve (1.3) for all sufficiently small E - E. This proves the theorem. q

and therefore

we can

3. Speed of convergence In this section we want to discuss the speed of convergence to steady state. We assume that the initial data g(x) of (1.3) are sufficiently close to the solution of the steady problem, so that we only have to discuss the behavior of the solutions of the linearized equation w, + (YW), = cwxx, w(x,

0) = g(x),

O
t) = w(1,

t>o, t) = 0.

(3.1)

G. Kreiss,

172

To determine

the speed of convergence

Kreiss

/ Convergence

to steady

state

we study the distribution

of eigenvalues

of

@(O) = +(1) = 0.

A+ + (Y@>X= E%xr Theorem

H.-O.

(3.2)

3.1. The eigenvalues of (3.2) are real and negative and their distribution

is given by (1.5).

Proof. Lemma 2.5 tells us that the eigenvalues are real and negative. First we consider the case f = 0, a > -b. We write (3.2) in the selfadjoint form (2.15) with p = y. Let h = A, be the largest eigenvalue. The corresponding eigenfunction $r does not change sign, and we can assume that $r > 0 for 0 < x < 1 and that max ) IJ~(x) 1 = 1. We assume that X, > -a2/8c. Then there is a constant K such that c(x) + X, > 0 for 0 z x >, 1 - Kc. Thus $I is monotone in the interval 0 < x 6 1 - Kc, and therefore It/r must have its maximum in the remaining interval, 1 - KE G x -c 1. By assumption max It/r(x) = 1 and therefore there must be a constant 6 > 0 such that +rX( 1) G -S/e. Now consider the corresponding eigenfunction @r(x) = exp(fr-‘~xy(S) Integrating

d6]+,(x),

+1,(l)

0 G G,(x)

= k,(l)~

G G4.

(3.2) gives us -8 2 e(+,,(l)

- @r,(O)) = XI/al+r dx 2 A,i’exp[fe-’

iXy do

dx = X,cd.

Thus h, < -min(

a2/8c,

a/de),

and the theorem is proven for this case. estimate When f+O, a=b=O, and Jdf(x)d x > 0 the corresponding way, since by Theorem 2.1 there are constants C, > 0 and K such that c(~)=f(h,(x)+O(~)+~~~‘(h(x)+0(\1;))~j>C~~O

1% I G II UII function

as trial function.

y(x)

Both +2 and ($r-‘y2

a = -b,

f = 0 or a = b = 0 and

f(x)

is

eel”).

We shall use the fact that for our selfadjoint smallest absolute value, A,, satisfies

for any smooth

in the same

forO,(x
We now consider the antisymmetric case when antisymmetric around x = :. We want to show that -A, = O(r-’

follows

2/u

is antisymmetric

II L+ II 2 = 2L1”(

problem

(2.15) the eigenvalue

@II23

+ + 0 satisfying

+ iy,)’

eigenvalue

the boundary

around

are symmetric

II + II 2= 2/o”‘exp (-cP1[‘2y

We chose

x = 4, and +(O) = +(l) = 0. Also

around

y2/4c + yX/2)2exp(

conditions.

x = 4. Therefore

-~~1~1’2y

d
d[)

dx,

d[]

- 112 dx

with the

173

G. Kreiss, H.-O. Kreiss / Convergence to steady state

and by (2.5) and Theorem 2.1 A2<

IIWl12

=

l14d12

l’ < .

~2~-2

l”(

y2/4c +~,/2)~

dx,/01’2( exp( &-l /b;

dr) - I)’ dx

e-2D/c 7

where C > 0, D > 0 are constants which do not depend on c. We shall now estimate the size of the second eigenvalue for the case with an interior boundary layer at x = 3. By assumption y(x) is antisymmetric around x = 3. Consider the eigenvalue problem (3.2) on half the interval, 0 G x G :, and denote its solutions by &(x),

Xi,

i=1,2

)... .

We know that ii has i - 1 sign changes, and we have already shown how the 1,‘s are bounded away from zero. The function forO
6iCx)

G2i(x)

=

i

-$i(~-+)

i, z=1’2’*“’

will satisfy (3.2) on the full interval, 0 Q x < 1 with X = X2, = xi. Also +2i changes sign 2( i - 1) + 1 times. Thus +2i is the 2i th eigenfunction and h2i is the 2ith eigenvalue. Therefore A, is bounded away from zero. This finishes the proof. 0

4. Numerical results We shall discuss difference approximations for the time dependent eigenvalue problem (3.2). We introduce gridpoints xi=ih,

i=O,l

tj = jk,

j=O,

,***,

1 ,***, N,

problem (1.2) and the

h= l/N,

where N is a natural number and k > 0 is the time step. We also introduce gridfunctions Ui’= U(Xi, ij). We approximate

(1.2) by the usual implicit method

(I-~kD_D_)uit’+fkDo(u:+‘)2=ui+k~,

i=l,2,...,N-1

(4.1)

with initial and boundary conditions

IA;=&

i=l,2

ui = a 7 +b.

,..., N>

j=l,

1, 2,...

.

Here h2D+D_ui

= ui+l - 22.4,+ ui_l

and

2hD,( ui)’ = ( u~+~)’ - ( u~__~)~

denote the usual centered difference operators. At every time step one has to solve a nonlinear system to determine u:+l. This is done by the iteration (I-ckD+D_)u;‘+‘)=

-+kDo(uj’))2+u:+kfi,

where u(O)is choosen by a predictor process.

I=O,

1, . . . .

(4.2)

174

G. Kreiss, H.-O. Kreiss / Convergence to steady state

1.c

.50

.25

0

.50

.25

.75

Fig. 7. Convergence when the shock is located at the boundary. N = 50, k = 0.1. Between each curve there are 5 time steps.

1.0 Here c = 0.04, f(x)

= asin(

U(X, 0) = isin(

In all our experiments the solution of (4.1) converges to a steady state solution. However, the speed of convergence depends on the location of the shock. If the shock is located at the boundary, corresponding to the first and third case of (1.5), then the convergence to steady state is quite rapid. See Fig. 7. If on the other hand the shock is located in the interior, corresponding to the other cases of (1.5), the convergence is, in general, very slow. When the shock is formed at an early stage it is in general in the ‘wrong’ place, depending on the initial data. From then on, the the shock moves slowly to the correct position. See Figs. 1 and 3. This process can be considered quasi-stationary, which makes it possible to use the same convergence acceleration as in [2]. Formally we can write our iteration (4.1) as H( un+i) = un+i - Un:= y”_

(4.3)

We can linearize the relation and obtain (I-

L)r”+’

= rn.

(4.4)

In our case Lr, = ddl+D_r,

- kD,( u:+lr,).

(4.5)

This is a discretization of the right hand side of the eigenvalue problem (2.14), with p = u”. If the process is quasi-stationary we can consider L to be independent of n. Then we have r n+j= (I_

L)-jrn

G. Kreiss, H.-O. Kreiss / Convergence to steady state

175

and P--l U

n+p=un+

c (I-L)-+.

j=O

If the eigenvalues hi, of L are negative the eigenvalues

K~,

of (I - L)-' satisfy

IK~ 1 < 1

lim un+P = u"+(I-(I-L)-l)-lr"=u"+(I-L-l)r".

and (4.6)

p-00

Instead of taking a large number of time steps we can take one large step, which we call an extrapolation step. We put u = u* + be,

(4.7)

where e is the solution of the equation Le=( L-l)r",

(4.8) and j3 is a stabilizing parameter. We choose /3 in such a way that H( un + j3e) has no component in the direction of e, i.e. (H(u”+

Be), e) = 0,

where ( . , .) denotes the usual inner product. There are other possible choices, for example choose p such that IIH(u” + be) II = $n II Hb” + Be) II. Of course (4.7) is not the steady solution we are seeking. We use the new u to restart the time iteration, and make a new extrapolation step once a new quasi-stationary state is reached. In our experiments we use an a priori fixed number of time steps between the extrapolation steps. Better strategies are under development. We have calculated the first eigenvalues and eigenvectors of the discrete linearized operator (4.5), provided MY+’ is the discrete steady state solution. The calculations show that the eigenvalues are negative and their distribution is of the same type as for the corresponding continuous case. See Table 1. In Figs. 8 and 9, the first few eigenvectors are plotted. Note that in

Table 1 Eigenvalues discretization Eigenvectors

of the eigenvalue problem (3.2) y is the solution of (1.3). Three different cases were treated. The is done according to (4.9, with N = 100 gridpoints. The eigenvalues were found using inverse iteration. corresponding to case (1) are plotted in Figs. 6(a) and (b) A,

A,

f(x) = sin(2ax)/2 a=b=O c = 0.04

- 8.64.10-3

- 4.34

- 5.32

f(x) = sin(2nx)/2 a=b=O e = 0.02

-4.62.10-6

- 5.617

- 5.622

- 1.24.10-9

- 12.8

A,

f(x)=0 a=l, b=-1 E = 0.02

- 13.5

176

G. Kreiss, H.-O. Kreiss / Convergence to steady state

.a) .275

1 0

I

I

I

I

.25

.50

-75

1.0

I .75

1 1.0

.275

-. 275

-. 55 I 0

I .25

I .50

Fig. 8. Eigenvector. (a) The first two eigenfunctions of problem (3.2), when y, the solution of (1.3), has a shock in the interior. In this case f = 0.04, f(x) = &r sin( 71x) cos( PX), a = b = 0, N = 100. (b) The third and fourth eigenfunctions of problem (3.2), when y, the solution of (1.3), has a shock in the interior. In this case c = 0.04, f(x) = in sin(nx) cos(nx), a = b = 0, N =lOO.

G. Kreiss, H.-O. Kreiss / Convergence to steady state

0

.25

.50

177

1.0

.75

Fig. 9. Eigenvectors. The first two eigenvectors, $1 and &, of problem (3.2), when y, the solution of (1.3), has a shock x=l.Inthiscase r=O.O8,f(x)=fasin(nx), a=b=O, iV=lOO.

the case of an interior shock the first eigenvector is exponentially small away from the shock region. Also, we have no doubt, and it is confirmed by the calculations, that the position of the shock does not change the nature of the eigenvalue distribution. In fact, in the proof of Theorem 3.1, y can be replaced by any function of the same structure. In our case, when the shock is located in the interior, (I - L)-’ has only one eigenvalue, IQ, close to zero. All other eigenvalues are small. Therefore, when we have reached the quasi-stationary state, r” is in the direction of the eigenvector corresponding to rci. See Fig. 10. Therefore we do not need to solve (4.8), and instead of (4.7) we use u = u” + pr”. (4.9) In Figs. 2 and 4 we have plotted accelerated. 5. A two-dimensional

u at different time stages to show how the convergence is

case

Consider the following problem u, + (+u2)x = &X u(O, y, t) = a, u(x,

0,

t) = u(x,

+ UyJ

O
u(I, y, 1) = -a, 1, t) = w(x),

O
t>,o,

a > 0,

u(x, Y, 0) =g(-%

(5 .I) Y),

178

G. Kreiss. H.-O. Kreiss / Convergence to steudj, stute

-.025 -

-. 0375 -

-.05/ 0

.50

.25

1.0

.75

Fig. 10. Differences between consecutive solutions at different time stages, when E = 0.04, f = in sin(ax) COS(VX), a = h = 0, u(x, 0) = i sin( TX). Between each curve there are 100 time steps. The calculation is made with time step k = 0.1 and N = 50 grid points.

where w(x) is the solution of the one-dimensional problem (1.3) with b = -u, and f(x) (2.3). A steady solution of (5.1) is u(x, y) = w(x). The speed of convergence can be studied by analyzing the corresponding eigenvalue ~~+(w~),=~(~*~+~_~“), We can solve (5.2) by separation

+=Oontheboundary. of variables.

(WX)’ - CX” = AX, Y” = -qY, with p=X)... j=2,3

x(0)

= x(1)

Let +(x,

(5 -2)

(5.3a)

= 0,

(5.3b)

Y(0) = Y(1) = 0,

4, =

(Jo)‘,

problem

y) = X(x)Y( y). Then

cq. We recognize (5.3a) as (3.2). Therefore . We can solve (5.3b). The solution is

Y,(y) = sin(jv),

= 0. See

j=l,2,

-h,

= O(e-I/‘)

and

-X, > 0(1/e),

... .

corresponding There is a whole sequence of eigenvalues, pt,, of order O(C). The eigenfunctions to this sequence, ~$r~,will be exponentially small away from the shock. All other eigenvalues will be of order 0(1/e). We expect that the time iteration will again lead to a quasi-stationary state, and that the residual will be composed of eigenfunctions corresponding to the eigenvalues of order O(C). Therefore e in (4.8) will be of the same form, and we can replace all components of e away from the shock by zero, thus obtaining a linear system of equations of order N instead of N*. More details will be given in another paper.

G. Kreiss, H.-O. Kreiss / Convergence to steady state

179

References [l] M.D. Salas, S. Abarbanel and D. Gottlieb, Multiple steady states for characteristic initial value problems, ICASE Report No. 84-57, NASA CR-172486, November 1984; also: Appl. Numer. Math. 2 (this volume) (1986) 193-210. [2] M. Hafez, E. Parlette and M. Salas, Convergence acceleration of iterative solutions for transonic flow computations, AIAA 85-1641. [3] J.D. Cole and J. Kevorkian, Perturbation Methocis in Applied Mathematics (Springer, Berlin, 1981). [4] E. Landau, Einige Ungleichungen fir zweimal differenzierbare Functionen, Proc. London Math. Sot. 13 (1913) 43-49.