Citation |

- Permanent Link:
- http://digital.auraria.edu/AA00002643/00001
## Material Information- Title:
- Eulerian-Lagrangian localized adjoint methods (ELLAM) discretization strategies for advection-dominated flows
- Creator:
- Trujillo, Rick Vance
- Publication Date:
- 1990
- Language:
- English
- Physical Description:
- 40 leaves : ; 29 cm
## Subjects- Subjects / Keywords:
- Fluid dynamics ( lcsh )
Lagrange equations ( lcsh ) Algorithms ( lcsh ) Boundary value problems ( lcsh ) Algorithms ( fast ) Boundary value problems ( fast ) Fluid dynamics ( fast ) Lagrange equations ( fast ) - Genre:
- bibliography ( marcgt )
theses ( marcgt ) non-fiction ( marcgt )
## Notes- Bibliography:
- Includes bibliographical references.
- General Note:
- Submitted in partial fulfillment of the requirements for the degree, Master of Science, Department of Mathematical and Statistical Sciences
- Statement of Responsibility:
- by Rick Vance Trujillo.
## Record Information- Source Institution:
- University of Colorado Denver
- Holding Location:
- Auraria Library
- Rights Management:
- All applicable rights reserved by the source institution and holding location.
- Resource Identifier:
- 25334816 ( OCLC )
ocm25334816 - Classification:
- LD1190.L622 1990m .T78 ( lcc )
## Auraria Membership |

Full Text |

