NASA Contnetor Report 178017

lCASE REPORT NO. 85-50

NASA-CR-178017

19860009593

leASE . .... ,"', : .. ".:-. " ." ", .

'......",' ".':v ....... ,..",- .......,•..

CONVERGENCE TO STEADY STATE OF SOLUTIONS OF BURGERS' EQUATION

Gunilla Kreiss

Heinz-Otto Kreiss

.: ': :.::~.~.~ -~_l " ;

:...:s....;J ;....: -'._,. '" "c:'.

;•• J

.

., .. ' .~~ -

Contract Nos. NAS1-17070

December 1985

Lfd<c.:...='" R~::;t~:\!?:,"! ::=~; 1::::1

L:cr:-.;~Y, :"~;;\':..'. \

,: .!;"\~,·l,r:T(Jr'J •. VIRGlj'fl/(

INSTITUTE FOR CO~WUTER APPLICATIONS IN SCIENCE AND ENGINEERING

NASA Langley Research Center, Hampton, Virginia 23665

Operated by the Universities Space Research Association

NI\SI\

Na~onal Aeronautics and

Space Administration

Langley Research center

Hampton, Virginia 23665

iii.6NJ.9064,._-# Ii_SSUE 9 PAGE o. ,,.s,,,.c:.:_-! CATEGORY 6"+ RPT_!!_: NASA-CR-i.?S01'7

!C:ASE-!-3_:._-L;:_O HAS ;1..F!6:1'7_:i301'7 L.','tq'F#: NASi-1'70'70 _'_':]:":]iO!X_.-iii!,_i_-!-::i-OZ,,,._!:f_.]-" NSi:::

DNS---i32-;[E'{!!d;,4 :......... ::_':;;:_._T.:iO0 _3:!. PAGES UNC:L_ASS!F]]TD DC!CUMENT

T'TL..;: {::::enw_:.rger_ce t:o s_;e--b,.:!y s't;at:@ of sek-".-.,,!:;:_.ons of gur_]ers:' _--_. q '-',-".i::i ,:::n, Ti.__]!P:,.

I:::'i rta ! [-;_ep,::,,,-t

UTH,- Aii<F:-;!!-.1;SS:, G. ii 2iI.(RE'];SS, H,, O. PP:A',: Ai(l:;.:c,,/s.! !n,_!!._s, of" Te,::'_, , Stc, cl.::hi::,].;-r,,,

_-_!,-,_ede_-}}; £{!i'i(C_::_;l.i!_(.:,i-;"]2._;_]i r; _:_!::. o .1: T,---,ch, , P,_:,._:s.-;'..dc.r,.::,)

o'r:_P" l...!_:_,..q:,l"Jai::!ic'_"_"__i.!::r,r,,, V.::.a,,A'.::"""°nauAvA..iJL t_,ci:!!- -:<,qd,,i..!T!;3!3Pa.Ce #idm'.ri:).EE.'{';',"'Tt1⁄2:[OT',,. L-":".',"','g!'z_'VF<.,..-::,_-i.,_:-.:,a[]e:,,,,,_;e:i,r_ rch

AP_, Hi[:; AO[-};iNF A,O:).

!O',: UNITED STATES Sub,,q_.tted 'for- pub!},.':::_:_:i;ion

,,..,,..,.'""-".:.:s,,: i"<'BC)!JN_],ARY ',.,."ALLJE PRO]ifiLEMf]i;ie_fURGER Et]_Ui_:,T]:DNi._:.-.C(])NVERi3ENC,..E/eF,ER.TUR:C_AT]]GN

THEORY / eSTE,,."-Dy STATE

INS- / E]:[-)ENVALUES- G2,,S )DYNAM._CS/ PRD:E_:LEI4 SD,"....ViNfi_i SIN!3ULAR]:Ty (HATI!EHAT:[GS)i

"f!!dE. DEPEN!)ENCE

i3A- Au L:h,:::,r

CONVERGENCE TO STEADY STATE OF SOLUTIONS OF BURGERS" EQUATION

Gunilla Kreiss

Royal Institute of Technology

Stockholm, Sweden

and

Helnz-Otto Kreiss

California Institute of Technology

Pasadena, California

Dedicated to Milton E. Rose

on occasion of his sixtieth birthday

Abstract

Consider the initial-boundary value problem for Burgers" equation. It is

shown that its solutions converge, in time, to a unique steady state. The

speed of the convergence depends on the boundary conditions and can be

exponentially slow. Methods to speed up the rate of convergence are also

discussed.

Research was partially supported by the Office of Naval Research under

