- Email: [email protected]
An approximate inertial manifold for computing Burgers' equation
Physica D 60 (1992) 175-184 North-Holland
An approximate inertial manifold for computing Burgers' equation L.G. Margolin a and D.A. Jones b a Los Alamos National Laboratory, Los Alamos, NM 87545, USA b University of California, Irvine, California 92717, USA
We present a numerical scheme for the approximation of nonlinear evolution equations over large time intervals. Our algorithm is motivated from the dynamical systems point of view. In particular, we adapt the methodology of approximate inertial manifolds to a t'mite difference scheme. This leads to a differential treatment in which the higher (i.e. unresolved) modes are expressed in terms of the lower modes. As a particular example we derive an approximate inertial manifold for Burgers' equation and develop a numerical algorithm suitable for computing We perform a parameter study in which we compare the accuracy of a standard scheme with our modified scheme. For all values of the parameters (which are the coefficient of viscosity and the cell size), we obtain a decrease in the numerical error by at least a factor of 2.0 with the modified scheme. The decrease in error is substantially greater over large regions of the parameter space.
1. Introduction
The methodology of approximate inertial manifolds (AIMs) generates numerical schemes for subgrid scale modeling based on the theory of dynamical systems. The fundamental idea of the AIM as applied to numerical schemes is that one can express the "higher" order modes in terms of the "lower" order modes. By higher order modes we mean, in general, the scales that are not resolved by the numerical discretization. The mechanism for the coupling between higher order and lower order modes, both mathematicaUy and physically, is the nonlinearity of the dynamical equation. To take advantage of this dependence numerically, we shall have to make some additional approximations. The crucial approximation is an algebraic closure, analogous to the closure made in algebraic stress models for turbulence. In effect we assume that the higher order modes are in a steady state and that their time derivatives can be ignored. This leads to an algebraic relation between the higher and lower modes that is especially well-suited
for numerical computation. In the case Fourier series or other spectral expansions are used, there is a natural choice for the lower and higher modes. One then looks for a global function that expresses, asymptotically in time, the higher Fourier modes in terms of the lower modes. The graphs in phase space of such functions are in general smooth finitedimensional manifolds that may be interpreted as approximating the global attractor. AIMs for this case have been constructed and studied for a variety of dynamical equations which include the 2D Navier-Stokes equations, the Kuramoto-Sivashinsky equation, and reactiondiffusion equations [2-5,8,9,11,14]. The purpose of this paper is to show how the ideas of AIMs can be applied to finite difference schemes (cf. [1,7,10,12,13]). The basic ideas are similar to those used for spectral algorithms. However, the development is complicated when one tries to be precise about defining the higher and lower order modes. We shall work with Burgers' equation to illustrate the AIM technique for finite differences.
0167-2789/92/$ 05.00 ~ 1992-Elsevier Science Publishers B.V. All rights reserved
176
L.G. Margolin, D.A. Jones / Inertial manifold for computing Burgers" equation
Ou Ou .OEu 0"~ + u-ff-~ - A-ff-~ = f ( x ).
(1.1)
Here u is the fluid velocity. 2 > 0 is the coefficient of viscosity, and f is a specified force. This simple equation contains the two necessary ingredients to apply AIMs, a nonlinearity in the convective term and a dissipative term. Extension of our ideas to other dynamical equations should be straightforward. Conceptually we begin by comparing an Nth order numerical calculation with a 2Nth order calculation. By Nth order we mean keeping N Fourier modes in a spectral calculation, or using N finite difference ceils. In a 2Nth spectral calculation, the modes {k = N + 1. . . . . 2N} are readily identified as the higher order modes; however, the evolution of the lower order modes {k = 1 , . . . , N } depends on the higher order modes through the nonlinearity. By finding an approximate relation between the higher order modes and the lower order modes, we can represent the effect of the higher order modes on the evolution of the lower order modes without resolving them - e.g., in an Nth order calculation. In a 2Nth order finite difference calculation, the 2N cells are essentially equivalent, and so the higher order modes cannot be simply identified as a subset of the cells. We consider instead an alternate basis in which the higher order modes are readily identified. W e then derive the evolution equations for the coefficients of this new basis set. In section 3 we derive the approximate inertial manifold that describes the slaving of the higher order modes to the lower order modes for our finite difference algorithm. In this section we make precise the rules by which we estimate the size of the various terms the appear. Two constants appear in the AIM that cannot be determined from the equation, but must be estimated from the boundary conditions. In section 4 we briefly discuss the issue of boundary conditions. As previously mentioned, approximation in the evolution equation for the
higher order modes leads to an algebraic relation that requires no additional boundary conditions. However, coefficients of the higher order modes are required in the boundary cells outside the mesh for the evolution of the lower order modes. We give a derivation of these terms, which are easily implemented. In section 5 we compare calculations with and without the additional terms derived from the AIM. We restrict ourselves to a steady state problem for which an analytic solution exists. When the convective term is approximated using centered differences (second-order accurate in space), the AIM method yields an improvement in accuracy of at least a factor of 2.0 over the original method, for most choices of the physical viscosity parameter and of the number of zones. The increase in accuracy is significantly greater in large regions of the parameter space. The increase in the amount of computation is negligible. The key to applying AIMs to finite differences is the construction of the higher order basis functions. We note that our choice is not unique. Also, while our particular choice determines the relation between the coefficients of the higherorder and the lower-order modes, the final algorithm (for the lower-order modes) is otherwise independent of that choice. Thus we are led to ask how sensitive our final algorithm is to the choice of the higher-order basis functions. In section 6 we describe an alternate choice for these higher-order basis functions. As we remarked earlier our results in this paper are not completely rigorous. More specifically, we will assume that we may neglect the time derivatives of the higher-order terms. In the spectral case this assumption has a rigorous mathematical foundation (see for example [4 ] ). We suppose that the physical basis for this assumption remains valid in the case of finite differences.
L.G. Margolin, D.A. Jones / Inertial manifoM for computing Burgers" equation
2. Finite difference basis sets
The calculation of the slaving of the higher modes to the lower modes for a finite difference method is somewhat more complicated than for a spectral method. In the spectral method of order N, the basis eigenvectors are the first N Fourier functions. When we go to a higherorder method, say 2N, all we do is to add the next N Fourier functions, i.e. of wave numbers N + 1 through 2N. Thus the higher order calculation for the spectral method has the following property: the first N basis eigenvectors do not change. Those first N basis elements constitute the lower order modes and it is important that they are a subset of basis elements for the 2N calculation. This does not happen naturally for finite difference methods. For example, if the (one dimensional) finite difference method is of order N, this means that we have divided the line segment representing the interval L into N disjoint segments, each of length L/N. If we want a higher-order method, say 2N, we would most naturally divide the interval into 2N disjoint segments - a refined mesh in which each cell has length L/2N. Note that none of the new (2N) basis elements belongs to the original (N) basis set. We get around this difficulty by defining a less natural basis set for order 2N. In the sense that the original N basis elements represent average values of some function within the cell, we can choose the extra N basis elements to represent first derivatives of the function within the cells. This means that when we go from order N to order 2N, the original basis elements remain unchanged. By doing this, we introduce a new difficulty. How do we define the difference scheme for the new "first derivative" basis elements? To address this problem, we note that we do know how to difference the equation for the order 2N calculation on the refined mesh, where all the cells are of length L / 2 N (the equality of the cell size is not crucial). In fact, it is the same difference scheme that we use for the order N
177
calculation. To use this information, we need to construct an approximate map from one basis set to the other. The following paragraphs describe this process. Suppose again that we are solving Burgers' eq. (1.1) on the segment [0,L]. We consider two sets of basis functions, each of 2N elements. These are shown graphically in fig. 2.1 and fig. 2.2. The basis set in fig. 2.1 consists of N elements which are the unit functions of a regular Nth-order calculation {hk}, plus N more elements, {Sk }. The element Sk is defined as minus one over the first half of the cell, and plus one over the second half of the cell. If we approximate an arbitrary function f (x) in this basis set, the coefficient of Sk is L
-L-/I,x>sk x 0
[s =
+
f
Of (kAx) Ax Ox 4 + O(Ax3)"
(2.1)
When we identify Sk as a higher order mode, we mean precisely that its coefficient is of O(Ax) whereas the coefficient of the hk modes, ak, is of O(1). The second basis set, shown in fig. 2.2, consists of 2N equivalent elements {Hi, Jk }. The elements labeled H k a r e one for the left half of the cell k (and zero everywhere else), whereas the elements labeled Jk are one for the right half of the cell k, We define the inner product notation (bra[ket} by L
(n , , nk ) -- f nk n,, dx.
(2.2)
0
With this inner product, the basis elements {hk, Sk} form an orthogonal set
(hk, hk,) = Ax Jkk,,
L.G. Margolin, D.A. Jones / Inertial manifold for computing Burgers" equation
178
hk.m
N
Sk
lu) ___~ ((,kl/-/k) + flklJr,,)). I
l~x
I
k-1
k+l
k-I
k
k+l
Fig. 2.1. The basis set {hk,Sk}.
Yk-I
I
l~X
I
I
~X
k-1
k
k+l
k-1
k
k+l
Fig. 2.2. The basis set {Ilk, Jk}.
(Sk, Sk') ~" ~
¢~kk',
(hk,Sk,) -~ O,
(2.8)
k=l
Note that eqs. (2.7) and (2.8) are two different approximations, and so are not necessarily consistent. We can use these approximations to generate two sets of approximate relations between the sets of expansion coefficients {a, b} and {a, fl}. By taking the inner product of eq. (2.7) with (HkJ, we find an expression in the coefficients {a, b}. To the extent that eq. (2.7) approximates eq. (2.8), we can equate the inner products to derive an approximate relation between the sets of coefficients,
(Hk, u) = ½Ax C~k ~-- ½Ax (ak -- bk).
(2.9)
(2.3)
as do ( Hk, Jk } ;
We next take the inner products of both equations with (Jkl and so derive
(H~,~/k,) = ½ax6~k,,
(Jk, U) = ½AXfk ~ ½Ax(ak + bk).
(Jk, Jk,) = ½~x ~kk,, (/'tk, Jk,) = O.
Thus we have the approximate relations (2.4) ~k ~ ak -- bk,
It is easy to find the inner products between elements of the two sets
(h,,,/-tk,) = ½ax '~kk', (ha, I~,) = ½aX,~kk,
(2.5)
and
(Sk, Hk, ) = --½~kX tJkk, ,
(Sk, A') = ½aX'~kk'.
(2.6)
Suppose now that we have some function u (x, t) that we want to approximate in our discrete bases. In the {h, s} basis we write N
lu) • ~_~ (aklhk) -4" bklSk) ). k~l
In the {H, J} basis, we write
(2.10)
(2.7)
flk ~ ak + bk.
(2.11)
Similarly, if we take the inner product of eqs. (2.7) and (2.8) with the vectors (hkl and (Ski, we derive the inverse relations
ak~_ ½(ak + flk),
bk~-- ½(flk--C~k).
(2.12)
Note that eq. (2.11) is an exact inverse to eq. (2.12). This is because of the way that we subdivided the line for the fine mesh, and in particular because
hk = Hk + Jk.
(2.13)
Other choices of the relation between the fine and coarse meshes would have led to eqs. (2.11 ). However, they would only be an approximate inverse to eqs. (2.12). The goal of this section is to write the evolution equations for the expansion coefficients
L.G. Margolin, D.A. Jones / Inertial manifold for computing Burgers" equation
{a, b}. We begin by writingthe {a,i} equations in the {Hk, Jk} basis: Oak kffil
2
IHk) \ Ot +
+
ISk> k=l
,__~ +
2Ax
(
?
Ax "--i
Oak . i ~ - - i 2 - 1 ~
]
421k+ik-I
Ax 2
= (fgB _ f~ ).
)
2ak (2.17)
2ak _ f ~ )
+
\or
2 --ak 2 42ak+l + a k _ m i k ~ ak+1
(Ocg~.~k
-- " - ~ "P
2~
+ tk-IAx 2-
_42ik
2
i k -- t k - I
179
Finally, we use eq. (2.11) to replace { a , t } by {a, b}. Further we make the approximations, as in eq. (2.1),
a / + 1 --
_42ak+i + ak -- 2ik _ fk#'~ Ax E
)
(2.14)
-- 0,
where fk%fk# are the coefficients of f ( x ) with respect to the IHk), IJk) basis respectively. In eq. (2.14) we show the semidiscrete form, in which the spatial derivatives are discretized, but the time derivatives are shown in analytic form. We are not concerned with the semidiscrete equations, but the time differencing is not relevant to our discussion. Furthermore, we have chosen a particular discretization for the spatial terms, both of which are centered in space. We will discuss the boundary conditions in detail in section 4. We convert the equations to the {h,s} basis in two steps. First, we change the basis elements themselves, by inserting a complete set of states
fka + fk'~ ~ 2 f (kAx) := 2j~, Ax d f (kAx)
dx
:= ½ fk'.
For the hk equation, Oak Ot
22ak+1 + a k - l - 2 a k - bk+l + bk-1 ~kX2 2
+
2
ak+ 1--ak_ 1 1 4A.X + ~(--2ak+lbk+l
+ 4akbk
--2ak-lbk-, +b~+, + b~_, - 2b~) = 3~,
(2.18)
and for the Sk equation Obk Ot
22ak+l -- ak-1 -- 6bk -- bk+l - bk-1 Ax 2 a2+l + a ~ - l - 2a2 1 + 4AX + -~--~(--2ak+lbk+l + 2 a k - l b k - , + b~+ 1 + b~_ 1 - 2 b ~ )
N
(2.15)
I = ~_~(lhk)(hk[ + [Sk)(Skl). I
Noting that each of the hk and Sk is independent, we derive two sets of equations, valid for each k. The hk equation, still in terms of {a, t } is + ak+l--ak
42ak+1
2Ax
(Oak
+ \ Og "1"
~
= (fk# + fk~). The Sk equation is
"k--t2-1
+ak-mil k
Ax 2 + t k _ l - - 2ak'~ 42 t k
ax2
}
(2.16)
1 = ~Ax f~.t
(2.19)
3. The approximate inertial manifold There is an important qualitative difference between eqs. (2.18) and (2.19) on the one hand, and their spectral analogs on the other. The difference is that in the spectral case, the linear term contains only terms in the lower wavenumbers. However in eq. (2.18 ), the linear term contains both ak and bk. The reason is that the coefficients ak and bk are not really independent in a
180
L.G. Margolin, D.A. Jones / Inertial manifold for computing Burgers" equation
"truncation error" sense. In fact, from eq. (2.12 ) we see Ou A x bk = ½ i l k - - O t k -- 19X 4 + O(Ax3)"
t~k ~-- -~
2/
kAx
u(x)dx
(k-1/2)Ax
(3.1) that
However, we can also write O~k ----u ( k~kx ) -- ux ( k A x ) ~ J x 1 ~(ak+l--ak-1)
O__.~uAx + O(Ax3)"
= OX 4
(3.2)
+ Uxx(kAx)~(Jx)
2
- Uxxx(kAx)13~,2(Jx)3 +
Since these two expressions can only be distinguished in the third order, and since the linear term is of the second order, it is not surprising that both expressions appear in the linear term ofeq. (2.18). In the next step we want to use eq. (2.19) to find an approximate relation between ak and bk that is based on the slaving of modes due to the nonlinearity. In particular, we have to be careful in estimating of the various terms. Combining eqs. (3.1) and (3.2) we have
O(Jx4).
Expanding the other terms in (3.4) in a similar fashion, we have that eq. (3.4) is approximately 22 / 1 0 3 u 3 :~.,,O(Axs)) ax2 ax If we think of bk as a continuous function, we have, using eq. (3.3), for the linear term in (2.19) _22ak+l -- a k - I -- bk+l - b k - i - 6bk /~X 2
bk ~ ak + l -- a k - I 8 -
cgU A x Ox 4 + O(Ax3)"
This relation is derived from the truncation analysis and in not due to the slaving of the higher order modes. It seems appropriate to use truncation analysis to estimate the magnitude of the other terms in the equations. That is, we will use Taylor expansions, and compare terms using the power of the cell size Ax that appears in the expansion. An approximate relation between ak and bk is derived by making some approximations in eq. (2.19). To begin, we ignore the time derivative of bk. This approximation is also made in the spectral case. For the linear term we go back to the definitions of a and b in terms of a and il. We obtain 22
AX 2 (Otk+l -- 3ilk + 3O~k -- i l k - l ) .
_O2b A-~.
(3.3)
(3.5)
The viscous term must be balanced by the nonlinear term. W e therefore have
02b ~-e f~ + ~--~(a2+l + a2_l - 2ak)
OX 2
2e A x 2 (ak+lbk+l -- a k - l b k - l ),
where e := Ax/42. Also, a similar analysis as 2 1 + b2_l - 2b 2 ,,~ (Ax4). above shows that bk+ This term is dearly small compared to the other terms in eq. (2.19) Thus, we have neglected this term in eq. (3.6). We continue by substituting the lowest order result of eq. (3.3) on the fight side of eq. (3.6) bk = ~(ak+l -- a k - l ) .
(3.4)
This appears to be a third derivative. Indeed, if we expand a function u about the spatial point k, we get, for instance, from
(3.6)
Using the approximations ~(ak+2 -- ak) ~-- ¼(ak-1 -- ak), ~(ak--ak-1)
~ ~(ak--ak-l),
L.G. Margolin, D.A. Jones / Inertial manifold for computing Burgers" equation
that the computational results are somewhat sensitive to the boundary conditions on b, and these choices are not appropriate. To find a more appropriate choice, we return to the partial differential eq. (1.1). On the left boundary, u = Uleft and so is constant in time. Thus on the left boundary
ak ~-- ½(ak-I q- ak+l),
we find ak+ lbk + 1 -- ak-lbk-1 1 2 I ~- ~(ak+
+
a2_i
--
181
2a~).
It follows that 02b e 02a 2 Ox 2 ~- 2 cox 2
OOU .OO2U Uleft~ - = /t~-~ + f ( 0 ) . eft.
(3.7)
The most general solution of eq. (3.7) has the form bk = ½e(a 2 + cl + c2x) - e F k ,
(3.8)
where Fk is F (kAx) and d F / d x = f . The conslants ct, c2 will be determined from the boundary conditions. This will be done in the next section. We note that it is no accident that both sides of eq. (3.7) are composed of second derivatives. eq. (2.19) is a differenced form of Burgers' equation after one additional differentiation in space, eq. (3.7) represents the steady state of this equation. Note however, that deriving the AIM (3.8) and the associated boundary conditions, we never made use of a knowledge of the analytic solution to Burgers' equation.
4.
(4.1)
Implementation
In this section we briefly discuss two aspects of the implementation of the AIM. The basic strategy is to compute with eq. (2.18). We use the approximate relation of eq. (3.8) to replace the bks where they appear. Now suppose that we specify the values of u on the left and the right boundaries. For an Nth order calculation, we might set al = Ule fl and a~v+1 = Uright, and then compute the interior values ak {k = 2 . . . . . N}. We note that to calculate a2 we will need bl. In eq. (3.8) these are not naturally defined. One could simply set b~ = blv+l = 0. However, it turns out
Returning to eq. (3.8), we can estimate the derivatives au 4bl - - - + 0(~), Ox Ax 02u 32 = ~ - - ~ (bE - b l ) Ox 2
+ O(f).
(4.2)
Note that the first and second derivatives in eq. (4.2) are centered at different spatial points, and neither are exactly at the left boundary. We neglect these corrections, and also terms of order e. Substituting (4.2) in (4.1) leads to
A• 5
leftbl = ( b 2 - b l ) +
5Ax2 ~-~ f ( 0 ) .
(4.3)
Finally, we note that the coefficient of the lefthand-side is proportional to e; further, the coefficient in front of f (0) is proportional to eAx. Thus for consistency these terms should be ignored, and on the boundary we have bl = b2,
bN+l = bN.
(4.4)
Next we approximate the constants cl and c2 defined in eq. (3.8). On the left boundary we have 8b12 Ax ~ u~eft + cl - 2 F ( 0 ) ,
(4.5)
where we have made the arbitrary choice x = 0 on the left boundary. On the right boundary 8bN+12 2 AX ~ Udght + Cl + c2L - 2 F ( L ) .
(4.6)
182
L.G. Margolin, D.A. Jones / Inertial manifoldfor computing Burgers"equation
We cannot apply the boundary conditions (4.4) since bl and bN+l a r e not yet known. Instead, we use the approximation of eq. (3.3) for b2 and bN and then solve eq. (4.5) and (4.6) for the constants c~ and c2. We mention one more point relating to the boundary conditions (4.4). Corrections due to the AIM in eq. (2.18) appear both in the dissipative term and in the nonlinear term. In the nonlinear term, the new terms have the flavor of higher order corrections. For example, if we set bk = 0, we regain the original differencing of this term. This is not true in the dissipative term. The first difference of bks is the some order as the second difference of aks. If we set bk to zero in this term, we do not reproduce the original differencing of this term. Thus in the dissipation, the bk terms do not represent higher order corrections. Now in deriving the boundary conditions (4.4), we have kept only lowest order terms. These are sufficiently accurate for the higher order corrections in the nonlinear term, but not for the new formulation of the dissipative term. One could try to find more accurate boundary conditions. However (4.4) is particularly simple to implement. Furthermore, we expect most of the benefit of the AIM to result from the nonlinear term. Hence, we replace the new form of the dissipative term with the old. The results in the next section were all computed with eq. (4.7), and using the boundary conditions (4.4) 2 _ ~2_1 Oak 2ak+t + ak-l -- 2ak ak+l Ot Ax 2 + 4Ax + (--2ak+lbk+l + 4akbk - 2ak-lbk-1
+b2+l + b2_1 - 2b 2) (4Ax) -l = fk.
(4.7)
5. Numerical results
In this section we present numerical results to illustrate the efficacy of the AIM as applied to
finite difference schemes. First we describe our model problem, and the measure we use for comparison. Then we compare the accuracy of calculations using our new scheme eq. (4.7) with calculations using a standard finite difference scheme, which is (4.7) with bk set to zero. We solve Burgers' eq. (1.1). We set Uleft = 10.0 and uright = 0.0 for boundary conditions. We choose an interval length L = 1.0 and set f = 0. We treat the total number of cells, N and hence Ax, and the viscosity 2, as parameters. In each calculation, we run the problem to steady state. There is a simple analytic solution for this problem, which we can use to calculate the accuracy of our solutions. In each cell we compare the computed value with the analytic solution U (x) evaluated at the cell center. We define the L 2 norm of the error by
L2(a, U) =
~_,[ak -- U(xk)122~
.(5.1)
k=2
We have also computed a n L 1 norm and an L °~ norm and find that all of our results are qualitatively unchanged. In each of the figures in this section, we plot a ratio, the L 2 norm of the error of the standard calculation divided by the L 2 norm of the error of the AIM calculation. Thus, a ratio greater than 1 means that the AIM calculation is more accurate in the global sense. In some cases this ratio may decrease as the number of cells is increased. However, we note that in every case the actual L 2 error of both the standard calculation and the AIM calculation monotonically decrease as the number of cells is increased. In figure 5.1 we show parameter studies where we hold fixed the viscous coefficient A and vary the number of cells. For 2 = 5.0, the ratio of errors is 3.3 (for N = 5), increases to a maxim u m value of 8.0 (for N = 15) and then slowly decreases. For larger values of A, the graphs are very similar, each showing a maximum at about N = 15. The shapes of the graphs tend to be-
L.G. Margolin, D.A. Jones/Inertial manifold for computing Burgers" equation
hk
183
Sk /
,
lamlxla=5.
k-I
lambda-20.
P
lambda=l.
o
== .
.
.
.
.
.
.
i
.
.
.
.
.
.
.
.
10
i
.
.
.
.
'.X
k+l
'
k-1
//k
'.X
k+l
Fig. 6.1. An alternate basis set for the {Sk}.
"6
.
,
k
.
.
100
.
.
1000
number of cells
Fig. 5. l. Ratio of errors in the L 2 norm for several values of ,l as the number of cells increases.
To summarize, the AIM leads to an improvement in accuracy of at least a factor of 2 essentially over the whole parameter space in cell size and viscosity. For large regions of the parameter space, the increase may be substantially greater.
100.
6. Alternate basis sets
0
~. 10
=
0
;
--
:
-
N =20
-
:
•
-
N=5
.
.
'3 0
1
.
.
.
.
.
. 110.
.
.
.
.
.
100
lambda
Fig. 5.2. Ratio of errors in the L 2 norm for several values of total cells as the viscosity increases.
come slightly narrower, and the maxima slightly higher as 2 increases. When 2 = 1.0, there is no maximum. All the graphs appear to asymptote to a value slightly greater than 3.0 as the number of cells increases. In figure 5.2 we show a complementary study in which we hold the number of cells fixed and plot the ratio of errors as the viscous coefficient A increases. Each of the graphs appears to reach a constant value as 2 increases beyond 3.0. For the plotted values, this constant value is largest at 7.8 when N = 20. Its actual maximum appears to be about 8.1 when N = 15 (not shown), which is consistent with the graphs shown in figure 5.1. For larger values of N, the constant value decreases slowly, being about 4.3 for N = 60 and about 3.9 for N = 100 (not shown).
The key to the derivation of the AIM is the basis set shown in figure 2.1. This particular choice led to approximate relation (3.8). However, given this relation, the form of the basis set plays no further role in the implementation. In this section we address the question of to what extent our results would change given a different basis set {s~}. To be specific, we choose the basis set shown in (6.1). Here the higher order basis functions are continuous lines. Note that both the original set {Sk} and the new set {s~} are antisymmetric with respect to the point k. If one follows the same procedure outlined above, one derives the simple result bk =
x/34AX (a~c + cl +
Ax
(6.1)
That is we derive the same relation except for an overall multiplicative constant. Moreover, the only essential ingredient in this derivation is the antisymmetry of the basis function Sk. Any antisymmetric basis function would lead to the same result, up to the overall multiplicative constant. The question naturally arises as to which is the best choice of basis function. The answer is that it depends on the particular function that must be represented (i.e the solution u). Among the various results of our parameter studies, the
184
L.G. Margolin, D.A. Jones / Inertial manifold for computingBurgers"equation
AIM of eq. (6.1) is some times more accurate than (3.8) and sometimes slightly less accurate. Other choices of the basis set {Sk} will lead to other multiplicative constants in (6.1). One finds experimentally that for any given problem, the most accurate result comes form choosing the overall multiplicative factor in a small range, about 1.0 ± 0.15.
Acknowledgements The authors appreciate the continued interest and advice of Edriss S. Titi and Darryl Holm. We gratefully acknowledge the support of the Institute for Geophysics and Planetary Physics (IGPP) and the Center for Nonlinear Studies (CNLS) at Los AJamos. This work was performed under the auspices of the U.S. Department of Energy at Los Alamos National Laboratory.
References [ 1] M. Chen and R. Temam, Incremental unknowns for solving partial differential equations, Num. Math., to appear.
[2] Ch. Dcvulder, M. Marion and E.S. Titi, On the rate of convergence of nonlinear Galerkin methods, submitted. [3] E. Fabcs, M. Luskin and G. Sell, J. Diff. Eq., 89 (1991) 355-387. [4] C. Foias, O.P. Manley and R. Temam, Math. Mod. and Num. Anal. M2AN 22 (1988) 93-I 14. [5] C. Foias, G. Sell and E.S. Titi,J. Dynam. Diff. Eq. l, No. 2, 199-243 (1989) [6] C. Foias and R. Temam, J.Math. Pures Appl., 58 (1979) 339-368. [7 ] C. Foias and E.S. Titi,Nonlinearity 4 (1991 ) 135-153. [8] M.S. Jolly,I.G. Kevrekidis and E.S. Titi,Physica D 44 (1990) 38-60. [9] M. Marion, J. Dynam. Diff. eqs. I (1989) 245-267. [10] M. Marion and R. Temam, 57 (1990) 205-226. [I I ] R. Temam, J.Fac. ofSci.,Univ. ofTokyo, IA, 36 (1989) 629-647. [12] R. Temam, SIAM J. Math. Anal. 21 (1990) 154-178. [13 ] R. Temam, Math. Comput., to appear. [14] E.S. Titi,J. Math. Anal. Appl. 149 (1990) 540-557.
An approximate inertial manifold for computing Burgers' equation L.G. Margolin a and D.A. Jones b a Los Alamos National Laboratory, Los Alamos, NM 87545, USA b University of California, Irvine, California 92717, USA
We present a numerical scheme for the approximation of nonlinear evolution equations over large time intervals. Our algorithm is motivated from the dynamical systems point of view. In particular, we adapt the methodology of approximate inertial manifolds to a t'mite difference scheme. This leads to a differential treatment in which the higher (i.e. unresolved) modes are expressed in terms of the lower modes. As a particular example we derive an approximate inertial manifold for Burgers' equation and develop a numerical algorithm suitable for computing We perform a parameter study in which we compare the accuracy of a standard scheme with our modified scheme. For all values of the parameters (which are the coefficient of viscosity and the cell size), we obtain a decrease in the numerical error by at least a factor of 2.0 with the modified scheme. The decrease in error is substantially greater over large regions of the parameter space.
1. Introduction
The methodology of approximate inertial manifolds (AIMs) generates numerical schemes for subgrid scale modeling based on the theory of dynamical systems. The fundamental idea of the AIM as applied to numerical schemes is that one can express the "higher" order modes in terms of the "lower" order modes. By higher order modes we mean, in general, the scales that are not resolved by the numerical discretization. The mechanism for the coupling between higher order and lower order modes, both mathematicaUy and physically, is the nonlinearity of the dynamical equation. To take advantage of this dependence numerically, we shall have to make some additional approximations. The crucial approximation is an algebraic closure, analogous to the closure made in algebraic stress models for turbulence. In effect we assume that the higher order modes are in a steady state and that their time derivatives can be ignored. This leads to an algebraic relation between the higher and lower modes that is especially well-suited
for numerical computation. In the case Fourier series or other spectral expansions are used, there is a natural choice for the lower and higher modes. One then looks for a global function that expresses, asymptotically in time, the higher Fourier modes in terms of the lower modes. The graphs in phase space of such functions are in general smooth finitedimensional manifolds that may be interpreted as approximating the global attractor. AIMs for this case have been constructed and studied for a variety of dynamical equations which include the 2D Navier-Stokes equations, the Kuramoto-Sivashinsky equation, and reactiondiffusion equations [2-5,8,9,11,14]. The purpose of this paper is to show how the ideas of AIMs can be applied to finite difference schemes (cf. [1,7,10,12,13]). The basic ideas are similar to those used for spectral algorithms. However, the development is complicated when one tries to be precise about defining the higher and lower order modes. We shall work with Burgers' equation to illustrate the AIM technique for finite differences.
0167-2789/92/$ 05.00 ~ 1992-Elsevier Science Publishers B.V. All rights reserved
176
L.G. Margolin, D.A. Jones / Inertial manifold for computing Burgers" equation
Ou Ou .OEu 0"~ + u-ff-~ - A-ff-~ = f ( x ).
(1.1)
Here u is the fluid velocity. 2 > 0 is the coefficient of viscosity, and f is a specified force. This simple equation contains the two necessary ingredients to apply AIMs, a nonlinearity in the convective term and a dissipative term. Extension of our ideas to other dynamical equations should be straightforward. Conceptually we begin by comparing an Nth order numerical calculation with a 2Nth order calculation. By Nth order we mean keeping N Fourier modes in a spectral calculation, or using N finite difference ceils. In a 2Nth spectral calculation, the modes {k = N + 1. . . . . 2N} are readily identified as the higher order modes; however, the evolution of the lower order modes {k = 1 , . . . , N } depends on the higher order modes through the nonlinearity. By finding an approximate relation between the higher order modes and the lower order modes, we can represent the effect of the higher order modes on the evolution of the lower order modes without resolving them - e.g., in an Nth order calculation. In a 2Nth order finite difference calculation, the 2N cells are essentially equivalent, and so the higher order modes cannot be simply identified as a subset of the cells. We consider instead an alternate basis in which the higher order modes are readily identified. W e then derive the evolution equations for the coefficients of this new basis set. In section 3 we derive the approximate inertial manifold that describes the slaving of the higher order modes to the lower order modes for our finite difference algorithm. In this section we make precise the rules by which we estimate the size of the various terms the appear. Two constants appear in the AIM that cannot be determined from the equation, but must be estimated from the boundary conditions. In section 4 we briefly discuss the issue of boundary conditions. As previously mentioned, approximation in the evolution equation for the
higher order modes leads to an algebraic relation that requires no additional boundary conditions. However, coefficients of the higher order modes are required in the boundary cells outside the mesh for the evolution of the lower order modes. We give a derivation of these terms, which are easily implemented. In section 5 we compare calculations with and without the additional terms derived from the AIM. We restrict ourselves to a steady state problem for which an analytic solution exists. When the convective term is approximated using centered differences (second-order accurate in space), the AIM method yields an improvement in accuracy of at least a factor of 2.0 over the original method, for most choices of the physical viscosity parameter and of the number of zones. The increase in accuracy is significantly greater in large regions of the parameter space. The increase in the amount of computation is negligible. The key to applying AIMs to finite differences is the construction of the higher order basis functions. We note that our choice is not unique. Also, while our particular choice determines the relation between the coefficients of the higherorder and the lower-order modes, the final algorithm (for the lower-order modes) is otherwise independent of that choice. Thus we are led to ask how sensitive our final algorithm is to the choice of the higher-order basis functions. In section 6 we describe an alternate choice for these higher-order basis functions. As we remarked earlier our results in this paper are not completely rigorous. More specifically, we will assume that we may neglect the time derivatives of the higher-order terms. In the spectral case this assumption has a rigorous mathematical foundation (see for example [4 ] ). We suppose that the physical basis for this assumption remains valid in the case of finite differences.
L.G. Margolin, D.A. Jones / Inertial manifoM for computing Burgers" equation
2. Finite difference basis sets
The calculation of the slaving of the higher modes to the lower modes for a finite difference method is somewhat more complicated than for a spectral method. In the spectral method of order N, the basis eigenvectors are the first N Fourier functions. When we go to a higherorder method, say 2N, all we do is to add the next N Fourier functions, i.e. of wave numbers N + 1 through 2N. Thus the higher order calculation for the spectral method has the following property: the first N basis eigenvectors do not change. Those first N basis elements constitute the lower order modes and it is important that they are a subset of basis elements for the 2N calculation. This does not happen naturally for finite difference methods. For example, if the (one dimensional) finite difference method is of order N, this means that we have divided the line segment representing the interval L into N disjoint segments, each of length L/N. If we want a higher-order method, say 2N, we would most naturally divide the interval into 2N disjoint segments - a refined mesh in which each cell has length L/2N. Note that none of the new (2N) basis elements belongs to the original (N) basis set. We get around this difficulty by defining a less natural basis set for order 2N. In the sense that the original N basis elements represent average values of some function within the cell, we can choose the extra N basis elements to represent first derivatives of the function within the cells. This means that when we go from order N to order 2N, the original basis elements remain unchanged. By doing this, we introduce a new difficulty. How do we define the difference scheme for the new "first derivative" basis elements? To address this problem, we note that we do know how to difference the equation for the order 2N calculation on the refined mesh, where all the cells are of length L / 2 N (the equality of the cell size is not crucial). In fact, it is the same difference scheme that we use for the order N
177
calculation. To use this information, we need to construct an approximate map from one basis set to the other. The following paragraphs describe this process. Suppose again that we are solving Burgers' eq. (1.1) on the segment [0,L]. We consider two sets of basis functions, each of 2N elements. These are shown graphically in fig. 2.1 and fig. 2.2. The basis set in fig. 2.1 consists of N elements which are the unit functions of a regular Nth-order calculation {hk}, plus N more elements, {Sk }. The element Sk is defined as minus one over the first half of the cell, and plus one over the second half of the cell. If we approximate an arbitrary function f (x) in this basis set, the coefficient of Sk is L
-L-/I,x>sk x 0
[s =
+
f
Of (kAx) Ax Ox 4 + O(Ax3)"
(2.1)
When we identify Sk as a higher order mode, we mean precisely that its coefficient is of O(Ax) whereas the coefficient of the hk modes, ak, is of O(1). The second basis set, shown in fig. 2.2, consists of 2N equivalent elements {Hi, Jk }. The elements labeled H k a r e one for the left half of the cell k (and zero everywhere else), whereas the elements labeled Jk are one for the right half of the cell k, We define the inner product notation (bra[ket} by L
(n , , nk ) -- f nk n,, dx.
(2.2)
0
With this inner product, the basis elements {hk, Sk} form an orthogonal set
(hk, hk,) = Ax Jkk,,
L.G. Margolin, D.A. Jones / Inertial manifold for computing Burgers" equation
178
hk.m
N
Sk
lu) ___~ ((,kl/-/k) + flklJr,,)). I
l~x
I
k-1
k+l
k-I
k
k+l
Fig. 2.1. The basis set {hk,Sk}.
Yk-I
I
l~X
I
I
~X
k-1
k
k+l
k-1
k
k+l
Fig. 2.2. The basis set {Ilk, Jk}.
(Sk, Sk') ~" ~
¢~kk',
(hk,Sk,) -~ O,
(2.8)
k=l
Note that eqs. (2.7) and (2.8) are two different approximations, and so are not necessarily consistent. We can use these approximations to generate two sets of approximate relations between the sets of expansion coefficients {a, b} and {a, fl}. By taking the inner product of eq. (2.7) with (HkJ, we find an expression in the coefficients {a, b}. To the extent that eq. (2.7) approximates eq. (2.8), we can equate the inner products to derive an approximate relation between the sets of coefficients,
(Hk, u) = ½Ax C~k ~-- ½Ax (ak -- bk).
(2.9)
(2.3)
as do ( Hk, Jk } ;
We next take the inner products of both equations with (Jkl and so derive
(H~,~/k,) = ½ax6~k,,
(Jk, U) = ½AXfk ~ ½Ax(ak + bk).
(Jk, Jk,) = ½~x ~kk,, (/'tk, Jk,) = O.
Thus we have the approximate relations (2.4) ~k ~ ak -- bk,
It is easy to find the inner products between elements of the two sets
(h,,,/-tk,) = ½ax '~kk', (ha, I~,) = ½aX,~kk,
(2.5)
and
(Sk, Hk, ) = --½~kX tJkk, ,
(Sk, A') = ½aX'~kk'.
(2.6)
Suppose now that we have some function u (x, t) that we want to approximate in our discrete bases. In the {h, s} basis we write N
lu) • ~_~ (aklhk) -4" bklSk) ). k~l
In the {H, J} basis, we write
(2.10)
(2.7)
flk ~ ak + bk.
(2.11)
Similarly, if we take the inner product of eqs. (2.7) and (2.8) with the vectors (hkl and (Ski, we derive the inverse relations
ak~_ ½(ak + flk),
bk~-- ½(flk--C~k).
(2.12)
Note that eq. (2.11) is an exact inverse to eq. (2.12). This is because of the way that we subdivided the line for the fine mesh, and in particular because
hk = Hk + Jk.
(2.13)
Other choices of the relation between the fine and coarse meshes would have led to eqs. (2.11 ). However, they would only be an approximate inverse to eqs. (2.12). The goal of this section is to write the evolution equations for the expansion coefficients
L.G. Margolin, D.A. Jones / Inertial manifold for computing Burgers" equation
{a, b}. We begin by writingthe {a,i} equations in the {Hk, Jk} basis: Oak kffil
2
IHk) \ Ot +
+
ISk> k=l
,__~ +
2Ax
(
?
Ax "--i
Oak . i ~ - - i 2 - 1 ~
]
421k+ik-I
Ax 2
= (fgB _ f~ ).
)
2ak (2.17)
2ak _ f ~ )
+
\or
2 --ak 2 42ak+l + a k _ m i k ~ ak+1
(Ocg~.~k
-- " - ~ "P
2~
+ tk-IAx 2-
_42ik
2
i k -- t k - I
179
Finally, we use eq. (2.11) to replace { a , t } by {a, b}. Further we make the approximations, as in eq. (2.1),
a / + 1 --
_42ak+i + ak -- 2ik _ fk#'~ Ax E
)
(2.14)
-- 0,
where fk%fk# are the coefficients of f ( x ) with respect to the IHk), IJk) basis respectively. In eq. (2.14) we show the semidiscrete form, in which the spatial derivatives are discretized, but the time derivatives are shown in analytic form. We are not concerned with the semidiscrete equations, but the time differencing is not relevant to our discussion. Furthermore, we have chosen a particular discretization for the spatial terms, both of which are centered in space. We will discuss the boundary conditions in detail in section 4. We convert the equations to the {h,s} basis in two steps. First, we change the basis elements themselves, by inserting a complete set of states
fka + fk'~ ~ 2 f (kAx) := 2j~, Ax d f (kAx)
dx
:= ½ fk'.
For the hk equation, Oak Ot
22ak+1 + a k - l - 2 a k - bk+l + bk-1 ~kX2 2
+
2
ak+ 1--ak_ 1 1 4A.X + ~(--2ak+lbk+l
+ 4akbk
--2ak-lbk-, +b~+, + b~_, - 2b~) = 3~,
(2.18)
and for the Sk equation Obk Ot
22ak+l -- ak-1 -- 6bk -- bk+l - bk-1 Ax 2 a2+l + a ~ - l - 2a2 1 + 4AX + -~--~(--2ak+lbk+l + 2 a k - l b k - , + b~+ 1 + b~_ 1 - 2 b ~ )
N
(2.15)
I = ~_~(lhk)(hk[ + [Sk)(Skl). I
Noting that each of the hk and Sk is independent, we derive two sets of equations, valid for each k. The hk equation, still in terms of {a, t } is + ak+l--ak
42ak+1
2Ax
(Oak
+ \ Og "1"
~
= (fk# + fk~). The Sk equation is
"k--t2-1
+ak-mil k
Ax 2 + t k _ l - - 2ak'~ 42 t k
ax2
}
(2.16)
1 = ~Ax f~.t
(2.19)
3. The approximate inertial manifold There is an important qualitative difference between eqs. (2.18) and (2.19) on the one hand, and their spectral analogs on the other. The difference is that in the spectral case, the linear term contains only terms in the lower wavenumbers. However in eq. (2.18 ), the linear term contains both ak and bk. The reason is that the coefficients ak and bk are not really independent in a
180
L.G. Margolin, D.A. Jones / Inertial manifold for computing Burgers" equation
"truncation error" sense. In fact, from eq. (2.12 ) we see Ou A x bk = ½ i l k - - O t k -- 19X 4 + O(Ax3)"
t~k ~-- -~
2/
kAx
u(x)dx
(k-1/2)Ax
(3.1) that
However, we can also write O~k ----u ( k~kx ) -- ux ( k A x ) ~ J x 1 ~(ak+l--ak-1)
O__.~uAx + O(Ax3)"
= OX 4
(3.2)
+ Uxx(kAx)~(Jx)
2
- Uxxx(kAx)13~,2(Jx)3 +
Since these two expressions can only be distinguished in the third order, and since the linear term is of the second order, it is not surprising that both expressions appear in the linear term ofeq. (2.18). In the next step we want to use eq. (2.19) to find an approximate relation between ak and bk that is based on the slaving of modes due to the nonlinearity. In particular, we have to be careful in estimating of the various terms. Combining eqs. (3.1) and (3.2) we have
O(Jx4).
Expanding the other terms in (3.4) in a similar fashion, we have that eq. (3.4) is approximately 22 / 1 0 3 u 3 :~.,,O(Axs)) ax2 ax If we think of bk as a continuous function, we have, using eq. (3.3), for the linear term in (2.19) _22ak+l -- a k - I -- bk+l - b k - i - 6bk /~X 2
bk ~ ak + l -- a k - I 8 -
cgU A x Ox 4 + O(Ax3)"
This relation is derived from the truncation analysis and in not due to the slaving of the higher order modes. It seems appropriate to use truncation analysis to estimate the magnitude of the other terms in the equations. That is, we will use Taylor expansions, and compare terms using the power of the cell size Ax that appears in the expansion. An approximate relation between ak and bk is derived by making some approximations in eq. (2.19). To begin, we ignore the time derivative of bk. This approximation is also made in the spectral case. For the linear term we go back to the definitions of a and b in terms of a and il. We obtain 22
AX 2 (Otk+l -- 3ilk + 3O~k -- i l k - l ) .
_O2b A-~.
(3.3)
(3.5)
The viscous term must be balanced by the nonlinear term. W e therefore have
02b ~-e f~ + ~--~(a2+l + a2_l - 2ak)
OX 2
2e A x 2 (ak+lbk+l -- a k - l b k - l ),
where e := Ax/42. Also, a similar analysis as 2 1 + b2_l - 2b 2 ,,~ (Ax4). above shows that bk+ This term is dearly small compared to the other terms in eq. (2.19) Thus, we have neglected this term in eq. (3.6). We continue by substituting the lowest order result of eq. (3.3) on the fight side of eq. (3.6) bk = ~(ak+l -- a k - l ) .
(3.4)
This appears to be a third derivative. Indeed, if we expand a function u about the spatial point k, we get, for instance, from
(3.6)
Using the approximations ~(ak+2 -- ak) ~-- ¼(ak-1 -- ak), ~(ak--ak-1)
~ ~(ak--ak-l),
L.G. Margolin, D.A. Jones / Inertial manifold for computing Burgers" equation
that the computational results are somewhat sensitive to the boundary conditions on b, and these choices are not appropriate. To find a more appropriate choice, we return to the partial differential eq. (1.1). On the left boundary, u = Uleft and so is constant in time. Thus on the left boundary
ak ~-- ½(ak-I q- ak+l),
we find ak+ lbk + 1 -- ak-lbk-1 1 2 I ~- ~(ak+
+
a2_i
--
181
2a~).
It follows that 02b e 02a 2 Ox 2 ~- 2 cox 2
OOU .OO2U Uleft~ - = /t~-~ + f ( 0 ) . eft.
(3.7)
The most general solution of eq. (3.7) has the form bk = ½e(a 2 + cl + c2x) - e F k ,
(3.8)
where Fk is F (kAx) and d F / d x = f . The conslants ct, c2 will be determined from the boundary conditions. This will be done in the next section. We note that it is no accident that both sides of eq. (3.7) are composed of second derivatives. eq. (2.19) is a differenced form of Burgers' equation after one additional differentiation in space, eq. (3.7) represents the steady state of this equation. Note however, that deriving the AIM (3.8) and the associated boundary conditions, we never made use of a knowledge of the analytic solution to Burgers' equation.
4.
(4.1)
Implementation
In this section we briefly discuss two aspects of the implementation of the AIM. The basic strategy is to compute with eq. (2.18). We use the approximate relation of eq. (3.8) to replace the bks where they appear. Now suppose that we specify the values of u on the left and the right boundaries. For an Nth order calculation, we might set al = Ule fl and a~v+1 = Uright, and then compute the interior values ak {k = 2 . . . . . N}. We note that to calculate a2 we will need bl. In eq. (3.8) these are not naturally defined. One could simply set b~ = blv+l = 0. However, it turns out
Returning to eq. (3.8), we can estimate the derivatives au 4bl - - - + 0(~), Ox Ax 02u 32 = ~ - - ~ (bE - b l ) Ox 2
+ O(f).
(4.2)
Note that the first and second derivatives in eq. (4.2) are centered at different spatial points, and neither are exactly at the left boundary. We neglect these corrections, and also terms of order e. Substituting (4.2) in (4.1) leads to
A• 5
leftbl = ( b 2 - b l ) +
5Ax2 ~-~ f ( 0 ) .
(4.3)
Finally, we note that the coefficient of the lefthand-side is proportional to e; further, the coefficient in front of f (0) is proportional to eAx. Thus for consistency these terms should be ignored, and on the boundary we have bl = b2,
bN+l = bN.
(4.4)
Next we approximate the constants cl and c2 defined in eq. (3.8). On the left boundary we have 8b12 Ax ~ u~eft + cl - 2 F ( 0 ) ,
(4.5)
where we have made the arbitrary choice x = 0 on the left boundary. On the right boundary 8bN+12 2 AX ~ Udght + Cl + c2L - 2 F ( L ) .
(4.6)
182
L.G. Margolin, D.A. Jones / Inertial manifoldfor computing Burgers"equation
We cannot apply the boundary conditions (4.4) since bl and bN+l a r e not yet known. Instead, we use the approximation of eq. (3.3) for b2 and bN and then solve eq. (4.5) and (4.6) for the constants c~ and c2. We mention one more point relating to the boundary conditions (4.4). Corrections due to the AIM in eq. (2.18) appear both in the dissipative term and in the nonlinear term. In the nonlinear term, the new terms have the flavor of higher order corrections. For example, if we set bk = 0, we regain the original differencing of this term. This is not true in the dissipative term. The first difference of bks is the some order as the second difference of aks. If we set bk to zero in this term, we do not reproduce the original differencing of this term. Thus in the dissipation, the bk terms do not represent higher order corrections. Now in deriving the boundary conditions (4.4), we have kept only lowest order terms. These are sufficiently accurate for the higher order corrections in the nonlinear term, but not for the new formulation of the dissipative term. One could try to find more accurate boundary conditions. However (4.4) is particularly simple to implement. Furthermore, we expect most of the benefit of the AIM to result from the nonlinear term. Hence, we replace the new form of the dissipative term with the old. The results in the next section were all computed with eq. (4.7), and using the boundary conditions (4.4) 2 _ ~2_1 Oak 2ak+t + ak-l -- 2ak ak+l Ot Ax 2 + 4Ax + (--2ak+lbk+l + 4akbk - 2ak-lbk-1
+b2+l + b2_1 - 2b 2) (4Ax) -l = fk.
(4.7)
5. Numerical results
In this section we present numerical results to illustrate the efficacy of the AIM as applied to
finite difference schemes. First we describe our model problem, and the measure we use for comparison. Then we compare the accuracy of calculations using our new scheme eq. (4.7) with calculations using a standard finite difference scheme, which is (4.7) with bk set to zero. We solve Burgers' eq. (1.1). We set Uleft = 10.0 and uright = 0.0 for boundary conditions. We choose an interval length L = 1.0 and set f = 0. We treat the total number of cells, N and hence Ax, and the viscosity 2, as parameters. In each calculation, we run the problem to steady state. There is a simple analytic solution for this problem, which we can use to calculate the accuracy of our solutions. In each cell we compare the computed value with the analytic solution U (x) evaluated at the cell center. We define the L 2 norm of the error by
L2(a, U) =
~_,[ak -- U(xk)122~
.(5.1)
k=2
We have also computed a n L 1 norm and an L °~ norm and find that all of our results are qualitatively unchanged. In each of the figures in this section, we plot a ratio, the L 2 norm of the error of the standard calculation divided by the L 2 norm of the error of the AIM calculation. Thus, a ratio greater than 1 means that the AIM calculation is more accurate in the global sense. In some cases this ratio may decrease as the number of cells is increased. However, we note that in every case the actual L 2 error of both the standard calculation and the AIM calculation monotonically decrease as the number of cells is increased. In figure 5.1 we show parameter studies where we hold fixed the viscous coefficient A and vary the number of cells. For 2 = 5.0, the ratio of errors is 3.3 (for N = 5), increases to a maxim u m value of 8.0 (for N = 15) and then slowly decreases. For larger values of A, the graphs are very similar, each showing a maximum at about N = 15. The shapes of the graphs tend to be-
L.G. Margolin, D.A. Jones/Inertial manifold for computing Burgers" equation
hk
183
Sk /
,
lamlxla=5.
k-I
lambda-20.
P
lambda=l.
o
== .
.
.
.
.
.
.
i
.
.
.
.
.
.
.
.
10
i
.
.
.
.
'.X
k+l
'
k-1
//k
'.X
k+l
Fig. 6.1. An alternate basis set for the {Sk}.
"6
.
,
k
.
.
100
.
.
1000
number of cells
Fig. 5. l. Ratio of errors in the L 2 norm for several values of ,l as the number of cells increases.
To summarize, the AIM leads to an improvement in accuracy of at least a factor of 2 essentially over the whole parameter space in cell size and viscosity. For large regions of the parameter space, the increase may be substantially greater.
100.
6. Alternate basis sets
0
~. 10
=
0
;
--
:
-
N =20
-
:
•
-
N=5
.
.
'3 0
1
.
.
.
.
.
. 110.
.
.
.
.
.
100
lambda
Fig. 5.2. Ratio of errors in the L 2 norm for several values of total cells as the viscosity increases.
come slightly narrower, and the maxima slightly higher as 2 increases. When 2 = 1.0, there is no maximum. All the graphs appear to asymptote to a value slightly greater than 3.0 as the number of cells increases. In figure 5.2 we show a complementary study in which we hold the number of cells fixed and plot the ratio of errors as the viscous coefficient A increases. Each of the graphs appears to reach a constant value as 2 increases beyond 3.0. For the plotted values, this constant value is largest at 7.8 when N = 20. Its actual maximum appears to be about 8.1 when N = 15 (not shown), which is consistent with the graphs shown in figure 5.1. For larger values of N, the constant value decreases slowly, being about 4.3 for N = 60 and about 3.9 for N = 100 (not shown).
The key to the derivation of the AIM is the basis set shown in figure 2.1. This particular choice led to approximate relation (3.8). However, given this relation, the form of the basis set plays no further role in the implementation. In this section we address the question of to what extent our results would change given a different basis set {s~}. To be specific, we choose the basis set shown in (6.1). Here the higher order basis functions are continuous lines. Note that both the original set {Sk} and the new set {s~} are antisymmetric with respect to the point k. If one follows the same procedure outlined above, one derives the simple result bk =
x/34AX (a~c + cl +
Ax
(6.1)
That is we derive the same relation except for an overall multiplicative constant. Moreover, the only essential ingredient in this derivation is the antisymmetry of the basis function Sk. Any antisymmetric basis function would lead to the same result, up to the overall multiplicative constant. The question naturally arises as to which is the best choice of basis function. The answer is that it depends on the particular function that must be represented (i.e the solution u). Among the various results of our parameter studies, the
184
L.G. Margolin, D.A. Jones / Inertial manifold for computingBurgers"equation
AIM of eq. (6.1) is some times more accurate than (3.8) and sometimes slightly less accurate. Other choices of the basis set {Sk} will lead to other multiplicative constants in (6.1). One finds experimentally that for any given problem, the most accurate result comes form choosing the overall multiplicative factor in a small range, about 1.0 ± 0.15.
Acknowledgements The authors appreciate the continued interest and advice of Edriss S. Titi and Darryl Holm. We gratefully acknowledge the support of the Institute for Geophysics and Planetary Physics (IGPP) and the Center for Nonlinear Studies (CNLS) at Los AJamos. This work was performed under the auspices of the U.S. Department of Energy at Los Alamos National Laboratory.
References [ 1] M. Chen and R. Temam, Incremental unknowns for solving partial differential equations, Num. Math., to appear.
[2] Ch. Dcvulder, M. Marion and E.S. Titi, On the rate of convergence of nonlinear Galerkin methods, submitted. [3] E. Fabcs, M. Luskin and G. Sell, J. Diff. Eq., 89 (1991) 355-387. [4] C. Foias, O.P. Manley and R. Temam, Math. Mod. and Num. Anal. M2AN 22 (1988) 93-I 14. [5] C. Foias, G. Sell and E.S. Titi,J. Dynam. Diff. Eq. l, No. 2, 199-243 (1989) [6] C. Foias and R. Temam, J.Math. Pures Appl., 58 (1979) 339-368. [7 ] C. Foias and E.S. Titi,Nonlinearity 4 (1991 ) 135-153. [8] M.S. Jolly,I.G. Kevrekidis and E.S. Titi,Physica D 44 (1990) 38-60. [9] M. Marion, J. Dynam. Diff. eqs. I (1989) 245-267. [10] M. Marion and R. Temam, 57 (1990) 205-226. [I I ] R. Temam, J.Fac. ofSci.,Univ. ofTokyo, IA, 36 (1989) 629-647. [12] R. Temam, SIAM J. Math. Anal. 21 (1990) 154-178. [13 ] R. Temam, Math. Comput., to appear. [14] E.S. Titi,J. Math. Anal. Appl. 149 (1990) 540-557.