R E G U L A R I Z A T I O N M E T H O D S F O R D I F F E R E N T I A L E Q U A T I O N S A N D T H E I R N U M E R I C A L S O L U T I O N B y Ping L i n B . Sc. (Mathematics), Nanjing University, Nanjing, P .R .Ch ina , 1984 M . Sc. (Appl ied Mathematics), Nanjing University, Nanjing, P .R .Ch ina , 1987 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE STUDIES INSTITUTE OF APPLIED MATHEMATICS DEPARTMENT OF MATHEMATICS We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA December 1995 © Ping L i n , 1995 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of H t x t t i e m J - i c s The University of British Columbia Vancouver, Canada Date T>ec -e^U-<- i"l , 19*1 ^ DE-6 (2/88) Abstract M a n y mathematical models arising in science and engineering, including circuit and device simulation in V L S I , constrained mechanical systems in robotics and vehicle simulation, certain models in early vision and incompressible fluid flow, lead to com-putationally challenging problems of differential equations with constraints, and more particularly to high-index, semi-explicit differential-algebraic equations ( D A E s ) . The direct discretization of such models in order to solve them numerically is typically fraught wi th difficulties. We thus need to reformulate the original problem into a bet-ter behaved problem before discretization. Index reduction with stabilization is one class of reformulations in the numerical solution of high index D A E s . Another class of reformulations is called regularization. The idea is to replace a D A E by a better behaved nearby system. This method reduces the size of the problem and avoids the derivatives of the algebraic constraints associated with the D A E . It is more suitable for problems wi th some sort of singularities in which the constraint Jacobian does not have full rank. Unfortunately, this method often results in very stiff systems, which accounts for its lack of popularity in practice. In this thesis we develop a method which overcomes this difficulty through a combination of stabilization and regularization in an iterative procedure. We call it the sequential regularization method ( S R M ) . Several variants of the S R M which work effectively for various circumstances are also developed. The S R M keeps the benefits of regularization methods and avoids the need for using a stiff solver for the regularized problem. Thus the method is an important improvement over usual regularization methods and can lead to improved numerical methods requiring only solutions to linear systems. The S R M also provides cheaper and more efficient methods than the n usual stabilization methods for some choices of parameters and stabilization matr ix. We propose the method first for linear index-2 D A E s . Then we extend the idea to nonlinear index-2 and index-3 problems. This is especially useful in applications such as constrained mult ibody systems which are of index-3. Theoretical analysis and numerical experiments show that the method is useful and efficient for problems wi th or without singularities. Whi l e a significant body of knowledge about the theory and numerical meth-ods for D A E s has been accumulated, almost none of it has been extended to partial differential-algebraic equations ( P D A E s ) . As a first attempt we provide a comparative study between stabilization and regularization (or pseudo-compressibility) methods for D A E s and P D A E s , using the incompressible Navier-Stokes equations as an in-stance of P D A E s . Compared with stabilization methods, we find that regularization methods can avoid imposing an artificial boundary condition for the pressure. This is a feature for P D A E s not shared with D A E s . Then we generalize the S R M to the nonstationary incompressible Navier-Stokes equations. Convergence is proved. Again nonstiff t ime discretization can be applied to the S R M iterations. Other interesting properties associated with discretization are discussed and demonstrated. The S R M idea is also applied to the problem of miscible displacement in porous media in reservoir simulation, specifically to the pressure-velocity equation. Advan-tages over mixed finite element methods are discussed. Error estimates are obtained and numerical experiments are presented. F ina l ly we discuss the numerical solution of several singular perturbation prob-lems which come from many applied areas and regularized problems. The problems we consider are nonlinear turning point problems, a linear elliptic turning point prob-lem and a second-order hyperbolic problem. Some uniformly convergent schemes wi th respect to the perturbation parameter are constructed and proved. A spurious solution phenomenon for the upwinding scheme is analyzed. i i i Table of Contents Abs t rac t i i L i s t of Tables v i i L i s t of Figures v i i i Acknowledgement ix 1 In t roduct ion 1 1.1 Regularization for Differential-Algebraic Equations ( D A E s ) 1 1.2 Regularization for the Incompressible Navier - Stokes equations . . . . 7 1.3 A Problem in Reservoir Simulation 8 1.4 Regularization for Differential Equations without Constraints 10 1.5 Contribution of This Thesis 12 2 Sequential Regular iza t ion Methods for Differential Algebra ic Equa-tions 15 2.1 Motivat ion of the S R M for General High Index D A E s 15 2.2 Linear Index-2 Problems 21 2.3 Problem Conditioning 23 2.4 Derivation of the S R M 27 2.5 Convergence Analysis of the S R M 30 2.6 Discretization and Implementation Issues 34 2.7 Numerical Experiments 41 2.8 More about the Proof of Theorem 2.1 46 iv 3 S R M for Nonl inear Problems 52 3.1. Nonlinear, Nonsingular Index-2 Problems 53 3.1.1 The case ct\ = 1 55 3.1.2 The case ctx = 0 58 3.2 Nonlinear, Singular Index-2 Problems 62 3.3 S R M for Nonlinear Higher-index Problems 64 3.3.1 The case of nonsingular GB 65 3.3.2 The case for constraint singularities 71 3.4 S R M for Constrained Mul t ibody Systems 72 3.5 Numerical Experiments 75 4 S R M for the Nonsta t ionary Incompressible Navier-Stokes Equat ions 83 4.1 D A E Methods for Navier-Stokes Equations 83 4.2 Preliminaries and the Properties of the Regularized Problems . . . . 88 4.3 Convergence of the S R M 94 4.3.1 Two linear auxiliary problems 94 4.3.2 The error estimate of S R M 97 4.4 Discretization Issues and Numerical Experiments 103 5 S R M for the Simula t ion of Misc ib l e Displacement in Porous M e d i a l l 2 5.1 Introduction 112 5.2 S R M Formulation 113 5.3 Convergence Analysis 116 5.4 The Galerkin Approximation and Its Error Estimates 120 5.5 Numerical Experiments 122 6 N u m e r i c a l Methods of Some Singular Per turba t ion Problems 127 v 6.1 One Dimensional Quasilinear Turning Point Problems 128 6.1.1 A repulsive turning point problem 128 6.1.2 A n attractive turning point problem 133 6.2 Notes about Spurious Solutions of Upwind Schemes 139 6.2.1 Inadequacy of Yavneh's argument 139 6.2.2 Our explanation 141 6.3 A Linear Hyperbolic-Hyperbolic Singularly Perturbed Initial-Boundary Value Problem 145 6.3.1 Construction of asymptotic solution and its remainder estimate 147 6.3.2 Difference scheme and its uniform convergence 150 7 Conclus ion and Future W o r k 154 7.1 Summary and conclusions 154 7.2 Discussion of future work 157 Bib l iography 160 vi'. Lis t of Tables 2.1 S R M errors for Example 2.2 using the midpoint scheme 43 2.2 S R M errors for Example 2.2 using the shooting-back technique . . . 44 2.3 S R M errors for Example 2.3 using backward Euler 44 2.4 S R M errors for Example 2.3 using forward Euler 45 2.5 Errors near singularity using modified formula (2.42) 45 2.6 Errors for problem without singularity using modified formula (2.42) 45 3.1 Errors for Example 3.3 using the explicit second order Runge-Kut ta scheme 76 3.2 Example 3.4 - bounded y and singularity at t* = .5 78 3.3 Example 3.5 - unbounded y and singularity at t* = .5 78 3.4 Errors for Example 3.6 using S R M ( 3 . 4 5 ) - ( 3 . 4 6 ) 80 3.5 Drifts of the S R M for the slider-crank problem 81 4.1 S R M errors for \x = 0.1 without upwinding 109 4.2 S R M errors for fj, = 0.001 with upwinding 110 4.3 S R M errors for /i = 0.001 with a pretty large time step k = h = 0.1 . 110 4.4 S R M errors for \x = 0.1 with a i = 0 I l l 5.1 Numerical results for Example 5.1 with grid size = 4^ 126 5.2 Numerical results for Example 5.1 with grid size = ^ 126 5.3 Numerical results for Example 5.2 126 v i i Lis t of Figures 2.1 planar slider-crank: ini t ia l state in solid line, subsequent states in dot-ted lines 17 3.1 Two-l ink planar robotic system 79 3.2 Solution for slider-crank problem with singularities 81 3.3 Acceleration of slider end 82 5.1 One element with velocity on each edge and pressure at the center . 123 v m Acknowledgement I am specially grateful to Dr . U r i Ascher, my thesis supervisor, who contributed significantly to the development of this thesis. That contribution ranged from fruitful suggestions and enlightening discussions to general guidelines and ideas about how to present the results. I am also grateful to Drs. John Heywood, T i m Salcudean, Michael Ward and Br ian Wetton, members of my supervisory committee, for their helpful comments and suggestions during this research and writing of the thesis. I would also like to thank the financial supports from K i l l a m pre-doctoral fellow-ship committee and from Department of Mathematics, University of Br i t i sh Columbia . Most of al l , I wish to acknowledge the tolerance displayed by my wife, Yuchun, and by my daughter, Xueshu (Susan) during this work. This thesis work often dominated the evenings, weekends and summers which could have been used in activities of greater interest to them. ix Chapter 1 Introduct ion Many mathematical models arising in science and engineering, including circuit and device simulation in V L S I , constrained mechanical systems in robotics and vehicle simulation, certain models in early vision and incompressible fluid flow, lead to com-putationally challenging problems of differential equations wi th constraints, and more particularly to high-index, semi-explicit differential-algebraic equations ( D A E s ) . The direct discretization of such models in order to solve them numerically is typically fraught wi th difficulties, and most methods proposed in the literature seek to circum-vent this by employing combinations of problem reformulation, regularization 1 and special discretization techniques. We wi l l consider the regularization of such mathematical models and the numer-ical solution of the resulting regularized formulations. These formulations are often singular perturbation problems because they typically depend on a small parameter which provides a measure of the closeness between the regularized and the original problems. We wi l l also apply our regularization method and idea to other relevant practical problems. 1.1 Regular iza t ion for Differential-Algebraic Equations ( D A E s ) D A E s are special implici t ordinary differential equations (ODEs) ^ h e concept of regularization was introduced by Tikhonov (see [109]). Its idea is that one solves a better behaved nearby problem instead of solving the original problem to circumvent some sort of difficulties. See the next section for more details. 1 Chapter 1. Introduction 2 where the partial Jacobian matrix fy(y,x,t) is singular for al l relevant values of its arguments. Here x' = A n extension to partial differential equations is considered in the next subsection. D A E s were motivated by applications like network analysis, circuit simulation and mechanical system simulation starting in the 1970's. They often arise as ordinary differential equations with additional variables and (equality) algebraic constraints. A n extensive list of applications is given in [92]. In the 1980's, D A E s have developed into a highly topical subject of computational and applied mathematics. Contributions devoted to D A E s have appeared in various fields, such as applied mathematics, scientific computation, mechanical engineering, chemical engineering, system theory, etc. Frequently, other names have been assigned to D A E s , e.g. semistate equations, descriptor systems, singular systems. Gear ([51], 1971) proposed to handle D A E s numerically by backward differentiation formulas ( B D F ) . For a long time D A E s had been considered not to differ essentially from regular implici t O D E s in general. Only since about 1980, because of computational results that could not be brought into line with the above supposition (e.g. Sincovec et al [103], 1981), the mathematical and particularly the numerical community have started investigating D A E s more thoroughly. W i t h their famous paper, Gear, Hsu and Petzold ([52], 1981) started a discussion on D A E s that wi l l surely be carried on for a while. The structure of D A E s is very much related to the concept of index, which is a measure of the amount of singularity of the system. There are several ways to define index. The most popular one is called differential index, which is defined as the min imal number of analytical differentiations in t such that (1.1) can be transformed by algebraic manipulations into an explicit ordinary differential system (in the original unknowns) x' '= (f>(x,t) Chapter 1. Introduction 3 which in turn is called the "underlying O D E ". There are several structural forms of D A E s which appear frequently in applications (see [92]). The differential index of these structural forms can be found by differentiating their algebraic constraints wi th respect to t (and substituting into the differential equations which complement the algebraic constraints). For instance, • Semi-explicit index-1 system x = f(x,y), (1.2a) 0 = g(x,y), (1.2b) if gy is invertible. • Hessenberg index-2 system x' = f(x,y), (1.3a) 0 = g(x), (1.3b) if gxfy is invertible. • Hessenberg index-3 system x' = f{x,y), (1.4a) y' = k(x,y,z), (1.4b) 0 = g(x), (1.4c) if gxfykz is invertible. Mechanical multibody systems with holonomic con-straints are examples.of Hessenberg index-3 D A E s . In [56] it is pointed out that higher index (> 2) D A E s , in the natural function space formulations, lead to ill-posed 2 problems because they do not have the usual 2 A problem is called ill-posed if it is not well-posed. A problem is called well-posed if it satisfies three conditions, i.e. the existence, uniqueness and stability of its solution. "Stability" means that the solution of the problem continuously depends on the "data", which may be initial data, boundary data, coefficients in the equation, values of the operator, etc.. Here high-index D A E s fail to satisfy a stability condition. Chapter 1. Introduction 4 stability property of differential equations in general. E x a m p l e 1.1 Consider (See [86]) x' = y,' (1.5a) 0 - x-p(t), (1.5b) where x(t) and y{t) are scalar functions and p(t) is a given function. This is a very simple index-2 DAE. The exact solution is x = p(t), y — p'(t). If we add a small perturbation 5smut, 8 << 1, to the right hand side of the second equation we have the exact solution x = p(t) + 8 sin uit, y = p'(t) + 8u> cos tot. Hence y — y = StocosLot could be very large if to 3> 1/8, i.e. the solution changes a lot under a small change in the right-hand side of the equation. • A numerical method which is directly applied to a complex, ill-posed problem may generally fail . Therefore, to solve D A E s , we have to stabilize such problems to bring about continuous dependence on the "data" (or stability). One such approach is to change the formulation of the problem but not its solution, e.g. in Example 1.1 we differentiate (1.5b) once and obtain x' = y (1.6a) 0 = y-p'{t), x ( 0 ) = p ( 0 ) . (1.6b) Now (1.6) is of index-1 and becomes a well-posed problem with the same solution as (1.5). We can solve (1.6) instead of (1.5) and gain well-posedness. However, such a direct index reduction procedure may cause the well-known drift difficulty (see [29]), i.e. the approximate solution of (1.6) may be far from satisfying the constraint (1.5b). Chapter 1. Introduction 5 Hence, methods have been designed to prevent moving away from the constraints. Baumgarte's stabilization [17] and projection invariant methods (cf. [16]) are popular among such methods. Most of these approaches treat in i t ia l value problems, and only a few apply to boundary value problems. See [29, 58] for various numerical methods for in i t ia l value problems; [93, 7, 15] for boundary value problems; and [16, 9, 35] for a survey of various stabilization reformulations. Another approach consists in adding some small perturbation terms (measured by a small positive parameter e) to the given D A E . The perturbed problem is close to the original problem (if e is small) and is well-posed. Such an approach is usually called regularization. This is a natural approach since the high-index D A E is i l l -posed; indeed in [43] a well-known Tikhonov regularization algorithm (see [109]) was applied to solve D A E s . However, such a method seems so general that it is not sufficiently related to the special structure of the D A E . There are two types of regularization methods which are probably more interesting for D A E researchers. One is called parameterization. One such possibility, the pencil regularization, was given independently by Boyarintsev [24] (or see his newer book [25] published in English) and Campbel l [32]. But the regularized problem is ensured to be well-posed only for constant coefficient cases. A further parameterization was proposed by Marz [85]. Her regularization is aimed at obtaining well-posed index-1 D A E s instead of obtaining well-posed O D E s . Heuristically, it seems evident that the D A E is less changed i f it is transformed into an index-1 D A E rather than an O D E . Marz 's regularization was proved to be well-posed for usual structural forms of D A E s . We refer to [59] for further results in this direction. Another class of regularization uses the penalty idea (see [84, 91]). It originates from penalty methods for constrained optimization problems. Note that an algebraic equation in a D A E can be viewed as a constraint. This method seems more natu-ral for D A E s . References [68, 70, 69] used the penalty regularization and singular Chapter 1. Introduction 6 perturbation theory to determine the solutions of D A E s when the in i t ia l or bound-ary values are given improperly (i.e. inconsistently). In Chapter 2 we wi l l mention these methods again with a bit more detail and indicate that Maxz's regularization is actually a kind of penalty method. Because the regularization method requires fewer differentiations of the constraints it is perhaps more suitable for D A E s which have singularities, i.e. whose constraints do not have full rank, e.g. when the matrix gxfy is singular at some isolated points in the index-2 system (1.3). These problems can be challenging for the methods that are usually employed and appear frequently in simulation of constrained mechanical systems. To our knowledge, there has not been a paper in the numerical analysis literature about this unt i l the recent two preprints [11] and [94] (although a number of relevant papers appear in the mechanical engineering literature). From a practical point of view, a number of codes which work well and efficiently (at least i f the regularization parameter e is not too small) are available for numeri-cally solving the regularized problems. We also note that the regularization method requires less smoothness of the coefficients of the differential-algebraic problem than other stabilization methods. These are the advantages of-the regularization method. The dominant disadvantage in the above regularization methods is that the parame-ter e must be small enough to maintain the accuracy of the numerical method we use for the regularized problem at an acceptable level. Hence, a stiff solver is necessary. Typica l ly in regularization methods, the parameter e must be chosen both "large enough" and "small enough": large so that the regularized problem would behave significantly better than the original, and small so that its solution wi l l not differ too much from that of the original problem. We wi l l present a class of new regularization methods, inspired by [19], which we call sequential regularization method ( S R M ) . The S R M can be viewed as a combina-tion of the penalty method with Baumgarte's stabilization in an iterative procedure; Chapter 1. Introduction 7 see §2.4 or §3.1 for specific instances. It is applicable for D A E s wi th constraint sin-gularities. Moreover, the regularization parameter e in the method is not necessarily small . Thus, a nonstiff solver can be used for solving the regularized problems. Some variants of the S R M are discussed for index-2 and index-3 D A E s wi th the goal of sim-plifying the computations. We wi l l apply the S R M to mechanical mult ibody systems as well. 1.2 Regular iza t ion for the Incompressible Navier— Stokes equations As noted before, D A E s have become a highly topical subject of applied and numerical mathematics. However, there seems to be st i l l a void in the literature about partial differential equations with constraints ( P D A E s ) . A typical instance of such problems is the well-known incompressible Navier-Stokes equations: Ui + (u • grad)u = /^Au — gradp + f, (1.7a) divu = 0, (1.7b) u\dn = b , u | < = 0 = a, (1.7c) in a bounded two- or three-dimensional domain fi and 0 < t < T. Here u(x,£) represents the velocity of a viscous incompressible fluid, p(x, t) the pressure, f the prescribed external force, a(x) the prescribed ini t ia l velocity, and b(t) the prescribed boundary values. The system (1.7) can be seen as a partial differential equation wi th constraint (1.7b) wi th respect to the time variable t. Hence, we call it a P D A E . It is easily verified that it is of index-2 without singularities since the operator divgrad = A is invertible (under appropriate boundary conditions). A huge number of methods have been designed to solve the nonstationary incom-pressible Navier-Stokes equations (1.7). Direct discretizations include finite difference Chapter 1. Introduction 8 and finite volume techniques on staggered- grids (e.g. [65, 26, 66]), finite element methods using conformal and nonconformal elements (e.g. [54, 110, 63, 64]) and spectral methods (e.g. [33]). Another approach yielding many methods has involved some ini t ia l reformulation and/or regularization of the equations, to be followed by a discretization of the (hopefully) simplified system of equations. Examples of such methods include pseudo-compressibility methods, projection and pressure-Poisson re-formulations (e.g. [36, 55, 72, 97, 102, 117]). The two types of regularizations we mentioned in §1.1 for D A E s were already proposed in the Navier-Stokes context quite a while ago (cf. [108, 72]). We are interested in the generalization of the S R M to this problem because the regularized problems can be made essentially nonstiff and then a more convenient difference scheme (e.g. an explicit scheme) in time is possi-ble. Moreover, the method retains the benefits of the penalty method. For example, computations for the velocity u and the pressure p are uncoupled and an artificial boundary condition for calculating the pressure p is not necessary. 1.3 A P r o b l e m in Reservoir Simulat ion The idea of the S R M can be applied to a reservoir-simulation problem — miscible displacement in porous media. Miscible displacement is an enhanced oil-recovery process that has attracted con-siderable attention in the petroleum industry. It involves injection of a solvent at certain wells in a petroleum reservoir, wi th the intention of displacing resident oi l to other wells for production. This oil may have been left behind after primary produc-tion by reservoir pressure and secondary production by waterflooding. The economics of the process can be precarious, because the chemicals it requires are expensive and the performance of the displacement is by no means guaranteed. Complex physical behavior in the reservoir wi l l determine whether enough additional oil is recovered Chapter 1. Introduction 9 to make the expense worthwhile. A numerical simulation of the complex process undoubtedly plays an important role. Mathematically, the process is described by a convection- dominated parabolic partial differential equation for each chemical component in the system. B y summing up the component equations, one can obtain an equation that determines the pressure in the system; this nonlinear equation is elliptic or parabolic, according to whether the system is incompressible or compressible. Thus, in this problem one encounters ell iptic, parabolic, and near-hyperbolic equations with complicated nonlinear behav-ior. For simplicity, we consider the miscible displacement of one incompressible fluid by another in a porous reservoir ft C R 2 over a time period [0,T]. Let p(x, i) and u ( x , t ) denote the pressure and Darcy velocity of the fluid mixture, and c(x, t) the concentration of the invading fluid. Then the mathematical model is a strongly coupled nonlinear system of partial differential equations (see [44, 82]): u = - a ( g r a d p - 7gradcf), (x, t) G ft x [0, T], (1.8a) divu = q(x,t), (x,<) € ft x [0 , r ] , (1.8b) dc 4>— - div(D(u)gradc) + u • grade = g(c), (x, t) € ft x [0, T ] , (1.8c) wi th the boundary conditions u n = 0, (x , t ) € T x [0,T], (1.9a) £ > ( u ) g r a d c • n = 0, (x, t) € T x [0, T], (1.9b) and in i t ia l condition c(x,0) = co(x), x e f t , (1.10) where a = a(x ,c) is the mobil i ty of the fluid mixture, 7 = 7(x,c) and c?(x) are the gravity and vertical coordinate, q is the imposed external rates of flow, <£(x) is Chapter 1. Introduction 10 the porosity of the rock, D is the coefficient of molecular diffusion and mechanical dispersion of one fluid into the other, g = g(x, t,c) is a known linear function of c representing sources, and n is the exterior normal to the boundary V = <9fi. The pressure-velocity equation (1.8a)-(1.8b) is elliptic (after eliminating u). The concentration equation (1.8c) is parabolic, but normally convection- dominated. It is derived from the conservation of mass which involves the Darcy velocity of the fluid mixture, but the pressure variable does not appear in i t . Thus a good approximation of the concentration equation requires accurate solution for the velocity variables. M i x e d finite element methods have been applied to the pressure-velocity equation, which can yield velocity solutions one order more accurate than those obtained using corresponding finite difference and usual finite element methods [40, 41, 42, 46, 47, 120]. However, the finite dimensional spaces for the velocity and pressure need to satisfy the Babuska-Brezzi condition, and the resulting linear system does not have a positive definite coefficient matrix. Moreover, the number of degrees of freedom in the linear system doubles that of finite difference or finite element methods. We are interested in designing an S R M for the pressure-velocity equation since the S R M formulation can produce as accurate a velocity approximation as the mixed finite element methods, and can avoid the above-mentioned problems in mixed finite element methods. 1.4 Regular iza t ion for Differential Equations without Constraints Regularization methods are also used to treat differential equations without con-straints when these equations have some sort of singularities or their solutions may have some sort of discontinuities. Examples are viscous solutions of hyperbolic con-servation laws [76], shape from shading problems with singularities [34] and transition Chapter 1. Introduction 11 phenomena in semi-conductor device simulation [13]. A frequently considered prob-lem is the following first-order partial differential equation du du a{u,x,t)— + y,bi(u,x,t)- h c(u,x,t) = 0 (1-H) at oxi with appropriate ini t ia l and boundary conditions. In general, (1.11) may have some kind of discontinuous solutions, e.g. a shock wave [76]. Or, i f we consider the steady-state case, (1.11) may have singularities, i.e. b{ may vanish at some points, as in the shape from shading problem [34]. Some regularization techniques have been designed for solving (1.11). A popular one is du 71 du a(u, x , + ed{u, x, t)Au + ^ bi{u, x,^)^~ + c(ui x,t) = 0. (1.12) The regularized problem often has a physical meaning by itself, e.g. a time-dependent advection-diffusion equation with a small diffusion term. Another choice could be d^u du ., du ed(u,x,t)( — - Au) + a(u,x,t)— + ^k(u,x,t)— + c(u,x,t) = 0. (1-13) i=l 1 This has physical meanings as well, e.g. a traffic flow problem [118] or so-called overdamped vibration problems [104]. Thus the approximate resolution of (1.11) becomes that of the singular pertur-bation problem (1.12) or (1.13). Unlike the S R M , the regularization parameter e has to be small in comparison with the mesh size to ensure that the solution of the regularized problem be a good approximation of that of the original problem. It is well-known that there are difficulties in solving these regularized problems numeri-cally wi th small e, e.g. the stability problem for the central difference scheme and the accuracy problem for the upwinding scheme in a boundary layer region in which the derivatives of the solution may be large (see [37, 61]). We wi l l consider some special cases of (1.12) or (1.13) and focus mainly on uniformly convergent methods. We wi l l discuss spurious solutions of a simple upwinding scheme as well. Chapter 1. Introduction 12 1.5 C o n t r i b u t i o n of T h i s Thesis Our objectives are to propose and to investigate regularization methods for various differential equations with or without constraints. Most attention is paid to ordinary and partial differential equations with constraints ( D A E s and P D A E s ) . We pro-pose and analyze a regularization method called the sequential regularization method ( S R M ) . A very important advantage of our regularization method ( S R M ) is that the problem after regularization need not be stiff. Hence explicit difference schemes can be used to avoid solving nonlinear systems and they make the computation much simpler. Improvements over stabilization methods and extra benefits for P D A E s are also achieved. In Chapter 2, the S R M is presented for linear index-2 D A E s with or without constraint singularities. A complete theoretical analysis is performed for both cases. It is proved that the difference between the exact solution of a D A E and the cor-responding iterate becomes 0(es) in magnitude at the 5th iteration, at least away from the starting value of the independent variable t. Hence, the regularization pa-rameter e need not be very small so the regularized problems are less stiff. B y some choice of parameters the regularized problems can be essentially nonstiff for any e. As an example, a simple difference scheme for solving the regularized problems is investigated. Implementation techniques are discussed to get an approximation in the whole region for boundary value problems and to economize storage for in i t ia l value problems. Numerical experiments support our theoretical results. Numerical examples also show that usual stabilization methods do not work for problems wi th constraint singularities. Most parts of this chapter are taken from the paper [11]. In Chapter 3, we extend the S R M to nonlinear problems and to D A E s wi th index higher than 2. Again , nonstiffness of the regularized problems is achieved. Rather Chapter 1. Introduction 13 than having one "winning" method, this is a class of methods from which a num-ber of variants are singled out as being particularly effective methods in certain cir-cumstances of practical interest. In the case of no constraint singularity we prove convergence results. The method is also applied to constrained mult ibody systems. Numerical experiments confirm our theoretical predictions and demonstrate the via-bi l i ty of the proposed methods. Most parts of this chapter are taken from the paper [12]-In Chapter 4, we generalize the S R M to P D A E s , in particular, to the nonsta-tionary Navier-Stokes equations. The convergence rate 0(es) at the 5th iteration is again proved for this P D A E case in appropriate norms. The S R M not only avoids the stiffness of the regularized problems which always occurs in pseudo-compressibility methods but also avoids providing an unphysical pressure boundary condition which has to be imposed in stabilization methods. Discretization and implementation issues of the S R M are considered as well. In particular, a simple explicit difference scheme is analyzed and its stability is proved under the usual step size condition (independent of the regularization parameter e) of explicit schemes. The stability result also indi-cates that the step size restriction can be relaxed as the viscosity becomes small . A numerical example is calculated to demonstrate these results. The S R M formulation is new in the Navier-Stokes context and it performs well. Most parts of the chapter are taken from the paper [79]. In Chapter 5, we apply the idea of the S R M to the simulation of miscible displace-ment in porous media. The problem is modeled by a nonlinear coupled system of two partial differential equations: the pressure-velocity equation and the concentration equation. Only the approximation of velocity is important for the approximation of concentration. The S R M idea is used for the pressure-velocity equation. A n 0(es) error estimate at the sth S R M iteration is also proved. A Galerkin finite element Chapter 1. Introduction 14 method is used for the discretization of the S R M formulation. It is capable of pro-ducing as accurate a velocity approximation as the mixed finite element method. But unlike the mixed finite element methodi ts stiffness matrix is symmetric positive definite and its finite element spaces need not satisfy the so-called Babuska-Brezzi condition. Most parts of this chapter are taken from the report [82]. Chapter 6 is devoted to numerical methods of some special cases of singular per-turbation equations in the form of (1.12) or (1.13). Sections §6.1 and §6.3 describe a collection of papers [79, 115, 107] which reflects the author's earlier research interests. We believe that, by considering these special cases, we make steps towards the general problems (1.12) or (1.13) which are undoubtedly very difficult. Also, these special cases have practical meaning in themselves, hence, are worthwhile to be considered independently. In §6.1, we consider the one dimensional steady-state case of (1.12). §6.2 covers a special two-dimensional steady-state instance of (1.12) given in [121] to show that upwind schemes may lead to spurious solutions even for problems wi th very smooth solutions. We indicate that this is actually an ill-posed problem when e is small . Hence, it is not strange that a direct discretization to the problem fails. In §6.3, we consider the linear one dimensional time-dependent case of (1.13). In this case, derivatives of the reduced problem may be discontinuous along the characteristic curves. Final ly, conclusions and possible future work are contained in Chapter 7. Chapter 2 Sequential Regular iza t ion Methods for Differential Algebra ic Equat ions 2.1 M o t i v a t i o n of the S R M for General H i g h Index D A E s The sequential regularization method ( S R M ) is motivated from the augmented La -grangian method applied to constrained multibody systems (index-3 D A E s in general) by Bayo and Avello [19] and an earlier paper [20]. So we start this chapter by con-sidering a mechanical system whose configuration is characterized by the generalized coordinates q. Let L be the system Lagrangian, defined by L = T-V, where T and V are the kinetic energy and the potential energy, respectively. Let Q represent non-conservative forces. Usually the Lagrangian coordinates are not independent, but rather are interre-lated through certain constraint conditions. When the connections between bodies are of holonomic type, these constraint conditions can be expressed mathematically in the following form: = 0. (2.1) Then Hamilton's principle leads to the Euler-Lagrange equations: where A is a vector function whose components are Lagrange multipliers. For general mult ibody systems, (2.2) becomes M(q)q" + $ [ A . = f(q, q') (2.3) 15 Chapter 2. Sequential Regularization Methods for Differential Algebraic Equations!6 wi th common ini t ia l conditions, where M is the mass matrix and $ g is the Jacobian of the constraint equations. (2.3) and (2.1) form the Euler-Lagrange equations for a constrained mult ibody system. This is an index-3 D A E i f $ g has full rank. We have indicated that direct discretization would not be good in general for such a higher index D A E . One usual way to treat this problem is index reduction (to an index-1 D A E or an O D E ) . The most straightforward transformation of the D A E (2.3), (2.1) to an index-1 D A E involves replacing the constraint (2.1) wi th its second derivative plus in i t ia l conditions: _ <P*(q(t),t)_ dt2 ' (2.4a) (2.4b) (cf. (1.6)). However, this causes well-known drift difficulties, i.e. the numerical solution of (2.3), (2.4) may drift away from the original constraints (2.1) as t ime proceeds. Hence we have to look for stabilized index reduction methods. A very popular method called Baumgarte's method proposed in 1972 [17] is a generalization of (2.4). It replaces (2.4a) with the equation $ " + a $ ' + 6$ = 0, where a and b are parameters chosen so that the roots of the polynomial < T ( T ) = T 2 + ar + b = 0, (2.5) (2.6) are both negative, i.e. the ini t ia l value problem for the differential equation (2.5) for $ is asymptotically stable (see [16]). The system (2.5), (2.3) can be written in the form " M ^ ' _ $ 9 0 A The matr ix 0 (2.8) Chapter 2. Sequential Regularization Methods for Differential Algebraic Equations!! Figure 2.1: planar slider-crank: ini t ia l state in solid line, subsequent states in dotted lines is nonsingular i f $ g has full row rank. Hence (2.7) can be integrated using standard numerical integrators. If $ 9 is rank-deficient, we have a potential difficulty in solving (2.7). Baumgarte's method may not work then. The problem is called singular i f $ ? does not have full rank. E x a m p l e 2.1 Consider two linked bars (see Fig. 2.1). One end of one bar is fixed at the origin, allowing only rotational motion in the plane. The other end of the other bar slides on the x-axis. The equations of motion form a nonlinear index-3 DAE p' = v M i / = / — GT X g(p) = 0 where Xj,yi,<pi are the coordinates of the center of mass of the ith bar, and Chapter 2. Sequential Regularization Methods for Differential Algebraic Equationsl8 P = (x1,yu4>ux2,y2,<h)T-If the left bar is strictly shorter than the right bar, then the Jacobian matrix G of the constraint functions of this problem has full rank. The problem is nonsingular. If the length of these two bars are the same, for example, each with length 2 and mass 1, then we have M = diag{l, 1,1/3,1,1,1/3} / = (0 , -9 .81 ,0 ,0 , - 9 . 81 , 0 ) T 9 = ( Xi — COS 0 1 \ t/i - s i n 0 i x2 — 2xi — COS 02 2/2 - 2y! - sin 0 2 V 2 sin 0 i + 2 sin 0 2 ) G - gv = ( 1 0 s in0 ! 0 0 0 0 1 - c o s 0 i 0 0 0 - 2 0 0 1 0 sin</>2 0 - 2 0 0 1 - c o s 02 V 0 0 2 cos 0! 0 0 2 cos 02 / Clearly, as the mechanism moves left through the point where both bars are upright ((f>i — I", 0 2 = the last row of G vanishes at this one point and a singularity is obtained. When arriving at this point with no momentum, this is actually a bifurcation point where two subsequent motion configurations are possible. We will consider only the case where the sliding bar continues to slide along the x-axis past the singularity, and note that the solution is smooth in the passage through the singularity. • In [19], Bayo and Avello proposed to solve the multibody system (2.3) using an augmented Lagrangian algorithm which is transplanted from the same method in the optimization context [5]. Their idea is to derive a modified formulation by adding to the expression of Hamilton's principle three terms: Chapter 2. Sequential Regularization Methods for Differential Algebraic Equationsl9 • a fictitious potential v* = Y, \«k<4*l • a set of Rayleigh dissipative forces Gk = -2akukuk$'k • a fictitious kinematic energy term T* = Y,\^t (2-11) k z where each ak is a very large real number (the penalty), and Lok and pk represent the natural frequency and the damping ratio of the penalized system (mass, dashpot and spring) corresponding to the constraint = 0. Then we get a modified Euler-Lagrange equation Mq" + $ J a ( $ " + 2 f y i $ ' + 0 2 $ ) + $ JA* = f(q, q') (2.12) or ( M + ^a$q)q" + $ f A* = f(q, q') - $?a((**)'q' + (*«)' + 2 0 / / $ ' + 0 2 $ ) , (2.13) where a, fl and p are diagonal matrices that contain the values of the penalties, the natural frequencies and the damping ratios of the fictitious penalty systems assigned to each constraint condition. Because A* is not given in advance Bayo, Jalon and Serna [20] propose an iteration to solve (2.12) or (2.13) by comparing (2.12) wi th (2.3): A: = A ; _ ! + a ( $ " + 2fy*$' + ft^U^, s = 1, 2, • • • (2.14) to get an approximation of A. In [19] the authors called the whole procedure an augmented Lagrangian algorithm and claimed that the algorithm works well, however without any theoretical analysis. (2.9) (2.10) Chapter 2. Sequential Regularization Methods for Differential Algebraic Equations20 In fact, the algorithm would certainly not work in general. In mult ibody motions q usually remains smooth even in passage through singular positions. However, as indi-cated in later sections, A may be unbounded at the singularity. So the iteration (2.14), as an approximate procedure to obtain A, is not appropriate in some cases. Also, should not be included in (2.12) in the singular case because unbounded coefficients may appear in it . From (2.14) the formulation corresponds to Baumgarte's formula-tion since when A* gets close to A * _ : (2.14) gets close to Baumgarte's stabilization. But Baumgarte's stabilization is not as good as some other stabilization techniques (see [8]). Moreover, in [19] the authors indicated that the iteration (2.14) is applied unt i l \\q" — q"_i\\ < 8, where 8 is a user-specified tolerance. This criterion is perhaps applied because they did not know the convergence rate of their iterative procedure. We do not recommend this criterion because it causes not only unnecessary extra iterations but also makes a storage-saving implementation difficult (cf. §2.6). Our a im is to construct a method for general D A E s which not only avoids the above shortcomings, but for which we also do not have a Lagrangian and Hamilton's principle. So another derivation of the algorithm is needed. In this and the next chapter, we propose a class of algorithms motivated by the augmented Lagrangian method for more general D A E s of order u, x(u) = f(x,x',...,x{,/-l\t) - B(x,t)y, (2.15a) 0 = g(x,t). (2.15b) The D A E (2.15) has index u + 1 if GB is nonsingular for a l H , 0 < t < tj, where G = gx. We are interested in the cases v = 1 or 2. The Euler-Lagrange equations (2.3) for mechanical systems with holonomic constraints are in this form with u = 2. The algorithm is derived by combining a modified penalty idea (a k ind of regularization) given in [91] wi th stabilization techniques such as Baumgarte's stabilization or the stabilization analyzed in [8, 35] in an iterative procedure. We call the method the Chapter 2. Sequential Regularization Methods for Differential Algebraic Equations21 sequential regularization method ( S R M ) (cf. [11, 12]). The method is applicable for more general higher index D A E s . More importantly, it works for both boundary and in i t ia l value problems and is justified by a theoretical analysis. The number of iterations can be determined beforehand depending on the penalty parameter a = 1/e, the mesh size h and the order of the method used. Since we specify the iteration number in advance we can design a procedure for the ini t ia l value case to perform the iteration without the need to store al l previous approximate values. The sequential regularization method is actually a functional iteration procedure in which the difference between the exact solution of a D A E and the corresponding iterate becomes 0(es) in magnitude at the 5th iteration, at least away from the starting value of the independent variable (which we shall call ' t ime') . Hence, unlike the usual regularization, the perturbation parameter e does not have to be chosen very small, so the regularized problems can be less stiff and/or more stable. Next we wi l l propose and analyze the S R M for the linear index-2 case wi th sin-gularities. Numerical experiments are given to verify our theoretical results. Some simple difference schemes for the regularized problems and implementation issues for the S R M are also discussed. 2.2 L inear Index-2 Problems We first write down the linear index-2 problem: where A(t), B(t) and G(t) are sufficiently smooth functions of t,0 < t < tj, A(t) £ R" l X n i , B(t) e R n* x n«, G{t) € R^x™*, and ny < nx. We consider the D A E (2.16) subject to nx — ny boundary conditions x 0 A(t)x - B{t)y + q(t) (2.16a) (2* 16b) Chapter 2. Sequential Regularization Methods for Differential Algebraic Equations22 Box{0) + B1x{tf) = (3. (2.17) These boundary conditions are assumed to be such that they yield a unique solution for the O D E (2.16a) on the manifold given by (2.16b). In particular, i f we were to replace (2.16b) by its differentiated form 0 = Gx' + G'x + r'= -^g(x,t), (2.18a) g(x(0), 0) = G(0)x(0) + r(0) = 0, (2.18b) and use (2.18a) in (2.16a) to eliminate?/ and obtain nx O D E s for x, then the boundary value problem for x wi th (2.17) and (2.18b) specified has a unique solution. In the ini t ia l value case B\ = 0, this means that (2.17) and (2.18b) can be solved uniquely for x(0). We wi l l give a more precise assumption in Lemma 2.1 below. The problem (2.16), (2.17) is of index-2 if GB is nonsingular for all t. However, here we allow GB to be singular at some isolated points of t. For simplicity of exposition, let us say that there is one singular point £*,0 <t*<tf. The inhomogeneities are q(t) £ R71* and r(t) £ Rny. We are only interested in the kind of singularities as in Example 2.1, where the solution x of (2.16), (2.17) passes through the singularity smoothly. Returning to Example 2.1 (where B = M~lGT), we can verify that, although the matr ix GM~lGT is singular at the singularity , the matrix M~1GT(GM~1GT)~lG is smooth for al l t. Also, two types of singular constraints (i.e., wi th vanishing rows or wi th some rows linearly dependent at some points) mentioned in [2] both have a similar property. Thus, for the linear model (2.16), we assume accordingly: Assumpt ion 2.1 The matrix function P — B(GB)~1G is smooth, or more precisely, P is continuous and P' is bounded near the singular point i * where we define P(U) = )im(B(GB)-1G){t). Chapter 2. Sequential Regularization Methods for Differential Algebraic Equations2d Because we are only interested in the case where (2.16) has a smooth solution for x (as is the case in Example 2.1), it is necessary to assume, in view of (2.16b): Assumpt ion 2.2 The inhomogeneity r{t) satisfies r £ S, where S = {w(t) £ Rny : Gz = w for a smooth function z{t). We note that Assumptions 2.1 and-2.2 are satisfied automatically i f GB is nonsin-gular for each t. O n the other hand, neither B(GB)~l nor (GB)~1G alone are smooth near a singularity in general. We also indicate here that to formulate the S R M (see §2.4) we only need the continuity of P. The further requirement in Assumption 2.1 on the derivative of P is needed for the regularity of the solution (cf. Lemma 2.1 and (2.24)) and the stability proof for the regularized problems (cf. §2.5). This require-ment can be avoided i f we make a more general assumption about the regularity and stability of the original problem (cf. §3.1). We consider both ini t ia l and boundary value problems. In §2.3 we briefly dis-cuss the conditioning of the problem (2.16) with singularities. In §2.4 we derive the sequential regularization method. In §2.5 we estimate the error of the S R M . In §2.6, we consider some discretization and implementation issues for both in i t ia l and boundary value problems. Final ly, in §2.7 several numerical examples demonstrate our theoretical results. 2.3 P r o b l e m Condi t ion ing Similar ly to [15] and to the method of pseudo upper triangular decomposition ( P U T D ) described in [2] (cf. §10.6; with the difference that we do pivoting to interchange the row with the singularity of the lowest order and the current row when al l the Chapter 2. Sequential Regularization Methods for Differential Algebraic Equations24 other rows vanish at some singular point), there exists a smooth matr ix function R(t) G RK-"»)x n * , which has full row rank and satisfies RB = 0, for each i, 0 < t < tf, where R can be taken to have orthonormal rows. As in [15, 16], define the new variable u = Rx, 0 < t <tf. • (2.19) Then, using (2.16b), the inverse transformation is given by x = Su- B{GB)~lr, (2.20) where S = (/ - B(GB)-1G)RT = (I- P)RT. B y the assumptions at the beginning of this chapter, this transformation is well-defined. Differentiating (2.19) and using (2.16a) and (2.20) we obtain the essential underlying O D E ( E U O D E ) : u' = (RA + R')Su - (RA + R')B(GB)-1r + Rq. (2.21) Hence the underlying problem of (2.21) is v! = (RA + R')Su + / , (2.22a) BoS(0)u(0) + 5i5(t/)w(</) = Pi. (2.22b) We make Assumpt ion 2.3 The boundary value problem (2.22) is stable, i.e. there exists a moderate-size constant K such that H < f f ( | | / | | + |)9i|), where \\u\\ = maxt{\u(t)\, 0 < t < tf}. Chapter 2. Sequential Regularization Methods for Differential Algebraic Equations25 Similar ly to Theorem 2.2 of [15], we have L e m m a 2.1 Let the DAE (2.16) have smooth coefficients, and assume that Assump-tions 2.1 and 2.2 hold. If the EUODE (2.21) with boundary conditions (2.22b) has a unique solution, then there exists a unique solution for the x of problem (2.16)-(2.17) which is smooth. This implies the unique existence of a smooth By as well. Furthermore, if Assumption 2.3 holds then there is a constant K such that \\x\\<K{\\q\\ + \\B{GB)-*r\\ + 1)91), < K(\\q\\ + \\B(GB)-lr\\ + | | ( J B ( G ' J B)- 1r) ,|| + \3\). R e m a r k 2.1 For problem (2.16) without singularities, we can get \\x\\ <K{\\q\\ + \\r\\ + \(3\), Hz'H < K(\\q\\ + ||r|| + ||r'|| + |/?|) as in [15, 16]. • The difference between the situation here and in the nonsingular case is that the perturbation inhomogeneities r yield reasonably bounded perturbations in the solution x only if they are (in general) from the subspace Range (G). From (2.16a) and (2.20), we can write y = -(GB)-xG(x' -Ax- q),t e [0,Q U (U,tf], (2.23) which could be unbounded at the singular point £* (whereas By is bounded). Note that G could be replaced in (2.23) by any appropriate matrix Q wi th the same size as G, e.g. Q can be BT. Chapter 2. Sequential Regularization Methods for Differential Algebraic Equations2d> R e m a r k 2.2 If B has full rank for each t, then (GB)~lG = {BTB)-lBTP. Hence, (GB)~lG is smooth. Hence, there exists a unique solution for the y of problem (2.16)-(2.17) which is smooth and can be expressed as (2.23) for each t. Furthermore, using Lemma 2.1, we have in this case \\y\\ < K(\\q\\ + \\B(GB)-lr\\ + | | (5 (C7i?) -V) ' | | + |/?|). (2.24) In the general case, however, we will have to consider By, rather than y alone, in the theorems of the next section. • A Baumgarte stabilization applied to (2.16) consists of eliminating y according to (2.18),(2.23), and stabilizing (see (2.30) below for the usual form). This gives the O D E x'=(I- B{GB)-lG){Ax + q)- B(GB)-l(G'x + r') - e"1 B(GB)-\Gx + r) (2.25) where e > 0 is a parameter (cf. [17, 8]). If there are no singularities then it fol-lows from the analysis in [16] that i f Assumption 2.3 holds then the boundary value problem (2.25),(2.17),(2.18b) is also stable. In other words, the "ini t ia l value stabi-l ization" works also for the boundary value case, because the new modes introduced by replacing (2.16b) with (2.18a) are separable and decaying, in agreement wi th the additional initial conditions (2.18b). However, in the singular case (2.25) may not work because the terms B(GB)~XG' and B(GB)~1r' are in general unbounded. Therefore, we develop an iterative method in the next section which builds up an approximation to By and x that avoids going through unbounded quantities. Chapter 2. Sequential Regularization Methods for Differential Algebraic Equations27 2.4 Der iva t ion of the S R M There are many discussions on regularization methods for D A E s . A direct regulariza-tion (cf. the pseudo-compressibility method in the Navier-Stokes context [108, 72, 97]) is: x' = A(t)x - B(t)y + q(t), " (2.26a) -ty' = G(t)x + r(t). (2.26b) This formulation is not popular because it requires conditions on A , B and G , for the purpose of stability of the system, and the existence of the first derivative of y, which is not necessarily true for the original problem (2.16) [86]. It may also change the properties of the original index-2 problem too much by jumping from index-2 to index-0 ( O D E ) . It seems evident that a regularization with fewer changes of the original problem (e.g. from index-2 to index-1) might be better. The penalty method [84, 91, 68, 70] is such a method. It reads x' = A(t)x - B(t)y + q(t), (2.27a) tE~ly = G{t)x + r(t), (2.27b) where E £ ~R,nyXnv is chosen such that BEG has non-negative eigenvalues. Hence, the system obtained by substituting (2.27b) into (2.27a) is generally stable. For example, we can choose, relying on Assumption 2.1, E — (GB)-1 (hence, BEG — P). Also , E = (GB)T could be a good choice in some circumstances. If B = M~lGT for some positive definite matrix M (as in the case of mechanical systems) then it is possible to choose E = I. Advantages of these choices of E w i l l be discussed in Chapter 3. For problems wi th singularities, we suggest using E = (GB)'1 to avoid a turning point problem. Another approach is parameterization [85, 59]: = A(t)x- B(t)y + q(t), (2.28a) Chapter 2. Sequential Regularization Methods for Differential Algebraic Equations28 0 = G{t)(x + ex') + r(t). (2.28b) Substituting (2.28a) to (2.28b), we get eGBy = Gx + r + eGAx. (2.29) This implies that parameterization can be seen as an instance of the penalty method wi th E = (GB)-1. Recently, [94] has reported a regularization for D A E s wi th sin-gularities based on a formulation obtained by the trust-region method in numerical optimization. A l l these regularizations require the regularization parameter e to be very small . Therefore a stiff solver is needed to solve the regularized problem. In this section, we derive a new regularization method which is called the sequential regular-ization method [11]. The S R M is an iterative procedure which combines the popular Baumgarte stabilization or other stabilizations with a modified penalty method. One purpose for doing so is that the regularized problems of the S R M can be non-stiff or at least less stiff. Hence, simple discrete schemes (e.g. explicit schemes) can be used. Other advantages of the method wi l l be discussed in Chapter 3. The Baumgarte stabilization of (2.16) reads (cf. (2.5)) a1^g(x,t) + a2g(x,t) = 0, g(x(0),0) = 0. (2.30) Apply ing the idea of the penalty method to equations (2.16a) with constraints (2.30), we obtain x' = A{t)x- B(t)y + q(t), (2.31a) y = y 0 + x0( x ><) + " 2 0 ( M ) ) - (2.31b) e dt where yo can be seen as an ini t ia l guess of the exact solution ye of problem (2.16), (2.17). If we take y 0 = He then the solution of problem (2.31), (2.17) is exactly the Chapter 2. Sequential Regularization Methods for Differential Algebraic Equations29 same as that of problem (2.16), (2.17). If we take yo = 0 then (2.31) coincides wi th the penalty method (2.27). Given any ini t ia l guess yo(t), the solution, say {a;i, y i } , of (2.31), (2.17) is an approximation of the exact solution {xe, ye} of (2.16), (2.17). Using this solution y\ as a new ini t ia l guess, we re-solve problem (2.31), (2.17). We expect that the solution obtained is a better approximation of the exact solution. Repeating the procedure, we invent the following iterative algorithm for solving (2.16): For 5 = 1 , 2 , . . . , solve the ODE problem xj = Axs -Bys + q (2.32) where ys = y . - i + ^E(cti^g(x„t) + a2g(xait)), (2.33) subject to the same boundary conditions (2.17). Note that yo(t) is a given ini t ia l iterate and that e > 0 is the regularization parameter. We call this algorithm a sequential regularization method ( S R M ) . Note that xs(t) and ys(t) are defined on the entire interval [0,tj] for each s. For the problem with singularities, the choice E = (GB)T generates turning point regularized problems which are complicated to solve and analyze. We thus choose E = (GB)-1 for the singular case. Noting that B(GB)~1G' may be unbounded at the singularity we then choose cti = 0 to avoid this term. Also, in practice we mult iply (2.33) by B and keep track only of the approximations Bys to the bounded function By, since y may be unbounded at the singularity. We thus have an S R M variant for the singular case: For 5 = 1 , 2 , . . . , solve the ODE problem xj = Axs -Bys + q (2.34) where Bys = Bys-i + -BEg(xs,t), (2.35) Chapter 2. Sequential Regularization Methods for Differential Algebraic Equations30 subject to the boundary conditions (2.17). If y is desired (at times other than at the singular point £*) then it can be easily retrieved from By in a post-processing step. 2.5 Convergence Analys is of the S R M We first prove a lemma which wi l l be used to discuss the convergence of the S R M . L e m m a 2.2 Let u, v be the solution of v! = (RA + R')Su + Siv + / i , (2.36a) Sv' + 7U = eS^u + eS3v + / 2 , (2.36b) BoS{0)u(0) + BiS(tf)u(tf) = '3 - 5 4 u ( 0 ) -S5v(tf), v(0) = v0, (2.36c) where all coefficients are bounded, 8 = 1 or S = e, 7 is a positive constant and Assumption 2.3 holds. Then, for e appropriately small or 7 appropriately large, we have the following stability inequality \\u\\ < ^ ( | | / i | | + | | / 2 | | + |/9| + K I ) , IHI</ir(£| | / i | | + | | / 2 | | + |)9| + |t;o|), where K is a positive constant. Proof: Let v = (ui, • • • ,vny)T. From (2.36b), we easily have \vi\ < -II^ IIHull + -||S,3||||'u|| + -H/2II + |u0|, i = 1, • • • ,ny. 7 7 7 Hence, taking the maximum of the left hand side for 1 < i < ny and choosing small e or large 7 appropriately such that e\\S3\\ < 7, we get \\v\\< j T ^ P z H + i , g I, • (2-37) 7 - e S 3 7 - e 6 3 Chapter 2. Sequential Regularization Methods for Differential Algebraic EquationsSl B y using Assumption 2.3, from (2.36a), there exists a positive constant K\ such that ||u|| < ^ idl^llUHI + HAH +1/?| + |5 4 |K0)| + I S B I K M ) < A-iail^H + iSsDiiuii + ii/iii + ^ i + ^ ^W) < A^edl^H + |gB|)||ga|| ^(115x11 + |56|)(||/a|| + 7kol) | K r N f ii , i f l . , . o M f > n ^IJH I H I + ^4s7\\ +A i (||/ i||+|0|+|5 4 || t , o | ) . Hence, by choosing smaller e or larger 7 such that A i e ^ l ^ | ^ ^ 5 2 ^ < 1, the stability inequality for u follows. Now the stability inequality for v follows from that for u and (2.37). • Now we estimate the error of the S R M (2.34), (2.35). Def ini t ion 2.1 J is an integer such that yo(0) = y e(0), y o(0) = y'e(0),..., y<J)(0) = y(J)(0), where yo(t) is the initial guess of the SRM iteration (2.34), (2.35) and ye(t) is the exact solution of the original problem (2.16). Set J = — 1 if yo(0) 7^ 2/e(0). • For in i t ia l value problems we may calculate y e ^(0) , i = 0 ,1 , . . . in advance by using the O D E and its derivatives. For boundary value problems we have J = — 1 in general since we don't know ye(0) beforehand. Theorem 2.1 Let the DAE (2.16) have sufficiently smooth coefficients, and assume that Assumptions 2.1, 2.2 and 2.3 hold. Then, for the solution of iteration (2.34), (2.35), we have the following error estimates (for J defined in Definition 2.1): xs(t)-xe(t) = 0(es) + 0(e J +V('AK i / e), Bys(t) - Bye(t) = 0(e°) + 0(eJ+1pa{t/e)e-^), for 0 < t < tf and s > 1. Here ps(r) = 0 if s < J +1; otherwise Ps(T) is a polynomial of degree s — J — 2 with generic positive coefficients and |p s(0)| = \(Byo)(J+1^(0) — (Bye)(J^(0)\. Chapter 2. Sequential Regularization Methods for Differential Algebraic Equations32 Proof: Let us = Rxs and ws = Pxs. Similarly to (2.20), we have xs = Sus + ws. (2.38) Furthermore, using (2.34) we obtain < = (RA + R')Sus + (RA + R')ws + Rq, (2.39a) ew's + ws = e(PA + P')Sus + e(PA + P')ws - e B y , _ i (2.39b) + 'ePq-B(GB)~1r, subject to BoS(0)us(0) + BiSitrfiisitf) = [3- Bows(0)- BlWs(tf), (2.40a) w3(0) = - J B(0)(C?(0)5(0)) - 1 r (0) . (2.40b) The iteration (2.35) for By becomes By, = Bys_1 + ]-(ws + B(GB)-lr). (2.41) The proof proceeds along familiar lines of singular perturbation analysis. Accord-ing to [111, 112] we can construct the asymptotic expansion of ws and us sequentially for s = 1,2, . . . , where we use Lemma 2.2 to estimate the remainders. Then, us-ing (2.41) and (2.38), we get the asymptotic expansion of Bys and xs respectively. Note that in these expansions the first terms are exactly xe and Bye. This process eventually yields the proof of the theorem. • To provide a better understanding about the sequential regularization method we give in §2.8 a detailed proof for the ini t ia l value case with no layers, s < J + 1. In that proof, the construction of the asymptotic expansion is directly for x and By. Moreover, the construction method we apply is somewhat different from [111, 112] and more relevant to the concept of D A E s . Chapter 2. Sequential Regularization Methods for Differential Algebraic Equations33 Next, we consider the S R M (2.32), (2.33). For the ini t ia l value problem with E = I and B = GT , this corresponds to Algor i thm (2.13), (2.14) of [19] for constrained mechanical systems (although they do it for the corresponding index-3 case) derived by a penalty-augmented Lagrangian formulation. Bayo and Avel lo have noted that under repetitive singular conditions this algorithm may lead to unstable behavior. For our index-2 case (2.16) with singularity, it appears to be impossible to choose a matr ix E such that problem (2.32),(2.33) is always stable, even i f we assume B = GT. A numerical example in §2.7 wi l l verify such instability phenomena even for the case of one singular point. However, for the case where constraints are without singularities, (2.33) (multiplied by B) may have a benefit over (2.35). That is, (2.33) yields an O D E problem for xs which is essentially not a stiff problem. Take E = (GB)~l as before and rewrite (2.33) as Bys = Bys-! + ^BE{al^g{x„,t) + a2g(xs,t)). (2.42) Then we give the following error estimation for (2.32), (2.42): Theorem 2.2 Let the DAE (2.16) have sufficiently smooth coefficients, and assume that G has full rank and that Assumptions 2.1, 2.2 and 2.3 hold. Then for the solution of the iteration procedure (2.32), (2.42) with a\ ^ 0, we have the following error estimates: xs - xe = 0(es), Bys - Bye = 0{es) for 0 < t < tf and s = 1 ,2, . . . . Note that no boundary layer terms appear here even for J = — 1 (See Definition 2.1)1 Proof: Denote us = Rxs and vs = Gxs. Hence xs = Sus + Fvs, (2.43) Chapter 2. Sequential Regularization Methods for Differential Algebraic Equations34 where S = (I - P)RT and F = B(GB)~1 = PGT(GGT)~1 are both sufficiently smooth. From (2.32),(2.42), we get u's = (RA + R')Sus + (RA + R')Fvs + Rq, (e + ax)v'a + a2vs = e(G' + GA)Sus + e(G' + GA)Fvs + eGBys-X + tGq - axr' - a2r, with the corresponding boundary conditions, and Bys = 5y s_! + -BiGBY^a^v. + r)' + a2(vs + r)). Repeating the procedure of the proof of Theorem 2.1 and using Lemma 2.2 again to estimate the remainder of the asymptotic expansion, we obtain us - ue = 0(es), vs-ve = 0 ( e s ) , Bys-Bye = 0(es), where ue = Rxe, ve = Gxe = —r. Hence, using (2.43) and xe = Sue + Fve, we obtain xs - xe = S(us - ue) + F(vs - ve) = 0(es). • R e m a r k 2.3 For the problem (2.32), (2.33) where GB is nonsingular, we have xs -xe = 0(es), ys-ye = 0(es). This estimate also holds for E = (GB)T. • 2.6 Discre t iza t ion and Implementat ion Issues The S R M iteration yields a sequence of O D E problems which are to be solved numer-ically. We only consider the most difficult case, i.e. (2.34), (2.35) wi th singularities. Chapter 2. Sequential Regularization Methods for Differential Algebraic Equations35 Inserting (2.35) into (2.34), the O D E problem to be solved at the 5th iteration is writ ten as the singular-singularly-perturbed problem (see [112, 89]) ex's + BE(Gxs + r) = eAx, - e ( £ y s _ i - q), (2.44a) Boxs{0) +B1xs(tf) = (3 , G ( 0 ) x s ( 0 ) + r ( 0 ) = 0. (2.44b) We consider finite difference (or collocation) discretizations of (2.44) on a mesh 7T : 0 = t0 < tx < ... < tN = tf hi = t{ — U-i, h = max hi l<i<N and denote by x*, y* the corresponding approximations of xs(ti), ys(ti), respectively. We now have essentially two small, positive parameters to choose: e and h. We assume that h is chosen small enough so that the E U O D E problem (2.22) may be considered as nonstiff and that the problem's coefficients are sufficiently smooth. In the B V P case the situation is the familiar one, much like the iterative solution of a nonlinear boundary value O D E using quasilinearization (see, e.g., [14]). Each of the linear boundary value O D E s (2.44) is discretized on a mesh n using, say, a symmetric finite difference scheme or some other method. We expect, as h —> 0, convergence to the solution of ( 2.44) and our theory then applies for the entire numerical algorithm. As an example, we give here a detailed analysis of the convergence of the backward Euler difference scheme for (2.44). A similar discussion and results can easily apply to the forward Euler difference scheme. The results for general higher order collocation schemes have been described in [11]. Now we write the backward Euler scheme of (2.44) as follows: e X l ~ h U +BMGix\ + n) = ^\-e{Biy9-X{U)-q^ (2.45a) B0xs0 + BlXsN = (3, G0xs0 + r 0 = 0, (2.45b) Biyl = BiVs-iiU) + -B%Ei{GiXst + n), (2.46) Chapter 2. Sequential Regularization Methods for Differential Algebraic Equations36 where we represent the value / (£ ; ) of a given function / at mesh point t{ by Mul t i p ly ing (2.45) by R{ and Pt, respectively, and denoting u\ = RlXst, w\ = Pix\ (then x* = SiU* + ws{ since (2.38) holds), we have ^rf^ = RiAiSiu'i + wt) + t ^ ^ i S ^ u U + wU) + R i q i , (2.47a) + w' = CRMS*! + < ) + c^rhiSi-iuU + wU) -e(5,-2/s_i(*,-) - qi) - BiEin, % = 1, • • • , A , (2.47b) Bo(S0u'0 + wa0) + B^SNU'N + wsN) = 8, ws0 = -B0E0r0 = P0xe(t0). (2.47c) A t first, we consider the stability of the following difference scheme corresponding to (2.47): rh _ Ui - Ui-1 r> A C Ri ~ R{-\ L/1Ui = : KiAibiUi : = UiAiWi hi h{ + * v > i - i + fi, (2.48a) . _ Wj - U>,-_i Pj - Pj-x L2wt — e 1- W{ = ePiAibiUi + e : i > i _ i u t _ i hi hi p. _ p. ePzAtwt + e 1 8 - 1 Wi-! + g{, i = 1, • • • , A , (2.48b) I hi B0S0u0 + BXSNUN = 8I = 3 - B0w0 - B\wN, w0 = B2. (2.48c) Using the discrete maximum principle for the difference operator L\, i.e. z0 > 0 and L2z{ > 0 zt > 0 , V i , (2.49) we easily get max \WJ\ < M(\32\ + max \LowA). (2.50) i < i < » Defining \\z\\i — maxi<j<; \ZJ\ and using (2.48b), we have H I . - ^ M c C i H i i + l H I O + ^ l + I M I i Chapter 2. Sequential Regularization Methods for Differential Algebraic Equations^! or \\w\\i < M(e\\u\\i + \32\ + \\g\\i). (2.51) Here M stands for a generic constant independent of i, e and h. O n the other hand, L\u{ is a one-step difference operator of (2.22a) — the underlying problem of (2.21). Using Assumption 2.3 and Theorem 5.38 of [14], we obtain H o o < Ki{\Pi\ + max \LhlUj\), (2.52) l<j<iV where ||^||oo = rnax0<j<iv \ZJ\ and K\ = K + 0(h). Here K is the stability constant defined in Assumption 2.3. Using (2.48a) and (2.48c) we have WuWcc^ M(\\w\\N + \p2\ + \p\ + \\f\\N). (2.53) Then, using (2.51), yields H o o < M(\\f\\N + \\g\\N + \\3\ + \fo\) (2.54a) I I H U < M(e\\f\\N + \\g\\N + \\P\ + \(32\). (2.54b) We thus obtain the stability inequalities (2.54) for the difference scheme (2.48) or (2.47). Now we discuss the convergence of (2.47). Using (2.54) , we have the estimates: K ( * , - ) - < l l o o < M ( | | r u | | i V + | | r w | | i v) (2.55a) H ( i , - ) - « > , ? l | o o < M(\\TU\\N + \\TW\\N), (2.55b) where r " and T-" are local truncation errors for the difference scheme (2.47) and they can be writ ten as r- = ^{<(e!) + i?"(e - ) [5(^-i)^^-i) + ^(^-i)] + m(ti)(Sua + wsy^t)} (2.56a) T? = e / i , - K ( ^ ) + P , ,(ni !)[5(t i_ 1)U a(< 1-_i) + U; 8(t <_ 1)] + P'(tt)(Sus + ws)'(rif)}, (2.56b) Chapter 2. Sequential Regularization Methods for Differential Algebraic Equations38 where £f,?7f £ ( t ,_i , t t ) , /^ = 1,2,3. To bound the truncation error, we need the derivative estimates of us and ws. From the asymptotic expansions of us and ius ( cf. Theorem 2.1) us = ue + eusi + e2(us2 + S s 2 ) H (2.57a) ws = we + e(wai + wsl) + e2(ws2 + ws2) H (2.57b) where ue = Rxe, we = Pxe , usj and wsj are functions of regular expansions, usj and w„j are boundary layer functions whose basic forms are p(t/e) exp(—t/e) (where p is some polynomial) and usj = wsj = 0 for j < s ( since S R M iteration cancels out the lower terms of regular expansions). We can expect that K | , K | , K | <M. (2.58) But Therefore K | = 0(e-> exp(-t/C) < MM_X i f ' 6 (2.59) M e otherwise . |r-|| W = 0 (*) , l |r - | |= ^ if « < * . = « . - * > ( 2 . 6 0 ) I C(ftJ otherwise . From (2.55), we thus have l k ( * , - ) - " . l o o < Mh . (2.61a) |^.(*t) - < l | o o < \ . (2.61b) Mft otherwise i.e. x,(*,-) - x'i = S(ti)(ua{ti) - <) + (wa(ti) - w'i) = 0(h). (2.62) However, we can not generally get a good approximation for Bye in the whole region if e is not very small compared with hi since in this case we generally have Bivi = Biyt-1(U) + -BiEi{Gix'i+ri) = Biy^t^ + ^BiEiiGiX^ + r^ + ^ w'i-Wsiti)) (2.63) = Biye(ti) + 0(e + exp(-U/e)) + 0{h/e). Chapter 2. Sequential Regularization Methods for Differential Algebraic Equations39 Fortunately, we can get 0(h) accuracy locally, i.e. in a smooth region or away from the layer region, say for ti0 <t<tj. Indeed, considering an equidistant mesh for simplicity, from (2.47b), we have L>Awt = e A w ° - ^ + A w l l = e P i A,-(5,-A«? + Aw?) + / ' " ^ " ' ( f l - i A t z ? + AwU) + 0 ( 7 f ) , (2.64) where Aws{ = ws(ti) — w°, A m - = us(ti) — u\ and we note that = 0(th) for io < i < N. Using the discrete maximum principle for L\ in £,•„. < t < tj on the barrier function Zi = [Aw" | A J - i o + max |<Jt-| ± Aw? 1 l°' i0<i<N 1 1 1 where <5;(= 0(h)) is the right-hand side of (2.64) and A = h/(e + h), and using (2.61), we get Zi > 0 or \Aw*\ < | A < | A ' - " ° + max \Si\ < M(hX-io + eh). (2.65) For e = h1+5, - 1 < 5 < 1, A 1 ' - 1 ' 0 = e x p ( - ( i - z 0 ) / i 5 ) = e x p ( - ( « , - - t^/h1-*). Taking i i such that e x p ( - ( « , - 1 -U^/h1'5) < Mh1+S = Me when h is sufficiently small , we get A'-*'0 < e x p ( - ( t i l - UA/h1'8) < Me for i > iu i.e. \Awf | < Mth, Vz > z'i, and any e satisfying e = h1+5, - 1 < <5 < 1. Combining this wi th (2.61b), we obtain |Atu?| < Meh, Vz > ti. (2.66) Then, using (2.63) and (2.66), we have the local error estimate \Biyl - BlVe(U)\ < M(h + e s), for i, < i < N. (2.67) Chapter 2. Sequential Regularization Methods for Differential Algebraic EquationsAO This means we can get good approximation in a region which is away from the in i t ia l layer. Once an accurate S R M solution, say {x*, B(t{)y*}, has been determined outside the in i t ia l layer it may be possible to obtain an accurate solution everywhere by ap-plying a few S R M iterations numerically solving (2.44a) ( changing BEG to —BEG) subject to the terminal value x{tf) = x*N,- (2.68) and choosing By0 satisfying B(tf)y0(tf) = B(tj)y*N. This procedure is feasible pro-vided that the terminal value problem (2.44a),(2.68) is well-conditioned (which holds if the terminal value problem for the E U O D E (2.22a) is well-conditioned). For the I V P case, where (2.44b) reduces to x(0) = x given, (2.69) we may, of course, proceed in the same way as for the B V P case. But now a few things are easier. First ly, for this case we can calculate Bye(0) and then choose Byo to be exact at t = 0. In fact, as indicated earlier we can also do this for higher derivatives of By at the ini t ia l value by repeated differentiation of (2.16). Such a preparation of the in i t ia l iterate Byo allows removing the layer error terms (or the condition t{ 3> e) in the error estimates (2.61). Secondly, one may use a more convenient difference scheme to integrate the I V P (2.44a),(2.69). If the E U O D E is sufficiently nonstiff to warrant use of a nonstiff integration method then this can be an attractive possibility here. Note, though, that —hi/e must be in the absolute stability region of the method (see (2.39b)). Thus, an explicit Runge-Kutta method of order p, for instance, may necessitate (at least) p S R M iterations in order for the error in the estimates of Theorem 2.1 to be of the same order as the error in the numerical approximation. Chapter 2. Sequential Regularization Methods for Differential Algebraic EquationsAl The most important difference between the I V P and B V P cases is that the iterative method described here does not appear to be necessarily optimal or even natural in the I V P context, certainly not from the storage requirement point of view: Note that the entire approximation of £>y s -i on [0,t/] needs to be stored. The situation here is similar to that encountered with other functional iteration methods like waveform methods. However, this difficulty can be resolved by rearranging the computation, assuming that the number of the S R M iterations, m , is chosen in advance. Thus, at each time step z, 1 < i < N, we calculate sequentially the quantities xj, Byj,x2, By2,..., a;™, By* To do this using a one-step scheme, say, we need only the corresponding quantities locally, over the mesh subinterval [£,-_!,£;), and By®. For the latter we may use, for instance, yf = y e (0), i.e. Byf = B(i,-)yo, 0 < i < N. The storage requirements are now independent of N and other typical I V P techniques such as local error control may be applied as well. 2.7 N u m e r i c a l Exper iments We now present a few very simple examples to demonstrate our claims in the previous sections. Throughout this section we use a constant step size h and set tf = 1. To make life difficult we choose h so that there is an i such that i ; = i * (if there is a singularity). In the implementation we monitor the size of the pivot in a Gauss elimination procedure for GB and slightly perturb ti away from t* when needed. A t a given time t, we use 'ex' to denote the maximum over al l components of the error in xs while 'ey' denotes the maximum over all components of the error in Bys. Similarly, 'drift' denotes the maximum residual in the algebraic equations. We first look at a boundary value problem. Chapter 2. Sequential Regularization Methods for Differential Algebraic Equations^! E x a m p l e 2.2 Consider the DAE (2.16) with - - ( - . ' : ) • • - ( . : . ) • - ( - : • ) G = ( 1 — 2* 1 - 2t) , r = - ( 1 - 2*)(e _* + smt)' subject to xi(l) + x2(0) = 1 /e . The exact solution is xe = ( e _ i s'mt), ye = 4^ singularity is located at t = 1/2, where ye becomes infinite while Bye stays bounded. We start computing with the iterate yo(t) = 0. In Table 2.1 we list errors when using the midpoint scheme = 4 - t * i L j - f l p J _ i + M By' , = By'-* +e-1B1_lEi_,(G,_LX-_l+ri_,) a % 2 J ~ 2 2 2 V * 2 1 2 2 7 where xs._L = x,+^'~1 (but no such relation is necessary for ys). We apply this scheme with h{ = h = .01 for various values of e. /£ is indicated in [11] that this scheme has 2nd order accuracy in ex and in ey, except for the case t <§C h when the error's order in By drops to 1. This is evident in the error column for t = 1.0. Note also the 0(e) improvement per SRM iteration when this term dominates the error (i.e. when es 3> h2). We note that the approximation for By at points within the initial layer is not accurate. To get a better approximation within the initial layer ( i.e. near the initial point t = 0), we solve a terminal value problem (2.44a) (changing BEG to —BEG), (2.68), as described in §2.6. Then we apply the SRM for the given problem with the improved values for Byo. In Table 2.2 we list the computed results after 3 iterations. They are obviously much better than the comparable ones in Table 2.1. • Chapter 2. Sequential Regularization Methods for Differential Algebraic Equations43 e iteration error at —) t= .Ul t = . l t = .3 t=.5 t=1.0 le-1 1 ex .38e-l .35e-l .56e-l .52e-l .39e-l ey .96 .40 .14 .34e-l .66e-l drift .87e-2 .49e-l .51e-l .0 .61e-l 2 ex .92e-2 .37e-l .89e-2 .65e-2 .72e-2 ey .91 .96e-2 .14 .34e-l .65e-2 drift .90e-2 .32e-l .61e-2 .0 .72e-2 3 ex .94e-2 .19e-l .12e-l .63e-2 .15e-2 ey .87 .20 .30e-l .43e-l .85e-3 drift .86e-2 .16e-l .43e-2 .0 .74e-3 le-2 1 ex .38e-2 .60e-2 .53e-2 .44e-2 .38e-2 ey .67 .30e-2 .40e-2 .48e-2 .62e-2 drift .65e-2 .80e-2 .38e-2 .0 .55e-2 2 ex .45e-2 .10e-3 .88e-4 .77e-4 .64e-4 ey .45 .32e-3 .70e-4 .44e-4 .23e-4 drift .44e-2 .15e-4 .14e-4 .0 .68e-4 3 ex .30e-2 .55e-5 .52e-5 .59e-5 . l l e - 4 ey .30 .18e-2 .26e-4 .18e-4 .67e-5 drift .29e-2 .17e-5 .26e-5 .0 .56e-5 le-3 1 ex .13e-2 .58e-3 .52e-3 .45e-3 .39e-3 ey .17 .47e-2 .41e-3 .49e-3 .62e-3 drift .17e-2 .79e-3 .38e-3 .0 .54e-3 2 ex .30e-3 .71e-4 .75e-5 .72e-5 .12e-4 ey .30e-l .17e-l ,34e-4 .13e-4 .54e-5 drift .30e-3 :51e-4 .20e-5 .0 .65e-5 3 ex .65e-4 .15e-3 .70e-5 .70e-5 .12e-4 ey .70e-2 .33e-l .12e-3 .14e-4 .56e-5 drift .69e-4 . l l e -3 .21e-5 .0 .59e-5 Table 2.1: S R M errors for Example 2.2 using the mi dpoint scheme Next we consider in i t ia l value problems. E x a m p l e 2.3 Consider the same DAE as for Example 2.2 with the same exact solu-tion but with initial values Xi(Q) = 1, x2(0) — 0 specified. From these initial conditions we can calculate y(0) = 1 in advance, and we choose the initial guess yo(t) = 1. Ta-bles 2.3 and 2.4 display error results for e = .1 and h = .001 using the backward Euler and the forward Euler schemes, respectively. As explained in §2.6 we calculate all iterates at each step before proceeding to the next. These tables show a significant improvement with each SRM iteration and no strong initial layer effect, as predicted by theory. Chapter 2. Sequential Regularization Methods for Differential Algebraic EquationsAA e error at —) • t=.01 t = . l t = .3 t=.5 le-1 ex .62e-2 .54e-2 .39e-2 .28e-2 ey .48e-2 .45e-2 .44e-2 .54e-2 5e-2 ex .76e-3 .58e-3 .32e-3 .17e-3 ey .22e-2 .18e-2 .10e-2 .50e-3 le-2 ex .57e-4 .49e-4 .35e-4 .26e-4 ey .10e-3 .85e-4 .52e-4 .32e-4 le-3 ex .49e-4 .42e-4 .31e-4 .23e-4 ey .56e-4 .46e-4 .30e-4 .19e-4 Table 2.2: S R M errors for Example 2.2 using the shooting-back technique iteration error at —> t=.001 t = . l t = .3 t=.5 t=1.0 I ex .2(Je-5 .72e-2 .37e-l .63e-l .11 ey .20e-2 .12 .15 .12 .59e-l drift .15e-5 .60e-2 .16e-l .76e-4 .15 2 ex .20e-5 .51e-2 .13e-l .10e-l .25e-2 ey .20e-2 .68e-l .45e-2 .20e-l .80e-2 drift .15e-5 .42e-2 .58e-2 .14e-4 .67e-2 3 ex .20e-5 .35e-2 .23e-2 .16e-2 .76e-3 ey .20e-2 .32e-l .26e-l .10e-l .37e-2 drift .15e-5 .29e-2 .12e-2 .10e-5 .12e-2 Table 2.3: S R M errors for Example 2.3 using backward Euler E x a m p l e 2.4 Here we investigate the use of the modified formula (2.42) instead of (2.35). First, we solve the previous example numerically using (2-42). In Table 2.5 we record error values at the singularity point t = .5 after 3 SRM iterations, starting with yo(t) = 1 and using as before e = .1 and h — .001 (cf. Tables 2.5). From these results it is clear that the SRM with (2.42) does not work well when a.\ 7^ 0: large errors in By are obtained near the singularity and these adversely affect the accuracy in x as well. However, the comparison changes when there is no singularity in the constraints: We now replace the algebraic constraint in Example 2.3 by xi + %2 — e - t — smt = 0 leaving everything else the same (including the singularity in B). In Table 2.6 we record maximum errors in x and By over all mesh points (denote those 'exg' and Chapter 2. Sequential Regularization Methods for Differential Algebraic Equations^ iteration error at —) t=.u(Jl t = . l t = .3 t=.5 t=1.0 1 ex .50e-6 .71e-2 .36e-l .63e-l .11 ey .20e-2- .12 .15 .12 .60e-l drift .50e-6 .60e-2 .16e-l .76e-4 .15 2 ex .50e-6 .51e-2 .12e-l .10e-l .44e-2 ey .20e-2 .68e-l .41e-2 .20e-l .70e-2 drift .50e-6 .42e-2 .58e-2 .14e-4 .67e-2 3 ex .50e-6 .35e-2 .43e-2 .18e-2 .98e-3 ey .20e-2 .32e-l .26e-l .97e-2 .46e-2 drift .50e-6 .29e-2 .12e-2 . l l e -5 .12e-2 Table 2.4: S R M errors for Example 2.3 using forward Euler ( a i , a 2 ) -> method (0,1 ex ) ey ( M ) ex ey ex (1,1) ey backward Euler forward Euler .16e-2 . .18e-2 . 10e-l 96e-2 .80e-3 .20e+l .25e-2 .18e+l .15 .64 .37e+3 .15e+4 Table 2.5: Errors near singularity using modified formula (2.42) 'eyg', respectively) for the starting iterates yo(t) = 1 and yo(t) — 0 (the latter does not agree with the exact ye(0)). ( a i , a2) -> yo method (o, 1) (h, 1) (1, 1) exg eyg exg eyg exg eyg = 1 backward Euler .46e-2 .44e-l .45e-2 .43e-l .22e-3 .28e-3 forward Euler .45e-2 .44e-l .44e-2 .43e-l .19e-3 .23e-3 = 0 backward Euler .22e-l .97 .22e-l .94 .12e-3 .75e-3 forward Euler .22e-l .97 .22e-l .94 .19e-3 .75e-3 Table 2.6: Errors for problem without singularity using modified formula (2.42) The modified method (corresponding to ( a i , a 2 ) = (1,1) in Table 2.6 is seen to work better for problems without singularities. The above calculations al l agree with our theoretical results described in §2.5 and §2.6. Chapter 2. Sequential Regularization Methods for Differential Algebraic Equations4Q 2.8 M o r e about the P r o o f of Theorem 2.1 To provide a better understanding about the sequential regularization method we now give a detailed proof of Theorem 2.1 for the ini t ia l value case with no layers, s < J + l. J is some positive integer defined in Definition 2.1. In this proof, the construction of the asymptotic expansion is directly for x and By. Moreover, the construction method we apply is somewhat different from [111, 112] and more relevant to the concept of D A E s . The same idea is applied to prove the convergence of the S R M for Navier-Stokes equations in Chapter 4. For s > J + l , additional ini t ia l layer expansions have to be developed. However, the construction of these layer expansions is precisely the same as in [111, 112] and so it is omitted here. In case that (2.17) are ini t ia l conditions (i.e. B\ = 0) our assumptions imply that (2.17) together wi th (2.18b) specify x(0), say x(0) = x (2.70) A t first, consider the case s = 1 of (2.34),(2.35): ex[ + B(GB)~l(Gxi +r) = eAxx - eBy0 + eq, with the ini t ia l conditions (2.70). This is a singular-singularly-perturbed problem (see [112, 89]). Let x\ = x10 + exn H + esxls H Comparing the coefficients of like powers of e, we thus have B(GBYlGxw = -B(GB)~lr (2.71a) B{GB)-lGxxl = -x'10 + Ax10-By0 + q, (2.71b) B{GB)-lGxu = - z i , - . ! + Axu-i, 2 < i < s + 1, (2.71c) where (2.71a) satisfies (2.70) and (2.71b) and (2.71c) satisfy homogeneous in i t ia l conditions corresponding to (2.70). Now, (2.71a) has infinitely many solutions in Chapter 2. Sequential Regularization Methods for Differential Algebraic Equations^! general. To realize the construction, we should choose x\o to satisfy (2.71a) and to ensure that the solution of (2.71b) exists. We choose xw to be the solution xe of problem (2.16)-(2.17), i.e. x[0 = Ax10 - Bye + g, (2.72a) 0 = Gxl0 + r, (2.72b) £ o *io (0) = 8. (2.72c) So xw — xe and (2.71b) has the following form B(GB)-1Gx11 = B(ye-y0). (2.73) Now we choose i n and a corresponding j/oi to satisfy x'n = Axn - By01 (2.74a) Gxn = GB{ye-y0), (2.74b) Boxn{0) = 0. (2.74c) Noting that Bye — —x'e + Axe + q is smooth, we have GB(ye — yo) £ <S\ Hence, using Lemma 2.1, there exists a smooth solution x n of (2.74), and xn satisfies (2.73). Indeed, using (2.74b) and Definition 2.1, we have G(0)a;n(0) = 0, so s n ( 0 ) = 0. A n d , from (2.74b) again, {GB)-xGxn = ye- y0, for each t £ [0,Q U {U,tf}. That is, B(GB)~1Gx11 = B{ye - y0),t £ [0,U) U (U,tf]. (2.75) Taking the l imi t of (2.75) at £*, we thus get that Xu satisfies (2.73) for each t £ [0, tj\. Moreover, from Definition 2.1, we have yoi(0) = y'ol(0) = ••• = yor 1 } (0 ) = 0, s < J + 1. Chapter 2. Sequential Regularization Methods for Differential Algebraic EquationsA8 Also we note that .Byoi is smooth. Generally, supposing we have got X i , - _ i , i ?yoi - i and yo,-i(0) = y O j _ 1 (0) = --- = y S n + 1 ) ( 0 ) = 0 for i > 2, we choose xu, yo; satisfying x'u = Axu - Byoi, Gxu = (GB)y0i-i, Boxu(0) = 0. B y the same argument as before, we obtain that xu satisfies (2.71c) for 2 < i < s + 1, and yoi(0) = y'oi(0) = ••• = y S r ° ( 0 ) = 0, s < J + 1. Also, Byoi is smooth. Next we denote the asymptotic solution xu+i = xio + ex i i H h e s x l s + e s + 1 a ; i s + i and ^ls+l = X i — Xls+l-Then tz'ls+1 + Pzis+i = eAzls+1 + e s + 2 ( - x ' l s + 1 + A x i s + i ) , *i s+i(0) = 0-Let u i s + i = Rzis+i and w\s+i = Pzis+i. Hence, we have (cf. (2.38)) . 2 l s + i = Suls+i + Wls+i and u'u+1 = (RA + R')Suls+1 + {RA + R')wls+1 + 0{es+1), Chapter 2. Sequential Regularization Methods for Differential Algebraic Equations49 ew'l3+1 + wla+1 = e{PA + P')Suls+1 + e(PA + P')wla+1 + 0{es+2), uls+1(0) = 0, wla+1(0) = 0. Using Lemma 2.2, we get wia+i = 0(es+2) and uis+i = 0(es+1), i.e. zla+1 = 0(es+1). Therefore, xj = x 1 0 + e.T 1 1 + :-- + £ s x l s + O ( e 5 + 1 ) . (2.76) Noting x 1 0 = xe, we thus obtain x l - x e = 0{t). (2.77) Then, by using (2.35),(2.76),(2.71),(2.72a) and (2.74a), it follows that Byi = By0 + -^B{GB)-\GXl + r) BVl = By0 + \{Pxw + B(GB)-1r + ePx11 + --- + esPxls + O(es+1)) = Bye + eByoi + '-' + t^Byos-i + O^) or Byx - Bye = 0(e). (2.79) Now we look at the second iteration s = 2 of (2.34), (2.35): ex'2 + B(GB)~1(Gx2 + r) = eAx2 - tByx + eq, with in i t ia l conditions (2.70). Let x2 = x20 + ex21 + e2x22 H . Noting that (2.78) gives us a series expansion for Byx we obtain, B{GB)~lGx2o = -B(GB)~lr, (2.80a) B{GBYlGx21 = -x'20 + Ax20 - Bye + q, (2.80b) B{GB)~lGx2i = -x'2l_x + Ax2i-X - By0l-U 2 < i < s + 1 (2.80c) Chapter 2. Sequential Regularization Methods for Differential Algebraic EquationsbO Again , (2.80a) satisfies ini t ia l conditions (2.70) and (2.80b) and (2.80c) satisfy the corresponding homogeneous ones. As the case of s = 1, we choose X2Q = xe. We thus have B{GB)-1Gx2l = 0. Then £21 is constructed to satisfy x'21 = Ax2i — Byu, (2.81a) Gx2l = 0, (2.81b) £0x21(0) = 0. (2.81c) Obviously x 2 i = 0 since (2.81) is uniquely solvable for x2i by Lemma 2.1. In general, similarly to the case of s = 1, we choose x2i satisfying. x'2l = Ax2l - Byu, (2.82a) Gx2l = - ( G f l ) ( y 0 , - i - yu-i), (2.82b) B0x2l(0) = 0. (2.82c) for 2 < i < s + 1. B y applying Lemma 2.2 and the same argument as in the case of s = 1 we get x2 = xe + tx2l + e2x22 + ••• + esx2s + 0(es+1) (2.83) or x 2 - x e = 0(e2). (2.84) Then, using (2.35),(2.80),(2.81a),(2.82a), (2.83) and (2.78), we conclude By2 = By^l-B{GB)-\Gx2-rr) = J B y e + e 2 5 y 1 2 + --- + e s - 1 5 y l s _ 1 + 0(e s ) (2.85) or By2 - Bye = 0{e2) (2.86) Chapter 2. Sequential Regularization Methods for Differential Algebraic Equations^! We can repeat this procedure, and, by induction, complete the proof for s < J + l. • Chapter 3 S R M for Nonl inear Problems In the previous chapter we derived the sequential regularization method and gave a detailed continuous and discrete analysis for linear index-2 D A E s . The method relates to a combination of stabilization and penalty-like methods. In this chapter we extend the method to nonlinear index-2 and index-3 problems (u = 1 and v — 2 in (2.15)), including constrained mult ibody systems. A number of variants are proposed, and particularly effective methods are singled out in certain circumstances. A l l results obtained here are certainly applicable to the linear case. The chapter is organized as follows: In §3.1 we consider problems without con-straint singularities. Two S R M variants are discussed. One variant involving ^ (corresponding to (2.32), (2.33)) leads to nonstiff problems. Taking E = I is particu-larly attractive. The other variant, corresponding to (2.32), (2.33) with a\ = 0, does not involve ^ . The choice E = possible (otherwise one can choose E = (GB)T), makes the computation particularly simple. Problems with constraint singularities are considered in §3.2. The S R M corresponding to (2.34), (2.35) is proposed for such problems. This variant works well in practice, but our proofs to date extend only to the linear case. In §3.3 we analyze and discuss various methods for index-3 problems. A number of S R M variants are possible, combining regularization with Baumgarte's stabiliza-tion or wi th invariant stabilization [8]. Of particular interest, in case of no constraint singularity, are the methods (3.47) and (3.33)—(3.35) which corresponds to invariant 52 Chapter 3. SRM for Nonlinear Problems 53 stabilization. The choice E = I leads to particularly simple iterations. A corre-sponding convergence result is given in Theorem 3.3. In case of a possible constraint singularity, the S R M (3.41) is recommended. These methods are reformulated in §3.4 for the special case of mult ibody systems with holonomic constraints. The "winning" methods are (3.44)-(3.45) with E = I for the nonsingular case and (3.46)-(3.47) for the case where the constraint Jacobian may have isolated rank deficiencies. In §3.5 we report the results of numerical exper-iments confirming our theoretical predictions and demonstrating the effectiveness of the proposed methods. 3.1 Nonl inear , Nonsingular Index-2 Problems The nonlinear index-2 D A E (y — 1 in (2.15)) reads: x' = f(x,t)-B(x,t)y, (3.1a) 0 = g(x,t), (3.1b) where / , B and g are sufficiently smooth functions of (x,t) G R™* x [0, £/], and y G Rnv. We consider this D A E subject to nx — ny boundary conditions b(x(0),x(tf)) = 8. (3.2) These boundary conditions are assumed to yield a unique 1 and bounded solution for the O D E (3.1a) on the manifold given by (3.1b). Concretely, i f we were to replace (3.1b) by its differentiated form (denoting G = |^) 0 = Gx' + gt {=dff) (3.3a) flf(x(0),0) = 0 (3.3b) 1 locally unique, or isolated solution in a sufficiently large neighborhood would suffice. Chapter 3. SRM for Nonlinear Problems 54 and use (3.3a) in (3.1a) to eliminate y and obtain nx O D E s for x, then the boundary value problem for x wi th (3.2) and (3.3b) specified has a unique solution. In the ini t ia l value case (i.e., when b is independent of x(tj)), this means that (3.2) and (3.3b) can be solved uniquely for x(0). In this section, we consider the case where GB is nonsingular. Generalizing the idea in §2.4, we have the following S R M formulation for the nonlinear index-2 D A E s (3.1), (3.2): for s = 1 ,2 , . . . , x's = f(x„t) -B(xa,t)ya, (3.4) where ys = ya-i + ^E(xs,t)(a1-^g(xa,t) + ot2g(xa,t)), (3.5) subject to the boundary conditions (3.2) and (3.3b). Note that yo(t) is a given ini t ia l iterate which we assume is sufficiently smooth and bounded and that e > 0 is the regularization parameter. The regularization matrix E is nonsingular and has a uniformly bounded condition number; possible choices are E = I, E = (GB)-1 and others (e.g. E = (GB)T, cf. [11, 94]). We note that if we take y 0 = y then x\ = x, where x and y are the solution of (3.1). If we take y0 = 0, then one S R M iteration is the usual penalty method (cf. [84, 91, 69]). As customary for the penalty method, we assume: Assumpt ion 3.1 The problem (3.4), (3.5),(3.2),(3.3b) has a unique solution and the solution is bounded if ys-\ is bounded. Assumption 3.1 is generally true for ini t ia l value problems. For general boundary value problems, we expect that it would hold for most practical cases since (3.4) (with (3.5) plugged in) may be seen as a perturbed problem of (3.1) according to the proof of Theorem 3.1 (see below), where the perturbation and its first derivative are both small i f e is small. Chapter 3. SRM for Nonlinear Problems 55 To analyze the S R M , we assume the following perturbation inequality: For 0 < t<tf, \\x{t) - x{t)\\ < M max (\8(T)\ + \8'(T)\), (3.6a) 0 < T < i f \W)-y{t)\\ < M max (\8(T)\ + |<T(T)|), (3.6b) 0<r<i / where || • || is some lp norm (say, the maximum norm), and x and y satisfy the following perturbed version of (3.1): x' = f(x,t)-B(x,t)y, (3.7a) 0 = g(x,t) + 8(t) ' (3.7b) wi th the same boundary conditions as (3.2). For in i t ia l value problems, (3.6) has been proved in [58], pp. 478-481. It is actually the definition of the perturbation index introduced in [58]. Furthermore, (3.6) also holds for boundary value problems if we impose some boundedness conditions on the corresponding Green's function (cf. [14]). The case ct\ ^ 0 in (3.5) is sufficiently different from the case oc\ = 0 to warrant a separate treatment. 3.1.1 The case a i = 1 Now we estimate the error of the sequential regularization method (3.4)-(3.5). We prove a theorem which says that the error after s S R M iterations is 0(cs) (i.e., each iteration improves the error by 0(e)) everywhere in t. This result coincides wi th that of the linear case. Theorem 3.1 Let all functions in the DAE (3.1) be sufficiently smooth and the above assumptions hold. Then, for the solution of iteration (3.4), (3.5) with Qi ^ 0, we have Chapter 3. SRM for Nonlinear Problems 56 the following error estimates: xs(t)-xe(t) = 0(es), y.(t)-ye(t) = 0(es), for 0 < t < tj and s > 1. Proof: Let vs = g(xs,t). Then, from (3.4), v's = G(xs, t)x's + gt(xs, t) = G(xs,t)f(xs, t) - G(xs,t)B(xs, t)ys + gt(xs,t). Using (3.5), we thus have (e(GBE)~1 + I)v's + a2vs = e(GBE)~\Gf + gt) - tE^y,^, (3.8a) v„(0) = 0. (3.8b) Therefore it is not difficult to get vs = g(xs,t) = 0(e), v's = g(xs,t)' = 0(e), (3.9) i f y s _ i is bounded (which implies that xs is bounded from Assumption 3.1). For s = 1, we have x'i = f(xi,t)-B(x1,t)y1 g(xut) = 0(e), g(xut)'= 0(e) since yo is chosen to be bounded. From assumption (3.6), we immediately get X! - xe = 0(e), j / i - ye = 0(e). (3.10) Then it is easy to see that j / i is bounded. So for s = 2, we obtain x'2 = f(x2,t)-B(x2,t)y2 g(x2,t) = 0(e), g'(x2,t) = 0(e) Chapter 3. SRM for Nonlinear Problems 57 B y using assumption (3.6) again, this yields x2 — xe = 0(c). Hence it can be verified, by substituting (3.3a),(3.la) for the exact solution, that the right hand side of (3.8a) becomes 0(e2). So, from (3.8), we can get g(x2,t) = 0(e2), g'(x2,t) = 0(e2). • Apply ing assumption (3.6), it follows that x2 - xe = 0(e2), y2-ye = 0(e2). (3.11) This also gives the boundedness of y2. We can repeat this procedure, and, by induction, conclude the results of the theorem. • From (3.8) it is clear that there is no stiffness here, so we can choose e > 0 very small, so small in fact that one S R M iteration would suffice for any desired accuracy, and discretize the regularized O D E , possibly using a nonstiff method like explicit Runge-Kutta . This gives a modified penalty method [/ + e-xBEG}x' = f - B y 0 - e~1BE(gt + a2g) (3.12) where B,E,g etc, all depend on x, with the subscript 5 = 1 suppressed. For the choice E = (GB)~\ let P = BEG = B(GB)~lG be the associated projection matr ix. Mul t ip ly ing (3.12) by *_ t P and by I — P, respectively, and then adding together, we have x' = f - rrr^Byo - -J-L-^BiGBr^Gf + 9 t + a2g) Thus the iteration obtained is similar to Baumgarte's stabilization x' = f - B(GB)~1[Gf + 9 i + a2g] (3.13) Chapter 3. SRM for Nonlinear Problems 58 In fact, the single S R M iteration tends to (3.13) in this case when e —> 0. Indeed, the parameter cx2 is the usual Baumgarte parameter, and choosing cx2 > 0 obviously makes equation (3.8a) asymptotically stable for the drift vs. As indicated in [12], for both of these methods we can apply post-stabilization instead, i.e. take a 2 = 0 but stabilize after each discretization step [8, 9]. For reasons of computational expense, it may be better to choose E = / in (3.12). The iteration obtained is simple, although a possibly large matrix (with a special structure) must be "inverted". E x a m p l e 3.1 The choice of E = I is utilized in Chapter 4 (see also [80]) for the time-dependent, incompressible Navier-Stokes equations governing fluid flow. The advantage gained is that no treatment of pressure boundary conditions is needed, unlike methods based on Baumgarte-type stabilizations which lead to the pressure-Poisson equation. • 3.1.2 The case ax = 0 For this case the drift equation (3.8) is clearly stiff for 0 < e <C 1. As in §2.5, we denote J such that yo(O) = y e(0), y'o(0) = y'e(0),y{0J)(0) = y( J)(0), (3.14) where J = —1 if yo(0) ^ y e (0), then we can prove the same result as Theorem 3.1 for s < J + 1. Note that we may choose y0 satisfying (3.14) for some m > 0 by expressing y in terms of x at t = 0 for ini t ia l • value problems. But this starting procedure generally does not work for boundary value problems. Hence we state and prove the theorem for ini t ial value problems and comment on the boundary value case following the proof. Chapter 3. SRM for Nonlinear Problems 59 Theorem 3.2 Let the assumptions of Theorem 3.1 plus (3.14) hold. In addition, suppose that the matrix function E(x,t) has been chosen so that GBE is positive definite. Then, for the solution of iteration (3.4),(3.5) with ot\ = 0, we have the following error estimates: xs(t)-xe(t) = 0(es), ys(t)-ye(t) = 0(es), for 1 < s < J + 1 andO <t <tf. Proof: We derive the result for the case s < J + 1 = 2. Following the proof, we wi l l comment on additional generalizations. The key is again the basic drift equation (3.8), which we rewrite here as ev's + (GBE)vs = e{Gf + gt-GBys_1), (3.15a) v3(0) = 0. (3.15b) where quantities are evaluated as before, at ( x s , t ) , unless otherwise noted. For 5 = 1, given the boundedness of yo we obtain as before >h - 0(e) To obtain a similar result for v[, however, a different procedure from that of Theorem 3.1 is needed. Note that at t = 0, the condition yo(0) — y(0) implies (Gf + gt-GByo)\t=o = 0 Hence from (3.15a), ui(0) = 0. Differentiating (3.15a) with respect to t and using v\ = 0(e), we get ev" + (GBE)v[ = 0(e), v[(0) = 0 and this yields vi = 0(e) Chapter 3. SRM for Nonlinear Problems 60 From assumption (3.6) we then get (3.10). Subtracting (3.1a) from (3.4) and using (3.10) gives also x[ = x' + 0(e) and boundedness of y[ is obtained from a differentiation of (3.5). For s = 2, given the boundedness of y i (by (3.10)) and j/i(0) = 2/(0), we get as for s = 1 v2 = 0(e), v'2 = 0(e) and hence also x2 — x + 0(e) This yields that the right hand side of (3.15a) is 0(e2), so v2 = 0(e2) Now comes the delicate part. To obtain an 0(e2) estimate also for v'2, so that the estimate (3.6) can be used to complete the proof, we differentiate the drift equation again, obtaining tv'2' + (GBE)v'2 = 0(e2) + e(Gf + gt - GBy,)' and v'2(0) = 0 obtained as for the 5 = 1 case. We are then left to show that F(x2, t) := (Gf + g t - GByi)' = 0(e) (3.16) For this purpose we must estimate v" first. Using the condition y'o(0) = y'(0), and also xi(0) = x(0), x[(0) = x'(0) (obtained from (3.4)), we can obtain (Gf(xut) +gt(x1,t) - GB(x1,t)y0)'\t=o = 0 Hence v"(0) = 0 from (3.15a) once differentiated. Differentiating (3.15a) twice we now obtain precisely as when estimating v[ above, v» = 0(e) Chapter 3. SRM for Nonlinear Problems 61 The boundedness of all needed quantities can also be obtained in the same way as before. Final ly , we note v[ = Gx[+gt(xut) = Gf(xut) + gt(xi,t) - GB(xut)yi Ready to show (3.16), we now write F(x2, t) = [(Gf(x2,t)-Gf(xut))+(gt(x2,t)-gt(xut)) + (GB(xut)-GB(x2, t))yi-v[}' Our previous estimates allow the conclusion that x2 = x\ + 0(e), x2 = x[ + 0(e), hence we can finally conclude the estimate (3.16) and obtain the result of the theorem for 5 = 2. The proof proceeds in a similar manner for larger J. Generally, one needs the estimate = 0(e), 1 < j < J + 1, and this necessitates (3.14). • R e m a r k 3.1 The convergence result holds for all s (i.e. also for s > J + l , assuming sufficient smoothness) away from an initial layer of size 0(e) in t. This is so because E is chosen so that we can express the solution for small e as a smooth outer solution which is bounded in terms of the right hand side as before, plus an initial layer of width 0(e) . Conditions (3.14) then ensure that the layer error is bounded by 0 ( e J + 1 ) for the first J + l iterations. • R e m a r k 3.2 For boundary value problems, there is no obvious technique to ensure J > — 1. For a given J, the results of Theorem 3.2 and Remark 3.1 can be extended as in Chapter 2. This requires a different proof technique, though. Basically, an asymptotic expansion for xs and ys is constructed, where the first term is the exact solution x, y. This latter proof technique follows more along traditional singular perturbation lines (see [112, 89, 69]), and is not as close to Theorem 3.1 and to DAE concepts. • Chapter 3. SRM for Nonlinear Problems 62 Taking a2 = I without loss of generality, we obtain the iteration x'„ = f - Bys-, - e-'BEgixs, t) (3.17) This is a singular, singularly perturbed problem (so e should not be taken extremely small compared to machine precision even if a stiff solver is being used). If GB is positive definite then we may choose E — I, and this yields a very simple iteration in (3.17) which avoids the inversion necessary in stabilization methods like Baumgarte's. However, i f an explicit discretization method of order p is contemplated then approx-imately p S R M iterations like (3.17) are needed, because one must choose e = 0(h), where h is the step size. 3.2 Nonl inear , Singular Index-2 Problems In this section we consider the nonlinear index-2 problem (3.1) with an isolated sin-gular point t*, i.e. GB is singular at t*. For simplicity, we assume that B and g are independent of t. Denote P(x) — B(GB)~1G. Motivated by constrained mult ibody systems (see Example 2.1), we assume P(x) to be differentiable in t, but f^ -(x) may be unbounded. For this reason, we consider only the case a\ = 0 in this section (cf. Chapter 2 for the linear case). In the drift equation (3.8) we then have essentially the singularly perturbed operator ev' + GBEv to consider. The choices of E = / or E = (GB)T yield a turning point problem (i.e., at least one of eigenvalues of the ma-tr ix GBE vanishes at the point £*), which complicates the analysis, even in the linear case , and degrades the numerical performance as well in our experience. Therefore, we choose E = (GB)~l. In the sequel we wi l l be careful to evaluate the effect of E only when its singularity l imi t is well-defined, as e.g. in P(x). A direct generalization of the linear case in Chapter 2 would give the S R M for-mulation (3.4), (3.5) where instead of updating y (because y may be unbounded at Chapter 3. SRM for Nonlinear Problems 63 t*) we update By by B{xs)ys = Bix^y,-! + l-B{xs){G{xs)B{xs))-1 g{xs). (3.18) However, (3.18) needs to be modified, since we may have RangeB(xs) ^ RangeB(xs-\). So we use the projection P(xs) to move from RangeB(xs-\) to RangeB(xs). Then we consider the following S R M formulation for singular problems: x's = f(x„t)-B(x,)ya, (3.19a) B(xs)ys = P{xB)B{xt-i)y,-X + -eB(xs)(G(xs)B{xs))-lg(xs), (3.19b) where xs satisfies the boundary condition (3.2). If the assumptions given at the beginning of §3.1 and in Theorem 3.2 remain valid, then the result of Theorem 3.2 st i l l holds. Unfortunately, for the singular problem, assumption (3.6) may not be true in general. To see this, consider one iteration, i.e. 5 = 1. The accuracy for the approximation of x depends on the extent that the bound (3.6a) holds. Numerical experiments show that we can get a good approximation of x near the singularity. But the situation for By is worse, and the bound (3.6b) often does not hold. Indeed, assume for the moment that we have a good, smooth approximation of Sciy xs —• x^ i.e. (3.7) holds with 8,8' = 0(e), and B(x)y is defined by (3.19b) for some B ( x s _ i ) y s _ i . From (3.7) we have B(x)y = P(x)f{x,t) + ri,' (3.20) where n = B(x)(G(x)B(x))~18'. It is not difficult to find that the exact B(x)y from (3.1) satisfies B(x)y = P(x)f(x,t). (3.21) Yet , even if n is small, B(x)y may not be a good approximation of By because | ^ may be unbounded at the singular point so that P(x) is not a good approximation of P(x). Chapter 3. SRM for Nonlinear Problems 64 E x a m p l e 3.2 In (3.1) let x = (xi,x2), g(x) = — cosa?i — cosx 2., and G = BT = ( s in 2 X\ s\nxism.x2\ j . sin X i sin £ 2 s in 2 x2 / Clearly, at a singular point x = (0,7r), the value of P depends on the direction from which it is approached. Thus, |^ is unbounded, even though P is a differentiable function oft. Further letting f = ( s inx 2 — s i n x i ) - 1 (s in £ 2 2 sin 2 2 — s i n x i ) r ; and given the initial conditions £ i ( 0 ) = — 7 r / 2 , £2 (0 ) = n/2, the exact solution is T 1 x(t) = (t — ir/2 t + ir/2) , y = (sinx 2 — sinXi)~ Thus, as t crosses t* = n/2, y(t) becomes unbounded, but By = (sin x2 — sin X i ) - 1 ( sin x, s i n x 2 ) T remains bounded. However, it is easy to perturb x(t) slightly and smoothly in such a way that the perturbed By becomes unbounded near t = t*, still satisfying (3.7) with a small 5. • Note that for the linear model problem (see Chapter 2), P = P(t) is independent of x. Hence we do not have the above difficulty in the linear case. For the nonlinear problem, the accuracy near the singular point is reduced and it no longer behaves like 0 (e s ) for more than one iteration. However, we do expect 0(es) accuracy away from the singular point, assuming that no bifurcation or impasse point is encountered by the approximate solution . 3.3 S R M for Nonl inear Higher- index Problems We now generalize the S R M to the more general problem (2.15). In particular, we consider the index-3 problem (y — 2). The Euler-Lagrange equations for mult ibody systems wi th holonomic constraints yield a practical instance of the problem. The Chapter 3. SRM for Nonlinear Problems 65 S R M formulations presented in this section are easy to generalize for more general problems (2.15) (index u + 1). The index-3 problem reads: x". = f(x,x',t)~ B(x,t)y, (3.22a) 0 ' = g(x,t), (3.22b) wi th given 2(nx — ny) boundary conditions, b(x(0),x(tf),x'(0),x'(tf)) = 0. (3.23) The meaning of G, B and the stabilization matrix E below remain the same as in the index-2 problems considered in previous sections. 3.3.1 The case of nonsingular GB We first use the idea from previous sections, viz. a combination of Baumgarte's stabi-l ization wi th a modified penalty method, to derive the S R M for the nonlinear index-3 problem (3.22). Then we apply a better stabilization [8] to generate a new S R M which is expected to have better constraint stability. Final ly, we seek variants which avoid evaluation of complicated terms in the second derivative of the constraints. First consider, instead of (3.22b), the Baumgarte's stabilization d2 d aidJ29(x,t^ + a2dt9^X,t') + ot39(x>t) = °>. (3.24a) 0(x(O),O) = O, !<7(*(0),0) = 0, (3.24b) where ctj,j = 1,2,3 are chosen so that the roots of the polynomial 3 a(r) = Y,a3r3~3 i=i are al l negative. Following the procedure of previous sections, we can write down an S R M for (3.22): for s = 1, 2 , . . . and y0 given, x': = f(xs,x's,t)-B(xs,t)ys, (3.25) Chapter 3. SRM for Nonlinear Problems 66 where xs satisfies boundary conditions (3.23) and (3.24b) and ys is given by It is not difficult to repeat the approach of §3.1 for the present case. Under assumptions similar to the index-2 case, i.e. (3.6) with a change to include 6~"(r) at the right hand side (cf. [58]) and Assumption 3.1 with the addition that the derivative of the solution is also bounded, we readily obtain extensions of Theorems 3.1 and 3.2 for the cases « i ^ 0 and ax = 0 (with a2 ^ 0), respectively. We do not allow « ! = a2 = 0 since in this case equations (3.25),(3.26) have different asymptotic properties. Note that the S R M (3.25),(3.26) with Qi = 0 avoids computing gxx; however, the iteration obtained now calls for solving problems which become stiff when e gets small , and to avoid gxx one should use a non-stiff discretization method. This formulation with E = I and a.\ ^ 0 is the same as that proposed in [20, 19] using the augmented Lagrangian method. Another way to generalize the S R M to higher index problems is based on invariant stabilization. Its advantages over Baumgarte's stabilization have been discussed in [8, 9]. We thus prefer the way based on this stabilization. Theoretical evidence is also mentioned in Remark 3.4. We first describe this new stabilization. B y two direct differentiations of the constraints (3.22b), we can eliminate y and get an O D E for which the original constraint (3.22b) together with its first derivative give an invariant. The idea of the method is to reformulate the higher index D A E (3.22) as ys = Vs-i + -E(xs,t)(a1—g{xa,t) + a2—-g(xs,t) + a3g(xs,t)). (3.26) e f(x,x',t), (3.27) a first order O D E (cf. (3.27)): (3.28) wi th an invariant 0 = h(z,t), (3.29) Chapter 3. SRM for Nonlinear Problems 67 where and to consider the stabilization families z' = f(z,t)-7F(z,t)h(z,t), (3.31) where F = DE for some appropriate matrix functions D and E such that E and HD are nonsingular and H = hz. The O D E (3.31) coincides wi th Baumgarte's stabilization for the index-2 problem (3.1) with D = B and E = E = (HD)-1. One choice for D is D = HT, but others wi l l be mentioned below. Note that (3.31) has the same solution as the original problem (3.22) for any parameter value 7. Al though the method has better constraint stabilization, both the evaluation of / and that of H involve gxx which may be complicated to calculate in practice. Next , we present an S R M method based on invariant stabilization which avoids the computation of / . In fact, we can avoid gxx altogether using the new stabilization. If we do not eliminate y by differentiations, f(z,t) in the stabilization (3.31) becomes ^ = ( , ( M ) - V 0 » ) - < 3 ' 3 2 ) Since y is not known in advance, we use an iterative S R M procedure to calculate y as in [20, 11]. The solutions of the iterative procedure no longer satisfy (3.22) precisely. Hence the iterative procedure has to be a regularization procedure and the parameter in (3.31) is changed to 7 = ^ to emphasize that it must be chosen sufficiently large. These lead to the following S R M formulation (for simplicity of notation, we only consider the special case where B and g are independent of t): < = ( ZlS) = ( f ( Z%( , ) - -F(zs)h(zs), (3.33) Z2s \ f{za,t) - B(zu)y,-i I e Chapter 3. SRM for Nonlinear Problems 68 where zs satisfies boundary conditions (3.23), (3.24b) and h = (g(zi), G(zi)z2)T. Thus the Jacobian of h is H = ( ^ ^ I , where L = z?qxx(zx). \ L(z) G(Zl) ) ' 2 i / V ^ We choose D and E so that F = B E ( ' °) = (BE ° ) (3.34) \0 I ) \ 0 BE ) v ; where, as in §3.1, E is chosen such that GBE has non-negative eigenvalues. Updat ing y by ys = ys-\ + -E(zls)G(zls)z2s (3.35) implies that the second part of the original index-3 system holds exactly, i.e. z2s = f(z*;t) - B(zls)ys. Next we analyze the convergence of (3.33)—(3.35). Again we assume that the solutions of (3.33), (3.23), (3.24b) exist uniquely and are bounded if y s _ i is bounded (see Assumption 3.1 ). Assumption (3.6) changes slightly: We first rewrite the system (3.22) as z[ = z2, . (3.36a) z'2 = f(z,t)-B(Zl)y, (3.36b) 0 = g(Zl). (3.36c) Then we assume the following perturbation bound, \\z(t) - z(t)\\ < M max (|<y(r)| + \S'(r)\ + + I W I + I W I ) , (3.37a) U S T S l f \\y(t) - y(t)\\ < M max (\S(T)\ + \S'(r)\ + \S"(r)\ + \0(r)\ + \9'(r)\), (3.37b) Chapter 3. SRM for Nonlinear Problems 69 where z and y satisfy a perturbed problem of (3.36), z[ = z2 + 0(t), (3.38a) z'2 = f(z,t)-B(Zl)y, (3.38b) 0 = flfO)+ (3.38c) wi th the same boundary conditions (3.23). Again , for ini t ia l value problems, (3.37) can be easily proved by following the technique presented in [58], and this can be extended for boundary value problems as well. Similar ly to the proof of Theorem 3.1, let h(zs) = I, where vs = g(zis), \wsJ ws — G(z\s)z2s. From (3.33), we get the drift equations (cf. (3.15)) ev's = -GBE(zs)vs + ews (3.39a) tv"s = -GBE(zs)v's - LBE{zs)vs + t[Lz2s + Gf(zs) - GB(zs)ys.1] (3.39b) wi th the in i t ia l conditions vs(0) = 0, ws(0) = 0. Apply ing (3.39a) and then (3.39b) for s = 1, we obtain = 0(e), wi = 0(e). Then (3.39a) further yields v, = 0(e2), v[ = 0(e) Comparing (3.38) with (3.33)-(3.35), we have to bound S = vu 0 = -e~lBEvi and their derivatives appearing in (3.37). We already-have that 8 = 0(e2), 8' = 0(e), '9 = 0(e) The procedure that follows continues to be similar to the one employed in the proof of Theorem 3.2, so we only sketch it here for s — 1. From (3.39a) we obtain v[(0) = 0. Using the condition yo(0) = y(0) gives [Lz21 + Gf(zl)-GB(z1)yo]\t=o = 0 Chapter 3. SRM for Nonlinear Problems 70 so, from (3.39b), also w[(0) = 0. Differentiating (3.39b) we get tw'l + GBE(zl)w'l = 0(e), ioi(0) = 0 hence w[ = 0(e). Differentiating (3.39a) we next have eo'l + GBE(z1)v[ = 0(e 2 ) , v[(0) = 0 so v[ = 0 (e 2 ) . This implies 8' = O(e), 8" = O(e) We can now use (3.37) and obtain the desired conclusion for 5 = 1, zl = z + 0(e), y i = y + 0(e), where {z,y} is the exact solution of the index-3 problem. Then, continuing to follow the proof procedure of Theorem 3.2, we obtain: Theorem 3.3 Let all functions in the DAE (3.22) be sufficiently smooth and assume the above assumptions (particularly (3.37)) hold. Assume in addition that yo satisfies (3.14)- Then, for the solution of iteration (3.33)-(3.35), the following error estimates hold: zs(t)-ze(t) = 0(es), (3.40a) ys(t)-ye(t) = 0(es) (3.40b) forl<s<J + l. R e m a r k 3.3 Extensions of this theorem to the boundary value case and to s > J +1 away from an initial layer are possible, similarly to the extensions for Theorem 3.2 contained in Remarks 3.1 and 3.2. • Chapter 3. SRM for Nonlinear Problems 71 R e m a r k 3.4 We note that, unlike Proposition 2.2 of [8], we do not assume \\H(z)f(z)\\2<^\\h(z)\\2 to discuss the stability and accuracy of the constraints. Also, from (3.39), we see the difference of the constraint stability or accuracy between SRM formulations based on Baumgarte's stabilization and the new stabilization. For the former, we only have v[ = G(zn)z'n = G(zu)z21 = Wi So if we obtain w\ = 0(e) then vx = 0(te). This can be much worse than what we get from (3.39a). • 3.3.2 The case for constraint singularities For the singular case GB may be singular at some isolated point t* as described in the previous sections. The situation here is similar to that for index-2 problems. A n examination of the drift equations (3.39) suggests that here, too, the choice E — (GB)-1 is preferable to E = I or E — (GB)T. The iteration for ys is modified as well. S t i l l assuming for simplicity that g and B do not depend explici t ly on t, this gives, in place of (3.33)—(3.35) the iteration z'u = z2s-l-B(GB)-1g(zs) (3.41a) 4s = f(*„t)-y. (3.41b) Vs = P{zs)ys-i + ^P(zs)z2s (3.41c) Also, as indicated in §3.2 for index-2 problems, we cannot expect 0(es) approxi-mation near the singular point any more. But we do expect that (3.40) holds away from the singular point, because the singularity is in the constraint and the drift manifold is asymptotically stable (following our stabilization). A numerical example Chapter 3. SRM for Nonlinear Problems 72 in §3.5 wi l l show that we do get improved results by using S R M iterations for the singular problem. 3.4 S R M for Constrained M u l t i b o d y Systems Constrained mult ibody systems provide an important family of applications of the form (3.22) and (3.1). We consider the system q' = v (3.42a) M(q)v' = f(q,v)-G(q)TX (3.42b) 0 = g(q) (3.42c) where q and v are the vectors of generalized coordinates and velocities, respectively; M is the mass matr ix which is symmetric positive definite; f(q, v) is the vector of external forces (other than constraint forces); g(q) is the vector of (holonomic) constraints; A is the vector of Lagrange multipliers; and G(q) = j^g. For notational simplicity, we have suppressed any explicit dependence of M , / or g on the time t. We first consider the problem without singularities. Corresponding to (3.22) in §3.3, we have B = M~lGr, so GB = GM~lGT. Other quantities like h and H retain their meaning from the previous section. In some applications it is particularly important-to avoid terms involving gxx, since its computation is somewhat complicated and may also easily result in mistakes and rugged terms. So [9] suggests post-stabilization using the stabilization matr ix F = M~1GT(GM~1GT)~1 ^ Q J j ( 3 - 4 3 ) twice, instead of involving H, at the end of each time step or as needed. They find that this F performs very well in many applications. However, while this stabilization avoids the gxx term m F, gxx is s t i l l involved in obtaining / , although only through Chapter 3. SRM for Nonlinear Problems 73 matrix-vector multiplications (see (3.27)). The S R M formulation (3.33)—(3.35) en-ables us to avoid the computation of / in the absence of constraint singularities. For the mult ibody system (3.42) we write the iteration as follows: For s = 1, 2 , . . . , find {qs, vs} by q's = vs-l-BE(qs)g(qs) (3.44a) v's = M-1f(qs,vs)-B{qs)Xs-1-^BEG(qs)vs (3.44b) Then update A by A s = A s _ j + l-EG{qs)vs. (3.45) It is easy to see that in this S R M formulation the gxx term is avoided com-pletely. Moreover, since GM~lGT is positive definite, we can choose E = I in (3.44),(3.45), obtaining a method for which Theorem 3.3 applies, which avoids com-puting (GM~1GT)~1. Although it requires an iterative procedure, a small number of iterations (p if an explicit discretization method of order p is used) typically provide sufficient accuracy. Numerical experiments wi l l show the 0(es) error estimate. Next we consider the singular problem, i.e. wi th the matrix GM_1GT being sin-gular at some isolated point t*, 0 < t* < tj. A typical example of singular mult ibody systems is the two-link slider-crank problem (see Example 2.1 and Figure 2.1) con-sisting of two linked bars of equal length, with one end of one bar fixed at the origin, allowing only rotational motion in the plane, with the other end of the other bar sliding along the x-axis. Various formulations of the equations of motion for this problem appear, e.g., in [60, 19, 11, 12, 94]. In our calculations we have used the formulation in Example 2.1 (see [11]), to make sure that the problem is not accidentally too easy. It consists of 6 O D E s and 5 constraints, wi th the last row of the Jacobian matr ix G vanishing when the mechanism moves left through the point where both bars are upright (0 i = f, 02 = ^f, where £ ; , y ; , 0 ; are the coordinates of the center of mass of the ith bar). Chapter 3. SRM for Nonlinear Problems 74 The last row of G vanishes at this point and a singularity is obtained. We note that the solution is smooth in the passage through the singularity wi th a nonzero velocity. When we attempt to integrate this system using a stabilization method like [8] which ignores the singularity, the results are unpredictable, depending on how close to the singular time point the integration process gets when attempting to cross it. In fact, radically different results may be obtained upon changing the value of the error tolerance. (Similar observations are made in [94].) In some instances a general purpose O D E code would simply be unable to "penetrate the singularity" and yield a solution which, after hovering around the upright (singular) position for a while, turns back towards the ini t ia l position (solid line in Figure 2.1). Such a pattern of motion may well look deceptively plausible. Methods which do not impose the constraints on the position level (e.g. methods consisting of differentiating the constraints once and solving the resulting index-2 problem numerically, or of projecting only on the velocity-level constraint manifold) perform particularly poorly (cf. numerical results in [94]). This is easy to explain: The position-level constraint corresponds to ensuring that the two bars have equal length. If this is not strictly imposed in the process of numerical solution, inevitable numerical errors due to discretization may yield a model where the lengths are not close enough to being equal, and this leads to the lock-up phenomena described e.g. in [60], which have a vastly different solution profile. We now wish to generalize the S R M to problem (3.42) with singularities since we have seen its success for the linear index-2 case. From the two-link slider crank prob-lem, we find that, although GM~lGT is singular at t*, P(q) = M~1GT(GM~1GT)~1G and M~1GT(GM~1GT)~1g are smooth functions of t for the exact solution or func-tions q satisfying the constraints, while M~lG7\GM~lGT)-\ M~1GT(GM~1GT)~1Gg and the derivative are not. Also, as indicated in [11], A is no longer smooth, while BX is since we assume the solution q to be sufficiently smooth. We only include Chapter 3. SRM for Nonlinear Problems 75 (3.46b) (3.46a) Xs = P(qs)Xs-1 + -P(qs)vs (3.47) As we indicated in §3.2, we do not expect 0(es) accuracy near the singular point. However, we do expect that the S R M iteration would improve the accuracy and that we st i l l expect to get 0(es) accuracy away from the singular point. Numerical experiments in §3.5 wi l l show such improvements. 3.5 N u m e r i c a l Exper iments We now present a few examples to demonstrate our claims in the previous sections. Throughout this section we use a constant step size h. To make life difficult we choose h when we can so that there is an i such that t{ = t*, namely, there is a mesh point hit t ing the singularity point t*, for singular test problems. A t a given time t, we use 'ex' to denote the max imum over all components of the error in xs. Similarly, 'drift' denotes the max imum residual in the algebraic equations. E x a m p l e 3.3 Consider the DAE (2.16), (2.17) with f 9 2^Xl x<1 e sin2 *). subject to a?i(0) = 1, £2(0) 0. The exact solution is xe = (e -t s'mt), ye = e*. This is a problem without singularities. Chapter 3. SRM for Nonlinear Problems 76 Using an explicit second order Runge-Kutta method with h = 0.001 we test various choices of E and ax (always taking a2 = 1) of the SRM formulation in §3.1. We list the computational results in Table 3.1. Observe that, for*ct\ / 0, the SRM works well methods e iteration error at —) • t = . l t=.5 t=1.0 « i = 1 le-8 1 ex . l l e -7 .94e-7 .19e-6 E = I drift .79e-8 .56e-7 .14e-6 ai = 1 le-8 1 ex . l l e -7 .92e-7 .18e-6 ' E = (GB)T drift .78e-8 .53e-7 .14e-6 OL\ = 1 le-8 1 ex . l l e -7 .95e-7 .19e-6 E = (GB)-1 drift .80e-8 .58e-7 .15e-6' Baumgarte ex .45e-6 .16e-6 .35e-6 drift .40e-6 .70e-7 .29e-6 a i = 0 5e-3 1 ex .60e-2 . l l e -1 . l l e -1 E — I drift .54e-2 .80e-2 .13e-l 2 ex . l l e -3 .26e-3 .22e-3 drift .96e-4 .20e-3 .27e-3 3 ex .32e-5 .65e-5 .46e-5 drift .29e-5 .47e-5 .54e-5 4 ex .26e-6 .23e-6 .28e-6 drift .13e-6 .51e-7 .12e-6 Cii = 0 5e-3 1 ex .70e-2 .12e-l .13e-l E = (GBf drift .64e-2 .13e-l .15e-l 2 ex .22e-3 .65e-3 .31e-3 drift .20e-3 .49e-3 .29e-3 3 ex . l l e -4 .16e-4 .69e-5 drift .10e-4 .10e-4 .52e-5 4 ex .85e-6 .91e-7 .29e-6 drift .75e-6 .77e-6 .14e-6 « i = 0 5e-3 1 ex .51e-2 .66e-2 .10e-l E = (GB)'1 drift .46e-2 .49e-2 .12e-l 2 ex .35e-4 . l l e -3 .21e-3 drift .30e-4 .79e-4 .24e-3 3 ex .86e-6 .23e-5 .47e-5 drift .77e-6 .17e-5 .53e-5 4 ex .26e-6 .18e-6 .26e-6 drift .26e-7 .31e-7 .13e-6 Table 3.1: Errors for Example 3.3 using the explicit second order Runge-Kut ta scheme for various choices of E. Its error is as good as Baumgarte's method whose parameter corresponds to the a2 of the SRM. For ai = 0, we see that the error improves at a rate of about 0(e) for various choices of E, including E = I. (Observe the errors att = 1; the error situation near t = .1 is different because of an initial layer.) Such an error Chapter 3. SRM for Nonlinear Problems 77 improvement continues until the accuracy of the second order explicit Runge-Kutta method, i.e. 0(h2), is reached. • The next two examples are problems with singularities. In the index-2 case of the Baumgarte stabilization the worst term is B(GB)~1gt for the type of the singularities in this paper. To show what happens when the Baumgarte method does not work well, we choose nonautonomous problems (i.e. gt / 0) as index-2 singular examples. Example 3.4 Consider the nonlinear DAE (2.16) with \2t + (t> - ±)e* J \x2) 9 = \{xl + x l - ( t - l - f - ( e - \ f ) subject to the initial condition Xi(0) = — |, £2(0) = — | . The exact solution is xe = (t — | , t2 — |), ye = e*. A singularity is located at t* = | . Using this example we test the SRM formulations of §3.2. We list the computational results in Table 3.2, where we take h = e = 0.001 for the case of a,\ = 0, and h = 0.001, e = 1 0 - 1 0 for the case of « i 7^ 0, and use the explicit second order Runge-Kutta scheme to easily see the iteration improvement (Ij stands for results of the jth iteration). From Table 3.2, we see the error's deterioration for the Baumgarte method and the SRM with a.i ^ 0. The SRM with a.\ = 0 performs better in the singular case. • Next we try an example in which y is unbounded at the singularity. Example 3.5 Consider the nonlinear DAE (3.1) with ^ ^ - z i + x 2 - s i n ( < ) - ( l + 2*)j ^ ^ 0 j g = x\-\- #i (#2 — sin(t) — 1 + 2t), Chapter 3. SRM for Nonlinear Problems 78 methods error at —•) t = . l t = .3 t = .5 t= .7 t= 1.0 a i = 1 ex .39e-6 .13e-5 .12e-3 .14e-3 .V6e-4 drift .24e-6 .16e-6 .10e-7 .39e-6 .75e-6 a i = 0 (11) ex .46e-3 .32e-3 .43e-4 .49e-3 .20e-2 drift .24e-3 .89e-4 .18e-8 .20e-3 .22e-2 a i = 0 (12) ex .81e-6 . l l e -5 .41e-5 .29e-5 .68e-5 drift .24e-6 .30e-6 .15e-10 .13e-5 .76e-5 a i = 0 (13) ex .23e-6 .26e-6 .34e-6 .29e-6 .29e-6 drift .90e-9 . l l e -8 .78e-13 .35e-8 .18e-7 a i = 0 (14) ex .23e-6 .26e-6 .36e-6 .27e-6 .29e-6 drift .47e- l l .33e- l l .10e-12 .29e- l l .28e-10 Baumgarte ex .43e-6 .45e-6 .34e-3 .39e-3 .21e-3 drift .24e-6 .16e-6 .61e-7 .24e-6 .75e-6 Table 3.2 Example 3.4 - boun ded y and singularity at t* = .5 subject to the initial condition xi(0) = 1, £2(0) = 0. The exact solution is xe = ( 1 — 2t, s'mt), ye = —cos£ / ( l — 2t). Taking the same parameters and using the same method as before, we get the results listed in Table 3.3. Clearly, the SRM with ct\ = 0 performs well for this situation, while methods error at —> t = . l t = .3 t = .5 t = .7 t = 1.0 S R M ( a i = 0) ex .40e-6 .25e-6 .14e-6 .46e-7 .60e-7 (13) drift .25e-8 .76e-9 .16e-15 .28e-9 .40e-9 Baumgarte ex .49e-7 .15e-6 .93e+l NaN NaN drift .39e-7 .59e-7 .52e+13 NaN NaN Table 3.3: Example 3.5 - unbounded y and singularity at t* = .5 Baumgarte method blows up upon hitting the singularity. • Our next example tests the formulation (3.33)-(3.35) or (3.44)-(3.45) for index-3 problems. Example 3.6 This example is made up from Example 2 in [9] (see Figure 3.1), which describes a two-link planar robotic system. We.use the notation of (1.11). Let Chapter 3. SRM for Nonlinear Problems 79 (*i> yi) fa, y?) q = (#1, 82)T and M = Figure 3.1: Two-link planar robotic system + m2(l\ + / 2 / 3 + kl2c2) m2(l22/3 + hl2c2/2) \ m2(lH3 + hl2c2l2) m2l\j3 ) where li = l2 = 1, m\ = m2 = 3 and c2 = cos 62. The constraint equation is g{q) = h sin 0i' + l2 sm(8t + 82) = 0. We choose the force term f = (Ii cos 9i + l2 cos(^! + #2)) cos t — 3 sin t l2 cos(#! + O2) c o s t-+ (1 — |c2) sin £ which yields the exact solution 8^ = s'mt, 82 = — 2 s in t and A = cost. Because M is (symmetric) positive definite and B = M~lGT -we can take E = I in the SRM formula (3.44)~(3-45)- Again we use the second- order explicit Runge-Kutta scheme, and set h = 0.001, e = 0.005. The results are listed in Table 3.4, where eq and ev stand for maximum errors in q and v = q', resp., and pdrift and vdrift stand for drifts at the Chapter 3. SRM for Nonlinear Problems 80 position and velocity level, resp. We see that the accuracy is improved significantly by the first two iterations. The third iteration is unnecessary here, because the error is already dominated by the Runge-Kutta discretization error. Qualitatively similar results are obtained for E = (GB)T and E = (GB)'1. More interestingly, though, for E = I we neither form nor invert GM~XGT, so a particularly inexpensive iteration is obtained. methods e iteration error at —>• t = . l t=.5 t=1.0 E = 1 5e-3 1 eq .41e-4 .66e-3 .26e-2 ev .75e-2 .74e-2 .69e-2 pdrift .22e-4 .28e-4 .'22e-4 vdrift .49e-2 .41e-2 .27e-2 2 eq .13e-6 .66e-6 .36e-6 ev .19e-5 .81e-6 .20e-4 pdrift .42e-9 .13e-7 .17e-6 vdrift .91e-7 .21e-5 .21e-4 3 eq .10e-6 .58e-6 .12e-5 ev .86e-6 .10e-5 .16e-5 pdrift .96e- l l .60e-9 .48e-8 vdrift .10e-8 .99e-7 .59e-6 Table 3.4: Errors for Example 3.6 using S R M (3.45)-(3.46) • Next we solve for the dynamics of the slider-crank mechanism described in E x a m -ple 2.1. this is a nonlinear index-3 D A E with isolated, "smooth" singularities. Example 3.7 We take h = e = 0.0001 and use the explicit second order Runge-Kutta method again. Singularities are located at (<f>i,<f>2) = (f>^f) (i-e-, each time the periodic solution passes this point). Corresponding to the case shown in [94], we choose </>i(0) = 7-f and </>i(0) = 0 and compute 6i = </>!- y , 02 = cf>2 + - , 9[ and 9'2. Using the formulation (3.47), (3-46), we calculate until t = 70 without any difficulty (see Figure 3.2). Chapter 3. SRM for Nonlinear Problems 81 0 10 20 30 40 50 60 70 Figure 3.2: Solution for slider-crank problem with singularities We also list the drift improvement as a function of the SRM iteration in Table 3.5. iteration number position drift at t=30 velocity drift at t = 30 I .669e-8 .b71e-4 2 .730e-ll .731e-7 Table 3.5: Drifts of the S R M for the slider-crank problem // we use the SRM formulations considered in §§5.5 and 3.4 for problems with-our singularities, or one of the usual stabilization methods with strict tolerances, the results become wildly different from the correct solution after several periods. Next we calculate the acceleration of the slider end in the horizontal direction under the initial data (j>i(0) — \ and (p[(0) = 2\/2. The same problem was discussed in [19]. The result shown in [19] is not perfect since the maximum and minimum values in each period appear to differ. Our result looks better (see Figure 3.3). • Chapter 3. SRM for Nonlinear Problems Figure 3.3: Acceleration of slider end Chapter 4 S R M for the Nonsta t ionary Incompressible Navier-Stokes Equat ions 4.1 D A E Methods for Navier-Stokes Equations While a significant body of knowledge about the theory and numerical methods for D A E s has been accumulated, not much has been extended to partial differential-algebraic equations ( P D A E s ) . The incompressible Navier-Stokes equations form, in fact, an example of a P D A E : to recall, these equations read ut + (u • grad)u = jiAu — gradp + f, (4.1a) divu = 0, (4.1b) u\da = b , u | t = 0 = a, (4.1c) in a bounded two- or three-dimensional domain f i and the time interval 0 < t < T. Here u(x, t) represents the velocity of a viscous incompressible fluid, p(x, t) the pres-sure, f the prescribed external force, a(x) the prescribed ini t ia l velocity, and b(t) the prescribed velocity boundary values. The system (4.1) can be seen as a partial differential equation with constraint (4.1b) with respect to the time variable t. C o m -paring wi th the D A E form (3.1) p corresponds to y, the grad operator corresponds to the matr ix B and the div operator corresponds to the Jacobian matr ix G. It is easily verified that (4.1) has index-2 since the operator diugrad = A (corresponding to the matr ix GB) is invertible (under appropriate boundary conditions). Indeed, the pressure-Poisson reformulation of (4.1) (see, e.g., [55]) corresponds to a direct index reduction of the P D A E , i.e. a differentiation of the constraint with respect to t followed by substitution of u* from the momentum equations. 83 Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 84 In this chapter we propose and analyze a sequential regularization method ( S R M ) for solving the incompressible Navier-Stokes equations. The method is defined as follows: wi th po(x,t) an in i t ia l guess, for s = 1,2, • • • , solve the problem e(us)t - g r a d ( a i ( ^ u u s ) i + a2divus) + e(us • grad)u s = epAus - egradp s_! + ef, u s | an = b, u s | t = 0 = a, ps - ps-i - ^(a1(divus)t + a2divus). This method is an extension of the S R M which was proposed and analyzed in previous chapters for ordinary D A E s , especially for the index-2 D A E s (3.4), (3.5). Here we can take E = I even for ct\ = 0 because (4.1) corresponds to (3.1) wi th B = GT. It is indicated in §3.1 that if we take a.\ ^ 0 then certain restrictions on choosing the ini t ia l iterate (cf. ( 3.14)) do not apply and, more importantly, the equation for xs is essentially not stiff if the original problem is not stiff. Hence, a non-stiff t ime integrator can be used for any regularization parameter e. For the Navier-Stokes application (4.2) we therefore choose ot.\ > 0 so that we can st i l l take e to be very small even when we use an explicit time discretization. So one S R M iteration is often good enough. However, we should not ignore the choice ct\ = 0 because §3.1 also indicates that with this choice the computation can be particularly simple. For (4.2), when a.\ > 0, although we use explicit time discretization, a symmetric positive definite system relevant to the discretization of the operator / + ^-grad div s t i l l needs to be inverted. If we take a.\ = 0, then we do not need to solve any system to obtain the discrete solution. In this case, (4.2) is not stiff only for relatively large e. So more than one S R M iterations are required generally. In the sequel, the convergence proof (4.2a) (4.2b) (4.2c) Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 85 and discretization stability analysis in §4.3 and §4.4 are mainly for the case of cti > 0. The discussion for the case of a\ = 0 can essentially be carried out in a similar way. We wi l l remark on this case in §4.3 and §4.4 and provide a numerical verification in §4.4. The importance of the treatment of the incompressibility constraint has long been recognized in the Navier-Stokes context. A classical approach is the projec-tion method of [36], where one has to solve a Poisson equation for the pressure p with the zero Neumann boundary conditions which is, however, non-physical. Re-cently, a re-interpretation of the projection method in the context of the so-called pressure stabilization methods, or more generally, "pseudo - compressibility methods" has been given in [97]. Some convergence estimates for the pressure can be obtained (cf. [101, 96]). In his review paper [97], Rannacher lists some well known examples of "pseudo-compressibility methods" (which are actually regularization methods): divu + ept = 0 , in ft x [0, T ) , p\t=o = Po, (artificial compressibility) divu + ep = 0, in ft x [0, T ) , (penalty method) divu — eAp = 0 , in ft x [0, T ) , f^|an = Po (pressure stabilization). If we generalize Baumgarte's stabilization to this P D A E example (4.1), we get u t + (u • grad)u = pAu — gradp + f, (4.3a) (divu)t + j divu = 0. (4.3b) El iminat ing ut from (4.3), we obtain an equation for p: —Ap + jdivu — div{(u • grad)u — pAu — f} = 0. We then find that this stabilization can be seen as a kind of pressure stabilization wi th 7 = e - 1 . Although it works, since we do not have a singularity here, it s t i l l sets up a non-physical boundary condition for the Poisson equation for p. Also, in this formulation, equations for u and p are not uncoupled. Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 86 In the S R M formulation (4.2) we do not need to set up boundary conditions for p. So it should be more natural than various pressure-Poisson formulations. This method relates to the idea of penalty methods but, unlike them, the regularized problems are not stiff for cti > 0 or less stiff for ctx — 0 since we can choose e to be relatively large. Hence, more convenient (nonstiff) methods can be used for time integration, and nonlinear terms can be treated easily. We wi l l indicate in §4.4 that e has l i t t le to do wi th the stability of the discretization there, i.e. the stability restriction is satisfied for a wide range of e for « i > 0. We also indicate there that, in the case of small viscosity, the usual time step restrictions for the explicit schemes can be loosened. A similar procedure following [5] (Uzawa's iterative algorithm ) in the framework of optimization theory and economics has actually appeared in the Navier-Stokes context for the stationary Stokes equations (i.e. without the nonlinear term and the time-dependent term in (4.1)) with c*i = 0 using the augmented Lagrangian idea, see Fort in and Glowinski [50]. (Also see [54] for a related discussion.) Note that, in their procedure, e _ 1 in (4.2c) is replaced by a parameter p. They prove that p = e _ 1 is approximately optimal. For the nonlinear case, they combine Uzawa's algorithm with a linearization iteration. They claim convergence but find it hard to analyze the convergence rate because their analysis depends on the spectrum of an operator which is non-symmetric in the nonlinear case. For the nonstationary case (4.1), the augmented Lagrangian method cannot be applied directly. Therefore, [50] first discretizes (4.1a) with respect to the time t (an implici t scheme is used). Then the problem becomes a stationary one in each time step. Hence, Uzawa's algorithm can be applied and converges in each time step. So, for the nonstationary case, their iterative procedure is, in essence, to provide a method to solve the time-discretized problem. Thus, their iterative procedure has litt le to do with time-discretization, or in other words, they st i l l do time-discretization directly for the problem (4.1). Consequently an implici t scheme is always appropriate because of the constraints Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 87 (4.1b), and a linearization is always needed to treat the nonlinear case. These properties are not shared by our method. We wi l l prove that the conver-gence results of the S R M in the previous chapter st i l l hold for the P D A E case (4.2). Hence, the solution sequence of (4.2) converges to the solution of (4.1) with the error estimate of 0(es) after the sth iteration. Therefore, roughly speaking, the rate is about 0(e). We prove the convergence results using the method of asymptotic ex-pansions which is independent of the optimization theory and is also applicable to the steady-state case. In addition, when the finite element method is used, the difficulty of constructing test functions in a divergence-free space to decouple the u,p system can be avoided by using the formulation of the S R M . We indicate here that, as many others do, we include the viscosity parameter u in the error estimates. So the estimates could deteriorate when u is very small . This is because we have an unresolved technical difficulty, associated with our inabil i ty to obtain an appropriate upper bound for the nonlinear term and wi th the weaker elliptic operator /J,AU (which is dissipative) as p —>• 0. In the S R M formulation, a supplementary dissipative term — ci2graddivus is introduced without perturbing the solution. As indicated in [50] for the stationary case, the relative advantage of such methods may therefore become more apparent for small values of the viscosity. The chapter is organized as follows: In §4.2 we define some preliminaries and discuss regularity properties of the solution of (4.2). The convergence of the S R M for Navier-Stokes equations is proved in §4.3. Final ly, in §4.4 a simple difference scheme is discussed and some numerical experiments are presented. These numerical experiments are only exploratory in nature. To summarize, our objective in this chapter is to present a method for the non-stationary Navier-Stokes equations from the viewpoint of D A E regularization, and to provide a way to apply a D A E method to P D A E s . It appears that such a formulation is new in the Navier-Stokes context and it is worthwhile because: Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 88 • Since e need not be very small, the regularized problems in the sequence (4.2) are more stable/less stiff. So more convenient difference schemes, e.g. explicit schemes in time, can be used under theoretical assurance. If we take cxi > 0, this is also true for small e. • The problem of additional boundary conditions, which arises in the pressure-Poisson formulation and projection methods, does not arise here. Fin i te element methods can be used easily and the elements do not have to conform to the incompressibility condition to separate the variables u and p. 4.2 Pre l iminar ies and the Propert ies of the Regular ized Prob lems Before we begin our analysis, we'first describe some notation and assumptions. As usual, we use L P ( H ) , or more simply L p , to denote the space of functions which are pth-power integrable in f i , and as its norm, where u = • • •, i t n ) . We denote the inner product in L 2 by (-, •) and let || • || = || • | | 2 . C°° is the space of functions continuously differentiable any number of times in $7, and C£° consists of those members of C°° with compact support in $7. H m is the completion in the norm We wi l l consider the boundary conditions to be homogeneous, i.e. b = 0 in (4.1c), to simplify the analysis. Nevertheless, through the inclusion of a general forcing term, the results may be generalized to the case of nonhomogeneous boundary values. We are interested in the case that (4.1) has a unique solution and the solution belongs to H 2 , where the arbitrary constant which the pressure p is up to is determined by U H m = ( E \\Dau\\')*. 0<|a|<m (4.4) Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 89 Hence, some basic compatibili ty condition is assumed (cf. [63]): a|an = 0, diva. = 0. (4.5) Furthermore, we assume sup | |f || < Mu \\a\\H2 < Mi *6[0,T] (4.6) where M i is a positive constant. We take po in (4.2) satisfying (4.4). Hence, it is easy to see that ps satisfies (4.4) for al l s. For simplicity, we only consider the two-dimensional case. We can treat the three-dimensional case in the same way, possibly with some more assumptions. Throughout the chapter M represents a generic constant which may depend on p as we have explained in the introduction. We wi l l also allow M to depend on the finite time-interval T since we are not going to discuss very long time behavior in this chapter. A t first, we write down some inequalities: • Poincare inequality: u|| < 7||grad u||, if u | n = 0. (4.7) More generally (see [87]), for u £ Hl(fl) • Young's inequality: abc < -ap + -bq + -cr p q r (4.9) if a, 6, c > 0, p, q, r > 1 and ^ i + i + i = l. p g r • Holder's inequality: (4.10) if p, q, r > 1 and ^ i + I + i = l. p g r Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 90 • Sobolev's inequality in two-dimensional space: | | u | | 4 < 7 * H * | | g r a d u | | * , (4.11) where 7 l = 2 if ft = R2. Suppose that w stands for the difference of two solutions of the S R M (4.2). Then w satisfies the homogeneous problem (4.27) (see next section). Hence, using the estimate in Lemma 4.2, uniqueness of the solution of the S R M (4.2) is easy to consider. The existence can be analyzed by following the standard existence argument for Navier-Stokes equations (e.g. [108, 62]) and that of penalized Navier-Stokes equations (e.g. [28]). Hence we assume the existence of the solution of the S R M and concentrate on the proof of the convergence of the method. Before we do so, we derive the following regularity results of the solution. L e m m a 4.1 For the solution {u s , p s } of (4-2) there exists a constant to such that when e < e0 we have the following estimates K I I H I + [T (^r\\(dlvus)t\\w + ^T\\divus\\2ul + | | (u s ) t | | 2 + | | A u s | | 2 + \\pa\\2Hi) dt Jo ez tL < M[\H\w + [T ( l l / l l 2 + l l g r a d ^ U 2 ) dt}. (4.12) Proof: For simplicity of notation we denote u s as v here. The proof for the case a i = 0 is just the same as that in [28]. So we only consider the case ax > 0. Hence, without loss of generality, we take « i = 1 and a2 = ot. We then write (4.2) as \ t — -grad((diuv)i + adivv) + (v • grad)v = JUAV - gradp s _i + f, (4.13a) v|an = 0, v | i = 0 = a, (4.13b) ps = ps-i —-((divv)t + adivv). (4.13c) Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 91 The proof follows the ideas in [28]. Mul t ip ly ing (4.13a) by v and integrating wi th respect to the space variables on the domain 0 , we get - — llvll2 H -||dzuv||2 -\ lldzuvll2 + jullgrad vII2 = - ( ( v • grad)v, v) - ( g r a d p s _ i , v ) + ( f , v ) < f l l ^ v l l 2 + ^ | | v | | 2 | | g r a d v | | 2 + ^ ||v||2 + ^ llgradp^H2 + %|| 2, le la lj n p where we use — ((v • grad)v, v) = |((rfiuv)v, v). Then let c = m i n ( " , ^ ) and Y = ||v||2 + ^ ||rfzuv||2. Using Poincare's inequality (4.7), we obtain ±Y + cY + \{u - ^F)||grad v||2 < l(\\f||2 + Hgradp^||2). (4.14) at l a u Note that Y(0) = ||a||2. Write (4.14) as ±(H - ^Y) - ^ | | g r a d v | | 3 ( > - —Y) > -^(Ugradp^ ||2 + ||f||2). at a la a ua Apply ing a standard technique for solving linear differential equations and taking e appropriately small (< eo) so that ^ _ 1HY(0) > £ and ?11Z fT (||ff + llgradp,.! a I aa Jo ^ dt < p Jo 4 we get P-— n < ) > 7 V t 6 [ 0 , T ] . (4.15) a 4 Then, using the same technique and (4.14), we have Y < ||a||2exp(-ct) + M e x p ( - c t ) / (||f ||2 + ||gradp s_i|| 2)exp(cz) dz Jo < M[||a||2 + f (||f||2 + llgradp.^H2) dz]. (4.16) Jo Thus, (4.12) holds for ||u||2. Integrating (4.14) directly and using (4.15) yields f ||gradu||2 dz < M[||a||2 + /* (||f||2 + Ugradp^||2) dz]. (4.17) Jo Jo Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 92 To prove other estimates for (4.12) we define the operator A w = —-grad ( (d iuw) t + adiww) — pAw = g, (4-18) where w satisfies w|gn = 0 and w | < = u = a. Let q = ((divw)t + adiww). Then we have (noting diww\t=o = 0) —pAw + gradg = g, (4.19a) divw = —e qexp(—a(t — z)) dz (4.19b) Jo This is a general nonhomogeneous Stokes problem. Using the results described in [28] (or cf. [108]), we get || A w | | + ||gradq|| < M[||g|| + e /* ||grad?|| exp(-a(< - z)) dz). (4.20) Jo Apply ing Gronwall inequality, it is easy to obtain ||gradg|| < M(||g|| + e f ||gradg|| dz) (4.21) Jo and /* ||grad?||2 dz<M f ||g||2 dz = M f* | | A w | | 2 dz. (4.22) Jo Jo Jo It thus follows that | A w | | < M(||g|| + e f ||g|| dz) = M(\\Aw\\ + t T | | A w | | dz), (4.23) Jo Jo and then / " | | A w | | 2 dz<M f ||g||2 dz = M f | | A w | | 2 dz. (4.24) Jo Jo Jo From (4.19b) and (4.22), we thus have /' ||graddww||2 dz = f ||gradg||2 dz < M f | | A w | | 2 dz. (4.25) JO Jo Jo 1 fl ~7i Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 93 Then 1 fl e2 / | |grad(>'uw) t|| 2 dz < M [ \\Aw\\2 dz (4.26) Jo Jo follows from (4.18). Now taking the scalar product of (4.13a) with Av, we have — \\(divv)t\\2 + — ^ -\\divv\\2 + -4-||grad vll 2 + | | A v | | 2 2e l l v ; 11 2edt" 11 2 dt"& 11 11 11 = -((v • grad)v, A v ) - (gradp,_i, A v ) + (f, A v ) . Note that i i i — ((v • grad)v, A v ) < ||v||4||grad v | | 4 | | A v | | < 7* ||v|| = ||grad v | | | | A v | | 2 | | A v | | < ^ ( | | A v | | 2 + | | A v | | 2 ) + jL(||v||2||gradv||2)||gradv||2 < <S[||Av||2 + M 2 | | A v | | 2 + M 2 e 2 ( f \\Av\\ dz)2} + -^(||v||2||grad v||2)||grad v||2, Jo loo'3 where we use (4.23) for the last inequality. Recall that we have estimates for | | v | | 2 and /0* ||grad V||2 dz. Therefore, taking 8(1 + M2) < | , it is not difficult to obtain ||gradv|| + - / \\(divv)t\\2 dz + | | A v | | 2 dz e Jo Jo < M[||grada||2+ /* (||f ||2 + Hgradp^||2) dz + e2 f \\Av\\2 dz] Jo Jo Taking e such that M e 2 < 1, we then get (4.12) for /0* | | A v | | 2 dz and ||gradv||. Noting (4.25), (4.26), from (4.2c), (4.12) also holds for /0* ||gradps||2 dz. App ly ing the inequality (4.8) and noting that ps satisfies (4.4), yields the bound for JQ \\ps\\2 dz. Hence, from (4.2c) we can obtain (4.12) for J0f \\divv\\2 dz and then / J | |(diuv)*| | 2 dz. We thus complete the proof. • From this lemma, we see that if we choose po such that J*Q ||gradpo||2 dz is bounded then, by induction, all terms in the left of (4.12) are bounded for any given s. Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 94 4.3 Convergence of the S R M In this section, we estimate the error of the S R M (4.2) in the solution of (4.1) by using the asymptotic expansion technique as in the proof of Theorem 2.1 of Chapter 2 (see §2.8). Note that, in the Navier Stokes context, the asymptotic expansion method was used in [54] to obtain a more precise estimate for a penalty method for the stationary Stokes equations and in [108] to calculate a slightly compressible steady-state flow. We wi l l mainly consider the case cx\ > 0. Hence, we take a.\ = 1 and a2 = ex for convenience. The result for a\ = 0 wi l l be described in Remark 4.3. A t first we discuss a couple of linear auxiliary problems. Then we go to the proof. 4.3.1 Two linear auxi l ia ry problems We discuss two linear problems in this section. One is ew< — grad(rfiuw)i — agraddivw + e(w • g rad)U +e(V • grad)w — epAw — egradg + ef, (4.27a) w | 9 a = 0, w | i = 0 = 0, (4.27b) where U , V and q are given functions. The other is w t + ( V • grad)w + (w • grad)V = pAw — gradp + / , (4.28a) (divw)t + adivw = g, (4.28b) w | 9 a = 0, w | t = 0 = a, (4.28c) where V , g and a are given functions, a satisfies the compatibil i ty conditions (4.5) and g satisfies (4.4). Now we show some properties of these two problems which wi l l be used later in the proof of the convergence of S R M . L e m m a 4.2 For the solution of problem (4-27), z ' /U and V satisfy I I • l l H 1 + fT I I • H H * dt < M , (4.29) Jo Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 95 then we have the following estimate e\\w\\2 + | | d ivw| | 2 < M e /* (||f | | 2 + | | ? | | 2) ds, (4.30a) Jo e | | g radw | | 2 + /* (e| |w,| | 2 + | |divw«| | 2 ) ds < Me f (||f|| 2 + | | g | | 2 ) ds. (4.30b) Jo Jo P r o o f : Mul t ip ly ing (4.27a) by w and then integrating on the domain ft yields 1 d 1 d 2 £ <f t" W " 2 + 2<ft " d * W W " 2 + aWdivwW2 + e^ l | g r adw | | 2 = —e((w • g r a d ) U , w) — e((V • g r ad )w , w) + e(q, divw) + e(f, w) < e | | g r a d U | | | | w | | 2 + | | | A W | | | | w | | 2 + e(q, divw) + e(f, w) (using (4.10)) i 1 < e 7 l 2 ( | | g r a d U | | + - | | d i u V | | ) | | w | | | | g r a d w | | + e(q,divw) + e(f ,w) (using (4.11)) 1 £ 1 < - e A i | | g r a d w | | 2 + ^ ( | | g r a d U | | + - | | ^ w | | ) 2 | | w | | 2 + e(q, divw) + e(f, w ) , where we have used — e((V • g r a d ) w , w) = ^((divV)w, w). Therefore, we have ^ ( e | | w | | 2 + \\divw\\2) - C ( i ) ( e | | w | | 2 + \\divw\\2) < -e / j | | g rad w | | 2 - (a + C(t))\\divw\\2 + 2e(q, divw) + 2e(f, w ) /f_,, M2 , _7 2 l l r , i 2 , . ! i/2 < - e ^ | | g r a d w | | 2 + e ^ | | w | | 2 + £ 4 - | | f | | 2 + e - | ' - 1 1 2 7Z \i a where < e ^ | | f | | 2 + e i | | g | | 2 , . (4.31) ji a 7l n i J T T I I , 1 n j , . . ^ n i N 2 C7(t) = - ^ ( | | g r a d U | | + - | | ^ V | | ) Not ing that w | t = 0 = 0 and divw\t=0 = 0, we thus get (4.30a). Now, mult iplying (4.27) by w t , then integrating with respect to x over ft, we get e | |w t | | 2 + H^uwil l 2 + ^Jt\\div™\? + e / " ^ l | g r a d w | | 2 = e((w • g r a d ) U , w t ) + e((V • g r ad )w , w t ) + e(q, divwt) + e(f, w t ) . (4.32) Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 96 We use the inequalities listed in the previous section to estimate the right-hand side of (4.32) and have the following: e ( (w-grad )U ,Wi ) < ^ | |w t | | 2 + Me( | |grad w | | 2 + | |w|| 2), e((V • grad)w, wt) < - | | w t | | 2 + M e s u p | V | 2 | | g r a d w | | 2 , . 4 n c(f,w«) < | | |w« | | 2 + Me| | f | | a , e(q,divwt) < -\\divwt\\2 + M e | | g | | 2 , where the bounds of s u p n | U | and s u p n | V | can be obtained by using the inequality sup | • | < M\\A • || a (see, e.g. [119]). Then, similarly to the procedure for obtaining (4.30a), we obtain (4.30b). • Next we consider problem (4.28). L e m m a 4.3 There exists a solution for problem (4-28). Moreover, for the solution of (4-28), we have the following estimate: rT W I H I + / ( H w H k + H w . l l 2 + \\p\\2m) dt<M (4.33) Jo i f So" | |f l | 2 dt and J 0 T Hfirll^i dt are bounded. Proof: First , we can solve divw from (4.28b) (noting that divw\t=Q = 0): divw ~ gi, (4.34) where gi = exp(—at) / gexp(as)ds (4.35) Jo satisfies (4.4) since g does. B y applying Corollary 2.4 in [54, p.23], the problem divw = p i , (4.36a) w\da = 0 (4.36b) Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 97 has many solutions. We pick one and denote it as w p . Then w := w — w p satisfies the linearized Navier-Stokes equations in the form of (4.28a) with a proper force term (denoted by f) and divw = 0, w|an = 0 and w | < = 0 = a - wp\t=0. B y noting that divwp\t=0 = gi\t=o = 0, the basic compatibili ty conditions like (4.5) for w are satisfied. From (4.35) and the assumption for g, J* 0 T (U^iH^i + ||(<?i)<||2) dt is bounded. Hence, based on the estimates for the solution of (4.36) (see [1] and [54]), it is not difficult to get I K | | H I + [T (||(wp),||2 + \\wp\\2n2) dt < M. (4.37) Jo Thus JQ ||f|| 2 dt is bounded. Simulating the regularity argument of [62] or [63] (multiplying the linearized Navier-Stokes equations by w, wt and P A w where P is a projection operator (cf. [62]), respectively), we can obtain | | w | | H i + [T + l l w i H 2 + \\p\\w) dt < M. (4.38) Jo Therefore, (4.33) follows from (4.37) and (4.38). Using the estimate (4.38) and fol-lowing a global existence argument (e.g. [62] or [108]), the existence of the solution for w can be obtained. We thus have the results of the lemma . • R e m a r k 4.1 The uniqueness of the solution of (4-28) follows from the standard ar-gument for Navier-Stokes equations (cf. [108]). • 4.3.2 The error estimate of S R M In this section we prove the convergence of iteration (4.2) based on the same procedure described in the proof of Theorem 2.1 of Chapter 2 (see §2.8) . We describe our results in the following theorem. Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 98 Theorem 4.1 Let u and p be the solution of problem (4-1), a n d u s and ps be the solution of problem (4-2) at the sth iteration. Then there exists a constant eo such that when e < e0 we have the following error estimates for all t £ [0, T]: | | u - u , | | H i < Mes, (4.39a) (/ \\p - Pa\\2 dt)* < Me% (4.39b) Jo where T is any given finite number and s = 1,2, • • •. Proof: A t first, consider the case s = 1 of (4.2). Let Ui = uio + e u n H h emulm H Comparing the coefficients of like powers of e, we thus have grad((<fo'uuio)t + adivuw) = 0, (4.40a) grad((efo'uun)t + adivuu) = [uw)t + (uw • grad)ui 0 -pAuw + gradpo - f, (4.40b) i - l grad((dzuui,-)t + adivuu) = (ui,-_i) t + E ( u i i ' grad)ui,-_i_j -pAuu-i , 2 < i < m + 1, (4.40c) where (4.40a) satisfies (4.2b) and (4.40b) and (4.40c) satisfy the homogeneous in i t ia l and boundary conditions corresponding to (4.2b). Now (4.40a) has infinitely many solutions in general. We should choose Ui 0 not only to satisfy (4.40a) but also to ensure that the solution of (4.40b) exists. A choice of Uio is the exact solution u of (4.1), i.e. (uio)t + (uio • g r ad )u 1 0 = ^ A u i o - gradp + f, (4.41a) (div\iw)t + adivuio = 0, (4.41b) uiolsn = 0, u 1 0 | < = 0 = a. (4.41c) Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 99 Note that divuw\t=0 = diva. = 0 and p is taken to satisfy (4.4). So u 1 0 = u and (4.40b) has the form grad((divun)t + adivun) — grad(p0 - ?)• (4.42) Now we choose U n and a corresponding p\\ to satisfy (un)i + (uio • grad)un + (un • grad)uio = ^ A u n - gradpu, (4.43a) (diuunjt + adivuw = p0 — p, (4.43b) unlsn = 0, u n | t = 0 = 0. (4.43c) Aga in we have divun\t-o = 0 and let pn satisfy (4.4). According to Lemma 4.3, Un and pn exist . Generally, supposing we have Ui;_i, pu-i for i > 2, choose Ui;, pu satisfying (u i i ) i + (mo • grad)uii + (u u - • grad)u 1 0 i-l = uAuu - gradp l t - ^ ( u j j • grad)ui,_i_ J , (4.44a) (divuu)t + adivuu = —pu-i, (4.44b) ui.-lan = 0, u i , - | t = 0 = 0, (4.44c) where we note that divuu\t=o = 0 and pu satisfies (4.4). App ly ing Lemma 4.3, al l uu and pu, i = 0,1, • • •, exist and satisfy (4.33). Next we estimate the remainder of the asymptotic expansion after the (m + l ) t h power of e. Denote u i m = uio + eun + • • • + e m + 1 u i m + i (4.45) (viim also satisfies (4.33)) and Wim = U i - U i m . (4.46) Then W i m satisfies e ( w l m ) t - g r a d ( d w w i m ) t - agraddivwlm Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 100 +e(w i m • grad)uj + e(u X m • g r a d ) w i m = epAwlm -m+1 e m + 2 { ( u i m + 1 ) i + [(ui,- • g r a d ) u i m + i _ , - ] - , u A u i m + i } , (4.47a) t=0 w i m | 9 n = 0, w i m | t = 0 = 0. (4.47b) Using regularity we have for Uu, i i i m and U i (see (4.12)) and Lemma 4.2 , we obtain | |w l m | | = 0 ( e m + 1 ) and | |gradw l m | | = 0 ( e m + 1 ) . Therefore U l = uio + c u u + • • • + e m u l m + 0(em+1). (4.48) in the H 1 - n o r m for the spatial variables. Noting U i 0 = u, we thus obtain U l - u = 0(e). (4.49) Furthermore, according to Lemma 4.2, we have | | ^ w l m | | = 0 ( e m + l ) , ( fT | | (dtt ;w l m )t | | 2 dt)* = 0 ( e m + l ) . Jo Then, by using (4.2c),(4.48),(4.41b), (4.43b), (4.44b) and the estimates for divwlm and ( & w l m ) ( , it follows that Pi = P + mi + • • • + e > i m + 0 ( e m + " ) (4.50) or Pi -p= 0(e). (4.51) in the sense of the L 2 - n o r m for both spatial and time variables, i.e. ( J 0 T || • | | 2 dt)*. Now we look at the second iteration s — 2 of (4.2). Let u 2 = u 2 0 + eu 2i H h e m u 2 m H . Noting that (4.50) gives us a series expansion for pi we obtain grad((cfem20)t + adivu2o) = 0, (4.52a) Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 101 grad((divu21)t + adivu21) = (u20)t + (u2o • grad)u20 -pAu20 + gradp-f, (4.52b) i-l grad((divu2i)t + adivu2t) = (u2,--i)t + J2(U^J " grad)u2 i_i —pAu2i-i - gradpi,_i,2 < i < m + 1. (4.52c) Again , (4.52a) is combined with the ini t ia l and boundary conditions (4.2b), and (4.52b) and (4.52c) are combined with the corresponding homogeneous ones. As in the case of s = 1, we again choose u2o = u. We thus have grad((fifo"uu2i)t + adivu2i) — 0. (4.53) Then u 2i is constructed to satisfy (u2i)i + (u 2 0 • grad)u2i + (u21 • grad)u20 = ^ ^ u2i - gradp21, (4.54a) (divvL2i)t + adivu21 — 0, (4.54b) u 2 i | 9 n = 0, u2i|*=o = 0. (4.54c) Obviously u 2 i = 0 and p2i = 0 is the solution of (4.54) and (4.4). In general, similarly to the case of s = 1, we choose u2;,p24- to satisfy (u2,-)t + (u 2 0 • grad)u2j + (u2i • grad)u20 = (4.55a) i-l pAu2l - gradp2; - X!( u 2J ' g r a d ) u 2 t - i - j , (4.55b) (divu2l)t + adivu2i = -p2i-i, (4.55c) u2i|3n = 0, u 2 8 ' | i =o = 0 (4.55d) for 2 < i < m + 1, where p2[ satisfies (4.4). B y the same procedure as for s = 1 we obtain error equations similar to (4.47) with the addition of a remainder term grad(pi — pirn) in the right-hand side, where pim stands for the asymptotic expansion (4.50) of pi. Apply ing Lemma 4.2 again, we get u 2 = u 2 0 + e u 2 1 + • • • + e mu 2 m + 0(em+1). (4.56) Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 102 Noting u2i = 0, u2 -u = 0(e2). (4.57) Then, using (4.2b),(4.56),(4.54b) and (4.55c), we conclude P2 = P + ep2i + • • • + e > 2 m + 0 ( e m + 2 ) (4.58) or p2-p=0(e2). (4.59) since p 2i = 0. We can repeat this procedure, and, by induction for s (choosing m larger than s), conclude the results of the theorem. • R e m a r k 4.2 Corresponding to Theorem 2.1, we expect that the error estimates (4-39) also hold for the SRM (4-2) with cxi = 0, at least, away from t — 0. • R e m a r k 4.3 In Theorem 4-1, we find that the result for p is in a weaker norm lo I " l | 2 dt. This is because we have difficulty in estimating the first order time-derivative of the right-hand side of(4-47), or concretely, the term /0T IKuim+i )«||2 ds. In [63] (Corollary 2.1) it is shown that /0T ||(uim+i)tt||2 ds may be unbounded ast —> 0 if we only assume the local compatibility conditions (4-5). In the case that this integral is bounded for 0 < t < T, we can get (4.60) Otherwise, we only can expect that (4-60) holds away from t = 0 by following the argument in [63]. • R e m a r k 4.4 Multiplying (4-27) by Aw, where A is the operator defined by (4-18), and following the later steps of the proof of Lemma 4-1, we can get Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 103 1 fT f (\\grad(divw)t\\2 + a||gradrfz'uw||2) dt Jo <M / T ( | | f | | 2 + | | g r adg | | 2 ) ^ . Jo Using this result to estimate the remainders of the asymptotic solutions in the proof of Theorem J^.l, we obtain rT (/ Wv-Psfm dt)* <Mes. (4.61) Jo • 4.4 Discre t iza t ion Issues and Numer i ca l Exper iments In previous sections, we have proposed the S R M and performed some basic analysis on i t . The S R M yields a sequence of P D E s which are to be solved numerically. The problem at the 5th iteration can be written as: e(us)t - grad(a 1 (^ 'uu s ) i + a2divus) + e(us • grad)u s = epAus + er s , (4.62a) us\dn = 0 , u s | < = 0 = a, (4.62b) where rs(t) is the known inhomogeneity rs = -gradp s _x + f. (4.63) A variational formulation of (4.62) gives: F i n d u s G HQ such that e-^-(us, 4>) + ^-^-(divUs, div<f>) + a2(divus, divcf>) +e/^(gradus,grad0) + b(us, us, <f>) = e(r s, 0), V0 € HQ, (4.64a) u s | i = o = a , divus\t=o = 0, (4.64b) Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 104 where the trilinear form 6(u,v,w) = (.(u • grad) v,w). From (4.64) we see that the finite element method in the spatial variables, com-bined wi th time discretizations, can be easily adopted. Note that we do not need to construct divergence-free test functions to separate the variables u and p. Neverthe-less, in this section, we are not going to discuss finite element methods further. Some discussions on using the S R M with the finite element method wi l l be given in the next chapter for a problem in reservoir simulations. Some numerical experiments are also given there. Here we only consider a very simple first-order difference scheme (forward Euler scheme in the time direction) in two dimensional space, as an ini t ia l attempt towards the discretization of the sequential regularization method for the P D A E . Concretely, we consider a rectangular domain such that an equidistant mesh can be used. Let (u,v)T stand for the approximation of us, and let k,hx,hy denote step sizes in time and spatial direction, respectively. Without loss of generality, we assume that hx = hy = h and that the domain is a unit square. Thus, mesh points can be expressed as Xi = ih, i = 0,1, • • •, / ; yj — jh,j = 0 ,1 , , • • •, J ; tn = rafc, n = 0,1, • • •, N, N = [T/k]. The difference scheme reads: + (4.65a) (xl{'U'xy + uyy)i (X2{Uxy + Uyy) (4.65b) u\an = 0, v\dn = 0, (4.65c) Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 105 where u Ui = h ity and My can be defined accordingly and the definitions for v are similar. Obviously, this is a first-order scheme explicit in time, where the nonlinear term is discretized somewhat arbitrarily. The scheme is easy to implement. Next we discuss its stability. For simplicity, we analyze the linear case (corresponding to the Stokes equations) first, and consider the full nonlinear equations (4.65) in Remark 4.7 below. We write the linear case of (4.65) as follows : a2(uxx + U y X ) + e/i(uxx + Uyy) + (4.66a) CV2(UXy + Uyy) + Cp(VXX + Vyy) + &T y , (4.66b) u\t=o = au, v\t=0 = av. (4.66c) Here we take a.i — 1 and a2 = a. The result for the case of a.\ = 0 wi l l be given in Remark 4.6. The following theorem gives the stability estimate for (4.66) in the sense of the discrete L 2 - n o r m : t'=0 j=0 where wh = (tOj,j), i = 0,1, • • • , / .— 1, i = 0,1, • • • , J — 1. Theorem 4.2 Let u anc? v be the solution of (4-66) and A = e(| |u| | 2 + \\v\\l) + \\us + vy\\l + ep(\\ux\\2h + \\uy\\2h + \\vx\\2h + \\vy\\2). (4.68) euj — cci(uxx + VyX)i = ev{ - ai(uxy + Uyy)i = u\aa = 0, u | 9 n = 0, Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 106 If f r < 1 — c> where c is any constant in (0,1), then A + kJ^ \\(us + vg)i\\l<Me maxJ\\ru\\2 h + \\rv\\l) (4.69) „ U < v l n < v i n = U - -where M is a generic constant dependent on p and c. Proof: We first write down the following identities and inequalities: • some difference identities [75]: (#)* = <ft/>* + < M £ V , (4.70a) (<f»p)i = Wi + feEfy, (4.70b) 2M>i = (Wi-kfay, . (4.70 c 4>4>xx = (dHf>i)s - (<fe)2, (4.70d) where the translation operator Elx(p(x,y,t) = (f>(x + ih,y,t). • an difference inequality [75]: Wx|U<2|HU (4-71) • a discrete version of the Poincare inequality (cf. [66]): U\\l<U4l + Uv\\l (4-72) if <p satisfies homogeneous boundary conditions. Mul t i p ly ing (4.66a) by au-\-bui and (4.66b) by av + bv^ and adding, then summing for al l ( i , j ) , i = 1, — 1, j = 1, • • •, J — 1 where we use the difference identities (4.70) (omitting lots of tedious algebraic manipulations), we obtain Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 107 ea(\\u\\l + \\v\\l)i + o ( « 6 + a)(\\us + vy\\l)i + -epb(\\ux\\2h + \\uy\\2h + \\vx\\2h + \\vy\\2h)i +e(b - aA:)(||M^||^ + \\vi\\2h) + aa\\ux + vv\\2h + (b - ^(ba + a)k)\\(us + Vy)i\\2h +eu.a(\\ux\\2h + \\uy\\2h + + - ^epbk(\\uxi\\2h + \\ugi\H + \\vxi\\2h + \\Vyi\\2h) < Me{\\ru\\2h + \\rv\\2h) + \epa(\\u\\2 + \\v\\2h) + ^eb5(\\Ui\\2h + | H | 2 ) , where 5 > 0 can be chosen less than c/b. Apply ing (4.71) and (4.72), we get h\{\\uxi\\2h + \\uyi\\2h + \\v-xi\\2h + \\VyiWl) < 4 ( | | n i | | 2 + \\vi\\l) and + \Ml < \\US\\1 + \\Uy\\2 + \\V42 + \\Vy\\l, respectively. Then we can choose a and b such that k 1 1 b- ak - p— - -bS > 0, b - -(ba + a)k > 0 n2 2 2 and obtain {A)i + depA + \\(ux + Vy)i\\l < Me(\\ru\\2h + \\rv\\2), where d is a constant independent of k, h, e and /J,. From this inequality, it is not difficult to see that (4.69) holds. • Remark 4.5 From (4-69) of Theorem 4-2, we find that the value of e will not affect the stability of the difference scheme. This means that the forward Euler scheme in the time direction works for any value of e. Also, the time step restriction k < (1 — c)h2/u. is actually loosened in the case of small viscosity (or large Reynolds number) which people are often interested in. This implies that the explicit scheme (4-66) to which an appropriate discretization of the nonlinear term (see the next remark) is added works very well. It enables us not only to avoid the complicated iteration procedure for nonlinear equations but also to choose the time step fairly widely in the case of small viscosity. • Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 108 R e m a r k 4.6 We have mentioned before that sometimes we may like to take a.\ = 0 to avoid solving any algebraic system. Following the same procedure as the proof of Theorem 4-2, we can get the stability condition for the case of cx\ = 0. That is, k < meh2, where m is a positive constant independent of e,h and p. We thus see that the stability of (4-66) with a.\ = 0 depends on the parameter e. This coincides with our experience with stiff problems discretized by explicit schemes. Fortunately, using the SRM, we do not need to take e very small. So the time step restriction is not much worse than the usual one corresponding to an explicit scheme applied to a non-stiff problem. • R e m a r k 4.7 For the nonlinear case (4-65), when the viscosity p is not small, we expect similar results since the nonlinear term can be dominated by the viscous term. When the viscosity is small, however, the scheme (4-65) is unstable. Although numer-ical computations indicate that we do get better stability if we increase a2, i.e. some kind of dissipation effect is obtained (we must note that such dissipation becomes small when the incompressibility condition is close to being satisfied), we suggest using spa-tial discretizations with better stability properties, e.g. upwinding schemes (cf. [102]), in the case of small viscosity. • R e m a r k 4.8 Applying corresponding difference identities for a nonuniform mesh (see e.g. [106]), the results of Theorem 4-2 may be generalized to difference schemes (4-66) on a nonuniform mesh. Hence, the difference scheme may be used for problems defined on more general domains. • Next we explain our theoretical results by calculating the solution of an artificial example. E x a m p l e 4.1 Consider the Navier-Stokes equations (4-1) with the exact solution Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 109 u = (u,v)T u = 5 0 £ 2 ( l - x ) 2 j / ( l - y ) ( l - 2 y ) [ l + e x p ( - t ) ] , v = - 5 0 y 2 ( l - y ) 2 x ( l - x ) ( l - 2 a ; ) [ l + e x p ( - t ) ] , P = [ - x ( | + 2 ) - y ( | - 2 ) + ^] [ l + e x P H ) ] 2 As indicated in §2.6 of Chapter 2, to carry out the SRM iterations, we do not need to store the entire approximation of ps-i on [0,T] for calculating us. Assuming that the number of the SRM iterations is chosen in advance, we can rearrange the computational order to make the storage requirements independent of N, where N represents the number of the mesh lines in the t direction. We first use constant steps k = 0.01 and h = 0.1. At a given time t, we use 'eu' to denote the absolute discrete L2-error in us while lep' denotes the absolute discrete L2-error in ps. Table 4-1 summarizes the computational results of the difference scheme (4-65) with Q i = Q 2 = 1 and viscosity f.i = 0.1. e iteration error at — • t = k t = 1.0 t = 2.0 t = 3.0 t = 4.0 t = 5.0 5e-l 1 eu 4.65e-3 2.69e-l 1 .55e- l 1.31e-l 1.15e-l l .U8e-l ep 2.49e-l 1.96e-l 1.57e-l 1.36e-l 1.25e-l 1.21e-l 2 eu 2.16e-3 2.53e-2 3.12e-2 3.18e-2 3.12e-2 3.06e-2 ep 1.80e-l 9.28e-2 7.37e-2 6.74e-2 6.48e-2 6.35e-2 3 eu 2.15e-3 1.77e-2 2.28e-2 2.48e-2 2.55e-2 2.57e-2 ep 1.80e-l 8.81e-2 6.69e-2 6.10e-2 5.91e-2 5.83e-2 le-3 1 eu 2.14e-3 1.73e-2 2.21e-2 2.41e-2 2.48e-2 2.50e-2 ep 1.80e-l 8.78e-2 6.61e-2 6.01e-2 5.82e-2 5.75e-2 Table 4.1: S R M errors for p = 0.1 without upwinding We notice that the errors improve as the iteration proceeds until es reaches the discretization accuracy 0(h), where s is the number of iterations. For small viscosity, say p = 0.001, the difference scheme (4-65) does not work. The errors blow up around t = 1. When we increase a2, say to 50, we do get pretty good results around t = 1; however, the errors still blow up at a later time. This Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 110 suggests that the scheme is not stable for small viscosity p. So we next discretize the nonlinear term using the upwinding scheme given in [102]. For the case of small viscosity, e.g. p = 0.001, we get good results (see Table \.2). e iteration error at —} • t = k t = 1.0 t = 2.0 t = 3.0 t = 4.0 t = 5.0 5e~I 1 eu 4.bbe-3 2.26e-I 2.50e-l 2.29e-l 2.13e-l 2.03e-l ep 2.58e-l 1.06e-l 6.67e-2 5.45e-2 5.12e-2 5.04e-2 2 eu 2.16e-3 7.74e-2 8.78e-2 9.13e-2 9.34e-2 9.53e-2 ep 1.84e-l 8.81e-2 6.22e-2 5.39e-2 5.11e-2 5.02e-2 3 eu 2.14e-3 7.69e-2 8.71e-2 9.06e-2 9.29e-2 9.48e-2 ep 1.83e-l 8.78e-2 6.21e-2 5.39e-2 5.11e-2 5.01e-2 le-3 1 eu 2.14e-3 7.69e-2 8.72e-2 9.07e-2 9.29e-2 9.49e-2 ep 1.83e-l 8.78e-2 6.21e-2 5.39e-2 5.11e-2 5.01e-2 Table 4.2: S R M errors for p = 0.001 with upwinding Recall that according to Remark 4-5, in the case of small viscosity, the time step size can be increased to some extent without adverse stability effects. To demonstrate this, we take k = h = 0.1, and p = 0.001. The numerical results in Table 4-3 support our claim. e iteration error at —5 • t = k t = 1.0 t = 2.0 t = 3.0 t = 4.0 t = 5.0 Ie-3 1 eu ep 2.18e-2 1.83e-l 8.61e-2 8.83e-2 9.43e-2 6.26e-2 9.70e-2 5.42e-2 9.86e-2 5.13e-2 9.99e-2 5.03e-2 Table 4.3: S R M errors for p = 0.001 with a pretty large time step k = h = 0.1 Although we use explicit schemes for SRM (4-2) with a\ > 0, we still have to solve a banded symmetric positive definite system. An alternative is to take cti = 0 to avoid solving any algebraic systems. Table 4-4 shows the computational results of the difference scheme (4-65) with cxx = 0 and a2 — 1. We take viscosity p = 0.1, h = 0.1 and k = 0.0005. Good results are obtained except for the pressure near t = 0 (cf. Remark 4-4) • Chapter 4. SRM for the Nonstationary Incompressible Navier-Stokes Equations 111 e iteration error at -> t = k t = 1.0 t = 2.0 t = 3.0 t = 4.0 t = 5.0 5e-l 2 eu 5.b4e-3 3.57e-2 2.94e-2 2.7Ie-2 2.62e-2 2.60e-2 ep 2.92e-0 9.70e-2 7.03e-2 6.17e-2 5.87e-2 5.77e-2 Table 4.4: S R M errors for p 0.1 with = 0 Chapter 5 S R M for the Simulat ion of Misc ib l e Displacement in Porous M e d i a 5.1 In t roduct ion As explained in Chapter 1, miscible displacement occurs in the tertiary oil-recovery process which can enhance hydrocarbon recovery in the petroleum reservoir. Numer-ical simulation plays an important role for this process because solvents (or chem-icals) are expensive and experiments are hardly possible. Mathematically, misci-ble displacement in porous media is modeled by a nonlinear coupled system of the pressure-velocity equation and the concentration equation with appropriate boundary and ini t ia l conditions. The pressure-velocity equation is elliptic, while the concen-tration equation is parabolic, but normally convection-dominated. Accuracy for ve-locity approximation is important to obtain a good approximation for concentration since the concentration equation only includes the velocity variable. M i x e d finite el-ement methods for the pressure-velocity equation have been applied for this purpose [40, 41, 42, 46, 47, 120]. We are interested in applying the idea of the S R M to this equation since it has benefits over mixed finite element methods. The S R M formulation for the pressure-velocity equation is a direct application of the sequential regularization method for time-dependent problems (see previous chapters). The velocity variable is involved in the linear systems only and the pres-sure variable is obtained by substitutions (without solving any linear systems) at each iteration level. We notice that the same formulation can also be obtained from the augmented Lagrangian method (originated from Uzawa's algorithm) [50, 31]. Unl ike 112 Chapter 5. SRM for the Simulation of Miscible Displacement in Porous Media 113 the augmented Lagrangian method using spectral analysis to discuss the convergence rate for discretized problems, we use asymptotic, methods directly for the differential problems. The asymptotic method is easier to use for more general and more com-plicated problems than spectral analysis because the latter is scarcely possible for non-symmetric operators. We wi l l later prove that our iterative schemes can improve the error to 0 (e s ) at the 5 t h iteration level, where e is a small positive number. In other words, the convergence rate of our iterative procedure is about 0(e). The-oretical convergence analysis and numerical experiments show that the number of iterations is extremely small, usually 2. The organization of this chapter is as follows. In §5.2, we describe our S R M iteration for the time-discretized problem (the differential problem in spatial vari-ables) and its advantages. In §5.3, we show its convergence . Then in §5.4 we give a fully-discretized scheme using the Galerkin finite element method for the problem formulation where the S R M is used for the pressure-velocity equation. Final ly , in §5.5, we present numerical examples to demonstrate the effectiveness and accuracy of our method. 5.2 S R M F o r m u l a t i o n To recall, we again write down the model problem. Consider the miscible displacement of one incompressible fluid by another in a porous reservoir (7 C R 2 over a time period [0, T ] . Let p(x,t) and u(x,t) denote the pressure and Darcy velocity of the fluid mixture, and let c be the concentration of the invading fluid. Then the mathematical model is a coupled nonlinear system of partial differential equations u = —a(gradp — 7gradd) , (x,t) 6 Q x [0,T], (5.1a) divu = q(x,t), (x,t)eilx [0,T], (5.1b) dc <f>a-- div(D(u)gradc) + u • grade = g(c), (x,t) <G tt x [0,T], (5.1c) Chapter 5. SRM for the Simulation of Miscible Displacement in Porous Media 114 with the boundary conditions u n = 0, ( x , i ) G T x [0,T], (5.2a) D(u)gradc-n = 0, (x,t) G T x [0,T], (5.2b) and the in i t ia l condition c(x, 0) = c 0 (x ) , x G ft, (5.3) where a = a(x, c) is the mobil i ty of the fluid mixture (we wi l l later denote it as a(c)), 7 = 7(x,c) and <i(x) are the gravity and vertical coordinate (we wi l l later denote 7 as 7(c)), q is the imposed external rates of flow, </>(x) is the porosity of the rock, D is the coefficient of molecular diffusion and mechanical dispersion of one fluid into the other, g = g(x, t, c) is a known linear function of c representing sources, and n is the exterior normal to the boundary F = <9ft. We assume that the mobil i ty is bounded below and above by positive constants and its gradient is bounded above by a positive constant. For existence of p, we assume that the mean value of q is zero and for uniqueness we suppose p has mean value zero. In recent years much attention has been devoted to the numerical simulation of this problem. In this chapter we are interested in solving the velocity-pressure equation (5.1a)-(5.1b) using the idea of the S R M for the time-discretized problem. The method considered herein is a direct application of the sequential regularization method (see previous chapters but without the sense of regularization) and is closely related to the augmented Lagrangian method [50, 31] (without the augmented La -grangian framework). We wi l l analyze the method using the technique presented in previous chapters, in particular Chapter 3. These analyses give the convergence of the iterative procedure and its convergence rate at the same time. 0 < m 0 < a(c) < M c 0, Chapter 5. SRM for the Simulation of Miscible Displacement in Porous Media 115 After a time discretization, we obtain the following system for u and p at the current t ime step: u = -a (c ) (g radp - 7(c)gradd), (x, t) € ft x [0, T], (5.4a) divu = <?(x, t), (x, ijeftx [0, T], (5.4b) where c is an approximation of c assumed to be known. Taking e to be a small positive number, we replace the system (5.1a), (5.1b), (5.2a) by the following iterative method: for s = 1, 2, • • •, find {us,ps} such that a(c)~1us grad(divus — q) = -(gradp8_x - 7(c)grado0, (x, t) G ft x [0, T ] , (5.5a) u s n = 0, ( x , i ) € r x [0 ,T], (5.5b) and ps = - ^(divua - q), (x,t) G ft x [0 ,T] , (5.6) where the in i t ia l guess po is required to satisfy the zero mean value property: / podx — 0. Thus;p s has mean value zero from (5.6) and (5.5b). We note that, by taking po = 0, \ each iteration is a k ind of penalty method. ' ti This iterative procedure has the following salient features: 1. We solve a small system (5.5) for the velocity u , and obtain the pressure p from (5.6) directly. We wi l l show that the accuracy of such a method is 0(es) at the 5th iteration level. Note that the system (5.5) is well-posed since, unlike the usual penalty method, we need not take e very small. Chapter 5. SRM for the Simulation of Miscible Displacement in Porous Media 116 2. The velocity-pressure equation was recently solved by the mixed finite element method [40, 41, 42, 46, 47, 120], in which the discrete spaces for u and p need to satisfy the Babuska-Brezzi condition, and the resulting linear system has a nonpositive definite coefficient matrix. In our method, u and p are obtained from equations (5.5) and (5.6) separately, and compatibili ty conditions between the discrete spaces of u and p are not needed. Moreover, system (5.5) leads to a symmetric positive definite coefficient matrix. 3. When the standard finite element method [38, 39, 44, 45, 48, 98, 99] is applied for the pressure equation, the velocity needs to be obtained by finite differencing the pressure variable, which gives less accuracy. Note that the accuracy of the approximate velocity is important, since the concentration equation involves the velocity only. The velocity in our method is obtained directly, without finite differencing. 4. The discrete version of our S R M formulation (5.5)-(5.6) gives the same accu-racy for the velocity as the mixed method and requires the solution of well-conditioned linear systems like Galerkin methods. Note that our numerical experiments wi l l show that a few (much less than. 10) iterations are usually enough for our iterative procedure. 5.3 Convergence Analys is Before we begin our analysis, we first describe some notation to be used throughout the rest of the paper. As in Chapter 4, we use Lp(fl), or simply Lp, to denote the space of functions whose pth power is integrable in f i , wi th the norm Chapter 5. SRM for the Simulation of Miscible Displacement in Porous Media 117 where u = • • •, u„ ) . We wi l l omit the subscript p in the norm notation when p = 2. H m is the closure of C£° (0 ) in the norm H H - = £ l l ^ l l 2 • \0<|o;|<m / In addition, we define the divergence space H(div) = { w £ L2(fl)2 : divw £ L2(fl)} with the following norm: | u | | H ( ^ ) - ( | | u | | 2 + | | ^ u | | 2 ) " . \g\\i^([a,p];B) = ( / \\g{-,t)\\2Bdt)? and | | f l f | | L o o ( T C T ^ ] ; B ) = sup \\g(-,t)\\B. We shall denote by (-, •) and (•, •) the inner products in Cl and on V, respectively. For a normed linear space B wi th norm || • \\B and a sufficiently regular function g : [a, (3] —> B, we define \\9(;t)\\B<lt)* a n d \\g\\L~(\a,l3];B) a<t<(3 If [a, (3] = [0, T] , we simplify the notation as ||<7||i,2(fl) and | |S ' | | L O O ( B ) J respectively. We shall also denote generic constants by M and K , which may be different at different occurrences. Before stating our convergence theorem, we first give a lemma. L e m m a 5.1 There exists a unique solution {u,p} to the problem u = -a(c)(gradp - f), (x, t) £ ft x [0, T ] , (5.7a) divu = q, (x,t) £ tt x [0,T], (5.7b) u - n = 0, (x,t) £ T x [0,T], (5.7c) where p and q have mean value zero. Furthermore, there exits a constant Mi such that the following estimates hold: | |u | | H (^) + \\p\\m < Mi[\\q\\ + ||f ||H(«K„)], Vt £ [0,T]. (5.8) Chapter 5. SRM for the Simulation of Miscible Displacement in Porous Media 118 Proof: Substituting (5.7a) into (5.7b) and (5.7c) we see that p satisfies a Poisson equation wi th Neumann boundary condition — div(agradp) = q — div(af), in ft, dP f P a— = at • n on 1. an Noting the zero mean value of p and using the standard results for the Poisson equa-tion, we obtain the uniqueness and existence of p and the estimate \\p\\m < M [\\q\\ + ||f ||H(«K„) + | |f ' n | | H - i / 2 ( r ) ] • A n application of the trace inequality [54] (p.28) ||f • n | | H - i / 2 ( R ) < ||f\\H(div), Vf G H(div) leads to the inequality (5.8) for p. Then the existence, uniqueness and the estimate (5.8) for u follow directly. • We are now ready to describe our convergence theorem and prove it using a similar technique as Chapter 3. T h e o r e m 5.1 Let {u ,p} be the solution of system (5.4a), (5.4b), and (5.2a) and {xisiPs} the solution of (5.5)-(5.6). Then we have ||u - u a | | H ( « H + \\p-Ps\W < ( 1 M M i £ ) 1 b o •5 = 1,2,---. (5.9) Here M i is from Lemma 5.1 and we assume M\e < | . Proof: We first consider the case s = 1 of (5.5)-(5.6). Write (5.5a) and (5.6) as a ( c) _ 1Ui + gradpx = ^(cjgradd div^ - q = e(p0 - px). Chapter 5. SRM for the Simulation of Miscible Displacement in Porous Media 119 Then subtracting equations (5.4a) and (5.4b), we have a ( c ) _ 1 ( u i - u) + grad ( p i - p) = 0 (5.10a) div(ui — u) = e(p0 p i ) . (5.10b) Using Lemma 5.1, we obtain ||ui - u\\H{div) + \\p! -p\\Hi < Mit\\po -Pi\\. (5.11) Wri t ing p0 — pi = po — p + p — pi, we immediately have | |u i - U | | H ( A « ) + Ibi -P\\H> < l e\\Po ~P\\- (5-12) Now we look at the second iteration s = 2. A t first, we can get equations similar to (5.10a) and (5.10b): a ( c ) _ 1 ( u 2 - u) + grad ( p 2 - p) = 0 (5.13a) div(u2 — u) = e(pi - p2). (5.13b) App ly ing Lemma 5.1 again, we have ||u2 - u||H(rf,-„) + ||P2 - p\\m < MitWpi - pad <M1e(\\p1-p\\ + \\p2-p\\). (5.14) Noting M\t < 1 and using the estimate of | |pi — p|| (see (5.12)), yield ||u2 - u | | H ( A » ) + ||P2 - p\\m < ( 1 MMie)2\\Po ~ P\\- (5-15) We can repeat this procedure, and by induction, conclude the results of the theorem. • F rom Theorem 5.1 we see that the convergence rate of our iterative scheme (5.5)-(5.6) is about 0(e). This implies that the number of iterations needed to achieve a prescribed accuracy is very small. The fast convergence of our method makes it dramatically different from penalty-like methods. Chapter 5. SRM for the Simulation of Miscible Displacement in Porous Media 120 5.4 T h e G a l e r k i n A p p r o x i m a t i o n a n d Its E r r o r E s t i m a t e s In this section, we approximate the velocity-pressure iterative scheme (5.5) and the concentration equation (5.1c) by using the standard Galerkin method. We only pro-vide a brief description about this approximation. More details are in [82]. Let W = {w £ H(div) : w • n = 0 on T}. The variational form of (5.5) can be written into the following: find u s : [0, T] —> W such that ( a _ 1 ( c ) u s , w) H—(divus, divw) = (ps-i H—q, divw) + (7(c)grad<i, w), Vw £ W . (5.16) The weak form of the concentration equation (5.1c) reads: find c : [0,T] —> H 1 ( f i ) such that (</»^,^) + ( J D(u)grad C ,g radz) + (u- grade, z) = (g(c),z), V z e ff^ft), (5.17) (c(x, 0), z) = (co(x),z), Vz £ H\n). (5.18) For hu > 0 and an integer k > 0, with respect to the velocity-pressure equation, we introduce finite element spaces W / i C W and Yh = {y : y = divw for w £ W^} associated wi th a quasi-regular subdivision of ft into triangles or rectangles of diameter less than hu. Similarly, we denote by Zh C i / 1 ( f t ) the finite-dimensional space for the concentration equation with grid size hc and approximation index /. Assume that the following approximation properties hold: inf | |w-w h \ \ < Khku+1\\w\\Hk+i, w e W , (5.19) inf \\div(w - wh)\\ < i^^*+ 1(||w||H*+i + \\divw\\Hk+i), w € W , (5.20) inf H* - zh\\ < Khlc+1\\z\\H,+l, z e H\n), (5.21) zh£Zh where K is a constant. The space W^, can be taken to be the vector part of the Raviart-Thomas [100] space of index fc, or Brezzi-Douglas-Marini [30] space of index k + 1. Chapter 5. SRM for the Simulation of Miscible Displacement in Porous Media 121 Given a parti t ion of [0,T], 0 = t0 < ti < ••• < £jv = T , we denote Jn = [tn,tn+1],Atn = tn+1 - £ n , and At = m a x { A t n } . Let {pn,un,cn} and { P n , U n , C n } be {p, u, c} and its approximation at time level tn. We define our approximation scheme at t ime tn(n = 0,1,2, • • •) by the following. Step 1: Given Cn, find { U „ , P n } £ ~Wh x ^ a s follows. Take the ini t ia l guess P0 = 0. For s — 1, 2, • • • , iteratively obtain t j s £ ~Wh and P , £ such that ( a _ 1 ( C n ) t j s , w) + - ( d i u t j s , divw) = ( P s _ ! + i g , divw) + (j(Cn)gradd, w ) , V w £ W f c , (5.22a) Ps = Ps-x - -e(div\Js - q). (5.22b) Let U r a = U s and Pn = Ps for some integer s. Step 2: When U n is known, find C„+i £ Zh such that C — C {(f> n+l -,z) + ( D ( U n ) g r a d C n + i , g r a d z ) + ( U n • g r a d C n + i , z) = (g(Cn+1),z), \/zeZh. (5.23) Note that in step 1 of the scheme, the ini t ia l guess could be more efficiently taken as Po = Pn-i for n > 1 and P0 = 0 for n = 0. Our numerical experiments w i l l show that the number of iterations can be generally taken to be s — 2 for the range of perturbation parameter e = 1 0 - 3 to 1 0 - 5 . The following error estimates of the above approximation are proved in [82]. Theorem 5.2 Let {p, u,c} be the solution to Problem (5.1c)-(5.1b), and { P , U , C } the solution to the scheme (5.22)-(5.23) with s iterations at each time step. Then there exists a constant K, for At sufficiently small, such that the following error estimates hold at time step tm ( m = 0,1,2, • • •, N): || Pm — Pm || + \\Um ~~ Um\\H(div) + \\Cm ~ Cm\\ Chapter 5. SRM for the Simulation of Miscible Displacement in Porous Media 122 < Ks[ts + / i ' + 1 | | c | | L o o ( [ 0 ) t m ] ; ^ + 1 ) + h l c + 1 \ \ c t \ \ L 2 m m ] . H i + i ) + (5.24) + ^ + 1 ( l l U I U ~ ( [ 0 , t m ] ; H * + i ) + | | d i u u | | L o o ( [ 0 i t m ] . H * + i ) ) + A t ( | | u i | | I , 2 ( [ 0 ] t m ] . L 2 ) + | | C t t | | L 2 ( [ o , t m ] i L 3 ) ) ] , m —1 ( A i n | | c „ + i — c n +i i i i ) 1 ^ 2 n=0 < Ks[es + hlc\\c\\Loo^tm].Hi+i) + hlc+1\\ct\\L2([0jm].Hi+i) + (5.25) + ^ + 1 ( I M U ~ ( [ 0 , t m ] ; H * + i ) + ll<fr'wu | | Loo( [ 0,t r a];tf*+3)) + A < ( | | U t | | L 2 ( [ 0 , t m ] ; L a ) + | | C „ \ \ L 2 ( [ 0 , t m ] ; Z 2 ) ) ] • This theorem tells us that for sufficiently small perturbation parameter e, the error estimates for the velocity, pressure and concentration are optimal. 5.5 N u m e r i c a l Exper iments In this section, we present some numerical examples to show how well our iterative scheme performs, and how the parameter e affects the number of iterations needed and accuracy required. For simplicity, we wi l l just consider the pressure-velocity equation, since the concentration equation has been analyzed previously [40, 41, 42, 45, 46, 47, 48, 98, 99, 120]. Consider the elliptic problem with Neumann boundary condition u = —a(gradp — f), x G ft, divu = g(x), x G ft, u • n = 0, x G T, where ft is a square and V its boundary. More general domains ft wi l l not present technical problems. The approximation scheme takes the form: F i n d I P G W^,, for s = 1, 2, • • • , ( a ^ U 3 , w) + \{divU% divw) = {P3'1 + \q, divw) + (f, w), Vw G W f c , (5.26) ps = pa-i _ k(divW -q). (5.27) Chapter 5. SRM for the Simulation of Miscible Displacement in Porous Media 123 Figure 5.1: One element with velocity on each edge and pressure at the center Par t i t ion the domain ft into a set of squares of side length h. We take the space Wh to be the vector part of the Raviart-Thomas [100] space of index 0. Thus W f c = [Vx ®V0) x (V0®Vi), where Vk is the set of one variable polynomials of order less than or equal to k. Consequently, the approximate pressure Ps lies in the space of piecewise constants. Part i t ioning the domain into triangles or rectangles or applying higher order approx-imation polynomials can be treated analogously. Let [/* denote the constant value of the flux in the positive x or y-direction on the edge a , a = L , R , B , T (representing left, right, bottom, and top, respectively), of each element. See Figure 5.1. Consider w to be the basis function (1 — ,x,0) (on the standard reference square). Apply ing the trapezoidal rule to (5.26) we have Chapter 5. SRM for the Simulation of Miscible Displacement in Porous Media 124 a ^ h 2 - ^ eh h pS-x qL + qR + qr + QB 4e (5.28) where qa is the value of q at the middle point of edge a. Similarly, letting w = (x, 0), (0, 1 — y), (0, y), and simplifying we have the following linear system for each element. eh2ailUsL-2{UsR-Ut + UT-UsB) = -2th thWUsR + WSR-UsL + U}-UsB) ps-i , QL + QR + qr + qs At = 2th pS-i qR + qr + qs At + eh2fL, + eh2fR, th2aBlUB-2{UR-UsL + U^-UB) -2th pS-i + QL + qR + qr + qs At + th2fB, th2aTlU^ + 2(USR -Ui + U}- USB) ps-i + 9L + qR + qr + qs = 2th At + ch2fT. (5.29c (5.29b) (5.29c) (5.29d) Note that equation (5.27) has the discrete version on each element: Ps = P s-1 1 e L UR-USL + U$- USB qL + qR + qr + qB (5.30) h 4 J y J From equations (5.29) we can easily form the element stiffness matr ix. Then assembling al l the element matrices and taking into account the boundary condition we obtain the stiffness matrix. The force vector can be obtained in an analogous way. A l l velocity and pressure errors are measured for iterates { U S , . P S } against the exact solution under the L°° norm, i.e. 11US - Ulloo ^ | | P S - p | | o o u \p\< Chapter 5. SRM for the Simulation of Miscible Displacement in Porous Media 125 The in i t ia l guesses are always chosen to be zero, so all . errors are 1.00 before the iterative procedure starts. E x a m p l e 5.1 Let the velocity u and the pressure p satisfy (x = (x; y)T) divu = 7T [cos(7ra;) + cos(7ry)], x £ ft, u gradp — ^ u • n = o, x e r, , X £ ft, where ft = [0,1] x [0,1], and T = 9ft. The true solutions for the velocity u and the pressure p are given by u = sin(7ra;) \ sin(Try) J ' p — —(cos(7rx) + cos(7ry) + x — y). 7T The pressure p and external flow rate q = divu are chosen in such a way that they both have mean value zero. E x a m p l e 5.2 Let the velocity u and the pressure p satisfy the nonhomogeneous prob-lem divu = ab3(eby - ebx), x £ ft, u = —a gradp — —x , X £ ft, u - n = #, x £ T , where ft = [0.5,1] x [0.5,1], F = <9ft, a = 0.05, and b — 10. The function g is chosen such that the true solutions for the velocity u and the pressure p are given by u = —a b2ebx + x -b2eby - v b (ehx - eby) Chapter 5. SRM for the Simulation of Miscible Displacement in Porous Media 126 e = l f r 1 e = lO"* e = 1 0 - J iteration velocity pressure velocity pressure velocity pressure 0 l.UU l.UU l.UU l.UU l.UU l.UU 1 1.14E-2 9.71E-3 1.53E-3 8.67E-4 2.69E-4 4.08E-5 2 2.69E-4 4.00E-5 1.29E-4 1.25E-4 1.28E-4 1.26E-4 3 1.29E-4 1.25E-4 1.28E-4 1.26E-4 - -4 1.28E-4 1.26E-4 - - - -Table 5.1: Numerical results for Example 5.1 with grid size = ^ e = 1(T 4 e = 10-" e = 10~ b iteration velocity pressure velocity pressure velocity pressure 0 l.UU l.UU l.UU l.UU l.UU l.UU 1 1.42E-4 1.16E-4 1.29E-4 1.25E-4 1.29E-4 1.26E-4 2 1.28E-4 1.26E-4 1.28E-4 1.26E-4 - -Table 5.2: Numerical results for Example 5.1 with grid size = ^ For Example 5.1, the results with (uniform) grid size are shown in Tables 5.1 and 5.2. The results of Example 5.2 with (uniform) grid size ^ and ^ are shown in Tables 5.3. More examples and results can be found in [82]. From Tables 5.1 through 5.3 we conclude that our iterative method performs as well as the theory predicts. In particular, it can achieve the same accuracy as mixed finite element methods (at least for velocity), while the linear systems to be solved are symmetric and positive definite. Also, the computational work of our method is much smaller than that of mixed methods, since the number of iterations required is usually very small . e = 10~ 3 grid size = ^ grid size = ± iteration velocity pressure velocity pressure U l.UU l.UU l.UU 1.00 1 1.18E-3 5.85E-3 3.23E-4 1.61E-3 2 1.02E-3 5.42E-3 2.00E-4 1.24E-3 Table 5.3: Numerical results for Example 5.2 Chapter 6 N u m e r i c a l Methods of Some Singular Per turba t ion Prob lems In this chapter, we discuss the numerical solutions of some singular perturbation problems which are all special cases of the regularizations (1.12) and (1.13), and have practical meanings in themselves as we indicated in §1.4. In §§6.1 and 6.3 we wi l l construct and analyze uniformly convergent methods for these singular perturbation problems. So we first describe the definition of the method ( cf. [90]). Def in i t ion 6.1 Ifu is the solution of a singular perturbation problem with a parame-ter e and uh is an approximation obtained using a uniformly convergent method, then there exist two constants eo and hQ independent of e and h such that when 0 < e < eo and 0 < h < ho we have an inequality of the form \\u-uh\\ < Chp, (6.1) or in a weaker version (cf. [79]) \\u-uh\\ <C(hp + er), (6.2) where C > 0, p > 0 and r > 0 are independent of e and of the mesh width h, and || • || is some appropriate norm. • We use M to represent generic (in the sense of 0(1)) positive constants indepen-dent of e and of the mesh width h. Some of these constants w i l l also be denoted by mo, m i , m2, M0 and M i , etc. 127 Chapter 6. Numerical Methods of Some Singular Perturbation Problems 128 6.1 One Dimens ional Quasilinear Turn ing Point Problems In this section, we consider the following two-point boundary value problem wi th Dirichlet data at the endpoint: — eu" — b(x,u)u' + c(x,u) = 0, i £ / = [ - l , l ] , (6.3a) « ( - l ) = E/_, u(l) = [/+, (6.3b) where e 1 usually and b(x,u) may be zero at some isolated points. Such prob-lems are usually called turning point problems. Constructing uniformly convergent methods for problem (6.3) is generally very hard. We wi l l consider a simpler case in which we know the point, say x*, at which b(x, u(x)) = 0 and bx(x, u) ^ 0 for al l u in the vicini ty of the solution. There are two types of turning point problems. They are called repulsive and attractive turning point problems corresponding to the negative and positive sign of j^b(x,u(x))\x=x*, respectively. The rest of the section is devoted to the two types of turning point problems and is actually a summary of two papers [79, 115] written by the author and his collaborator. We assume that the problem which we consider (in the form of (6.3)) has a unique solution. We also assume that the coefficients of problem (6.3) are sufficiently smooth (usually C2(I x R ) is enough). 6.1.1 A repulsive tu rn ing point problem In this section, we assume cu(x,u) > c0 > 0 on [ - l , l ] x R (6.4) Hence, a max imum principle holds for (6.3). B y constructing a barrier function it is not difficult to get max lul < r, (6.5) -1<X<1 v ' Chapter 6. Numerical Methods of Some Singular Perturbation Problems 129 where r = max_ 1 < a ; < 1 \c(x,0)\/cQ + max(\U-\, \U+\). Furthermore we make the fol-lowing assumptions 6(0, u) = 0 A ( 0 , u) < 0 . for \u\ < r, (6.6a) b(x,u)^0 for |u | < r and x / 0. (6.6b) According the condition (6.6) there exists a small positive number 8 such that bx(x,u) < —bo < 0 for \x\ < 8 (6.7a) b(x,u)>b_1>0 for - 1 < x < -8 (6.7b) b(x,u) < -bx < 0 f o r < J < x < l , (6.7c) where 8, bo, 6-i and bi are positive constants independent of e. We first give the bounds on the derivatives of the solution which is useful in the proof of uniform convergence: \u^(x)\ < M ( l + e - i e x p ( - m 0 ( x + l ) /e) for a; G [-1,<J], (6.8a) |«^| < M for \x\ < 8, (6.8b) \u^(x)\ < M ( l + e - i e x p ( - m 0 ( l - x)/e) for x G [5,1], (6.8c) where 8 is sufficiently small. (6.8a) and (6.8c) are a direct application of the result for problems without turning points which is derived by [113] and [77] independently. (6.8b) is proved in [79] by examining the Newton's sequence. This result means that the repulsive turning point problem does not have any interior layers. Next we consider how the solution v(x) of the reduced problem: b(x,v)v'- c(x,v) = 0 , - 1 < x < 1, (6.9a) c(0,u(0)) = 0. (6.9b) approaches the solution u(x). As in [21] (Remark 2.11), applying (6.8b) and results in [83], we have \u(x) - v(x)\ < M ( e + e x p ( - m 0 ( x + l) /e)) for x G [—1, —5], (6.10a) Chapter 6. Numerical Methods of Some Singular Perturbation Problems 130 \u(x) - v(x)\ < Me for x 6 [-6,6], (6.10b) \u(x) - v(x)\ < M ( e + e x p ( - m 0 ( l - x)/e)) for x e [8,1], (6.10c) Now we start to construct an almost uniformly convergent algorithm. In the construction, we combine the initial-value technique [67] (a modification is given in [81]) with the idea in [23]. We want to indicate here that [67, 81, 23] are al l for problems without turning points. We first rewrite (6.3) as eu"+{f(x,u))'- g(x,u) = 0, x € / = [-1,1], (6.11a) u(-l) = £/_, u(l) = U+, (6.11b) where f(x,u) — J0ub(x,s)ds and g(x,u) = c(x, u) — fx(x, u). Integrating (6.11a), we get eu' + f(x,u) = [Xg(t,u(t))dt + K for - 1 < x < 1, (6.12) Jo where the integration constant is K = eu'(0). Let E(x)= fX g{t,u(t))dt. Jo Then problem (6.3) is reduced to the following equivalent nonlinear in i t ia l value prob-lems: eu'l + f{x,ul) = E(x) + K, -1 < x < 0, (6.13a) m(-l) = U-; (6.13b) and eu'2 + f{x, u2) = E(x) + K, 0 < x < 1, (6.14a) u2(l) = U+. (6.14b) Chapter 6. Numerical Methods of Some Singular Perturbation Problems 131 Replacing E(x) by E(x)= f g(t,v{t))dt, Jo where v{t) is the solution of the reduced problem (6.9), and neglecting K in (6.13a) and (6.14a), we obtain the approximate problems: ey[ + f{x,yi) = E{x),-l<x<0, (6.15a) y i ( - i ) = tf-; (6.15b) and ey'2 + f{x,y2) = E(x),0<x<l, (6.16a) y2(l) = U+. (6.16b) It is proved in [79] that \Ui(x) - yt(x)\ < Me for 0 < |x | < 1 and * = 1,2. (6.17) F ina l ly we consider the numerical solution of (6.15) and (6.16). Note that by using the change of variable x = —x the numerical method for problem (6.15) wi l l follow from that for problem (6.16), so we proceed considering only problem (6.16). For convenience we write (6.16) as ey2 — m\xy2 = —xe(x,y2), 0 < x < 1, (6.18a) y 2 ( l ) = U+, ' . • (6.18b) where e(x,y) = d(x,y) + mxy. In the expression d(x,y) = (f(x,y) + E(x))/x is bounded and mj is a positive constant to be determined below. O n the interval [0,1], introduce an arbitrary mesh {xi,i = 1,2, • • • , 7Y, wi th XN = 1}. Integrating problem (6.18) on the subinterval [ 2 ; , x ; + i ] , and replacing function Chapter 6. Numerical Methods of Some Singular Perturbation Problems 132 e(x,y2(x)) by e (x ; + 1 , y 2 (^ t+i ) ) > suggests a difference scheme: S/2,.- = kiV2,i+l + C 1 ~ ki)e(xi+i>y!*,i+i)/mii (6.19a) y 2 V = <7+, i = l , 2 , . - - - , i V - 1, . _ (6.19b) where k{ = e x p ( — | m 1 e _ 1 ( x 2 + 1 — x2)). According to the idea of [23], we can choose a suitable mesh _ J 1 + (e/m) l n ( l - A T - 1 ) ^ ( 1 - e ) ) , i = N - Nu • • •, A , x N - N l = 1 - ht [ Xi,maxi(xi - X J _ I ) = h2, i = 1, • • • , N - Nx - 1, (6.20) where ht = e\ l n e | / m and hi = 1/iVj. We have the following error estimate (see [79]): Theorem 6.1 Let y2(xi) and y2i be the solutions of problems (6.18) and (6.19), respectively, i = 1, • • • , N. Under the mesh (6.20), taking ™-i > -fu(x,u) = -b(x,u), then we have max \u2(xi) - UoA < Mh, (6.211 Ki<N 1 V ' A" ~ V ' where h = max(hi,h2). Apply ing the same method, we can also get a numerical solution yxi,i — — A , • • •, — 1, of problem (6.15) such that _ max_ i \yi(Xi) - yhu\ < Mh. (6.22) Let f for l = -N, y\ = I u(0) for i = 0 (6.23) I y$ti for i = l,---,N be an approximation of the solution u(x) of problem (6.3), where v(0) can be solved from c(0,u(0)) = 0. Apply ing (6.17), (6.10b) , (6.21) and (6.22), we have _max N \u(Xi) - yf | < M(h + e). (6.24) Chapter 6. Numerical Methods of Some Singular Perturbation Problems 133 Therefore, the numerical method we propose is an almost uniformly convergent method. In [79], a technique is given to modify the method to achieve uniform convergence. A numerical example is also given in [79]. 6.1.2 A n attractive turn ing point problem In this section, we shall extend the results from [114] and [78]. We make the following assumptions: Moreover, we shall assume that e is sufficiently small. Note that, unlike the previous section, the max imum principle may not hold since we do not assume (6.4) here. The corresponding reduced problem has a discontinuous solution consisting of two smooth curves, u + and u _ , which satisfy: uniformly convergent methods, we shall estimate the solution u(x) and its derivatives and the quantities: v± = (u — u±)^k\x), for k = 0,1,2 and x £ I±, bi(x,u±)u'± — ci(x,u±) = 0, u ± ( ± l ) = U±. denote the maximum norm in C(I). To consider where /_ = [ - l , 0 ] , J + = [0,1]. We collect these estimates in the following theorem. Theorem 6.2 \(xv±(x))'\ < M(pV{x)\ x G Chapter 6. Numerical Methods of Some Singular Perturbation Problems 134 where \{xv±{x))"\ < Mip+.p^Vix)), x € / ± , e|u"(x)| < M ( e + y ( x ) ) , x e I, e\u"'(x)\ < M(p + iTlV(x)), x e l , V(x) = e x p ( - — | x | ) , t1 with an positive constant mo independent of e. • ' From the estimates we see that the attractive turning point problem has an interior layer but no boundary layers. This theorem is built on several lemmas whose proofs are very technical. To show how technical these proofs are we prove one lemma below. We refer to [115] for complete details. L e m m a 6.1 |w(x)| < M. Proof: Let J |x | , x e i\[-p,fi\. " l g + f, *e[-p,p}. We can easily verify that p G Cl{I), max( |x | , p) > p{x) > - m a x ( | x | , p). (6.25) Next, consider the following Riccat i ini t ia l value problem P(a) := col + xKa + M0p{x) + ec? = 0, (6.26) a(0) = 0. (6.27) It has a uniformly bounded solution. Indeed, by applying Newton's method to (6.26), (6.27) (cf.[88] or [79]), with the ini t ia l guess: oo(x) = - r ^ p ( t ) e x p ( - ^ ( x 2 - f))dt, Jo e Ze Chapter 6. Numerical Methods of Some Singular Perturbation Problems 135 the conditions of the Newton-Kantorovich theorem are satisfied since: | |ao(x) | |oo < M , ll^'M" 1!! <M^ | |Ol - « o | | o o < Mp, \\P"(a)\\ < 2. Here, || • || is the operator norm corresponding to || - | | o o - Hence, there exists a solution a(x) to (6.26), (6.27), such that ||a - ao | |oo < Mu, thus a is bounded uniformly in e. Furthermore, we have xct(x) < 0 for x e I. (6.28) Indeed, because of the max imum principle and ea' + xKa = -Mop(x) - to? < 0, a(0) = 0, a(x) < 0 for x > 0 and a(x) > 0 for x < 0. Let <p(x) — e x p [ / a(t) dt]. Jo It holds that 0 < m 0 < ip(x) < M , ip'(x) = a(x)ip(x) = 0(1). (6.29) Note that <p is a solution to the following equation: e<p" + xKip' + M0p(x)tp = 0. Let us now consider an auxiliary problem: eu" + xbi(x, u(x))u' — c(x, u) = 0, (6.30) u ( - l ) = U-, u(l) = U+, (6.31) Chapter 6. Numerical Methods of Some Singular Perturbation Problems 136 and let us make the transformation u(x) = z(x)ip(x). Then z(x) satisfies : ez" + [2ect(x) + xb^x, u(x))]z' - c(x, z) = 0, (6.32) *<-» = £ry = %rv (6'33) where c(x,z) = —[e(p"(x)-\-xbi(x,u(x))(p'(x)](f(x)~xz-\-(p(x)~1c(x,z(p(x)) — M0p(x)z — x[bi(x, u(x)) — b*]a(x)z + ip(x)~1c(x, zip(x)). Choosing M 0 sufficiently large and using (6.25) and (6.28) we have cz(x,z) = M0p(x) — xa(x)(b(x, u(x)) — &*) + cu(x, z<p(x)) > m 2 max( |x | , p). Hence, the problem (6.32), (6.33) satisfies the maximum principle and therefore has a unique solution z, such that w , ) l < w _ i ) l + W i ) l + ^ - Q ^ L < M . This shows that (6.30), (6.31) has a unique solution which is equal to ut. Then (6.29) completes the proof. • The numerical method is closely related to that from [114]. The same special non-equidistant mesh Ih is used. It has the following mesh points: 2i Xi = X(ti), ti = - 1 + — , i = 0 ( l )n , n = 2no, no G N , where X(t) = [u{t):=^v *G[0 ,ao] 7r(i) := S(t - a 0 ) 3 + \u"(a0){t - a0)2 +u>'(a0)(t - a 0 ) + u(a0), t G [a0,1] I -X(-t), * e [ - i , 0 ] Chapter 6. Numerical Methods of Some Singular Perturbation Problems 137 a 0 is an arbitrary parameter from (0,1), i 5 is determined from 7r(l) = 1, S O that A G C2(I±) and A G Cl(I), and the parameter f3 is chosen from ( 0 , 7 _ 1 ( 1 — a 0 )~ 2 ] . It may look as i f A is artificial, but its part a; is a suitable rational approximation to the logarithmic function representing the inverse of the interior layer function V(x) for x > 0. Then TT is just a smooth extension of u. Let hi = Xi - x^, i = l ( l ) n , 1 hi = -(hi + hi+1), and let wh denote a mesh function on Ih \ { — 1,1}, which wi l l be identified wi th the R n _ 1 - v e c t o r : wh = [u>i,iu2, • • • , ^ n - i ] T , (wi := wf). Moreover, let us introduce the following standard finite-difference operators: D'±Wi = ±(wi±1 - Wi)/hi, D"wi = - Wi)/hi + (wi+i - Wi)/hi+1]/hi. We shall use the following discrete L 1 - n o r m : 71-1 i = l For al l this cf. [114]. Final ly, the constants M w i l l now be independent of Ih as well . Before discretizing the problem (6.3), as in the previous section, we rewrite (6.3a) in the following conservation form: Tu:=-eu"-'f{x,u)f + g(x,u) = 0, x G / , (6.34) Chapter 6. Numerical Methods of Some Singular Perturbation Problems 138 where f+(x,u), x e 1+ g+{x,u), x e /+ pu f±(x,u) = / xb1(x,s)ds, J u± (x) g±(x,u) = c(x, u) — xbi(x, u±(x))u'±(x) + / (xbi(x,s))xds J u± (x) = c(x,u) — xci(x,u±(x))-{- I (xbi(x,s))xds. Ju4-!x) Then the discrete problem corresponding to (6.34), (6.3b) is given by: Thwi = 0, i = l ( l ) n - 1, (6.35) where rph \ Tllwu i = l ( l ) n 0 [ Tfttfi, i = n0 + l{l)n - I T±Wi = -eD"wi - D±f±(xi,Wi) + g±(xi,Wi), and where w0 and wn should be replaced by U- and U+ respectively. The discretiza-tion is a generalization of that for the mi ld ly nonlinear case considered in [114] and [78]. Let us introduce the following assumption in addition to i) — iv): v) gu(x,u) = (xbi(x,u))x -f cu(x,u) > g* > 0, x £ /, u G R . The following error estimate is proved in [115]. Theorem 6.3 The discrete problem (6.35) has a unique solution wh and the following estimate holds: \\wh - u f c | | * < M-[u + exp(-n)], n where uh = [u(xi), u(x2),.. • , u ( x „ _ i ) ] T . • A numerical example is given in [115]. Chapter 6. Numerical Methods of Some Singular Perturbation Problems 139 6.2 Notes about Spurious Solutions of U p w i n d Schemes 6.2.1 Inadequacy of Yavneh's argument First order upwind schemes have been found by several researchers to yield spurious solutions. Such behavior has been attributed to the excessive artificial viscosity intro-duced by the numerical scheme, to multiple solutions of the nonlinear set of algebraic equations that is obtained from the discretization, and to poor resolution by grids that are too coarse (cf. [121, 27]). A simple example in [61] also shows that an upwind scheme may lead to a spurious solution near a boundary layer. Yavneh in his P h . D . thesis [121] (also see [27]) claimed that the spurious solution may occur even in cases where none of the above apply, and that such behavior is seen to occur even in the linear scalar advection-diffusion equation, despite the fact that the truncation-error tends to zero with the mesh-size throughout the entire domain, the solution being smooth everywhere. However, his proof is incomplete. Yavneh considers a linear advection-diffusion equation —eAu + -^rue = 0, (6.36a) r 2 u (a,9) = Ui , u(b,6) = u0 (6.36b) over the circular disk 0 < # < 2 7 r , 0 < a < r < 6 . Here 1 1 Au = urr + -ur + -^uee is the Laplacian in polar coordinates and e > 0. The unique solution (the uniqueness comes from the maximum principle) is _ Marco-'] a In order to see what might go wrong with the numerical solution of this problem, we rewrite the equation in Cartesian coordinates Lu = - e A u —-ux + —r^-—-uy = 0, (6.37a) Chapter 6. Numerical Methods of Some Singular Perturbation Problems 140 u \ x 2 + y 2 = a = m, u\x2+y2=b = u0, (6.37b) where A u = uxx + uyy. As is well known, discretization of this equation by central finite differences loses its stability when e is small compared to the product of the mesh-size and the ab-solute value of either coefficient of the first derivatives in the equation. A common practice is to retain stability by increasing the absolute values of the coefficients of the discretized second derivatives. A particular method of this type is the first order upwind difference scheme. This sort of discretization gives a second-order central difference scheme to the following equation: -c\uxx - e2uyy V „ux + X uy = 0 (6.38) xz + y xl + y1 with the same boundary conditions as (6.37b),where h\y\ h\x\ e i = e + Wt o I—oT> £2 = e + 2(x2 + y2) 2(x2 + y2) We plot the relationship of these continuous and discrete equations as follows: Problem (6.37) <— Upwind scheme (first order) t II Problem (6.38) <— Corresponding central scheme Yavneh considers the case e < h/2r ( e > h/2r is not interesting since in this case the central difference scheme is stable) and shows that the difference between the solutions of problems (6.37) and (6.38) has a lower positive 0(1) bound throughout the domain except near the boundaries as the mesh-size tends to zero. Based on this result, he then claims that the first order upwind scheme yields a spurious solution of problem (6.37). From the above diagram, we see that the argument is incomplete since we do not know if the central difference scheme is a good approximation for problem (6.38). The proof of this fact is not that obvious since we no longer know if Chapter 6. Numerical Methods of Some Singular Perturbation Problems 141 the solution of (6.38) is smooth. (6.38) is a partial turning point problem whose layer property could be complicated. 6.2.2 Our explanation In this section, we show that the stability constant of problem (6.36) depends on e. The smaller the parameter e, the closer the problem is to being ill-posed. Hence, we can expect that any direct discretization (of course, including upwind schemes) on the problem would be unstable when e is sufficiently small . Therefore, it is not strange that the first order upwind scheme fails for this problem. Let's consider (6.36) with a forcing term f(x, y). Suppose that u(x, y) is a solution of the problem, i.e. Lu = f(x,y) (6.39) and u satisfies boundary conditions (6.36b). We construct u*(x, y) = u(x, y) + ( r 2 - a2){b2 - r 2 ) , (6.40) which satisfies the boundary conditions and Lu* = f(x, y) + e(12(x2 + y2) - 4(a 2 + b2)). (6.41) In other words, i f we make a small perturbation e(12(x 2 + y2) — 4(a 2 + 6 2)) to the right-hand side of the equation, the solution changes a lot (by ( r 2 — a2)(b2 — r 2 ) ) . That means the problem (6.39),(6.36b) is not well-posed as e is small. Then we may not expect its approximation (direct discretization) to be stable. Indeed we can prove that any difference scheme wi l l not be stable (as e is small). Suppose we have a discretization (consistent with order hr) LhUh = f(x,y), Bhuh = 0 (discrete B C ) also Lhu*h = f(x, y) + e(12(z 2 + y2) - 4(a 2 + b2)) Chapter 6. Numerical Methods of Some Singular Perturbation Problems 142 Bhu*h = 0 Then Lh(uh - u%) = e(12(x2 + y2) - 4(a 2 4- b2)) Bh(uh - u*h) = .0 If (Lh,,Bh) is stable, or more precisely, i f the coefficient matr ix Ah of the linear algebraic system corresponding to the discretization satisfies Halloo < M , (6.42) where M is a generic positive constant independent of e and h, we get Uh — u = 0(hr) u*h-u* = 0(hr) uh-u*h = 0(e) Hence u-u* = 0(e + hr) This is a contradiction since u — u* = ( r 2 — a2)(b2 — r2). So (Lh,Bh) can not be stable. In fact, we can show C l a i m 1 \KX\\>0( J - j - r ) . (6-43) max(e, hT) Proof: If e > hr and WA^W^ = 0(l/es), S < 1, then U h - U = 0('^)<0(e1-5), u*k-u* < 0(e1~5), uh - u\ = 0(e1'5). Chapter 6. Numerical Methods of Some Singular Perturbation Problems 143 Therefore u — u 0(e l-8\ can be arbitrarily small. This is a contradiction in (6.40). So we must at least have IIA^Hoo = 0(l/e). On the other hand, if e < hr, we can prove that II Vlloo > 0(1). Thus (6.43) is proved. • For the first order upwind scheme, we actually can prove C l a i m 2 A-h1\\ = 0{- (6.44) 'max(e, h)' Before we prove this claim, we first write down the scheme and then prove a discrete max imum principle. The scheme is Uj+lJ-UjJ "t,J-"t-l,J h-xi h-xi — 1 -e : ; e 2(^3:2 "F hxi—\) 2^yi h'yj—1) (forx > 0,y > 0) (forx > 0,y < 0) (forx < 0, y > 0) (fonc < 0, y < 0) x2i+y2j + 1,J "••,.) h-xi + I 1 yJ ,7 ""<J ^ hyj — l y.i U J - U > - 1 J + X, U j - u i , j - l h x i - i x2+y? yj u + Xi u xj+yj h-xi xj+y* yj _ 2 1 . .2 u - 1 , 7 + ~ 2 1 -.2 u 0. (6.45) L e m m a 6.2 TVie difference scheme (6-45) has a discrete maximum principle, that is, LhUh > 0 and •u^lr^ > 0 ==> > 0, where Th represents all the mesh points on the boundaries x2+y2 = a2 and x2+y2 = b2. Proof: Suppose that uh > 0 is not true. Then there exists an interior point (xs,yt) of the domain such that ust < 0. Furthermore we can choose u S j t such that uSft = Chapter 6. Numerical Methods of Some Singular Perturbation Problems 144 min ; j Uij and at least at one adjacent point u 8 j of us,t the rigorous inequality Uij > uStt holds. Hence (for simplicity, we only consider x > 0,y > 0 case), « 5 , t - « g , t us,t—Us,t Us,t—Us,t « j , | - " a , l Lhuh\ < ^ hxs hxa-i ^ hyt hyt-i \{hxs + h yt us,t - uStt xs us>t - usj _ ^ 5 + y* hxa X2S + This is a contradiction to LhUh > 0, which completes the proof. • Now we prove Claim 2. Proof: We construct the barrier function _ x2 _ y2 Wh - \u0\ + Wi \ + Mi - 7 — T V m ^ x \LhUh\, (6.46) max(e,n) « , j where h = max(hx, hy), while hx = max; hx{y hy = maxj hyj. Then it is easy to verify that Lh(wh ± u f c) > 0, to/j ± uh\rh > 0 when M i is sufficiently large. Using the discrete maximum principle we obtain wh±uh> 0. We thus have . . . i b 2 - x 2 - y 2 \uh\ <wh = \u0\ + \ui\ + Mi —p— m a x \ L h u h \ . max{6,n) h3 This means \A~hl\U < M l—- (6.47) max(e: a) So (6.44) follows from (6.47) and Claim 1. • . Hence, the first order upwind scheme for problem (6.39),(6.36b) is not convergent when e is small compared with mesh-size h. Of course, it produces a spurious solution in general. For Yavneh's example (f(x, y) = 0 in (6.39)), numerical calculation verifies the appearance of the spurious solution. From the result of Claim 2 we expect that a Chapter 6. Numerical Methods of Some Singular Perturbation Problems 145 good second-order scheme would converge i f t is smaller than h but much larger than h2. Numerical calculation in [121] verifies this too. Note that, in polar coordinates, the first order upwind scheme converges since its truncation error has a factor e (because the derivatives of u wi th respect to 9 are all zero). In Cartesian coordinates, no such behavior exists. 6.3 A Linear Hyperbo l i c -Hyperbo l i c Singular ly Per tu rbed In i t i a l -Boundary Value P r o b l e m Previous to this work, we discussed several singularly perturbed problems of hyper-bolic type in joint papers [105, 106], constructed some difference schemes according to the properties of the problems, established discrete energy inequalities for the solu-tions of difference problems, and, based on the inequalities, proved that the difference schemes are uniformly convergent in the sense of the discrete energy norm. Bu t in those papers the equations considered did not include a first derivative term wi th respect to the space variable x. Here we discuss a more complete initial-boundary value problem: Lcu = e(utt - uxx) + a(x,t)ut + b(x,t)ux + c(x,t)u = f{x,t), (x,t) G G = {0 < x < /, 0 < t < T} (6.48a) u(s, 0) = <f>(x), ut(x, 0) = ip(x), (0 < x < I) (6.48b) u(0,t) = 0, u(l,t) = 0, (0 < t < T) (6.48c) where a(x,t),b(x,t),c(x,t), f(x,t),(/)(x) and ip(x) are sufficiently smooth functions and a(x,t) > a0 > 0 for all (x,t) G G. Moreover, b(x, t), f(x, t), 4>(x) and if>(x) satisfy the following compatibili ty conditions: C I : (f)(0) = 0, V(0) = 0, (f)(1) = 0, rj>(l) = 0; C 2 : 4>"(0) = 0, 6(0,0)^(0) = / (0 ,0 ) , </>"(0 = 0, 6(/ ,0)^(/) = / ( / , 0 ) . Chapter 6. Numerical Methods of Some Singular Perturbation Problems 146 The subcharacteristics of the reduced operator Low — a(x, t)wt + b(x, i)wx + c{x, t)w (6.49) are timelike (cf.[53]) with respect to the characteristic directions of equation (6.48a), that is, \b(x,t)\/a(x,t) < 1 for (x , t ) G G. (6.50) For simplici ty we assume that b(x,t) >bo>0. Hence, (6.50) becomes b(x,t) < a(x,t) for (x,t)'e G. (6.51) The reduced problem of (6.48) is LOUQ = f(x,t), (6.52a) uo(x,0) = <f>(x), u o (0,f) = 0, (6.52b) where LQ is defined as (6.49). Therefore, the solution of problem (6.48) has boundary layer at t = 0 and x — 1. The reduced problem (6.52) is a first order hyperbolic initial-boundary problem in a semi-bounded region D = {(x,t), x > 0, t > 0}. According to [22], under the conditions Cl and C 2 , the first partial derivatives of the solution Uo(x, t) of (6.52) are continuous, but the second derivatives are discontinuous along the characteristic line. In order to construct the asymptotic solution, it is necessary that uo(x,t) G C2(D). We wi l l give a method to overcome the difficulty without additional compatibil i ty conditions. This section is a summary of the joint paper [107]. In what follows we first give an energy estimate of the solution of problem (6.48). Then, the asymptotic solution is constructed under the compatibility conditions C l and C2 and uniform validity is proved in the sense of the energy norm. Thirdly, an exponentially fitted difference scheme for problem (6.48) is proposed and a discrete energy inequality is established. Chapter 6. Numerical Methods of Some Singular Perturbation Problems 147 Final ly , we show that the discrete problem is uniformly convergent in the sense of discrete energy norm. 6.3.1 Cons t ruc t ion of asymptot ic solution and its remainder estimate We first give the following energy inequality which wi l l be used in estimating the remainder of the asymptotic expansion. L e m m a 6.3 Let u(x,t) be the solution of (6:48), and let the conditions we described previously be satisfied. Then as e is sufficiently small, we have Proof: Mul t ip ly ing equation (6.48a) by 2ea~lut + «, then integrating the obtained equation in the region Gt = {(x,s) |0 < x < I, 0 < s < £}, and performing the standard argument for the energy method, we thus have (6.53). • Next we construct the asymptotic solution. As we indicated before, under the compatibil i ty conditions Cl and C 2 , the solution Uo(x,t) of the reduced problem belongs to C1(D), but not C2(D). In order to realize the iterative procedure of the asymptotic solution, the continuity of the second derivatives of uo(x,t) is needed. However, this usually gives rise to the increase of the compatibil i ty conditions. In this section, we construct the asymptotic solution under the conditions Cl and C2 without adding other ones. B y making some transformations on (6.52) so the ini t ia l and boundary conditions become homogeneous and the right-hand function becomes equal to zero in the neighborhood of t = 0, we obtain a problem in which al l compat-ibi l i ty conditions we need are satisfied. Then based on the transformed problem, an «|| + e\\ut\\ + e\\ux\\ < MK(G, e) (6.53) where Chapter 6. Numerical Methods of Some Singular Perturbation Problems 148 approximate problem for (6.48) is constructed and the solution of its reduced problem belongs to C2(D). Then an asymptotic expansion can be constructed. Now we deal with problem (6.52). Let Uo(x,t) = Wo(x,t) + (f>(x), noting that <f>(0) = 0, then L0wo = F(x,t), w0{x,0) = 0, wo(0,t) = 0, (6.54) where operator LQ is defined as (6.49) and F(x,t) = f(x,t) — Lo4>- Introduce a function u:(y) G C°° satisfying o «><»M) w 1 i (s > 1) and 0 < u(y) < 1. Defining F = u)(t/5)F(x,t), we have n M ) - F ( ^ ) = ( i - . ( > ^ ) = { f M ) | ° ^ / 2 ) Therefore, F(x,t) and F(x,t) have difference only in the region {0 < t < 5, 0 < x < 1}. Replacing F(x,t) by F{x,t), problem (6.54) is changed into L0w0 = F(x,t) = u(t-)(f{x,t)- L0</>), (6.55a) w0(x, 0) = 0, u>o(0, t) = 0. (6.55b) Problem (6.55) for any S > 0 satisfies all compatibili ty conditions we need. Thus wo is sufficiently smooth. Transforming back, i.e. letting UQ = WQ + 4>, we obtain an approximate problem for (6.52) L0uo = / ( M ) . (6.56a) uo(x,0) = <f>(x) , u0(0,<) = 0, (6.56b) where f(x,t) = u(t-)f(x,t) + (l-co(t-))L0cf>. Chapter 6. Numerical Methods of Some Singular Perturbation Problems 149 Here Uo(x, t) is sufficiently smooth since WQ and <f> are. From the appendix of [22], we have u 0 - u 0 = 0(8). (6.57) Now we can get an approximate problem of (6.48) Ltu = f(x,t), (6.58a) «(0 , x) = <p{x), ut(0, x) = ip(x) , u(0, t) = u(l, t) = 0. (6.58b) Using the energy inequality (6.53), we can get \\u - u\\ = 0(8*). (0 <t <T) (6.59) Note that problem (6.49) is the reduced problem of (6.58). Hence, the reduced prob-lem of (6.58) has enough smoothness such that the usual procedure to construct the asymptotic expansion can be performed. For example, we can construct the asymp-totic expansion of (6.58) in the form of u(x,t) = u(x,t) + z, where u(x, t) = « 0 ( x , t) + evj}0)(x, r) + 4°(£, 0 + t), r = t/c, £ = (I - x)/e, u0 is the solution of (6.49) and v^jV^ and v[1^ satisfy the following equations (vi%T + a{x,0)(viO))T = 0, ( ^ O ) ) T ( x , 0 ) + (u o ) t ( .T ,0) = tb(x), l i m vi0) = 0, T—>-00 (v{0\i + b(l,t)(v{0\ = 0, vH\0,t) + uo(l,t) = 0, l i m vH\U) = 0 <->oo Chapter 6. Numerical Methods of Some Singular Perturbation Problems 150 and (v[% + b(l,t)(v{\ = v[l)(0,t) = 0, l i m v[l) = 0, £ - » o o respectively. Using the energy inequality (6.53) and estimating carefully, we have the following bound | |z | | = \\u(x,t) - u(x,t)\\ < M(5^ + + e8~l + e28~2). (6.60) Denote the asymptotic expansion of problem (6.48) as Ul(x,t) = u0(x,t) + evi°\x,r) + v$\{,t) + ev[l){C,t). 2 Apply ing (6.57),(6.59) and (6.60) and taking 8 = es, we finally obtain \\u(x,t)-Ul(x,t)\\< MeK (6.61) R e m a r k 6.1 If the coefficients, right-hand side and initial values of problem (6.4-8) satisfy enough compatibility conditions such that the solution u(x,t) £ C3(G), then we can construct the asymptotic solution without needing the function u(t/5). The asymptotic solution is still Ui (replacing UQ by UQ in the equations for and v$). Moreover, we have the following remainder estimate \u(x,t)-Ul(x,t)\< MeK (6.62) • 6.3.2 Difference scheme and its uniform convergence Taking uniform meshes in the directions of x and t, we obtain a discrete region Gd = {(x{,tj), i — 0, • • • , N, j = 0, • • • , [T/k], Xi = ih, tj — jk], where h and k Chapter 6. Numerical Methods of Some Singular Perturbation Problems 151 are mesh sizes and Nh = I. Denoting ud(x,i) as the approximate value of u(x,t), we establish the difference scheme of problem (6.48) by using the exponential fitting idea: where L{h^ud(x, t) = 7 l ( x , t, k)u$i(x, t) - 7 2 ( x , t, h)udxx(x, t) + a(x, t)uf +b(x, t)ud(x, t) + c(x, t)ud(x, t) = f(x, t), (x, t) £ Gd (6.63a) ud(x, 0) = cf>(x), ud(x, k) - ud(x, 0) = kxb(x), (6.63b) ud(0,t) = ud(l,t) = 0, (6.63c) (X,t) = ( X i , t j ) , udxS = {u%, u | = (u$)h a(x, t)kexr)(—a(x,t)k/e) 1 — exp(—a(x, t)k/e) ' b(x,t)hexp(—b(x,t)h/e) 1 — exp(—b(x, t)h/e). ' u r f(a;,t) - u r f(x - h,t) h ud(x,t) - u" (x,t- k) k ud(x + h,t) - ud(x,t) 7 i (s , t , fc) •y2(x,t,h) uds(x,t) = ul(x,t) For this difference scheme we can establish the following discrete energy inequality. L e m m a 6.4 Let ud(x,t) be the solution of (6.63) and let mesh sizes h and k satisfy inequality b(x,t)h < a(x,t)k for all (x,t) G Gd- Then when e,h and k are sufficiently small we have \\A\l + Il7i«rll2 + I IVl^4 \ \ 2 s < MK(h, fc, e), (6.64) where K(h,k,e) = ^ f c ^ ^ ^ : ^ 2 -I- H-yxTx^ H? H- II^Ti^^ll? -I- II^ /O^ TWM^ II? -|-j=lz = l N h^2v(ih,sk)2, s = 1,- • •, J, J — [T/k], p = k/e. i=i \v\\2. Chapter 6. Numerical Methods of Some Singular Perturbation Problems 152 Proof: Essentially, the proof is a simulation of the continuous case but much more complicated because of the fitting factors 71 and 72. A sketch of the proof is in [107]. • R e m a r k 6.2 If fitting factors are not used, then the condition b(x,t)h < a(x,t)k can be removed in the proof. • Next we prove the uniform convergence. We assume that enough compatibil i ty conditions are satisfied so that the solution u(x, t) of problem (6.48) belongs to C3(G). Based on the asymptotic expansion we may assert that the following estimates | | ^ | < M(e- + 0 < * < 3, 0 < i < k (6.65) hold. For convenience we also assume c\k < h < c2k, ci > 0, a(x,t)/b(x,t) > c 2 > 0. (6.66) Using derivative estimate (6.65), it is not difficult to verify L{h'k\u(x,t)-ud(x,t)) = 0(^2 + -e), k2 (u - ud)\t=0 = 0, (u - ud)\t=k = min( — ,k), {u - u%=0 = 0, (u - u%=l = 0. Apply ing discrete energy inequality (6.64), yields \u-u%<M(- + -). (6.67) e2 e O n the other hand, we can obtain 4 M ) ( f i i - «") = L{h>k)(u0 + e^ 0 ) + + - / / - x = 0(e + h + k + exp(—m 0—-—), (m0 > 0) Chapter 6. Numerical Methods of Some Singular Perturbation Problems 153 and («i - ud)\t=0 = 0(e), (uu-uDl^k = 0(1) , («i - ud)\x=0 = 0(e"), (u! - ud)\x=l = 0, where n is an arbitrary positive number. First making a change of variable such that the boundary value conditions become homogeneous, then using the discrete energy inequality (6.64), we have - ud\\s < M(^ /max(e , /0 + k + en/h). (6.68) Hence, from the remainder estimate (6.62) of the asymptotic solution u\, we obtain another error estimate | | M -U%< M(^/max(e, h) + k + en/h). (6.69) Combining (6.67) with (6.69) we achieve our uniformly convergent estimate. Theorem 6.4 Supposing that (6.66) is satisfied and the solution of (6.^8) u(x,t) £ C3(G). Then the solution ud(x,t) of the discrete problem (6.63) uniformly converges to u(x,t) in the sense of the discrete energy norm, i.e. I l u - u ^ l . < Mh*, s = 0,1, (6.70) Proof: Using (6.69), (6.66) as e 5/ 2 < h and using (6.67), (6.66) as e 5 / 2 > h, we. have (6.70) immediately. • R e m a r k 6.3 If we apply the remainder estimate (6.61) and corresponding procedure there we may discuss the uniform convergence under the compatibility conditions Cl and 0 2 only. • Chapter 7 Conclusion and Future W o r k 7.1 S u m m a r y and conclusions The main focus of the thesis is on constrained ordinary differential equations ( D A E s ) , constrained partial differential equations ( P D A E s ) , and their applications. A new class of methods for solving high index D A E s has been developed, which we call the Sequential Regularization Method or S R M for short. These methods offer significant advantages over some known solution techniques, such as regularization and stabiliza-tion methods, and are applied to the nonstationary Navier-Stokes equations governing incompressible fluid flow and to a mathematical model of reservoir simulation. High index D A E s (index > 2) are usually difficult to discretize directly [29, 86]. We thus need to reformulate the original problem as a better behaved problem before discretization. Index reduction with stabilization is a popular reformulation for the numerical solution of semi-explicit high index D A E s . Another class of reformulations is regularization, where the D A E is replaced by a better behaved nearby problem. Such a method reduces the size of the system to be solved and avoids the deriva-tives of the algebraic constraints associated with the D A E problem. Regularization is particularly suitable for problems with certain singularities where the constraint Jacobian does not have full rank. Unfortunately, this approach often yields very stiff problems, which accounts for its lack of popularity in practice. The S R M is proposed to overcome this difficulty. It keeps the benefits of regularization methods and avoids the need to use stiff solvers for the regularized problems, because the regularization 154 Chapter 7. Conclusion and Future Work 155 parameter does not need to be very small. Thus, we obtain an important improvement over usual regularization methods which leads to easier numerical methods (explicit t ime discretization for regularized problems). The S R M also provides cheaper and more efficient alternatives to the usual stabilization methods for some choices of pa-rameters and stabilization matrix. We first propose and analyze the method for linear index-2 D A E s . Then we extend it to nonlinear index-2 and index-3 D A E s . This is especially useful in applications such as constrained mult ibody systems which are of index-3. Numerical experiments show that the method is useful and efficient for problems with and without singularities. Whi le a significant body of knowledge about the theory and numerical methods for D A E s has been accumulated, almost none has been extended to partial differential-algebraic equations ( P D A E s ) . As a first attempt we provide a comparative study be-tween stabilization and regularization (or pseudo-compressibility) methods for D A E s and P D A E s , using the Navier-Stokes equations as an instance of P D A E s . Compared with stabilization methods, regularization methods can avoid imposing an artificial boundary condition for the pressure p. This is a feature for P D A E s not shared wi th D A E s . We generalize the S R M to the nonstationary incompressible Navier-Stokes equations. Similar to D A E s , explicit schemes in the time direction can be used for the P D A E because of the reduced stiffness (taking the regularization parameter rel-atively large) or even essential nonstiffhess obtained for some choice of parameters. Unl ike usual regularization methods, the time step restriction for the explicit scheme can be independent of the regularization parameter e. The time step restriction is further loosened for the case of small viscosity. A simple discretization (such as the forward Euler difference in time and a first-order scheme in space) is analyzed and implemented. Numerical results support our theoretical results. The method works for both two- and three-dimensional problems. In recent years considerable attention has been devoted to numerical reservoir Chapter 7. Conclusion and Future Work 156 simulation, e.g. miscible displacement in porous media. We have applied the S R M idea to the pressure-velocity equation in its 2-dimensional mathematical model equa-tions. This procedure is first analyzed at the differential level and then discretized by finite element methods. Theoretical analysis and numerical experiments show that this procedure converges at the rate of 0(e) per iteration, where e is a small positive number. The fast convergence rate makes our iterative method dramatically different from penalty methods . In addition, the perturbation parameter e does not have to be carefully chosen, unlike the case for other iterative methods. Indeed, our numerical experiments show that two iterations are usually enough for a variety of problems. Compared with mixed finite element methods, the discrete version of our scheme can provide the same accurate approximations for velocity and pressure, which is crucial in reservoir problems since velocity is intimately involved in the concentration equa-tion. However, in contrast to mixed finite element methods, our scheme requires only the solution of symmetric positive definite linear systems which have a smaller number of degrees of freedom corresponding to the velocity variable. Since our method com-pletely decoupled the velocity and pressure variables, the so-called Babuska-Brezzi condition is not needed in constructing the finite dimensional spaces for velocity and pressure. The method is easily applied to three-dimensional problems. Another topic of this thesis is singular perturbation problems, which come from many applied areas and regularized problems. We discuss numerical solutions of several singular perturbation problems. Uniformly convergent schemes wi th respect to the perturbation parameter e are constructed and analyzed for nonlinear repulsive and attractive turning point problems and a second-order hyperbolic problem. We are the first to be able to construct uniformly convergent schemes for these problems. Also, an interesting spurious solution phenomenon from an upwinding scheme is analyzed for an elliptic turning point problem. We find that the spurious solution is caused by a m i l d instability of the problem (the constant for the stability inequality is of 0 (^ ) ) . Chapter 7. Conclusion and Future Work 157 This type of instabili ty is not as serious as supersensitivity [73, 74]. It can be handled by using higher order upwinding schemes as long as their accuracy hr <C e when e is not too small . 7.2 Discussion of future work The results of this thesis can be extended in a number of directions. Efficient s imulat ion of kinemat ic chains w i th closed loops The kinematic chain problem is an example of mult ibody systems, such as robot manipulators. Consider a chain consisting of n links. The numerical simulation of the problem is usually treated as two separate problems: (i) the forward dynamics problem for computing system accelerations, and (n) the numerical integration prob-lem for advancing the state in time. In recent years, many different algorithms have been proposed for solving the forward dynamics problem wi th tree structure (without closed loops), ranging in computational complexity from 0(n3) (e.g. the composite rigid body method ( C R B M ) [116]) to 0{n) (e.g. the articulated-body method ( A B M ) [49]). However, it has never been possible to construct an 0(n) algorithm for the chain problem wi th many closed loops. The S R M (plus explicit discretization) opens up a way to do this. According to the idea of the S R M , we can remove the extra constraints caused by closed loops and incorporate them into external force terms of the system. The reformulated problem has the same structure as that of the problem without closed loops. We thus expect to develop and test an 0(n) algorithm for closed-loop chains. Fu l ly nonlinear D A E s Consider a fully nonlinear i ndex - f D A E ^ f(t,x,x,---,x^ \y)? Chapter 7. Conclusion and Future Work 158 0 = g(t,x), (i.e. where algebraic variables y appear nonlinearly together wi th the differential variables x.) Some mechanical multibody systems are in this form with v = 2, e.g. the motion of a point mass on a parabolic orbit with the forces of gravitation and Coulomb friction (see [4]). The S R M can be applied to this problem. For instance, corresponding to the index-2 problem (with v = 1) we have: for s = 1,2, • • •, x s — fifi xsiys)i ys = Vs-i + ^e(g(t,xs)), where the function e should be chosen such that the matrix fyeggx has positive eigen-values. Efficient s imulat ion of real fluid problems using S R M Because the computations for u and p are uncoupled and explicit time discretiza-tion can be used, we expect that the S R M incorporated with a finite difference or finite element discretization would provide a cheap and efficient method for simulating more realistic fluid problems. For long time simulations, a reinitialization technique (e.g. projecting back to the divergence free space after a number of time steps) may be useful. A comparative study between the S R M and other methods would be interesting. Solving the system result ing from the discret ization of the operator / + ^-graddzu This operator comes from using S R M to solve the nonstationary Navier-Stokes equations with cti ^ 0. The system is easily made to be banded symmetric positive definite. Hence a direct method can be used to solve it. A n interesting observation is that the usual iterative methods do not work well. This is probably due to the lack of ell ipticity of the system. Some research on solving this problem using mult igr id Chapter 7. Conclusion and Future Work 159 and domain decomposition techniques (at least for e not too small) is about to be completed by Arno ld , Falk and Winther [3]. Based on a technique described in [57], iterative methods (including multigrid) would be feasible under some pre-processing of the system (to increase the ell ipticity). This was also suggested by W . Hackbusch in a private communication. Bibl iography C. Amirouche and V . Giraul t , Decomposition of vector spaces and application to the Stokes problem in arbitrary dimension, Czechoslovak Mathematical J . 44(1994), No. 119, pp. 109-140. F . Amirouche, Computational Methods in Multibody dynamics, Prentice-Hall , 1992. D . N . Arno ld , Private communication, 1995. M . Arno ld , A perturbation analysis for the dynamical simulation of mechani-cal multibody systems, Preprint 94/24, Universitat Rostock, Fachbereich Mathe-matik, Germany, 1994. K . Arrow, L . Hurwicz and H . Uzawa, Studies in Nonlinear Programming, Stan-ford University Press, 1968. U . Ascher, On some difference schemes for singular singularly-perturbed bound-ary value problems, Numer. Ma th . 46 (1985), 1-30. U . Ascher, On numerical differential algebraic problems with application to semi-conductor device simulation, S I A M J . Numer. A n a l . 26(1989), 517-538. U . Ascher, H . Ch in and S. Reich, Stabilization of DAEs and invariant manifolds, Numer. M a t h . 67 (1994), 131-149. U . Ascher, H . Ch in , L . Petzold and S. Reich, Stabilization of constrained me-chanical systems with DAEs and invariant manifolds, Journal of Mechanics of Structures and Machines 23 (1995), 135-157. U . Ascher, J . Christiansen and R . D . Russell, Collocation software for boundary value ODEs, Trans. Ma th . Software 7(1981), 209-222. U . Ascher and Ping L i n , Sequential regularization methods for higher index DAEs with constraint singularities: Linear index-2 case, S I A M J . Numer. A n a l . , to appear. U . Ascher and Ping L i n , Sequential Regularization Methods for Nonlinear Higher Index DAEs, S I A M J . Sci. Comput. , to appear. U . Ascher, P . A . Markowich, P. Pietra, C . Schmeiser, A phase plane analysis of transonic solutions for the hydrodynamic semiconductor model, M a t h . Models & Methods in Appl ied Sciences 1(1991), 347-376. 160 Bibliography 161 U . Ascher, R . Mattheij and R. Russell, Numerical Solution of Boundary Value Problems for Ordinary Differential Equation, Prentice-Hall , 1988. U . Ascher and L . Petzold, Projected implicit Runge-Kutta methods for differential-algebraic equations, S I A M J . Numer. Ana l . 28 (1991), 1097-1120. U . Ascher and L . Petzold, Stability of Computational Methods for Constrained Dynamics Systems, S I A M J . Sci. Comput. 14 (1993), 95-120. J . Baumgarte, Stabilization of constraints and integrals of motion in dynamical systems, Comp. Ma th . A p p l . Mech. Eng. 1 (1972), 1-16. G . Bader and U . Ascher, A new basis implementation for a mixed order boundary value ODE solver, S I A M J . Scient. Stat. Comput. 8 (1987), 483-500. E . Bayo and A . Avel lo , Singularity-free augmented Lagrangian algorithms for constrained multibody dynamics, J . Nonlinear Dynamics 5 (1994), 209-231. E . Bayo, J . G . de Jalon and M . A . Serna, A modified Lagrangian formulation for the dynamic analysis of constrained mechanical systems, Computer Methods in Appl ied Mechanics and Engineering 71(1988), 183-195. A . E . Berger, H . Han and R . B . Kellogg, A priori estimates and analysis of a numerical method for a turning point problem, M a t h . Comp. 42(1984, 465-491. L . E . Bobisud, The second initial-boundary value problem for a linear parabolic equation with a small parameter, M i c h . Ma th . J . 15(1969), 495-504. L P . Boglaev, On numerical integration of a singular- perturbed initial-value prob-lem for ordinary differential equation, Zh . Vychis l . mat. i mat. fiziki (Comp. Ma th . & Math . Physics) 25(1985), 1009-1022. Y . Boyarintsev, Reguljarnyje i Singularnyje Sistemy Linejnych Obyknovennych Differencial'nych Uravnenij, Nauka (Sibirskoje otdelenije), Novosibirsk, 1980. Y . Boyarintsev, Methods of Solving Singular Systems of Ordinary Differential Equations, John Wiley &s Sons, 1992 (originally published in 1988 in Russian). A . Brandt , Multigrid techniques: 1984 guide with applications to fluid dynamics, The Weizmann Institute of Science, Rehovot, Israel, 1984. A . Brandt and I. Yavneh, Inadequacy of first order upwind difference schemes for some recirculation flows, J . Comp. Phys. 93(1991), 128-143. B . Brefort, J . M . Ghidaglia and R. Temam, Attractors for the penalized Navier-Stokes equations, S I A M J . Ma th . Ana l . 19 (1988), 1-21. Bibliography 162 K . Brenan, S. Campbel l and L . Petzold, The numerical solution of higher in-dex differential/algebraic equations by implicit Runge-Kutta methods, S I A M J . Numer. A n a l . 26 (1989). F . Brezz i , J . Douglas, Jr. , and L. D . Mar in i , Two families of mixed finite elements for second order elliptic problems, Numer. Math . , 47(1985), pp. 217-235. F . Brezz i , and M . Fortin, Mixed and Hybrid Finite Element Methods, Springer-Verlag, New York, 1991. S.L. Campbel l , Regularizations of linear time varying singular systems, Auto-matica 20(1984), 365-370. C . Canuto, M . Y . Hussaini, A . Quarteroni and T . A . Zang, Spectral Methods in F l u i d Dynamics, Springer-Verlag, 1987. P. Carter, Computational methods for the shape from shading problem, P h . D thesis, University of Br i t i sh Columbia, 1993. H . S. Ch in , Stabilization methods for simulations of constrained multibody dy-namics, P h . D thesis, University of Br i t i sh Columbia , 1995. A . J . Chorin , Numerical solution of the Navier-Stokes equations, M a t h . Comp. 22 (1968), 745-762. E . P . Doolan, J . J . H . Mi l le r and W . H . A . Schilders, Uniform Numerical Methods for Problems with Initial and Boundary Layers, Dubl in , Boole Press, 1980. J . Douglas, Jr. , The numerical solution of miscible displacement in porous media, in : Computational Methods in Nonlinear Mechanics (J . T . Oden, Ed . ) , Amster-dam, North Holland, 1980. J . Douglas, Jr. , Simulation of miscible displacement in porous media by a modified method of characteristics procedure, in: Lecture Notes in M a t h . 912, Springer-Verlag, (1982). J . Douglas, Jr. , R . E . Ewing and M . F . Wheeler, The approximation of the pressure by a mixed method in the simulation of miscible displacement, R A I R O Numer. A n a l . 17(1983), pp. 17-33. J . Douglas, Jr , R . E . Ewing and M . F . Wheeler, A time-discretization proce-dure for a mixed finite element approximation of miscible displacement in porous media. R A I R O Numer. A n a l . 17(1983), pp. 249-265. R . C . Duran, On the approximation of miscible displacement in porous media by a method of characteristics combined with a mixed method, S I A M J . Numer. A n a l . 25(1988), pp. 989-1001. Bibliography 163 H . W . Engl , M . Hanke and A . Neubauer, Tikhonov regularization of nonlin-ear differential-algebraic equations, Institutsbericht Nr . 385, Johannes-Kepler-Universitat, Institut fur Mathematik, L inz , 1989. R. E . Ewing (Ed.) , The Mathematics of Reservoir Simulation, S I A M Philadel-phia, 1983. R. E . Ewing and T. F . Russell, Efficient time-stepping methods for miscible displacement problems in porous media, S I A M J . Numer. A n a l . 19(1982), pp. 1-67. R . E . Ewing, T . F . Russell and M . F . Wheeler, Simulation of miscible displace-ment using mixed methods and a modified method of characteristics, S P E 12241, 7th S P E Symposium on Reservoir Simulation, San Francisco, 1983. R. E . Ewing , T . F . Russell and M . F . Wheeler, Convergence analysis of an ap-proximation of miscible displacement in porous media by mixed finite elements and a modified method of characteristics, Comp. Meth . A p p l . Mech. Engrg. 47(1984), pp. 73-92. R . E . Ewing and M . F . Wheeler, Galerkin Methods for miscible displacement problems in porous media, S I A M J . Numer. A n a l . 17(1980), pp. 351-365. R . Featherstone, Robot Dynamics Algorithms, Kluwer Academic Publishers, Nor-well , M A , 1987. M . Fort in and R . Glowinski , Augmented Lagrangian Methods: Applications to the Numerical Solution of Boundary-Value Problems, Studies in Mathematics and Its Applications 15, North-Holland, 1983. C . W . Gear, The simultaneous numerical solution of differential-algebraic equa-tions, I E E E Trans. Circuit Theory, CT-18(1971), 89-95. C . W . Gear, H . H . Hsu and L . Petzold, Differential- algebraic equations revisited, Proc. O D E Meeting, Oberwolfach, Germany, 1981. R . Geel, Singular Perturbations of Hyperbolic Type, Mathematical Center Tracts 98, Amsterdam, 1978. -V . Giraul t and P . A . Raviart , Finite Element Methods for Navier-Stokes Equa-tions, Springer, Berlin-Heidelberg 1986. P . M . Gresho, Some current issues relevant to the incompressible Navier-Stokes equations, Comput. Methods A p p l . Mech. Engrg. 87 (1987), 201-252. E . Griepentrog and R. Marz , Differential-algebraic equations and their numerical treatment, Teubner Texte zur Mathematik 88, Leipzig, 1986. Bibliography 164 W . Hackbusch, Multigrid Methods and Applications, Springer-Verlag, Ber l in , Hei -delberg 1985. E . Hairer and G . Wanner, Solving Ordinary Differential Equations II, Springer-Verlag, 1991. M . Hanke, Regularization of differential -algebraic equations revisited, Preprint Nr . 92-19, Humboldt- Univ . Ber l in , Fachbereich Mathematik, 1992. E . J . Haug, Computer-Aided Kinematics and Dynamics of Mechanical Systems Volume I: Basic Methods, A l l y n and Bacon, 1989. P . W . Hemker, A Numerical Study of Stiff Two-point Boundary Problems, M a t h -ematical Center Tracts 80, Amsterdam 1977. J . G . Heywood, The Navier-Stokes equations: On the existence, regularity and decay of solutions, Indiana Univ . Ma th . J . , 29 (1980), 639-681. J . G . Heywood and R. Rannacher, Finite element approximation of the nonsta-tionary Navier-Stokes problem, I. Regularity of solutions and second-order error estimates for spatial discretization, S I A M J . Numer. A n a l . 19(1982), 275-311. J . G . Heywood and R . Rannacher, Finite element approximation of the nonsta-tionary Navier-Stokes problem, IV: Error analysis for second-order time dis-cretization, S I A M J . Numer. A n a l . 27 (1990), 353-384. C. Hirsch, A general analysis of two-dimensional convection schemes, V K I Lec-ture Series 1991-02 on Computational F l u i d Dynamics, Von Karman Institute, Brussels, Belgium, 1991. T . Y . Hou and B . T . R . Wetton, Convergence of a finite difference scheme for the Navier-Stokes equations using vorticity boundary conditions, S I A M J . Numer. A n a l . , 29 (1992), 615-639. M . K . Kadalbajoo and Y . N . Reddy, Initial-value technique for a class of nonlinear singular perturbation problems, J . Op t im. Theory and Appls . 53(1987), 395-406. L . Kalachev and R . O'Malley, Jr. , The regularization of linear differential-algebraic equations, S I A M J . Ma th . Ana l . , 1996, to appear. L . Kalachev and R. O'Malley, Jr. , Regularization of nonlinear differential-algebraic equations, S I A M J . Ma th . Ana l . 25 (1994), 615-629. L . Kalachev and R. O'Malley, Jr., Boundary value problems for differential-algebraic equations, N u m . Func. Ana l . A p p l . 16 (1995), 363-378. Bibliography 165 [71] M . Knorrenschild, Differential/algebraic equations as stiff ordinary differential equations, S I A M J . Numer. A n a l . 29 (1992), no. 6, 1694-1715. [72] H . - O . Kreiss and J . Lorenz, Initial-Boundary Value Problems and the Navier-Stokes Equations, Pure and Appl ied Mathematics, V o l . 136, Academic Press, 1989. [73] J . Laforgue and R . E . O'Malley, On the motion of viscous shocks and the su-persensitivity of their steady-state limits, Methods and Applications of Analysis 2(1994), 465-487. [74] J . Y . Lee and M . J . Ward, On the asymptotic and numerical analyses of expo-nentially ill-conditioned singularly perturbed boundary value problems, Studies in Appl ied M a t h . 94 (1995), 271-326. [75] M . Lees, A priori estimates for the solutions of difference approximations to parabolic partial differential equations, Duke M a t h J . 27 (1960), 297-311. [76] R . J . LeVeque, Numerical Methods for Conservation Laws, Lectures in Mathe-matics E T H Zurich, Birkhauser Verlag, Basel, 1990. [77] P. L i n , Numerical solutions of some singular perturbation problems, M . Sc. thesis, Nanjing University, Nanjing, China, 1987. (in Chinese) [78] P. L i n , Remarks on a mildly nonlinear turning point problem, App l i ed Ma th . Letters 6(1993), 29-32. [79] P. L i n , A numerical method for quasilinear singular perturbation problems with turning points, Computing 46(1991),155-164. [80] P. L i n , A sequential regularization method for time-dependent incompressible Navier-Stokes equations, S I A M J . Numer. Ana l . , to appear. [81] P. L i n and Y . Su, Numerical solution of quasilinear singularly perturbed ordinary differential equations without turning points, A p p l . Ma th . Mech. (English Ed.) 10(1989). [82] P. L i n and D . Q. Yang, An iterative perturbation method for the pressure equa-tion in the simulation of miscible displacement in porous media, S I A M J . Sci . Comput. , submitted. [83] J . Lorenz, Combinations of initial and boundary value methods for a class of singular perturbation problems, In: P . W . Hemker and J . J . H . Mi l l e r (eds.), N u -merical Analysis of Singular Perturbation Problems, London: Academic Press, 1979, 295-315. Bibliography 166 [84] P. Lotstedt, On a penalty function method for the simulation of mechanical sys-tems subject to constraints, Royal Inst, of Tech. T R I T A - N A - 7 9 1 9 , Stockholm, 1979. [85] R. Marz , On tractability with index 2, Preprint 109, Humboldt -Univ . Ber l in , Sektion Mathematik, 1986. [86] R . Marz , Numerical methods for differential- algebraic equations, Part I: Char-acterizing DAEs, Preprint Nr . 91-32/1, Humboldt-Universitat zu Ber l in , 1991. [87] J . Necas, Les methodes directes en theorie des equations elliptiques, Academia, Prague, 1967. [88] K . Ni i j ima , On the behavior of solutions of a singularly perturbed boundary value problem with a turning point, S I A M J . Ma th . A n a l . 9(1978), 298-311. [89] R . O'Mal ley , Jr. , Singular Perturbation Methods for Ordinary Differential Equa-tions, Springer, New York, 1991. [90] E . O 'Riordan and M . Stynes, A globally uniformly convergent finite element method for a singularly perturbed elliptic problem in two dimensions, M a t h . Comp. 57(1991), 47-62. [91] K . Park and J . Chiou, Stabilization of computational procedures for constrained dynamical systems, J . of Guidance 11 (1988), 365-370. [92] L . R . Petzold, A brief history of numerical methods for differential-algebraic equa-tions, Lecture notes in the meeting of the 50th anniversary of the journal " Mathematics of Computations", 1993. [93] L . R . Petzold and P. Lotstedt, Numerical solution of nonlinear differential equa-tions with algebraic constraints II: Practical implications, S I A M J . Sci. Stat. Comput . 7(1986), 720-733. [94] L . R . Petzold, Yuhe Ren and T . Maly , Numerical solution of differential-algebraic equations with ill-conditioned constraints, Technical Report 93-59, University of Minnesota, 1993 [95] P. Rabier and W . Rheinboldt, On the computation of impasse points of quasilin-ear differential algebraic equations, Ma th . Comp. 62 (1994), no. 205, 133-154. [96] R . Rannacher, On Chorin's projection for the incompressible Navier-Stokes equa-tions, In: Rautmann, et al. (eds.): Navier-Stokes equations: Theory and numeri-cal methods. Proc. Oberwolfach Conf., 19.-23.8.1991, Springer, Heidelberg 1992. [97] R. Rannacher, On the numerical solution of the incompressible Navier-Stokes equations, Z. angew. Ma th . Mech. 73(1993), 203-216. Bibliography 167 [98] T . F . Russell, Finite elements with characteristic finite element method for a mis-cible displacement problem, S P E 10500, Proc. 6th S P E Symposium on Reservoir Simulation, Dallas, 1982. [99] T . F . Russell, Time stepping along characteristics with incomplete iteration for a Galerkin approximation of miscible displacement in porous media, S I A M J . Numer. A n a l . 22 (1985), pp. 970-1003. [100] P. A . Raviart and J . M . Thomas, A mixed finite element method for second order elliptic problems, in : I. Galligani and E . Magenes, Eds., Mathematical Aspects of the Fini te Element Method, Lecture Notes in Mathematics 606, Springer-Verlag, Ber l in and New York, 1977, pp. 292-315. [101] J . Shen, On error estimates of projection methods for the Navier-Stokes equa-tions: First order schemes, S I A M J . Numer. A n a l . 29 (1992), 57-77. [102] D . Sidilkover and U . M . Ascher, A multigrid solver for the steady state Navier-Stokes equations using the pressure-Poisson formulation, Technical Report 94-3, Dept. of Computer Science, University of Br i t i sh Columbia, 1994. [103] R . F . Sincovec, A . M . Erisman, E . L . Y i p and M . A . Epton, Analysis of descriptor systems using numerical algorithms, I E E E Trans. Au t . Control , AC-26(1981), 139-147. [104] D . R . Smi th and J .T . Palmer, On the behavior of the solution of the telegraphist's equation for large absorption, Arch . Rat . Mech. A n a l . 39(1970), 146-157. [105] Y . Su and P. L i n , Numerical solution of singular perturbation problem for hy-perbolic differential equation with characteristic boundaries, Proceedings of the B A I L V Conference, Shanghai, Boole Press, Dublin(1988), 20-24. [106] Y . Su and P. L i n , Uniform difference scheme for a singularly perturbed linear 2nd order hyperbolic problem with zero-th order reduced equation, Appl ied M a t h . Mech . 11 (English Ed.)(1990), No. 4. [107] Y . Su and P. L i n , An exponentially fitted difference scheme for the hyperbolic-hyperbolic singularly perturbed initial-boundary problem, Appl ied M a t h . Mech. 12 (English Ed.)(1991), 237-245. [108] R . Temam, Navier-Stokes Equations, North-Holland, Amsterdam, 1977. [109] A . N . Tikhonov and V . Y a . Arsenin, Methods for Solving Ill-posed Problems, Nauka, Moscow, 1979. [110] S. Turek, Tools for simulating nonstationary incompressible flow via discretely divergence-free finite element models, Preprint, Universitat Heidelberg, M a y 1992. Bibliography 168 [111] A . Vasil 'eva and V . Butuzov, Asymptotic Expansion of Solutions of Singularly Perturbed Equations, Nauka, Moscow, 1973. [112] A . Vasil 'eva and V . Butuzov, Singularly Perturbed Equations in the Critical Case, M R C Technical Summary Report # 2039, 1980. [113] R . Vulanovic, A uniform numerical method for quasilinear singular perturbation problems without turning points, Computing 41(1989), 97-106. [114] R. Vulanovic, On numerical solution of a mildly nonlinear turning point prob-lem, R A I R O M a t h . Model . Numer. A n a l . 24(1990), 765-784 . [115] R. Vulanovic and P. L i n , Numerical solution of quasilinear attractive turning point problems, Computers Ma th . Appl ic . 23(1992), 75-82. [116] M . W . Walker and D . E . Or in , Efficient dynamic computer simulation of robotic mechanisms, Journal of Dynamic Systems, Measurement, and Control 104 (1982), 205-211. [117] B . Wetton, Error analysis for Chorin's original fully discrete projection method and regularizations in space and time, Preprint, Institute of Appl ied M a t h , U n i -versity of Br i t i sh Columbia, 1995. [118] G . B . Whi tham, Linear and Nonlinear Waves, Wiley, New York, 1974. [119] W . X i e , A sharp pointwise bound for functions with L2-Laplacians and zero boundary values on arbitrary three-dimensional domains, Indiana University M a t h . J . 40 (1991), 1185-1192. [120] D . Q. Yang, Mixed methods with dynamic finite element spaces for miscible displacement in porous media, J . Comput. A p p l . M a t h . 30(1990), pp. 313-328. [121] I. Yavneh, Multigrid techniques for incompressible flows, P h . D . thesis, The Weizmann Institute of Science, Israel, 1991.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Regularization methods for differential equations and...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Regularization methods for differential equations and their numerical solution Lin, Ping 1995
pdf
Page Metadata
Item Metadata
Title | Regularization methods for differential equations and their numerical solution |
Creator |
Lin, Ping |
Date Issued | 1995 |
Description | Many mathematical models arising in science and engineering, including circuit and device simulation in VLSI, constrained mechanical systems in robotics and vehicle simulation, certain models in early vision and incompressible fluid flow, lead to computationally challenging problems of differential equations with constraints, and more particularly to high-index, semi-explicit differential-algebraic equations (DAEs). The direct discretization of such models in order to solve them numerically is typically fraught with difficulties. We thus need to reformulate the original problem into a better behaved problem before discretization. Index reduction with stabilization is one class of reformulations in the numerical solution of high index DAEs. Another class of reformulations is called regularization. The idea is to replace a D A E by a better behaved nearby system. This method reduces the size of the problem and avoids the derivatives of the algebraic constraints associated with the D A E . It is more suitable for problems with some sort of singularities in which the constraint Jacobian does not have full rank. Unfortunately, this method often results in very stiff systems, which accounts for its lack of popularity in practice. In this thesis we develop a method which overcomes this difficulty through a combination of stabilization and regularization in an iterative procedure. We call it the sequential regularization method (SRM). Several variants of the S RM which work effectively for various circumstances are also developed. The S RM keeps the benefits of regularization methods and avoids the need for using a stiff solver for the regularized problem. Thus the method is an important improvement over usual regularization methods and can lead to improved numerical methods requiring only solutions to linear systems. The S RM also provides cheaper and more efficient methods than the usual stabilization methods for some choices of parameters and stabilization matrix. We propose the method first for linear index-2 DAEs. Then we extend the idea to nonlinear index-2 and index-3 problems. This is especially useful in applications such as constrained multibody systems which are of index-3. Theoretical analysis and numerical experiments show that the method is useful and efficient for problems with or without singularities. While a significant body of knowledge about the theory and numerical methods for DAEs has been accumulated, almost none of it has been extended to partial differential-algebraic equations (PDAEs). As a first attempt we provide a comparative study between stabilization and regularization (or pseudo-compressibility) methods for DAEs and PDAEs, using the incompressible Navier-Stokes equations as an instance of PDAEs. Compared with stabilization methods, we find that regularization methods can avoid imposing an artificial boundary condition for the pressure. This is a feature for PDAEs not shared with DAEs. Then we generalize the S R M to the nonstationary incompressible Navier-Stokes equations. Convergence is proved. Again nonstiff time discretization can be applied to the S RM iterations. Other interesting properties associated with discretization are discussed and demonstrated. The S RM idea is also applied to the problem of miscible displacement in porous media in reservoir simulation, specifically to the pressure-velocity equation. Advantages over mixed finite element methods are discussed. Error estimates are obtained and numerical experiments are presented. Finally we discuss the numerical solution of several singular perturbation problems which come from many applied areas and regularized problems. The problems we consider are nonlinear turning point problems, a linear elliptic turning point problem and a second-order hyperbolic problem. Some uniformly convergent schemes with respect to the perturbation parameter are constructed and proved. A spurious solution phenomenon for the upwinding scheme is analyzed. |
Extent | 7038853 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-02-18 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0079749 |
URI | http://hdl.handle.net/2429/4782 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 1996-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_1996-09118X.pdf [ 6.71MB ]
- Metadata
- JSON: 831-1.0079749.json
- JSON-LD: 831-1.0079749-ld.json
- RDF/XML (Pretty): 831-1.0079749-rdf.xml
- RDF/JSON: 831-1.0079749-rdf.json
- Turtle: 831-1.0079749-turtle.txt
- N-Triples: 831-1.0079749-rdf-ntriples.txt
- Original Record: 831-1.0079749-source.json
- Full Text
- 831-1.0079749-fulltext.txt
- Citation
- 831-1.0079749.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0079749/manifest