NO0014-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. NASI-17070 while the authors were in

residence at the Institute for Computer Applications in Science and

Engineering, NASA Langley Research Center, Hampton, VA 23665-5225.

i

1. Introduction. In many gasdynamical problems one tries to calculate tile steady state

solution by solving the corresponding time dependant problem. One hopes that for t --* oo the

solution converges to a unique steady state. Recently, M. D. Salas, S. Abaxbanel and D. Gottlieb

[1]considered the initial-boundary wlue problem

ut -t- _(u2)x = f(x), t __O, O< x < Tr,

(1.1)

_(x,0)=_(x).

They used

f(_)=sin_cos_, g(_)=bsinx, O<b,

and showed that the solution u(x, t) of the above problem converges to a steady state v(x), as

t --_ c_, but that v(x) depends on the initial data.

In this paper we consider the viscous problem

1

ut+-'(u2)x=6uxx+f(x), t>0, 0<x< 1, e>0, (1.2a)

2

with initial and boundary conditions

u(x,O) =g(x), (1.2b)

u(0, t)=a, 4(1, t)----b,

and the corresponding steady state problem

l(y2)x=_y..+f(x), 0<_<1, _>0,

(1.3)

y(0)=a, y(1)=b.

For simplicity we restrict ourselves to two cases:

1) a>O>_b, a>_-b, f(x)-O,

2) a = b = 0, f is such that there exists an (_ with 0 < a < 1 such that f(x) > 0 for 0 < x < a,

f(x) < 0 for o_< x < 1, f(O) = f(1) = O, f:_(O) > fo > 0 and fx(1) > .to.

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 y(x) as _ --* 0 exists. Thus, if

nm,_(zt,) = y(:_)

t_Oo

1
exists, we obtain a unique steady state solution of the inviscid equation (1.1) if we first let

t _ c_ and then _ _ 0. This is in contrast to the procedure in [ 1], where the two limit

procedures are taken in the reverse order.

We shall prove that the eigenvaluesof the eigenvalueproblem

),p = -(y_)_ +6p_, p(0) = p(1) = 0, (1.4)

are all negative. Therefore, the solution of (1.2) converges to the solution of (1.3) provided

u(x, O) = g(x) is sufficiently close to y(x). In another paper we shall prove that u(x, t) converges

to y(x) as t --_ oo for arbitrary initial data. The speed of convergence is determined by the

eigenvalues, )_j, of (1.4). We shall show that the eigenvalue distribution depends on f(x) and

on a, b in the following way:

There is a constant c > 0 which does not depend on _ such that

(1) if a>-b, f--0 then 0>-c/_>.XI>.X2>...

(2) if a=-b, f-=0 then -Al=O(e -1/_)>0, -e/e>)_2>A3>...

(1.5)

(3) if a=b--O, I( )dx#O, then

_o 1

(4) if a ----b = O, f(x)dx = O, then - ),I = O(€-1/€) > O, --e > ,_2 _> _3 _>''"

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 -A1 = O(e-1/_). This

is confirmed by numerical experiments. We see that at first u(x, t) quite rapidly approches the

same hmit as the inviscid equation (1.1), which consists of solutions of the stationary equation

connectedby a shock. Oncethe viscousshockhas beenformed,the solutionof (1.2)becomes

quasi-stationary and the shock creeps extremely slowly to the "right" position. We can explain the behavior,becauseby linearizingaroundthe quasistationarysolutionwefind that the

eigenvaluesof the correspondingeigenvalueproblemhave a similardistributionas earher.