EULERIAN-LAGRANGIAN LOCALIZED ADJOINT METHODS
(ELLAM) DISCRETIZATION STRATEGIES FOR ADVECTION-DOMINATED FLOWS by Rick Vance Trujillo B.S., Metropolitan State College of Denver, 1988 A thesis submitted to the Faculty of the Graduate School of the University of Colorado in partial fulfillment of the requirements for the degree of Master of Science Department of Mathematics This thesis for the Master of Science degree by Rick Vance Trujillo has been approved for the Department of Mathematics by Thomas F. Russell Stephen F. McCormick Contents 1 INTRODUCTION 1 2 BASIC CONCEPTS 6 3 EXACT WEAK FORMULATION 10 4 VANISHING ADJOINT AND TRACKING ALGORITHM 15 5 APPROXIMATE WEAK FORMULATION 20 6 INLET BOUNDARY DISCRETIZATION 23 7 OUTLET BOUNDARY DISCRETIZATION 28 8 MEASURE ADJUSTMENT AND IMPLICIT ELLAM 31 9 FORWARD TRACKING 36 10 SUMMARY 38 References 40 Rick Vance Trujillo (M.S., Applied Mathematics) Eulerian-Lagrangian Localized Adjoint Methods (ELLAM): Discretization Strategies for Advection-Dominated Flows Thesis directed by Professor Thomas F. Russell Eulerian-Lagrangian Localized Adjoint Methods (ELLAM) are generalizations of the successful Eulerian-Lagrangian Method (ELM) discretization strategy for the Finite Element Method (FEM) solution of the equations that govern advection- dominated flows. In addition to the non-oscillatory and large-Courant-number properties of ELM solutions, ELLAM offers the features of numerical conserva- tion and systematic treatment of boundary conditions. The ELLAM discretization strategy presented here is guided by the analytic solution of advection equations by the Method of Characteristics (MOC) and is designed to treat advection-reaction- diffusion equations with variable coefficients and in multiple spatial dimensions. In this paper, we focus on the theoretical aspects of ELLAM that bear directly upon numerical implementation. We explore exact and approximate weak formulations, vanishing adjoints, tracking algorithms, measure adjustments, and boundary dis- cretizations. The form and content of this abstract are approved. I recommend its publication. 1 INTRODUCTION Advection-reaction-diffusion relations arise in many fluid dynamics problems (e.g. groundwater flows and petroleum reservoir simulations), or in related prob- lems (e.g. semiconductor modeling). These relations pose stiff tests for numerical methods, due to the advection or hyperbolic component of the problem. Tradi- tional Eulerian (static-grid) techniques generally suffer from numerical diffusion and/or spurious oscillations when confronted by advection-dominated problems; for a collection of perspectives on this subject, see [5]. These difficulties can be ameliorated by using very fine computational grids or, preferably, by using opti- mal test functions (OTF) that satisfy a local adjoint relation. The OTF strategy ( see, e.g., [2] ) allows the test functions to adapt in space in order to minimize spatial discretization error; however, they are static in time. As a result, the range of Courant numbers that can treated by OTF methods is usually restricted to be less than unity. For advection-dominated flows, physics dictates that numerical methods with Lagrangian, or moving-grid, properties that can adapt to changing domains of influence would be more appropriate. In particular, the Eulerian-Lagrangian Methods (ELM) [1, 4] that track the advecting flow in order to find the proper domain of influence for the hyperbolic subproblem represent an approach that is closer to the underlying physics of the problem. On grids whose fineness cor- responds to an accurate representation of the moving solution front, ELM can obtain non-oscillatory solutions that are not numerically diffused. However, ELM suffer from two major difficulties: they fail to conserve mass and to handle general boundary conditions. The Eulerian-Lagrangian localized adjoint method (ELLAM) of Celia, Ewing, Herrera, and Russell [3, 6] combines ideas from OTF and ELM by using space- time test functions that satisfy a local adjoint relation. In this setting, the test functions are defined on space-time finite elements that track the evolution of an advecting flow between static computational and data grids at different time levels. Furthermore, ELLAM is able to conserve mass and systematically treat general boundary conditions. Previous papers [3, 6] have presented these ideas for advection-diffusion equations with constant coefficients in one spatial dimen- sion. For the case of advection-diffusion equations with variable coefficients and in multiple spatial dimensions, Russell and the author have recently presented what we shall refer to as the implicit ELLAM formulation [7]. This new formu- lation represents a significant break from ELM and earlier ELLAM formulations in style and substance. The implicit ELLAM formulation views the advection- diffusion problem in Eulerian coordinates and only applies Lagrangian concepts to the formal adjoint of the advection term. In this way, the incorporation of boundary conditions and mass conservation, which distinguishes ELLAM bom ELM, is straightforward. However, the treatment of the adjoint implictly induces a Lagrangian coordinate system which is at the heart of the entire method. Thus, 2 we are compelled to find the appropriate ways to introduce Lagrangian concepts for ELLAM in a variable coefficient and multiple spatial dimension setting. The goal of this thesis is to develop methods for incorporating Lagrangian ideas in an ELLAM setting. Our approach will explicitly use Lagrangian co- ordinates and view the entire ELLAM process from the perspective of a La- grangian reference frame. Finally, we will develop techniques for converting from Lagrangian to Eulerian coordinates, in order to link our explicit ELLAM formu- lation to the implicit ELLAM formulation of [7]. The explicit ELLAM discretiza- tion strategy that will be presented here is designed to treat advection-reaction- diffusion equations with variable coefficients and in multiple spatial dimensions. In Section 2: Basic Concepts, we will pose the non-divergence and divergence forms of the model problems for advection-reaction-diffusion processes. Also, we will see how the analytic solution of the advection-reaction equation by the Method of Characteristics (MOC) introduces us to Lagrangian concepts in a very natural way. Furthermore, we will use our knowledge of the MOC to guide us through the entire ELLAM formulation. In Section 3: Exact Weak Formulation, we will discover how to imbed our knowledge of the analytic solution of the advection-reaction equation into the weak forms of our model problems. We define a space-time, Lagrangian domain in terms of backtracking, from a computational grid to a data grid, along the characteristics of the analytic solution. In this manner, we can view ELLAM 3 as a space-time, Lagrangian variant of the Finite Element Method (FEM) for non-advective parabolic equations. In Section 4: Vanishing Adjoint and Tracking Algorithm, we directly en- counter the inner workings of the Lagrangian coordinate system that is created by the vanishing adjoint. With an eye towards numerical implementation, we de- velop a simple, predicted tracking algorithm that is based on the exact solution of an approximate adjoint relation. And with an eye toward solution accuracy, we develop a higher-order, corrected tracking algorithm that is based on the approximate solution of the exact adjoint relation. In Section 5: Approximate Weak Formulation, we modify the results of Sec- tion 3, in order to incorporate the developments of Section 4. Although the pre- dicted tracking algorithm defines an approximate Lagrangian coordinate system for our equation framework, we can incorporate the corrected tracking algorithm, which is more faithful to the exact Lagrangian coordinate system, to characterize the solution. In Sections 6 ( Inlet Boundary Discretization ) and 7 ( Outlet Boundary Dis- cretization ), we apply the framework of Section 5 and systematically incorporate boundary conditions into our explicit ELLAM formulation. Since we are working in a space-time domain, we are led to incorporating time as a boundary variable. The resulting measure adjustment ( from space to time ) to the integrals of our weak formulation is an Eulerian concept and represents a significant development. 4 In Section 8: Measure Adjustment and Implicit ELLAM, we formally de- velop measure adjustment of the integrals of our weak formulation as a means to convert from Lagrangian to Eulerian coordinates. With measure adjustment, the explicit ELLAM formulation conserves mass and, coupled with the results of Sections 6 and 7, distinguishes itself from its ELM predecessors. Furthermore, measure adjustment leads us to developing the analogous concept of reference frame adjustment to Lagrangian coordinates for implicit ELLAM. Thus, we are able to construct a gateway between implicit and explicit ELLAM formulations. After performing measure adjustment, we find that the implicit and explicit for- mulations are equivalent to within O(At)2 discretization error and lead us to the same discrete equations. In Section 9: Forward Tracking, we identify the forward tracking relation that is derived from the predicted backtracking relation of Section 4 as the mech- anism behind measure and reference frame adjustment. Employing the often-used piecewise linear Raviart-Thomas representation of the velocity field, we outline the numerical implementation of forward tracking. As a result, we find that ex- act forward tracking and mass conservation is possible, whenever we avoid shocks and the tracking of characteristics is well-defined. 5 2 BASIC CONCEPTS The ELLAM framework is designed to treat the two canonical forms of advection-reaction-diffusion equations in m spatial dimensions: Non-Divergence Form Cal = Au + Hu + Vu /(x, t) (1) where Au = ut + t/(x,i) Vu Hu = k(x,t)u Du = V (D(x,f)Vu) are the advection, reaction, and diffusion components of (1) respectively, and v(x, t) = (i(x, Â£),..., um(x, Â£)), fc(x, Â£), and D{x, t) are the velocity field, reaction and diffusion coefficients, respectively, where x = (xi,... ,xTO); and Divergence Form Qu = Au + H*u + Du = /(x, Â£) (2) where H*u = Hu + (V u(x, t))u is a modified reaction term that contains an artificial reaction induced by the di- vergence of the velocity field, v(x, Â£), in the Lagrangian reference frame. In a later section on measure adjustment, we will see that the velocity dependent reaction is canceled when we convert to Eulerian coordinates. Thus, our formulation for the Divergence ( or Conservation law ) form will conserve solution mass and is consistent with physical law. In relation to previous work [7], we note that the usual statement of the di- vergence form is (.4* + 71 + V)u = /, where A* = ut + V (u). However, methods that apply Lagrangian concepts explictly need to work with the La- grangian derivative .A(-) instead of Therefore, explicit ELLAM splits A* into A and TV, in order to obtain a consistent formulation for (1) and (2). In the Lagrangian setting of explicit ELLAM, the adjoint of the advection term for (1) and (2), A, is A, so Lagrangian ideas can be used throughout the formu- lation. In the Eulerian setting of implicit ELLAM, the adjoint of the divergence advection term, A*, is A, so Lagrangian ideas can be applied to the adjoint. However, in Eulerian coordinates, the adjoint of the non-divergence term, A, is A* and there is some ambiguity in the treatment of the adjoint. To complete our specification of the problem, we note that equations (1) and (2) hold for x Â£ E C Rm, t J = (t0,T] and with Dirichlet or flux boundary conditions, u(z,Â£) = h(x,t) or (yu DVu)(z,t) n(x,t) = h(x,t) (3) for x Â£ 3E, t Â£ J and an initial condition, u(x,t0) = w0(s), for x G E. The ELLAM formulation for (1) and (2) is guided by a knowledge of the 7 analytic solution for the numerically problematic advection component of the equations. For the time being, we will focus on the non-divergence form (1) of the advection-reaction-diffusion equation. Taking D(x, t) = 0 and using the Method of Characteristics (MOC), the solution to the non-diffusive case of (1) is given by the following system of ODEs: dt(s,r) dr t(s,tn+1) = dx;(s,T) -- = Xi(s,tn+1) = Si, for i = l,..., m (5) ^dr ^ + r) r)U = r)r) u(s,f*+1) = un+1(s) (6) where tn+1 is the new time level, s = (si,..., sm) 6 E is a parametrization of the grid points at the new time level, and un+1(s) is the solution at the new time level. The idea is to take a point (*i,..., xm, Â£n+1), called the head of the characteristic, from the computational grid and track it backwards in time along the characteristics defined by (5) to its image (sj,... ,a^,Â£*), called the foot of characteristic, on a data grid at a previous time level. For the present case, we assume that t* = tn and our backtracked image remains in the interior of 1, tn+1 (4) 8 Â£ and not on the boundary. With a knowledge of the solution at the old time level, u(xj,..., x*m,tn) = un(x*), the solution to (6) gives us an expression for the solution at the new time level, u(zi,... xm,Â£n+1) = ttn+1(x). Given that (4) tells us that t = r, equation (6) tells us that / -*+1 un+1(x) = un(x*)exp I y fc(x(s,r),r) dr ftn+1 ( ft**1 \ + Jtn exp[~JT k{x{s,T,),T,)dr'J /(x(s,r),r)dr. (7) Now, let At = tn+1 tn. After making 1-point forward and backward Euler approximations to the first and second integrals of (7) and 0(AÂ£)2 approximations of the exponential terms, we obtain un+1(x) = un(x) (1 fcn(x*) AÂ£) + /n+1() At + 0(AÂ£)2 (8) as a second order approximation to the analytic solution (7). The formulation for the divergence form (2) proceeds along the same line with k(x, t) + V v(x, t) replacing k(x,t) in (6), (7), and (8). Symbolically, the key to our discretiza- tion strategy will be to obtain the weak form of (8) for the advection-reaction component of (1) and (2). 9 3 EXACT WEAK FORMULATION For clarity and to establish, a context for our discussion, we will consider the case of a two-dimensional rectangular computational grid, En+1, that is dis- cretized into rectangular subpatches, E?*1. As in standard finite element methods (FEM), we pose the weak form of (1), namely /ni+i w Cu = /n+i / w. However, in contrast to standard methods, fl?*1 is a space-time finite element that is formed by backtracking the computational patch, E?*1, along the exact characteristics given by (5) to the data grid, En. O?*1 defines a Lagrangian coordinate system where the spatial measure on the /n?+1 integral is tied to the head of the char- acteristic and the measure of E-1"1. Also, fly1"1 is used to denote the process rewriting the advection operator from a transport term (ut + v Vu ) in Eulerian coordinates to a derivative ( utau ) in Lagrangian coordinates. The test function w refers to those local space-time test functions, Wij, whose values at tn+1 are non-zero on E-*"1. Additionally, w is standard in space on the computational domain, E?*1, at tn+1, but is defined in ft-1"1 as the solution of a local adjoint relation. With the proviso that the subtleties of this Lagrangian, space-time domain will detailed in later sections, we will concentrate on the mechanics of creating a weak form of (8). For simplicity, we assume that we are working on the portion of the com- putational grid whose backtracked images remain in the interior we consider boundary interactions in a later section. Considering each component of (1) separately, we start with, the advection operator and integrate by parts to obtain / wAu = / A(uw) / uAw. (9) Jo"*1 v ynn+i ^ J J *1 Since .A(-) denotes differentiation along characteristics and integration over fl?*1 can be viewed as integration along these same characteristics, the first integral on the right hand side of (9) is straightforward. For the second integral on the right hand side of (9), we define the test function, -u;, to be the solution of the adjoint relation, Aw = 0. Therefore, fapi wAu = /En+1 ttn+1(a) w X-+1 w (10) since the adjoint relation tells us that w = ton+1(z) = wn(x*). Moving on to the reaction term, we apply a forward Euler discretization to obtain J^iw'R,u = At J ^ kn(x*)un(x*)w + O(At)2. (11) Finally, the diffusion term receives the standard FEM treatment for elliptic equations as we invoke Greens theorem to obtain f wVu= [ f V (w DVu) + f DVu'Vw. (12) Jn$+I v yn~+ v ' Then, we apply the divergence theorem and backward Euler discretization re- spectively, to the right hand side integrals of (12) to obtain wDu = -L Jn+l w D Vu n +At j Dn+1 (z) Vun+l(x) Vw + 0(At)2, 11 where n(x, t) is formed by backtracking the unit outward normal on How- ever, the dE?*1 boundary integral vanishes, since w = 0 on 5E-1"1. In practice, E?*1 is subdivided and w is supported on more than one element and the bound- ary integral does not vanish on those elements. However, the boundary terms from the sudivisions do not need to be computed or stored, since they will sum to zero when combined to form the equations for all of E?*1. Therefore, we have ^n+i wVu = At jf ^ Dn+1(x) Vun+1(s) Vw + 0{At)2. (13) In relation to eariler work in streamline coordinate techniques [5], the treatment of diffusion, even in advection-dominated flows, is crucial. The backward Euler integration puts the diffusion on the new time level where there is no distinction between Eulerian and Lagrangian coordinates. Therefore, the formal incorpora- tion of diffusion is straightforward. However, if the flow is rapidly converging or diverging, as determined by the divergence of the velocity field, backward Euler integration can make large errors unless the time step is small relative to the time scale of the physical phenomena. Thus, the time evolution of the solution can be used to determine the time step that is required for the non-destructive incorporation of diffusion. To close the system, we apply a backward Euler rule to the source term to obtain /E. r+\*)w + 0(Atf. (14) 12 Combining (9), (11), (13), and (14), the weak form of (1) becomes [ un+1(x)w + At f Dn+1Vun+1-Vw Jv?+i v ' 7e"+1 U J = L, "(*') (1 *"(**) Ai) + Af J r" + 0(At). (15) - Lii Lii In relation to standard FEM, we can view (15) as a conversion of (1) and (2) into a'non-advective parabolic problem in the characteristic reference frame, Setting; D = 0 in (15), we obtain the desired weak form of (8) /.,(*)= r,K(i<)(l-*:"(x')At) + rA()+0(Ai)!. (16) "u L*J The development for the divergence form (2) proceeds along the same line with fc(,Â£)+ V v(x,t) taking the place of k(x,t). To establish a basis for comparision, the implicit ELLAM treatment of (1) ( ignoring boundary terms ) yields [ un+1(x) w + Atf Dn+1 Vun+1 Vw / / (V v) u J'Lij JVij Jj+WEy = [ un(x) (1 kn(x) At)w + Af f /Ti+1 w + O(At)2. tÂ£7 For the divergence form, the (V v)u integral is missing in the implicit ELLAM formulation; however, the term reappears in the data integral of the explicit ELLAM formulation. In Section 8, we will find that the additional (V v) u and the are equivalent to the Â£-+1 integration of the data terms on the right hand side of (15). Therefore, we have formally captured the essence of the MOC solution of the non-diffusive case of (1) and (2) in our weak formulation. However, the true 13 essence of the ELLAM discretization strategy can be found in the vanishing of the adjoint relation, A.w, which in turn invokes the solution of (5), the definition of the test function, w, and the definition of the Lagrangian reference frame, Q?.+1. y 14 4 VANISHING ADJOINT AND TRACKING ALGORITHM The key to the ELLAM formulation is the solution of the vanishing adjoint relation, Aw = 0, which is given by the MOC and the system of ODEs defined by (4), (5) and dw(s,r) ^ = tun+1(s) (17) dr where twn+1(s) is usually a standard piecewise polynomial test function that is defined at the mesh points of E?*1. The most familiar examples of iun+1(s) are bilinear test functions on a 2-D rectangular mesh where u;|j+1(si,s2) = 1 at (Xi,yj) and is linearly interpolated to 0 at the adjacent mesh points and is identically zero on the remainder of En+1. Equation (17) says that the test function w is constant along the characteristics defined by (5) and that the ELLAM scheme is structurally conservative. That is, when we define the test functions to sum to 1 on E (e.g. = 1, V(x,y) Â£ En+1), the test functions sum to 1 over the entire space-time domain fin+1. It is clear that the solution to (5) is not crucial for conservation; that is, the solution of (17) for w = 1 is independent of ii(s, t). However, the vanishing of the adjoint relation reduces to the solution of (5)- Our solution to (5) for the adjoint relation is called the tracking algorithm. This tracking algorithm is crucial to the entire ELLAM scheme, because it defines the domain and orientation of integration, Q, it defines w in the fi?*1 space- time domain, and it defines the backtracked image x* and hence all data integrals. What we will find is that our tracking algorithm is best described by a Heun predictor-corrector scheme for (5). In this setting, the space-time elements and the test functions are defined by a predictor that causes an approximate adjoint relation to vanish. However, a better approximation, the corrector, is used to track the solution ( trial functions ), so that our method remains faithful to the exact adjoint relation and to the spirit of the MOC. The guidelines for our tracking algorithm are that it be 1) consistent, 2) simple, and 3) accurate. To see the reason for condition 1), we convert to the integral equation representation of (5), rtn+1 ?(s) = ?+1(.s) / -Ui((s,r),r)dr for z = l,...,m. (18) Jtn Given the way in which the solution appears on both sides of (18), ;(s,t) has a dual role as a) a solution of (18) and b) defining an integration path for Vi(z(s,r),r). Thinking in terms of a Picard iteration scheme (Method of Succes- sive Approximations), it is clear that for general Vt((s,r),r) it is only possible to acheive an internal consistency in the dual roles of a) and b) in a limiting sense. However for special v,-, e.g. = constant, it is possible to directly acheive consistency. The reason for a simple tracking algorithm is one of practicality. Tracking will be the numerical workhorse for evaluating the crucial data integrals in (15) and 16 (16). A complex scheme can grind a numerical implementation to a halt given that Vi depends on x(s, r), not just x;(s,r), and (18) represents a coupled system for i(s, r). Also, our knowledge of the velocity field is generally restricted to the current and previous time levels. Therefore, dynamic time-stepping algorithms, such as Runge-Kutta, which require knowledge of V{ at intermediate times may be ill-suited, if not cost-prohibitive, for our task. The simplest tracking rule for x* in a variable coefficient setting is based on a backward Euler approximation of (18) x*(s) = x-^(s) + O(At)2 where X;^(s) = x"+1(s) v"+1(xn+1(s))AÂ£ for i = (19) If we use the predicted characteristic given by (19) in a predictor-corrector scheme for (18), as suggested in (22) below, we find that we can acheive the consistency criterion when the vt- are constant along the characteristics defined by (19). That is, the corrected value for x\ will coincide with the predicted value. The concept of an accurate tracking algorithm needs some clarification. For our purposes, accuracy is defined with respect to the adjoint relation and not to the integral equation (18). We need to construct x*(p) so that a modified form of the original adjoint relation is identically zero; however, x*(p) can be a coarse approximation to the true backtracked image x*. Given the discussion of the previous two criteria for a tracking algorithm, we define u(x, t) to be the 17 velocity field defined by a pointwise extrapolation of nTl+1() as a constant along the approximate characteristics of (19) v(x,t) = backtrack(vn+1(x)) (20) i.e. un+1(x) = Un+1(a;) = Then, we can solve (18) exactly with U* replacing and hence, we can make Aw = 0, where Aw = wt + v Vw, (21) by using the coarse sE^hcharacteristics of (19). Although the accuracy of the approximation of (18) is not crucial for the modified adjoint relation, Aw = 0, accuracy is necessarily all-important for the solution un+1(). With respect to the solution, we must be able to build a higher-order approxi- mation, x*(c\ that is invisible with respect to the discretization errors committed in the construction of (15). That is, x*(c) must direct the solution un+1(x) to the correct (with respect to the discretization errors) history value un(x*(c)). If *(c) fails to characterize the correct domain for un+1(x), we expect some Courant number limitation (Cu < 1) for the accuracy, but not the stability, of the solution technique. To define x*(c\ we return to (18) in residual velocity form xt*(s) = aÂ£+1(s) f u,(x*(p)(s, t), r) dr J wn+l f (v u)i(a*^(s,r),T) dr + 0(At)3 where the 0(A<)3 truncation error comes from the expansion of the integrands about the point x*(p\ Then, we employ a 2-point trapezoidal approximation to 18 obtain *;w = *:(0)+o(At)J where *?c\s) = - xAt(u ti)(i'(p)(s)) (22) = *?+1(*)-5AÂ£(f+I(a."+*(.)) + .?(x^'>W)). (23) Despite the appearance of (22) and (23), x*(p) is actually an O(At)2, and not an 0(AÂ£) approximation to x*(c\ To see this fact, consider (u v)"(*^) = At u?+1(x) v?(x ) v?+l(x) v^(x
At ) + 0(AÂ£)2, )) w + 0{At)3. (30) ) + 0(AÂ£)2 ]AÂ£ ) At + O(At)3. ) kn(x ) At = [ n+i un(x |