If -A1 = O(e-1/*) then the speed of convergenceis so slowthat the abovemethod to

calculatethe steadystate isimpractical,see figures(1)and (3). However,wecan usethe same

technique as Ha:fez, Parlette and Salas in [2]to speed up the convergence. See figures (2) and

(4).

Unfortunatly, not only the speed of convergence but also the condition number of the

stationary problem deteriorates. We have to calculate with O(el/_) decimals to obtain correct

2
results. To avoid an excessive number of decimals we have used a quite large z 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

--_1i = O(j2e), j = 1, 2,...,

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 D. Brown, W. Kath, H. O. Kreiss and W. Henshaw, M. Nanghton (private communication).

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

with uniqueness, which can be proven by standard techniques.

Lermna 2.1. If the steady equation (1.3) has a solution, then it is unique.

Proof. Let u, v be two solutions. Then w = u - v is the solution of

= =

=ew_, P u+v, w(O) w(1)=0. (2.1)

If w _ 0 then the zeros of w are isolated. Let _ with 0 < _ <_ 1 be the first zero to the right

of x =0. Without restriction we can assume that w > 0 for0 < x <7, i.e. wz(0) > 0 and

w_(_) _<0. Integration of (2.1) gives us

-*(Iwx(_)l + Iw_(O)l)= _[_=]g = _[pw]o = o.

Thus wz(0) = w=(_) = 0. We can consider (2.1) as an initial value problem with initial data

w(0) = w_(0) = 0 whose solution is w(x) = O, and the lemma is proved.

We shall now discuss the properties of the solution. Let us start with the case f(x) - O,

a > 0 > b, a > -b. Integrating (1.3) gives us

1 2

eyx=_y --c, O<x<l,

(2.2)

v(0)= a.

3
The constant c has to be determined so that y(1) = b. We necessarily have e = d2/2 > a2/2,

because with c <_a2/2, Yx >- 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)

2e f _2d_-_d2 -fd ,

a 0

i.e.

- d = edz/_.

Therefore y(1) = b implies d = a + O(e-'/'), and

1 - re-a(x-x)/_ a - b

y(x) = al + re-a(1-x)/_' with r = --a+b" (2.3)

Away from the boundary layer at x=l we have y(x) = a +O(e-a(X-z)/_). Thus, for e --+0, y(x)

converges to a for 0 < x < 1.

If a = -b we consider (2.2) on the interval 0 < z < 1⁄2, with boundary conditions y(0) =

a, y(1⁄2) = 0 and obtain a solution yt(x) of the form (2.3). The solution on the whole interval

is given by

Yl(X), if0<x< 1

--2'

1 <x<l. y(x)= --yx(1-- x), if__ _

In figures (9) and (10) we have plotted y(x) for two different sets of boundary values.

Consider case 2, where f only vanishes at x = O,c_,1 and a = b = O. Without restrictions

we can assume that

1

f(x) > O. (2.4)

0

If this is not true, we transform the problem by introducing new variables,

_=l-x, /=_jr, Y=-VThe new problem satisfies (2.4).

Lemma 2.2. Let y(x) be the solution of (1.3), F(x) = fo f(_)d_ and h(x) = _.

Then

yx(1) < yx(0) _<If,, gl = oma., u z<l {Is h_(z)l}+ Ih=do)l.

4
Proofi Integration of (1.3) gives

1 2

_(vx- vx(0))= _v - F,

(2.5)

v(0)=o,

where y_(0) is determined by y(1) = 0. If u = y - h, then u is the solution of

.1-x2

u_ = yx(O)- h_+ _-luh _ -_ u ,

=o.

Assume that yx(0) > K1. It follows that yx(0) -hz(x) is positive and thus u and ux are positive

for all x > 0. In particular u(1) > 0 and y(1) = u(1) + h(1) > 0, which contradicts y(1) = 0.

Thus y_(0) _<KI. Also

ey_(1) = eyx(O) -- F(1) _<eyx(O).

This proves the lemma.

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 and y(x) has exactly one maximum. If F(1) = 0 then there exists

anbwith0<5<lsuchthaty(x)>0for0<x<5, and y(x)<Oforb<x<l. Alsoy(x)

has exactly one minimum and one _,._a.ximum.In both cases [y(x)l <_max [F(x)l.

Proof. At extrema yx = 0 and

<0 for 0 < x < o_,

Yxx-----e-lf= ----0 forx----_, (2.6)

>0 for c_< x < 1.

Thus y cannot have a minimum to the left of a maximum. Since y(0) = y(1) ---0 there are only

three possibilities, namely

y > 0 for 0 <: x < 1, y has exactly one maximum, (2.7a)

y<0 for 0<x<l, y has exactly one minimum, (2.7b)

y>0 for 0<x<5, 0<5<1, (9.7c_

y < 0 for 5 < x < 1, y has exactly one maximum and one minimum.

5
We shall prove that if F(1) > 0 then (b) and (c) are not possible, and that if F(1) = 0

then (a) and (b) are not possible.

Let F(1) > 0. Suppose (2.7b) holds. Then

y=(o)_<0, v=(1)> 0.

By (2.5)

0 < e(yx(1) -- yx(0)) = -F(1) < 0. (2.8)

This is a contradiction, so (2.7b) cannot hold. Now supposse (2.7e) is valid. Then yx(0) > 0

and by (2.8)

yx(0) _>e-iF(l).

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) yz(0) -- yz(1), which is

only possible if y_(0) = y_(1) = 0. Differentiating (1.3) gives us

eYzzx = YYxx + (Yx)2 -- fx. (2.9)

Thus

y(0)= y_(0)=y=(0)=0, w=(0)<0,

y(1) = yz(1) = yx_(1) = 0, Yxzx(1) < 0.

This implies that y must change sign at least once, which contradicts the assumption, and

therefore (2.7c) must hold.

It remains to show that ly(x)l is bounded by maxlF(z)l. Since y(0) = y(1) = 0, the

maximum absolute value of y is found at a local extrema, where yz = 0. Thus, from (2.5) it

follows that

[y(x)[ < max IF(x)-eyz(0)l < max [F(x)l.

--o_<x<1 --o_<z_<1

This finishes the proof.

Wecan use the usual singular perturbation methods to discuss the behavior of the solution

in detail, see for ex. [3].

Theorem 2.1. Let F(1) > 0, assume that (1.3) has a solution and that e is sufficiently

small. Then y(x) has a boundary layer at x = 1. For 1- O(e[log(e)]) < z < 1, y(x) is close to

w(z) which is the solution of

1 2 F(1), -_<x< 1, w(1) 0. (2.9)

In any interval 0 < x0 _<x _<1 -- O(e Ilog(e)[)

y(x)=h(x)+e_,(x,e),h(x)=_=: _(_), (2.10)

6
where ul and its derivatives are bounded independantly of e. For 0 < z < x0 < a we have

y(x) = h(x) + $u(5:), 5:= x/v_, (2.11)

where u and the derivatives dUu/d_Y are bounded independantly of e. Thus, for e --*0, y(x)

converges to h(x) for 0 < x < 1.

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

denote the intervals 0 < _ < 1, 1 < _ _ xo/_rg and 0 < Yc< xo/x/'g, respectively. We shall also

use

Ilfll,:= m_tIf(_)l,

where I is an interval.

We introduce a new variable in (1.3),

V(x) = h(z) + eu(z/Vq).

This gives us

_ - (_g(_)+ v_)_ - hx_= -h=_, 0_<_____0/vq, _,(0)= 0, _(_0/v_)=_o,(2.12)

where uo = ux(xo,e) is bounded independantly of e. From xo < c_and the assumtion fx(0) _>

fo > 0 it followsthat h_(x) > h0 > 0 for 0 < x < x0. Therefore we can use the maximum

principle. The maximum of u is found either on the boundary or at a local extrema, where

ue = 0. At local extrema

I_1<-[h=_l<-h__Uh,=(x}Ut=: c_.

Thus

II_11_ -<max(uo, c_). (2.13)

Next we want to estimate Ilu_llt. First we consider the interval Ix = [0,1]. By (2.12) and

(2.13) there are constants Cx and C2 such that

II_i_llr,-<c, Iluilb,+ c2.

It is well known, see Landau [4],that one can estimate Ill,lit,in terms of II_llt,,andII,,_llt,,

i.e. for every constant 6 there is a constant C(6) such that

I1_11--- ,, 611_llt,+ c(6)ll_lb,.

Thusfor6 = 1⁄2(C-_l) weobtaina boundforII_llt,, whichgivesus a bound for II_lb,.

Especially, lu_(1)lisbounded.

In the remaining interval I2 = [1,Xo/v/'g], we have

F _>F(x/_ = efx(0)(1 + O(vq)).

Thus

sg+vT_=,/_/v7 +,/7__3 +o(,/z),

7
i.e. for sufficiently small x/_

1

_g+ vq,__>_d-_.

At local extrema of u_, u_ = 0 and we have, by (2.12),

Thus

x0

II_+II__, _ max(lu+I_+(_)I, (l)l, #),

and u_ is bounded independantly of e in the whole interval. By differentiating (2.12) bounds

for higher derivatives of u can be obtained.

It is also clear that as e -+ 0, y(x) convergesto h(z). This finishesthe proof.

If F(1) = 0 then the solution switches at _ from x/_+ O(_) to -v/_ + O(e). In each

subinterval 0 < z < _ and _ < z < 1 the local behavior of the solution is of the same type as

in the first case. As _ --*0, y(z) converges to h(x) for 0 _<x < _ and to -h(z) for • < x _<1.

In general, tlte position of _ can only be obtained by detailed calculation. However,if f(z) is

antisymmetric around z = _I then -x = 3"1 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 large €the steady state equation (1.2) has a solution.

Proof. By integrating (1.3) twice, we can write the equation in the form

F(+)+++,+++o, ,+:,,°,

0 0

I I

1⁄2.]"++(+),++- ]_,,-'(+),++ ++o=O,

o o

or after the change of variable y = r/_

z X

+,(+:++,,+1+,+,+)++-/+(+)+++:+

0 0

I I

++o:.

o o

8
For r/ = 0 the above equations have a unique solution• Therefore the same is true for all

sufficiently small 77.This proves the lemma.

Lemma 2.5. Let p(x) be a smooth function. Consider the eigenvalue problem

>,to-----(p!a)_ + _tox_, !a(0) = VP(1)= 0. (2.14)

The eigenvalues are real and negative.

ProoL We introduce a new variable ¢(x) by

1⁄2_-'f v(_)d_

,

and obtain

,k¢ -- _¢=z - c¢ --: L¢, c(x) ----1 1

7pz(x) + _o (p(x))2,

(2.15)

¢(0) = ¢(1) = 0.

(2.15) is selfadjoint and therefore the eigenvalues are real. Let _o_ 0, ), be a solution of (2.14),

and let :_be the first zero of _oto the right of x = 0. We can assume that _ > 0 for 0 < x < ._.

Thus _x(0) _>0 and _z($) _<0, and integration of (2.14) gives us

0

It follows that _ _<0. If ), = 0, the only possible solution of (2.14) would be _(z) - 0. Thus

)_< 0, which proves the lemma.

Now we can prove

Theorem 2.2. The equation (1.3) has a unique solution for all e > 0.

Proof. We have already shown that (1.3) has a solution for sufficiently large e. We will

now employ continuation in e to prove existance for all e > 0. Assume we have shown existance

for e > g. We want to show that there is a solution for e = g. By lemma 2.3 the solutions of

(1.3) are uniformly bounded for g < e < g + 1. Therefore the same is true for the first three

derivatives. Thus we can select a sequence of solutions

y(x, ev), u = 1,2,..., lim eu ='g, b*-"* CO

such that

• di di

]]moo-_xyy(x, ev) = -_xiY(X,'d), j = O, 1, 2 ,

9
and y(x,-g) is the desired solution.Linearizing the equation around y(x,'g) gives us

(y(x,_)6y)z = e(by)zx + (_ -_)y(x, _), by(O) = by(i) = 0. ,.

By the previous lemma ), = 0 is not an eigenvalue of the above equation and therefore we can

solve (1.3) for all sufficiently small _ - _. This proves the theorem.

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

wt + (yw)x = ew_, O< z < l, t >_O,

0)= (3.1)

w(0, t) = w(1, t) = 0.

To determine the speed of convergence we study the distribution of eigenvalues of

_ + (W)_ = _x_, _(0) = _(1) = 0. (3.2)

Theorem 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 A = A1

be the largest eigenvalue. The corresponding eigenfunetion _bl does not change sign, and we

can assume that _bl > 0 for 0 < x < 1 and that max I_bl(x)l = 1. We assume that A1> -a2/8e.

Then there is a constant K such that c(x) + A1> 0 for 0 > x > 1 - Ke. Thus _blis monotone

in the interval 0 < x <_ 1 - K_, and therefore _bl must have its maximum in the remaining

interval, 1- ff_ _<x < 1. By assumption max _bl(x) = 1 and therefore there must be a constant

3 > 0 such that _bxx(1)< -6/_. Now consider the corresponding eigenfunction

_l(X)=e ' _l(X), _1_(1)=_(1), 0 _<_, (x) _<q_,(z).

Integrating (3.2) gives us

fo

-b > e(Wl,(1) - WI_(0)) = ,kl _dx >

_01 = d

> AI e1⁄2*-_f, v _dz= Aled.

lO
Thus

_,x< - man• a2(=,£1,

de

and the theoremis provenfor this case.

Whenf _ O,a = b= O,and f_ f(x)dx > 0 the correspondingestimatefollowsin the same

way,sinceby theorem2.1there are constantsCo > 0 andK suchthat

c(x)=2(h_(z)+O(v/e')+le-'(h(z)+O(v_)2)->G0>O for O<x<l-Ke.

We now consider the antisymmetric case when a = -b, f = 0 or a = b = 0 and f(x) is

We want to show that antisymmetric around x = 3"

-_,1 = O(e-le-1/_).

We shall use the fact that for our selfadjoint eigenvalue problem (2.15) the eigenvalue with the

smallest absolute value, ),x, satisfies

for any smooth function € _ 0 satisfying the boundary conditions. We chose

_ i -_ lf/2 _(_)a_

€(_)=_ 1/_ __ o

as trial function, y(x) is antisymmetric around x = 1⁄2, and €(0) = €(1) = 0. Also

!12

15.--1

y2 y_, -_ of_(_)d_

L6 = (_ + -yJe

1 Therefore Both ¢2 and (1⁄4e-ly2 + 1⁄2yz)2 are symmetric around x = 3"

1/2 2

IIL¢I2J= 2f (_ + u_2e-_-12f_"12Ud_dx'

o

112

f _ /-1/2 i -t

II€1--122 €-_1"o ud_(e_ f°_d_-- 1)2dx,

0

11
and by (2.5) and theorem 2.1

,4, + q) 2d_

)`_< IILCll 2_ o _<c2_=e_2w,,

I1€112,/_ ",-,_,,d_

f(e o -1)2d_ 0

where C > 0, D > 0 are constants which do not depend on e.

We shall now estimate the size of the second eigenvalue for the case with an interior

boundary layer at x = _.1 By assumption y(x) is antisymmetrie around x = 3"IConsider tile

eigenvalue problem (3.2) on half the interval, 0 < x < I and denote its solutions by

_Si(x), ),i, i=1,2 ....

We know that _5ihas i- 1 sign changes, and we have already shown how the _i's axe bounded

away from zero. The function

1

W2i(x)= _Si(x) for0<x_<_ i----1,2

--#i(x -1) for1⁄2<x_<l .... '

will satisfy (3.2) on the full interval, 0 _<x _< 1 with ), = )`2i = _i. Also _i changes sign

2(i - 1) + 1 times. Thus _2i is the 2ith eigenfunction and ),2i is the 2ith eigenvalue. Therefore

),2 is bounded away from zero. This finishes the proof.

4. Numerical results. We shall discuss difference approximations for the time dependant problem (1.2) and the eigenvalue problem (3.2). We introduce gridpoints

1

(xi=ih, ti=jk), i=0,1,.., j=0,1,...,N, h=_,

where N is a natural number and k > 0 is the time step. We also introduce gridfunctions

_ =_(x,,t;).

We approximate (1.2) by the usual implicit method

(I-ekD+D-)ui+'+lkDo(_+x)2----_+kfl, i-- 1,2,...,N- 1 (4.1)

12
with initial and boundary conditions

u°=9i, i=l,2,...,N-1,

 N=b, S=1,2,...

Here

hZD+D-=_= =_+1- 2t_+ __, and 2hOo(=_)2= (=_+,)2_ (__,)2

denote the usual centered difference operators. At every time step one has to solve a nonlinear

system to determine _+1 . This is done by the iteration

(I- ekD+D_)u} t+l) - lkn tu(0'2 (4.2) --5 _,0t_ j +_+kf_, 1=0,1... ,

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

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 figure (5). 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 figures (1),(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(u n+l) = u"+1 - un := rn. (4.3)

We can linearize the realation and obtain

(I - L)r n+l = rn. (4.4)

In our case

Lri = ekD+O_ri - kho(up+lri). (4.5)

This is a discretization of the right hand side of the eigenvalue problem (2.14), with p = un. If

the process is quasi-stationary we can consider L to be independant of n. Then we have

rn+j = (I- L)-ir '*

and

p--1

un+v = un + _'_(I- L)-Jr n.

j;----O

If the eigenvalues ),i, of L axe negative the eigenvalues tci, of (I - L) -1 satisfy [xi] < 1 and

lim u"+v = u" + (I - (I - L)-l)-lr" = u" + (I - L-1)r ". (4.6)

10_co

13
Instead of taking a large number of time steps we can take one large step, which we call an

extrapolation step. We put

u = un + fie, (4.7)

wheree isthesolutionoftheequation

Le = (L - I)r", (4.8)

and fl is a stabilizing parameter. We choose/3 in such a way that H(u n +/ge) has no component

in the direction of e, i.e.

(H(u"+ e)= 0,

where (., .) denotes the usual inner product. There are other possible choices, for example

choose fl such that

IIS( ,"+ fle)ll= IIn(, "+ fle)ll.

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 un+l 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

continous case. See table (1). In figures (6),(7) the first few eigenvectors are plotted. Note that

in 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) -1 has only one eigenvalue,

xx, close to zero. All other eigenvalues are small. Therefore, when we have reached the quasistationary state, rn is in the direction of the eigenvector corresponding to xl. See figure (8).

Therefore we do not need to solve (4.8), and instead of (4.7) we use

u = un + fir n. (4.9)

In figures (2),(4) we have plotted u at different time stages to show how the convergence is

accelerated.

14
5. A twodimenslonal case. Consider the followingproblem

1 k _ _

0<x<l, 0<y<l, t>0,

.....

u(O,y,t)=a, u(1,v,t)=-a, a>0, (5.1)

u(x,0, t) = u(x,1, t) = w(x),

_(_, v, 0)= a(_, v),

where W(x) is the solution of the one dimensional problem (1.3) with b = -a, and f(x) -- O.

See (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 problem

ItW+ (w_)z = e(Wxx+ Wry), W= 0 on the boundary. (5.2)

We can solve (5.2) by separation of variables. Let _(x, y) = X(x)Y(y). Then

(wX)' - eX" = XX, X(O) = X(1) = O, (5.3a)

Y" = -qY, Y(0) = Y(1) = 0, (5.3b)

with/z = ), - eq. We recognize (5.3a) as (3.2). Therefore -),t = O(e-t/*) and -),j > O(1]e),

j = 2, 3, .... We can solve (5.3b). The solution is

Yi(Y) =sin(flry), qi = (Jzr)2, J= 1,2 ....

There is a whole sequence of eigenvalues,/qi, of order O(_). The eigenfunctions corresponding

to this sequence, Wli, will be exponentially small away from the shock. All other eigenvalues

will be of order O(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(e).

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 2. More

details will be given in another paper.

15
REFERENCES

[1] M. D. Salas, S. Abarbanel, D. Gottlieb, Multiple steady states for characteristic initial value

problems, Icase rcport No 84-57 , NASA CR-172486, November 1984.

[2] M. Hafez, E. Parlette, M. Salas, Convergence acceleration of iterative 8olutions for transonic

flow computations, AIAA 85-1641.

[3] J. D. Cole, J. Kevorkian, Perturbation methods in Applied Mathematics, Springer 1981.

[4] E. Landau, Einige Ungleichungen _r zweimal differenzierbare Funktionen, Proc. London

Math. Soc. 13(1913) 43-49.

16
Table I.

Eigenvalues of the eigenvalueproblem (3.2), y is the solution of (1.3). Three different cases were treated.

The discretization is done according to (4.5), with N = 100 gridpoints. The eigem:alues were found using

}averse iteration. E_genvectors corresponding to case (1) are plotted in figure (6a, b).

AI A2 A 3

f(_=)_(2_)/2

a ----b ----0 -8.64- 10-s -4.34 -5.32

= 0.04

:(_=)_n(2_)/2

a = b = 0 -4.62.10 -6 -5.617 -5.622

= 0.02

f(x)=o

a = 1, b -- -1 -1.24.10 -° -12.8 -13.5

e = 0.02

17
LJ)

•

o

o

•

o

LJ)

•

o

I

o

.....-+---------r----------r---..:::::~!!!!!!I---....,r_--------=:::::..,

b.oo 0.25 0.50 0.15 1.00

Figure 1. Convergence in time without convergence acceleration. Numerical solutions at

different time stages for the case e = 0.05, f == 0, a = 1, b = -1, u(x, 0) = 1+ 2(e-2Z -1)/(1- e-2 ).

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.

18
o

.:.. "'~---IIIIIl:;~~::::~-=::::::::~

lJ)

•

o

o

•

o

lJ)

•

o

I

o

· -I--------,r-------r---=~~~~.......---....=::::::::.... --b.oo 0.25 0.50 0.15 1.00

Figure 2. Convergence in time with convergence acceleration. Numerical solutions at different

time stages for the same case as in figure 1. Between each curve there are 15 time steps and one

extrapolation step. The same time step, k=O.2, and number of grid points, N=50, are used.

19
o

• .....

LJ)

•

o

o

•

o

LJ)

•

o

I

o

.--t-----------r-------------,r----------""T-----------, .....

b.oo 0.25 0.50 0.15 i.OO

Figure 3. Convergence in time without convergence acceleration. Numerical solutions at

different time stages for the case e = 0.04, f = fsin(1rx)cos(1rx), a = b = 0, u(x,O) = ksin(1rx).

Between each curve there are 100 time steps. The calculation is made with time step k = 0.1 and N=50

grid points.

20
Q

LD

O

E:)

Q

LD

O

I

O

-- i i I i

{3.00 0.25 0 °50 0 .-f5 t .00

Figure 4. Convergence in time with convergence acceleration. Numerical solutions at different

time stages for the same case as in figure 3. Between each curve there are 20 time steps and one

extrapolation step. The sume time step, k=0.1, and number of grid points, N--50, are used.

21
O

O_

° 1

O

O

O

0

O,

• I I I

_.00 0.2S O.SO O._S 4.00

Figure 5. Convergence when the shock is located at the boundary. Here _ = 0.04, f(x) =

sin(rfx), u(x, 0) = 1⁄2sin(zrx), N = 50, k ----0.1. Between each curve there are 5 time steps.

22
0.25 0.50 0.15 1.00

lJ)

rN

•

o

o

o

O~~ _

•

o

o

lJ)

lJ)

•

o

I

0.00

lJ)

rN

•

o

I

Figure 6a. Eigenvectors. The first two eigenfunctions of problem (3.2), when y, the solution of (1.3),

has a shock in the interior. In this case e =0.04, f(x) =f sin(i'l"X) cos(lTx), a =b=0, N =100.

23
LJ)

rN

•

a

a

a

a

•

a

LJ)

rN

•

a

I

a

LJ)

LJ)

•

a

I

0.00 0.25 0.50 0.15 1.00

Figure 6b. Eigenvectors. The third and fourth eigenfunctions of problem (3.2), when y, the solution

of (1.3), has a shock in the interior. In ~his case g =0.04, j(x) = sine1rx) cos(1rx), a =b=0, N =

100.

24
o

o

o

•

o

l()

N

•

o

I

o

l()

N

•

o

I

l()

rN')

•

o

I

o

ol()-+ .- -. --.~------_____.

•

CfJ .00 0.25 0.50 0.15 1.00

Figure 1. Eigenvectors. The first two eigenvectors, <PI and <Pz, of problem (3.2), when y, the solution

of (1.3), has a shock x = 1. In this case e = 0.08, f(x) = f sin(7l"x), a = b = 0, N = 100.

25
o

o

qT------=::::::::::::::-------~~~~~!""'"' o

Lf)

N

--o

•

o

I

o

Lf)

N

o

•

o

I

Lf)

rN)

o

•

o

I

0.25 0.50 0.15 i.OO

o

o

~-1--------.------------r-------.,_------____, •

cB .00

I

Figure 8. Differences between consecutive solutions at different time stages, when e = 0.04,

f = fsin(rrx)cos(rrx), a = b = 0, u(x, 0) = tsin(rrx). Between each curve there are 100 time steps.

The calculation is made with time step k = 0.1 and N=50 grid points.

26
O

LD

O

O

I

°O •0 0 •S _.0

Figure 9. The solution of (1.2) when f = 0, a = 1, b = 0 and e = 0.05.

27
O

c;-

o

C3--

LN

c;-

i

o

lOo0 0 o5 I .0

Figure 10. The solution of (1.2) when f _--0, a -- 1, b - -1 and e = 0.05.

28

1. Reporl No. NASA CR-178017 I2. GovtrnmMl Ace_,an No. 3. RtciPtenl's UUlog No.

ICASE Renort No. 85-50

4. Tille .nd Sublille 5. R.port 0..I.

CONVERGENCE TO STEADY STATE OF SOLUTIONS December 1985

OF BURGERS' EQUATION 6. Ptrlorm."9 Orpanllllion Code

7. AUlhorhl B. Performing Orpanililion Report No.

Gunilla Kreiss and Heinz-Otto Kreiss 85-50

10. WorkUnil No.

9. Performing Or~nizalion Name .nd Addrtll

Institute for Computer Applications in Science

and Engineering 11. Contracl or Granl No.

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

Hampton, VA 23665-5225 13. Typt 01 Repon .nd Period Covtred

12. Sponsoring Agency Name and AddrelS Contractor Report

National Aeronautics and Space Administration

Washington, D.C. 20546 14. Sponsori"9 Agency Code

505-31-83-01

15. Supplementarv NOleS

Langley Technical Monitor: Submitted to SIAM J. Appl. Math.

J. C. South Jr.

Final Report

16. AbstraCI

Consider the initial-boundary 'va1ue problem for Burgers' equation. It is

shown that its solutions converge, in time, to a unique steady state. The

speed of the convergence depends on the boundary conditions and can be

exponentially slow. Methods to speed up the rate of convergence are also

discussed.

17. Key Words ISuggesled by Authorhl I 1B.' Distribution Slatement

singular perturbation 64 - Numerical Analysis

eigenvalue problem

convergence acce1~ration

Unclassified - Unlimited

19. SecuritY Oassi!. (01 this reportl 20. Security Classif. lofthis page) 21. No..of Pages 22. Price

Unclassified Unclassified 30 A03

For sale by the National TechnicallnfOlmation Service, Springfield. Virginia 22161

NASA-Langley, 1985