A SUPRA-CONVERGENT SCHEME FOR T H E SOLUTION OF DIFFERENTIAL EQUATIONS ON AN ARBITRARY MESH by KAREN SCOTT B.Sc, The University of Calgary, 1986 A THESIS SUBMITTED IN PARTIAL FULLFILMENT OF T H E REQUIREMENTS FOR T H E DEGREE OF MASTER OF SCIENCE in T H E FACULTY OF GRADUATE STUDIES Department of Mathematics Institute of Applied Mathematics We accept this thesis as conforming to the required standard T H E UNIVERSITY OF BRITISH COLUMBIA July 1988 ©Karen Scott, 1988 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 The University of British Columbia 1956 Main Mall Vancouver, Canada V6T 1Y3 DE-6G/81) A B S T R A C T T h i s thesis descr ibes a c o m p a c t n u m e r i c a l scheme to solve second-order l i nea r d i f ferent ia l equa t ions w i t h D i r i c h l e t b o u n d a r y c o n d i t i o n s o n a r b i t r a r y meshes. T h e scheme uses a s taggered g r i d to achieve second-order a c c u r a c y o n a n o n u n i f o r m m e s h . R e g u l a r a n d s i n g u l a r l y p e r t u r b e d p r o b l e m s i n one a n d h igher d i m e n s i o n s are c o n -s ide red . N u m e r i c a l e x p e r i m e n t s are p resen ted w h i c h s u p p o r t the t h e o r e t i c a l resu l t s . T h e m e t h o d s p resen ted are also a p p l i c a b l e to equa t ions w i t h o ther t ypes o f b o u n d a r y c o n d i t i o n s , a n d to n o n l i n e a r a n d h ighe r -o rde r d i f ferent ia l equa t ions . 11 Table of Contents 1. Introduction 1 1.1 Definition of Terms 2 1.2 Previous and Related Work 4 2. A Second-Order Three-Point Scheme 8 2.1 Deriving a Second-Order Three-Point Scheme 8 2.2 Determining the Order of the Error 16 2.3 Numerical Results 23 3. A n Extension to Two Dimensions 36 3.1 Deriving a Two-Dimensional Midpoint Scheme 38 3.2 Deriving a Second-Order Nine-Point Formula 40 3.3 Numerical Results 47 4. Singularly Perturbed Problems 53 4.1 Exponential Fitting 53 4.2 Discretization Errors 59 4.3 Numerical Results 61 5. Other Differential Equations 67 6. Conclusion 70 7. References 72 iii L i s t o f F i g u r e s 1. Standard Evaluation Points on a Mesh 9 2. Evaluation Points on a Staggered Mesh 10 3. Determining a Sine-Type Mesh 24 4. Usual Three-Point Formula : Results for Problems 1 to 6 27 5. Three-point Midpoint Formula : Results for Problem 1 28 6. Three-point Midpoint Formula : Results for Problem 2 29 7. Three-point Midpoint Formula : Results for Problem 3 30 8. Three-point Midpoint Formula : Results for Problem 4 31 9. Three-point Midpoint Formula : Results for Problem 5 32 10. Three-point Midpoint Formula : Results for Problem 6 33 11. Three-point Midpoint Combined Results for Problems 1 to 6 34 12. Three-point Trapezoidal Combined Results for Problems 1 to 6 35 13. The Two-dimensional Mesh 41 14. Results for Two-dimensional Problem 1 49 15. Results for Two-dimensional Problem 2 50 16. Results for Two-dimensional Problem 3 51 17. Results for Two-dimensional Problem 4 52 18. Results for Singular Perturbation Problems 1,2 on Uniform Meshes 63 19. Results for Singular Perturbation Problems 1,2,3 on Nonuniform Meshes 64 20. Results for Singular Perturbation Problem 4 on Nonuniform Meshes 65 21. Results for Singular Perturbation Problems 1,2,3,4 with h « e 66 iv Acknowledgement Thanks to Uri for all his patience. V 1. INTRODUCTION One of the major problems encountered in the numerical solution of boundary-value problems by finite difference methods is the selection of a mesh on which the approximate solutions can be represented with the desired accuracy. Nonuniform meshes typically are used for problems which have solutions with large variations in gradient. Mesh points are chosen close together in regions where the gradient is large, and further apart where less detail is needed for an accurate approximation. If only uniform meshes are used for such problems, the efficiency and stability of the difference scheme can be greatly affected. However, the accuracy of the solution (as determined by the usual error bounds) is sometimes of lower order on a nonuniform mesh than on a uniform mesh. Thus, for the cases where a nonuniform mesh is more desirable, it is beneficial to derive schemes that have the higher "uniform" order of accuracy on the nonuniform meshes. Such schemes are sometimes called supra-convergent schemes be-cause their order of convergence (that is, the order of the error in the solution) is higher than the apparent order of their local truncation error. In Chapter 2 of this thesis a three-point method is derived to solve a second-order scalar ordinary differential equa-tion on a nonuniform mesh. The discretization method is shown to be second-order and stable. Stability bounds for the method are also found. The stability bounds are very important for predicting the size and growth of the error. The other chapters of this thesis extend the method to a wider range of differential equations. In Chapter 3 1 the method is extended to solve a two-dimensional elliptic differential equation with Dirichlet boundary conditions. Singularly perturbed problems are discussed in chapter 4. Chapter 5 suggests methods that can be used for higher order problems, nonlinear problems, and more complicated boundary conditions. 1 . 1 Definition of Terms Consider a linear ordinary differential equation y" = {r{x)yy + q{x)y + f{x) (1.1a) on the interval a < x < b with Dirichlet boundary conditions y{a) = 0, y{b) = 0 (1.16) where y is a continuous function of x; r is a given continuous function; and q and / are given piecewise continuous functions. The independent variable, x, can be discretized over the interval [a, b] to create a mesh irn with n sub-intervals (or n+1 mesh points including the end points of the interval) so that x\ = a, x n + i = b TTn = {xi,X2,X3, ...Xn,Xn+1} where X{+i = X{ + hi, hi > 0, i=l,..,n. hi is called the step size of the i t h sub-interval. The mesh is called uniform if the mesh points are equally distributed throughout the interval, i.e. hi = h = , i = 1,.., n Once the mesh has been chosen, a mesh function y. Vn = {yi, V2, V3, yn+i}T is found to approximate the solution y(x) at each of the mesh points i = 1,.., n + 1. Discretization schemes, supplemented by the obvious discretization of the Dirichlet boundary conditions (1 .1c) , are used to determine the values of the mesh function. The error, et = y,-— y(xi), «' = l,..,n + 1, in these calculated values has two components : roundoff error and discretization error. The roundoff error is a product of using floating point arithmetic with a finite machine precision. In the following we shall assume that this error is negligible, so et- is the discretization error, that is, the error produced by the discretization scheme assuming that there is no roundoff error. The local truncation error, is the maximum amount by which the exact solution fails to satisfy the discretization scheme. It depends on the spacing of points in the mesh; thus the local truncation error (and consequently the discretization error) are measured in terms of h, the maximum step size on the mesh. The discretization scheme has first-order accuracy if Tn = O(h), and second-order accuracy if Th = 0(h2). The discretization scheme is called stable if small perturbations in the data produce only small changes in the calculated solutions. For a stable discretization, the usual bounds on the error in the approximate solutions are determined by the local truncation error, Th, and the stability constant, K The constant K is independant of h for all meshes with h sufficiently small. Thus, for a stable scheme, the orde of accuracy is also the order of convergence. Many centered second-order difference schemes will give second-order approximations to the solutions on a uniform mesh, but only first-order approximations on nonuniform meshes. i = 1,.., n + 1. 3 Often differential equations are converted to systems of first-order equations before they are discretized. (This is done in the derivation of the schemes discussed in this thesis.) On a uniform mesh discretizations of such systems can frequently be rewritten as three-point discretization formulas a; y t-_j + bi yi + ct- y t - + 1 = #t, i = 2,.., n, where the a^'s, 6 t 's, c;'s, g^s are constants determined by the discretization; from the three-point formula an n x n system of equations can be set up to solve for the mesh function yw Ay, = g. (1.2) The tridiagonal matrix A has the 6t- 's forming the diagonal, the a; 's the lower diagonal, and the Ci 's the upper diagonal. The right hand side g is the vector g = {0, 92, gs,-,gn, #}T-Because of the high efficiency available for solving (1.2), the three-point form is highly desirable. However, on a nonuniform mesh it is generally not possible to transform the discretized first-order system to a three-point formula. One goal of this thesis is to establish a procedure for producing second-order three-point discretization formulas for linear second-order differential equations on a nonuniform mesh. 1 . 2 P r e v i o u s a n d R e l a t e d W o r k Many attempts have been made to derive schemes that are second-order accurate on nonuniform meshes. Hoffman [3] tried to achieve second-order local truncation errors by using coordinate transformations to map nonuniform meshes onto uniform 4 meshes. He found, though, that if the local truncation errors of the centered finite difference schemes are second-order on the uniform mesh, the discretization error on the nonuniform mesh is second-order only if the coordinate transformation to the nonuniform mesh is a polynomial of second degree or less. This is a very restrictive condition. Chong [2] and others have used a variable mesh system to produce second-order accuracy on a nonuniform mesh. However, these methods result in complex finite-difference schemes to be applied on a very restricted type of mesh; the methods do not readily extend to higher dimensions. Tikhonov and Samarskii [12] tried a different approach. They applied certain second-order difference schemes to a class of second-order differential equations (r(x) y')' - q{x) y + f{x) = 0 with Dirichlet boundary conditions. By looking at the Green's functions for the prob-lem, they were able to show second-order convergence of the difference schemes on nonuniform meshes. They also found stability bounds for the differential equation. But their results were only valid under the conditions that r(x) > 0, q(x) > 0. In 1969, Keller [4] proposed stable methods for solving first-order systems of equations through discretizations on nonuniform meshes. This method can be used on equation (1.1) if the equation is converted to a first-order system, y' = Q(x) y + f (x). In this method the midpoint scheme is applied to each equation in the first-order system. The discretized system is used to create the banded matrix system BIYJI = b! (1.3) with y n ={yi, y 2 , " , y n + i } T b i ={6, Ku K2,...,Kn, #}T 5 and the matrix Bi has the form (Ba Si Ri Bi S2 i?2 V Sn Rn BbJ where the SYs and R^s are 2 x 2 matrices that are determined by the equation Si Yi + Ri Y i+ i = K{, i = 1,.., n, and Ba and Bb represent the boundary conditions. On an arbitrary nonuniform mesh the system (1.3) generally yields solutions that are second-order accurate. However, Keller's work does not readily extend to singularly perturbed problems, and (if the system is a conversion of (1.1a) to first order) the resulting linear system is not tridi-agonal. Manteuffel and White [8,9] solved the latter problem. They manipulated the system (1.3) for equation (1.1a) to give a compact-as-possible discretization scheme by reducing the system to the tridiagonal one A yw = g. They assumed stability to show that this scheme retains the second-order accuracy. On the other hand, Keller's results apply to much more general problems than those resulting from converting (1.1a) to a first-order system. First-order systems representing multipoint boundary-value prob-lems and initial-value problems with linear and integral constraints can all be solved using this simple scheme. He also claims, but does not show, that his results extend to problems with nonlinear boundary conditions. In a later paper with White [5], they show that one-step centered difference schemes besides the midpoint scheme can be used to give similar results. They also show that with linearization (Newton's method in particular) this method can be used to solve for isolated solutions of nonlinear prob-lems. This nonlinear work was expanded on by Lentini and Pereyra [7], who applied deferred corrections to Keller's scheme to produce a self-adapting, variable-order finite difference scheme for nonlinear multipoint boundary-value problems. 6 A l l o f these schemes were deve loped as ve r t ex -cen te red schemes, m e a n i n g t ha t u n k n o w n s are eva lua t ed o n l y at m e s h p o i n t s . M a n t e u f f e l a n d W h i t e [8] d i scuss c e l l -cen te red schemes for (1.1a), whe re u n k n o w n s are o n l y eva lua t ed at the p o i n t s m i d - w a y be tween the m e s h p o i n t s . T h e y show t h a t , even t h o u g h t h e i r scheme is incons i s t en t ( tha t i s , the a p p r o x i m a t i o n to y" is i ncons i s t en t ) , i t s t i l l y i e ld s a t r i d i a g o n a l s y s t e m w i t h second-orde r a c c u r a c y o n a n o n u n i f o r m m e s h for l i nea r b o u n d a r y - v a l u e p r o b l e m s w i t h s epa ra t ed l inear b o u n d a r y c o n d i t i o n s . A g a i n t hey assume bu t do no t show s ta -b i l i t y . T h i s thesis gives a p rocedu re to p r o d u c e s table compac t - a s -poss ib l e difference schemes w i t h second-order accu racy o n n o n u n i f o r m meshes w i t h o u t the l i m i t a t i o n s present i n some of the ear l ier w o r k . T h e schemes are ne i the r ve r t ex -cen te red or c e l l -cen te red , b u t a c o m b i n a t i o n of b o t h . T h e dif ferent ia l e q u a t i o n (1.1a) does no t need to be c o n v e r t e d t o a first-order s y s t e m i n order for these schemes to be a p p l i e d , even t h o u g h s u c h a conver s ion is used for d e r i v i n g the schemes. 7 2 . A S E C O N D - O R D E R T H R E E - P O I N T S C H E M E 2 . 1 D e r i v i n g a S e c o n d - O r d e r T h r e e - P o i n t S c h e m e T h e e q u a t i o n u n d e r c o n s i d e r a t i o n i n one d i m e n s i o n is the e q u a t i o n (1.1) y" =(r{x)y)' + q{x)y + f{x), y(o) = 9, y{b) = 0 , (1.1) o n the i n t e rva l a < x < b . T h i s e q u a t i o n is r e w r i t t e n as the f i r s t -order s y s t e m of O D E ' s y' = r{x)y + z, (2.1a) z' = q{x)y + f{x). (2.16) T h e t r a n s f o r m a t i o n to the s y s t e m (2.1) has some advantages over the m o r e c o m m o n t r a n s f o r m a t i o n ')'•(•0 MCW0)- <"> y'J \q(x) + r'(x) r(x)J \y'J \f(x)J F o r the s y s t e m (2.1) , r'(x) does no t need to be eva lua t ed , whereas to solve the s y s t e m (2.2) i t is necessary to evalua te r'(x). T h i s f u n c t i o n is no t a lways e x p l i c i t l y ava i l ab le . 8 Also for some applications, (r(x) y)' may be smooth while y' is not. Since (2.1) does not separate r{x) y and does not require evaluation of the potentially non-smooth function y', the system (2.1) can be used in such cases but the system (2.2) cannot. Here the transformation to (2.1) is done partly to facilitate the discussion of singularly perturbed problems in Chapter 4. The method described in this section will work equally well for any transformation of the equation (1.1) to a first-order system of equations. Typically the system (2.1) is solved by discretizing both equations on a mesh nn. Difference schemes are applied separately to these two equations to find approximate values of y and z at the mesh points. (See Figure 1) Figure 1 Standard Evaluation Points on a Mesh Zi-l Z{ Zi+i Vi-i Vi J/t+i 1 : 1 * Xi—\ X{ %i+l As Keller shows in his paper [4], the system B i y n = b x (1.3) results. This system is not tridiagonal in general. Now let us introduce the concept of a staggered mesh. "Staggered mesh" means that equations (2.1a) and (2.1b) are not discretized at the same points on the mesh. The values of y are approximated at the meshpoints , yi ~ y{xi) i = 1, ..,n+ 1, 9 but the values of z are approximated at the points half-way between the mesh points , Zi+l/2 ~ z{Xi+l/2) Xi+1/2 = %i + y i= l , . . , n (See Figure 2 below). F i g u r e 2 Evaluation Points on a Staggered Mesh Zi-l/2 z t '+ l /2 3 / t - i y t + i ^ i - i X i + i On a staggered mesh, the discretization of the system of first-order ordinary dif-ferential equations (2.1) by one-step centered difference schemes can be represented by the banded matrix system Bvn = b (2.3) where V , r ={2/1, ^1 + 1/2, 2/2, 2 2 + 1 / 2 , — , 2 n + l / 2 , Vn+l}T b={6, 0, / 2 , 0, / 3 , . . . , / n , 0, #} T and the matrix B has the form 10 B = ( 1 Si - 1 *1 " 2 1 ~h2 - 1 t2 U3 -1 ^3 V3 _1_ h3 Vn 7 - Wn V 1 tn 1 The s t 's and t\s are determined by the equation * i Vi + ti y i + i - zi+1/2 =0 , i = 1, . . ,n, the U i ' s , v,-'s, and w^s by - 1 1 , Ui yi_! + ^ Z i _ i / 2 + v{ y,- + ^ - 2 l + 1 / 2 + Wi y i + 1 = i = 2,..,n hi hi and hi = \(hi + / i ; _ i ) . Because the discretization was performed on a staggered mesh, the system (2.3) can be reduced to a tridiagonal system A yn = g. (1.2) The i row of the tridiagonal system, i = 2,.., n, is created by adding from (2.3) J-row[2t] + row[2i - 1] - J - row[2^ - 2 ] . (2.4) h2i h2i This operation eliminates the auxiliary unknowns by elementary row operations (thus the operation is stable) and leaves an equation in only the principal unknowns t/,-_i, y t-, yi+i- The boundary conditions, as represented in (2.3) form the 1 s t and n + 1st rows of the system (1.2). The tridiagonal system is a matrix representation of a three-point formula a ; Vi-i + bi yl + c t y l + i = gi i = 2,..,n. 11 In order for the resulting three-point formula to have second-order accuracy on a staggered mesh with the development shown here, any difference scheme used must be one-step and centered, (and hence second-order accurate). The reasons for this will be discussed in the next section. For example if the midpoint scheme is used, the discretization for equation (2.1a) on [x{, xi+i] , i = l , . . ,n is y%+1h. V% = rixi+i/2) y{xi+i/2) + z{xi+1/2) (2.5a) where Z t + i / 2 = Xi + a s a D o v e - Since r(x) is a given function, r(xi+i/2) c a n D e evaluated exactly. The value for y ( £ t + i / 2 ) will be approximated by |(j/i+i + yi) , and the value for z(xi+1/2) will D e approximated by Zi+i/2- So (2.5a) becomes y t + \ . V l = \ r ( 3 ; t + i / 2 ) {Vi+i + yi) + Zi+i/2- (2.6a) Equation (2.1b) can also be discretized using the midpoint scheme. But in-stead of being discretized on [x , - ,Xi+i ] , equation (2.1b) is discretized on the interval \xi-i/2> xi+i/2], i = 2,..,ra . This interval has length hi, and midpoint, called xt-, Xi = Xi + - (hi - hi-i). The midpoint scheme applied to (2.1b) therefore yields 'i+l/2 ~ Zi-l/2 = q { i i ) y [ i i ) + / ( - ,} { 2 5 b ) hi The known functions q(x) and f(x) are easily evaluated at x;. A n approximation y t for y(x t) is obtained by quadratic interpolation on yi, y;+i. (For details on this calculation, see, for example [1, 2, 14].) This yields y{ = a.i y t _ ! + & y{ + q{ yi+1, where (3/i I + / i i _ 1 ) ( / i i _ 1 -hi) Oil = Pi IQhi^^hi + / i ; _ i ) (3hj + hj-i)(3hj-i + hj) IShi-yhi (3hi_i + hi)(hi — hi_i) 16hi(h{ + K^i) 12 Note that quadratic interpolation has second-order accuracy. Equation (2.5b) then becomes, 2 t + 1 / 2 ~ — = q{xi) (a» Vi-i + Pi Vi + li Vi+i) + f{*i)- (2-6b) hi We can use the difference formulas (2.6a) and (2.6b) directly to develop a three-point formula. Then the system (1.2) can be formed so that, ultimately, approximate values of y at the mesh points can be determined. We use equation (2.6a) as a definition Of zi+l/2 Zi+1/2 = V l + \ . V l ~ 2 r(x*'+i/2) (JM-i + Vi)- (2.7) Substituting this definition in for Zi+1/2, Zi-\/2 m equation (2.6b) to eliminate all occurrences of z leaves an equation in only y. h j - i j y j + i - yj) - h j ( y j - y , - i ) | l/2(hi + h i ~ i ) h i h i ^ i r { x i - i / 2 ) (yi + Vi-i) ~ r(xj+1/2)(yi+1 + yt-) _ hi + hi-1 q ( i i ) (cti y t-_! + p\ y{ + n yi+1) + f(x), (2.8) i = 2,.., n. Equation (2.8) can be rearranged to give the three-point formula ai y t _ i + bi y{ + c, y i + 1 - gu i = 2,.., n, a . = 2 + r(xj_1/2) _ ^ h i _ i ( h i + h i - i ) hi + h i - i i - 2 r(xj_1/2) -r(xi+1/2) b i = ~h~h— + i m — q(Xi> tiihi-x hi + hi-1 2 r ( z i + 1 / 2 ) C l " htihi: + hi.1) + ht+hi^ q { X l ) 1 " 9i = f{xi). (2.9) For a uniform mesh, the interval [xi_i/ 2 ,2^+1/2] has length h and midpoint X{. So then discretization of (2.1b) by the midpoint scheme gives zi+l/2 _ zi-l/2 { \ , ft \ (o ir\\ 7 — = q{Xi) yi + f ( X i ) . (2.10) h 13 S u b s t i t u t i n g (2.7) i n to (2.10) y i e l d s the " u s u a l " t h ree -po in t f o r m u l a Q>i Vi-i + bi yi + c» y , - + i = gi, i = 2 , . . , n, where _ 1 r ( g t - _ 1 / 2 ) a i ~ h * + 2h ' , - 2 , . , r ( g t - _ 1 / 2 ) - r ( z t - + 1 / 2 ) 6i = ^ - « ( * 0 + Y h , _ J L _ R (*•••+1/2) C * /* 2 2/* ' N o t e t h a t e q u a t i o n (2.11) is the same as e q u a t i o n (2.9) w i t h on = 7; = 0, a n d Bi = 1. T h e t h ree -po in t fo rmulas (2.9) a n d (2.11) were d e r i v e d u s i n g the m i d p o i n t scheme. O t h e r difference schemes c a n also be used . If the t r a p e z o i d a l scheme were u sed o n (2 .1b) , t h e n i n s t ead of e q u a t i o n (2 .5b) , the d i s c r e t i z a t o n o f (2.1b) w o u l d be zi+l/2 ~ zi-l/2 _ 1/2(^ + ^ -1) ~ - 9(^-1/2) (yi + yi-i) + g / ( x i - i / 2 ) + ^ 9(^+1/2) {Vi+i +Vi) + \ f{xi+1/2). (2.12) S u b s t i t u t i n g (2.7) i n to (2.12) t o e l i m i n a t e z f r o m the e q u a t i o n leads t o the th ree -po in t f o r m u l a where 4 2r (x t -_! / 2 ) 1 a; + T — T — - 0 9 ( ^ i - i/2), ^ t -_i(/ i , : + hi-i) hi + h^! 2 _ - 4 4 2 r ( a t - _ 1 / 2 ) - 2 r ( x t + 1 / 2 ) hi-i(hi + hi-i) hi(hi + hi-i) / i ,-+/i t-_i - ^ 9(^-1/2) - ^ 9(^+1/2), 4 2r{xi+1/2) 1 " - 9 ( ^ i + l / 2 j , hi(hi + hi-i) hi + hi-i 2 9i = /(xt+1/2) + /(^i-1/2) (2.13), 14 i = 2,.., n. Finite element and quadrature approaches can also be used. For example, an implementation of the Ritz-Galerkin method using roof functions as basis functions yields a scheme that is related to the trapezoidal scheme : / ' [q{x) y{x)+f{x)) dx+ f ^\q{x) y{x) + f{x)) dx. (2.14) When the two integrals in (2.14) are replaced by simple midpoint quadrature formulas, the equation becomes z i + l/2 — zi-l/2 — ^hi-i q{xi_l/2) {yi+ y t-_i) + ^/i t_i /(x t-_i/ 2) + hii q{xi+l/2) (yt- + yi+i) + f(xi+1/2). (2.15) This yields the three-point formula where 1 r ( g i - i/2) hi-i 9(^-1/2) a 1 2 4 7 - 1 1 r(xj_1/2) - r(xi+1/2) hj-i g{xj_1/2) + h{ g{xi+1/2) f"~fc,--i ht+ 2 4 - _ L _L r ( Z « + l/2) _ fet g ( ^ t - l/2) C i ~ f c i 2 4 ffi = ^ ( / ( ^ ' + l/2) + / (^ i-1/2)), i = 2,..,n. (2.16) Note that if the three-point formulae (2.13) or (2.16) are applied on a uniform mesh, they do not reduce to the usual three-point formula (2.11), but (2.13) and (2.16) are still expected to be second-order accurate. 15 2 . 2 D e t e r m i n i n g t h e O r d e r o f t h e E r r o r How is the order of the discretization error of these three-point formulae deter-mined? First the local truncation error over the ith sub-interval, r^., is found for equation (2.7). This is done using the Taylor expansions for j/(x t+i) and y(x t) about the point X j + 1 / 2 . (These expansions hold for all i — 1, ..,n ). h h2 h3 y{xi+1) =y{xi+1/2) + -^y'{xi+1/2) + -£-y"{xi+l/2) + -£y'"{xi+1/2) + 0{h{4), h- h2 h3 y{Xi) =y{xi+1/2) - jy'{xi+1/2) + ^ -y"{xi+1/2) - ^ y " ' ( x i + 1 / 2 ) + 0{h{4). Then it is easily determined that h3 y{xi+l) - y(Xi) =ht y'{xi+1/2) + ^ Ly»'(Xi+1/2) + 0{hi4), (2.17) h2 y{xi+1) + y{xi) =2 y{xi+1/2) + -^y"{xi+l/2) + 0{h{4), (2.18) i = 1,.., n. Substituting the expansions (2.17) and (2.18) into equation (2.7) gives an expression for the local truncation error in equation (2.7). h3 rhi = - z i + 1 / 2 + y'{xi+1/2) + ^ y"'{xi+1/2) + 0(h4)-1 h2 - r{xi+1/2) [2y{xi+1/2) + ~^y"{xi+1/2) + 0(/i t- 4)), (2.19) % = 1,.., n. But y' — r y + z according to equation (2.1a), or y'{xi+1/2) = r(xi+1/2) y{xi+1/2) + z i + 1 / 2 . (2.20) By replacing y ' ( x t + l / 2 ) with its equivalent (2.20), the local truncation error over the ith sub-interval for equation (2.7) can be simplified to rht = -£y"(zi+1/2) r(xl+y2) + f^y'''(xi+1/2) + 0(hi4).y''(xl+1/2) r ( x i + 1 / 2 ) + 0(/ i l - 4 ). (2.21) 16 If h is denned to be the maximum step size on the interval [a,b], h = max (h{), l<i<n then, over the interval [a,b], the maximum local truncation error of equation (2.7) obviously satisfies rh = max \rh. \ - 0(h2). l < i < n The local truncation error for (2.6b) can be found using the Taylor expansions for zi+\/2-> Zi-i/2 about the point ii. As before, let hi = \(hi + / i t-_i), the length of the interval [xi^1/2,Xi+1/2], and let 1; = z(xi). (Note that h = max hi = 0(h) .) 2 < i < n z (x l + 1 / 2 ) = ~zt+h±z>+1- ( £ ) % + \(%) %' + 0{hi), (2.22) 2 / - \ 3 2 1 2 \ 2 1 6 2 <*-V2) = « - J % + \ (j) ~z"-l[j) + ( 2- 23) Now substitute the exact expansions into (2.6b) and include the 0(h2) error for quadratic interpolation to get rhi = ~z[ + ±-A h2 ~z'l' + 0(hi4) - i q(ii) & + 0(ht2)) - f(ii). (2.24) But if this expression is simplified by using equation (2.1b) as a definition for Zi, the local truncation error for (2.6b) is found to be rhi = ^ h2 z'l' + 0(hi4) + q(ii) 0(h{2). (2.25) Therefore, both equations (2.6a) and (2.6b) have second-order local truncation error. The three-point formula (2.9) derived using these equations (2.6a) and (2.6b) will have a second-order discretization error on an arbitrary mesh even though the local truncation error of the formula is only first order. (The local truncation error is found by substituting Taylor expansions for t/t-+i, y^ , y t -_ i , and y^ about the point x t + 1 / 2 into equation (2.8).) The same will hold for any three-point formula developed as outlined 17 above using one-step, centered, second-order difference schemes on a staggered mesh : the discretization errors will be second-order. A more generalized look at the errors in the system will show why. Recall that the discretization can be represented by the system B yw = b. (2.3) Since each of the difference equations represented in (2.3) has second-order local trun-cation errors, the system (2.3) as a whole will also have second-order local truncation errors. To show stability of (2.3) we first demonstrate that, after appropriate scaling, B resembles a multiple-shooting matrix (-Y{x2,xx) I \ -Y(x3,x2) I M V - y ( x n , x n _ i ) / BbJ (2.26) Ba = where the Y(XJ,X{) are fundamental solutions of the system (2.1) satisfying F ( x t , x t ) = J. Thus these are 2x2 matrices that solve the homogeneous form of (2.1) fr(x) 1\ y ' ( x , *,•)=• )Y{x,xi). \q(x) OJ The matrices Ba and BD represent the Dirichlet boundary conditions 1 0 \ / 0 0 \ Bb=( 0 0 / \ i 0 / First note that, for the multiple shooting matrix M , both dependent variables in the system must be evaluated at the mesh points. Presently in (2.3) y is evaluated at the mesh points but z is not; however the evaluation of z at the midpoints of the mesh in (2.3) can be viewed as an approximation to the values of z at the mesh points. By Taylor series we know that z ( x i + 1 / 2 ) = z(x{) + 0(h). z i + i / 2 , an approximation to z(xi+1/2), can therefore be considered a first-order approximation to z(Xi), and 18 consequently the system (2.3) can be seen to evaluate both y and z at the mesh points with first-order accuracy. Next we see that the elemental structure of B is different from that of M . M has (2n+2) rows and columns whereas B has only (2n+l). We can extend B to the "proper" size by adding the equation zn+3/2 ~ zn + l/2 \{hn+i + hn) q{xn+i) yn+i + f(xn+i) to the system. zn+3/2 w i U g i y e a first-order approximation to z(b). If we also move the first row of the system (2.3) so that the boundary conditions are represented in the last two rows of the system, and then multiply on the left by f1 \ 1 t2 I -w3 1 t3 V l j then the correct zero structure of M is attained. The system that results from these transformations is By, = Dib (2.27) where we denote V«- = { j / l , Zl, J/2, Z2,.., yn, zn, Vn+l, zn + l } T ' , with b = {o, i2, o, / 3,-, in, o, f{xn+1), e, i?}T, 19 and B = - 1 U2 T -«2 ^2 fco ^ *2 - 1 -1 7 1 + 1 <l[xn+l) 'n+1 0 o that f. 0 V o o i o J Note that multiplication by Di is the same as solving for in equation (2.6a) and then substitutimg into equation (2.6b). A look at the local truncation errors shows O(hi), so Di = I + 0(h). Thus, for h small enough, the inverses of B and B can be bounded by similar constants. From here we proceed as in Theorem 5.38 in [1]. To simplify the rest of the discussion we write the system (2.27) in block form using 2x2 matrices S; and Ri (Sx R1 \ S2 R2 B = Rn \Ba BbJ B can be equilibrated by multiplying by the matrix D2 = diag(R-\ R2\..,R-\ I) to form the matrix (R-1S1 B = D2B R2 S2 Rn1Sn I \ Ba BbJ We now show that the block entries R^~1Si approximate the fundamental solution matrices Y(x;+i>xi)- As discussed earlier, when the scheme represented by (2.3) 20 is a p p l i e d to the f u n d a m e n t a l s o l u t i o n s , it c a n be cons ide red a first-order difference scheme a n d so y i e l d s Si Y{xi,Xi) + Ri Y{xi+Uxi) = 0{hi) (2.28) i — 1, ..n. Since H-R^ Hoo = 0(h{), the m u l t i p l i c a t i o n of (2.27) by D2, w h i c h co r responds to m u l t i p l y i n g (2.28) by i 2 t - 1 , gives R^Si + Y{xi+1,xi) = 0{h2). (2.29) Seen i n th i s l igh t i t becomes o b v i o u s t h a t the m a t r i x B is r o w - w i s e a second-order a p p r o x i m a t i o n to the f u n d a m e n t a l s o l u t i o n m a t r i x , M , a n d the e r rors i n (2.29) c a n be represen ted i n b l o c k d i a g o n a l f o r m by E = diag[E1, E2,..., En, 0) where the last e n t r y co r re sponds to the er rors i n the a p p r o x i m a t i o n o f the D i r i c h l e t b o u n d a r y c o n d i t i o n s . N o t e t h a t \Ei\\ = 0{h 2). T h u s we w r i t e B = M + E, or [M + E)\w = /J2£>ib. (2.30) N o w t o d e m o n s t r a t e s t a b i l i t y o f (2.3) , we find a su i t ab l e b o u n d o n B 1 (a b o u n d o n B"1). F i r s t r eca l l t h a t the inverse o f M is g i v e n by / G{x1,x2) ... G(xuxn+l) Y(x1,x1)Q~l \ G(X2,X2) ... G(X2,Xn+l) Y[x2,x1)Q~1 M = \G(xn+1,x2) ... G(xn+1,xn+1) Y(xn+l,xl)Q~l J where Q is a n o n s i n g u l a r m a t r i x r epresen t ing the f u n d a m e n t a l so lu t ions at the b o u n d -a r y c o n d i t i o n s Q = BaY(x1,x1) + BbY(xn+1,xi) 21 and G(xi,Xj) is the Green's function of the problem (see for example [l]). Hence a bound on the norm of M~l is given by H M - 1 ^ < n K where /c, the conditioning constant of the problem, gives a measure of the sensitivity of the boundary-value problem to perturbations in the data. Then it is easy to see that H M ^ D a l l o o < K + 0(h) and | |M - 1 . f? | | oo < c/cfo, for some constant c. If c/c/i < 1 then Wil + M^E)-1^ < l + 0{h). Now B-1 = (I + M~1E)~1M~1D2. So for any mesh with h sufficiently small, l l ^ - 1 ! ! = K + o{h) and a similar bound holds for B~x. Therefore, (2.3) represents a stable difference scheme with a stability constant K = K + 0(h). We already know that (2.3) has second-order local truncation errors when applied to the problem (1.1) on a staggered mesh. Thus, by stability and consistency arguments, the solutions obtained using the difference methods derived in Section 2.1 and represented by the system (2.3) are second-order approximations to the exact solutions. Since the system (1.2) is just another representation of the same method, it also yields solutions with second-order discretization errors. 22 2 . 3 N u m e r i c a l R e s u l t s Numerical results were calculated for six test boundary-value problems with Dirich-let boundary conditions. The exact solution is known for each boundary-value problem. 1) v" = \ y ' - ^ - y - l e " 3 1 + 4e~2* +1 ( i + * ) e - * , y(0) = 2, y(2) = e - 4 + e- 2 , solution : y = e~ 2x + e~ x. 2) x2y" = -xy' + 4y, y{l) = 2, y(2)=4.25, 2 1 solution : y = x H — - . x 3) x2y" = xy' + 3y - 9x 2 - 8x - 3, y(l) = 7, y(5)=211, solution : y = xz + 3x 2 + 2x + 1. 4) y" = A 2 y + (A 2 + 4TT2) cos 2(7rx) - 2TT 2, y(0) = 0, y(l) = 0, . . exp(A 2x — A 2 ) + exp(—A2x) 2 , * solution : y = ; - r : COS ( 7 T X ) . 1 + exp(-A^) °) y" — — 7 r tan(7rx) y', y(0) = 0, y(0.25) = y/2/2, solution : y = sin(7rx). y ( l ) = 0 , y (2)=2 ln(1.75), 7 solution : y = 2 ln( - ) . 8 — x z The solution to each problem was approximated 30 times on each of eight different types of meshes. Each time a different number of mesh points was used. A random 23 number generator was used to pick the number of mesh points between 9 and 300. The points on the meshes were calculated using formulas relating to transcendental functions. The first mesh type was created according to a sine curve. Points are taken at even sub-intervals along the sine curve over the interval [0,7r/2] : • / n t \ xi = sm\Wi—rT\ J ' t = 0 , . . ,n + l \2{n + l)J where n+1 is the number of mesh sub-intervals. (See Figure 3 below.) The mesh is then scaled and translated to fit the interval of the boundary-value problem 24 under consideration. The resulting formula is Xi = a + sin ( — ] ( 0 _ a ) « = 0,... n + 1 V2(n + l ) y K h where a and b are the endpoints of the interval. Meshes of other types were similarly created according to the cos curve, the ln curve, the l o g 1 0 curve, and linear combinations of the above. The meshes were created in this way to ensure non-uniformity in the meshes, and to allow calculated variations in the amount of non-uniformity. It also simplified the creation of a mesh with a random number of mesh points. For each problem on each mesh, the approximate solution at each mesh point was compared to the exact solution at the mesh points. The pointwise error was recorded. From these local errors an average error and a maximum error were calculated. The average error is calculated to ensure that the maximum error is not unrepresentative of the local errors. The traditional "sub-divide the mesh" method for experimentally determining the order of accuracy of discretization schemes involves first the repeated sub-division of each mesh subinterval, starting from an arbitrary first mesh, and then a comparison of the esulting solutions. This method is inappropriate for our purposes, however, because the refined meshes are not arbitrary anymore. So, in the experiments reported in this thesis, the order of accuracy of the approximate solution is determined using the more statistical method described in Manteuffel and-White [8]. A log-log plot is created of the maximum local error against the maximum step size for each problem on each mesh. Figure 4 shows the first-order errors obtained when the usual three-point formula (2.11) is applied to problems 1) to 6) on the nonuniform meshes. Figures 5 to 10 show the results for problems 1) to,6), respectively, when the three-point formula (2.9) derived from the midpoint scheme is used. The solid line appearing is the linear least-squares fit to the plotted points. The slope of this line gives the order of accuracy of the method used. Figure 11 shows the combined results for problems l ) to 6). The accuracy of the three-point formula (2.9) is determined to be second-order. Results are also given for the three-point trapezoidal formula (2.13). Figure 12 shows 25 that this method is also second-order accurate. These results are for problems 1) to 6) together. 26 Figure 4 Usual Three-Point Formula : Results for Problems 1 to 6 slope = 0.98 o 27 F i g u r e 5 T h r e e - p o i n t M i d p o i n t F o r m u l a : R e s u l t s for P r o b l e m 1 28 Figure 6 Three-point Midpoint Formula : Results for Problem 2 Figure 7 Three-point Midpoint Formula : Results for Problem 3 slope = 2.01 30 Figure 8 Three-point Midpoint Formula : Results for Problem 3] Figure 9 Three-point Midpoint Formula : Results for Problem 5 32 Figure 10 Three-point Midpoint Formula : Results for Problem 6 slope = 2.43 33 Figure 11 Three-point Midpoint Combined Results for Problems 1 to 6 34 Figure 12 Three-point Trapezoidal Combined Results for Problems 1 to 6 3 . A N E X T E N S I O N T O T W O D I M E N S I O N S T h e goa l here is t o deve lop a n i n e - p o i n t f o r m u l a w i t h second-order a c c u r a c y to solve a d i f fe ren t ia l e q u a t i o n posed i n t w o d i m e n s i o n s , o n a n a r b i t r a r i l y spaced r e c t a n g u l a r m e s h by u s i n g the same p rocedures as were used for the o n e - d i m e n s i o n a l case i n C h a p t e r 2. A l t h o u g h a five-point f o r m u l a w o u l d have been preferrable f r o m some pe r spec t ives , the d e r i v a t i o n of such a f o r m u l a does n o t seem poss ib le u s i n g the m e t h o d s d e s c r i b e d here. A l l of the p rocedures u sed to der ive a second-order th ree -po in t f o r m u l a i n one d i m e n s i o n c a n be ex t ended t o the d e r i v a t i o n o f a second-order n i n e - p o i n t f o r m u l a i n t w o d i m e n s i o n s . T h e first p r o b l e m at h a n d is t o define w h a t is m e a n t by second-order accuracy . A d e f i n i t i o n ana logous to the o n e - d i m e n s i o n a l case is used . Second-o rde r a c c u r a c y means t h a t the e r ro r has o rde r 0(h2), where h is the m a x i m u m step size o n the g r i d 7 r n j m i n e i ther the x - d i r e c t i o n o r the y - d i r e c t i o n . h = max(Hi, H2), where Hi is the m a x i m u m step size i n the x - d i r e c t i o n , a n d H2 is the m a x i m u m step size i n the y - d i r e c t i o n N o w , cons ide r the t w o - d i m e n s i o n a l e l l i p t i c d i f ferent ia l e q u a t i o n Uxx + Uyy = {A{x,y)U)x + (B{x,y)U)y + C(x,y)U + D{x,y) (3.1) h < x < Q u 0 2 < y < © 2 36 where A , B , C are given smooth functions, and D is a given piecewise smooth function. This equation is discretized over an arbitrary (non-uniform) rectangular grid TTn,m where 7Tn,m = {{xi,yj),i = l , . . , n + l , i = l , . . , m + 1} xi+i = %i + hi, yj+\ = yj+kj. h{,kj > 0 are the step sizes in the x- and y-directions, respectively. Note that any dicontinuous points in D(x,y) must be included as points in the grid 7Tn,m- On this type of grid, Uij is written to denote an approximate solution for U{xi,yj). Let us define R(x,y) and Q(x, y) by R = UX- A{x,y) U (3.2a) Q = Uy-B{x,y) U. (3.26) Then, equation (3.1) can be rewritten as Rx + Qy = C(x,y)U+ D[x,y). (3.2c) Together equations (3.2) form a first-order system of equations that is equivalent to (3.1). In one dimension, a second-order three-point formula was derived by using second-order difference schemes on a staggered mesh to solve both equations in the system of first-order ordinary differential equations. Similarly in two dimensions a nine-point formula is derived by solving the equations (3.2a), (3.2b), and (3.2c) by second-order difference schemes on a staggered mesh. Let the term 6ox[(x;, yy), (xt-+i,yj+i)} refer to the rectangular region on the grid with the four corners (xi,yy), (x ; , yy + i ) , ( x t + i , y y ) , ( x t + i , y y + i ) . In two-dimensions we create a staggered mesh by discretizing (3.2a) and (3.2b) over box[(x»,yy), ( x i + i , y y + i)], and (3.2c) over 37 box[(xi_ 1 / 2 ) 2/?-i/2)> {xi+i/2i Vj+i/v)]- The resulting formula will then be expected to have second-order accuracy. Hence to derive a second-order nine-point formula, the first step is to derive a second-order difference scheme to solve equations (3.2a), (3.2b), and (3.2c). A two-dimensional midpoint scheme is used here, although other schemes are possible. 3 . 1 D e r i v i n g a T w o D i m e n s i o n a l M i d p o i n t S c h e m e We start by considering Taylor expansions of an arbitrary continuous function f(x,y) over box[(zi, yy), (z^+i, yy+i)] . The midpoint scheme calls for evaluations of the equations under consideration at points midway between the mesh points. To derive a midpoint scheme, Taylor expansions are taken about the point ( x j + 1 / 2 , y y + i / 2 ) h' f{xi+l,Vj + l) = f(Xi+i/2,yj + l/2) + y fx{Zi+i/2,yj+l/2) + 1 / h \ / y(z» + i/2 ,yy+ i/2) + 2 ( y ) /«(z . -+ l/2 ,yy+l/2)+ h{ k- 1 f k \ 2 y y fxy{xi+i/2,yj+1/2) + ^ ( y ) fyy{xi+i/2,yy+1/2) + 0{h3) (3.3) h f{xi,yj+i) = f{xi+i/2,yj+1/2) - y fx(xi+1/2, yy+i/2) + y / y («T+i/2»yj+1/2) + 2 ( y ^ fxx{xi+1/2,yj+l/2) -h k 1 /k\2 y y / x S / ( x l + 1 / 2 , y y + 1 / 2 ) + - ( y ) f yy{xi+i/2,yj+1/2) + 0{h3) (3.4) f[xi,yj) = /(z ;+ i/2,yi+i/2) - y /i(a:i+i/2>yy+i/2) -38 y / y(Zt + i/2 , 3 / j + l/2) + g ("^) ^ ^ ( X i + l / 2 ' y ; + l/2) + / i i A; 1 (k \ 2 •f -f fxy{Xi+i/2,yj + l/2) + 2 ) / y y f o + i / 2 , yy+1/2) + 0 ( / i 3 ) (3.5) f{zi+i,yj) = / ( z t + i/2, yy+1/2) + y fx{xi+i/2,yj+i/2) kj . . l / / i i \ 2 , y / y ( z t + i / 2 , yy+1/2) + - I y I / z * ( z t + i / 2 , y ; + i / 2 ) -/ i - A; 1 (k\2 Y y / * y ( * i + i/2, yy+1/2) + 2 ( ) / y y ( * t +1/2, yy+1/2) + C*(/ i 3 ) (3.6) N o w s u b t r a c t i n g (3.4) f r o m (3.6) gives hi / x ( z i+ i/2 , yy+ i/2) - ^ / y ( x i + i / 2 , y y + i / 2 ) = f{xi+uyj) - f{xi,yj+1) + 0{h3) (3.7) a n d s u b t r a c t i n g (3.5) f r o m (3.3) gives hi fx{xi+i/2,yj+i/2) + kj fy{xi+1/2,yj+i/2) = f { x i + u y i + 1 ) - f{xl,yj) + 0{h3) (3.8) A d d i n g equa t ions (3.7) a n d (3.8) w i l l g ive an expres s ion for / x (£{+i/2,yy+1/2) fx{Xi+i/2,yj + l/2) = - j - [ / ( z i + i , y y ) + f(xi+1,yj+1) - f(xi,yj+1) - f(xi,yj)} + 0(h2) (3.9) 2hi S u b t r a c t i n g equa t ions (3.7) a n d (3.8) y i e ld s an expres s ion for fy(xi+1/2, yy+1/2) fy{xi+i/2, yy+1/2) = ^ - [ / ( x i + 1 , y y + 1 ) - f{xi+i,yj) + f{xi,yj+1) - / ( z i , y y ) ] + 0 ( / i 2 ) (3.10) E x p r e s s i o n s (3.9) a n d (3.10) c a n n o w be used as a second-order m i d p o i n t scheme i n t w o d i m e n s i o n s . 39 3 . 2 D e r i v i n g a S e c o n d O r d e r N i n e - P o i n t F o r m u l a Analogously to the one-dimensional case, equations (3.2a) and (3.2b) are dis-cretized at the mesh points, but equation (3.2c) is discretized at the points midway between the mesh points. This gives the staggered grid that is necessary for the nine-point formula to have second-order accuracy. So for the midpoint scheme over box[(x;,yy), (xi+i, yy+i)] equations (3.2a) and (3.2b) become •R(s»+i/2.yj+i/2) = Ux(xi+1/2,yj+i/2) ~ Mxi+i/2,yj+i/2)U{xi+l/2,yj+i/2) (3.11a) Q{Xi+l/2,yj+l/2) = Uy(xi+1/2,yj + l/2) - B(xi+1/2,yj + i/2)U(xi+1/2,yj + l/2) (3.116) i=:l,..,n j = l,..,m Equation (3.2c) is to be evaluated over box[(x t _ 1 / 2 ) yy-1/2)) {xi+i/2i yy+i/2)]- (See Fig-ure 13). In the resulting equation Rx{xi,yj) + Qy{xi,yj) = C(x;,yy) C/(x;, yy) -+- D(x;,yy) (3.11c) i = l,..,n j = l,..,m C(x,-, yy) and -D(x t, yy) can be evaluated exactly because C and D are known functions. For a uniform mesh (xi,yy) will be the point (x t-,yy), however for an arbitrary mesh ixi,yj) = {xi + hi/4 - hi-i/4 , yy + fcy/4 - fcy-i/4). 40 Figure 13 T h e T w o - D i m e n s i o n a l M e s h 1 ! i z i - i hi. xi - i h Xi+l i Vi-i U(xi,yj) is a p p r o x i m a t e d by Ui,j, f o u n d by b i q u a d r a t i c i n t e r p o l a t i o n . (See [14], for e x a m p l e , for de ta i l s o f the ca l cu l a t i ons . ) Uitj =ctj [ai Ui-ij-i + Pi Uij-i + 7i Ui+ij-i] + (3j [ai Ui-ltj + Pi Uij + n Ui+ij] + lj [cti Ui-ij+i + Pi Uij+i + 7t C/i+ij+i] where li (Shj + hi-i){hj_i — hi) 16hi-i(hi + hi-i) (Shj + /t1-_1)(3/tt-_i + /ij) 16/ii-j/ii (3/lj_! + /lj)(/lj - /^t-l) I6hi(hi + /ij_i) 41 _ (3kj + kj-i)(kj-i - kj) a j ~ 16kj-i[kj + fcy_i) ' 0. = (3fej + fcj-i)(3fc,-_i + kj) 3 Wkj-ikj ' _ (3/cy_i + kj)(kj - kj-i) l j ~~ 16kj{kj + ) ' t r ( x t + 1 / 2 , yy+1/2) is a p p r o x i m a t e d u s i n g the c o m m o n four -corner f o r m u l a U[zi+i/2,yy+1/2) ~ - [ ^ r ' + i , y + i + C^i.y+i + Ui+i,j + Ui,j]-W h e n U is used i n p lace of the a r b i t r a r y f u n c t i o n / i n express ion (3.9) , a f o r m u l a for Ux is o b t a i n e d , a n d e q u a t i o n (3.11a) is d i s c r e t i z e d as D . _Ui+i,j + Ui+i,j+i - Uj,j+i ~ Ui,j Ri+1/2,3 + 1/2 - ^ ^Mxi+i/2,yj+i/2){Ui,j + Ui+ij + Uij+i + U i + i i j + i ) (3.12a) i = l,..,n j ' = l , . . , m . S i m i l a r l y a d i s c r e t i z a t i o n is o b t a i n e d for (3.11b) by u s ing U i n p lace of / i n exp re s s ion (3.10) ^ _Uj,j+i + UJ+IJ+I - UJ+IJ - Ujj Vt+i/2,y+i/2 - ^ ^B{xi+1/2,yj+1/2){Ui>j + Ui+1>j'+Ui!j+1 + U i + l t j + i ) (3.126) i = l,..,n j = l,..,m. T o d i sc re t i ze (3.11c) , R is used as the a r b i t r a r y func t i on i n (3.9) a n d Q i n (3.10). E q u a t i o n (3.11c) is to be d i s c r e t i z e d over b o x [ ( x I _ 1 / 2 , y y - i/2)> (^+1/2, J /y+i/2)]? b u t express ions (3.9) a n d (3.10) were d e r i v e d for b o x [ ( x ; , y y ) , ( x t + i , y y + i ) ] - The re fo re the express ions (3.9) a n d (3.10) m u s t be sh i f ted to th i s new b o x to be used i n the d i s -c r e t i z a t i o n of (3.11c) fx(xi,yj)= (3.9a) [ / ( S j + l / 2 , y j - l / 2 ) + /(St+ l/2 , i/j + l/2) ~ /( S t - 1 / 2 . yy + 1/2) ~ / ( ^ - l / 2 , y y - l / 2 ) ] Q , F E 2 ^ 2{l/2hi + l/2hi-1) [ ' 42 fy{xi,yj)= (3.10a) + 0{h2). [f{Xi-l/2,Vj + l/2) + f{Xi+l/2,yj + l/2) ~ f{Xj+l/2,yj-l/2) ~ fjXj-1/2, Vj-l/2)} , 2(l/2fcj + l/2fcy_i) This gives the discretization {Rj+\/2,j + \/2 - Rj-l/2,j-l/2 - Rj-l/2,j + l/2 + Rj+l/2,j + l/2) h{ + hi-\ (<?t+l/2,j + l/2 ~ Qt-l/2,j-l/2 + Q i - l / 2 J + l/2 - Qt+l/2,y+l/2) _ &y ~h fcy — i C(£ t-,yy) ^ y + D ^ y y ) (3.12c) t' = 2,. . ,n, j ' = 2,. . ,m. Now (3.12a) can be used as a definition for -R;+i/2,y+i/2; a n d (3.12b) as a definition for Qi+i/ 2 )y+i/2- When these equations are substituted into equation (3.12c), the derivation of the second-order nine-point formula is complete. 43 Second Order 9-point Formula in Two Dimensions (3.13) -i.y-1 2 +hi - iA(x i _ 1 / 2,yi_i / 2 ) 2 +/c J _iB(x t _ 1 / 2,yi_i/ 2 ) 4hi-i(h{ + /ii_i) 4fcy_1(fcy + /cy-i) C(x;,yy) cti a3 + "4 + hi-iA(xi^1/2iyi-i/2) + hi-iA(xi_1/2,yi+i/2) 4/ii_i(/it-+/ij-_i) 2fcy + 2fcy_i - fcyfcy_i/J(xt_1/2,yl-i/2) + feyfcy_17J(xt_ 1 / 2 , yl + X/2) 4kjkj_i(kj + fcy_x) C(x;,yy) a t /?y] + 2 +/ i i _ 1 ^(x i _ 1 / 2 ,y i + 1 / 2 ) 2 - kjB{xi_1/2,yi+1/2) 4h{-i(hi + /j-i-i) C(xi,yy) a; 7y] 4kj(kj + A;y_i) + 44 Ui,j-1 -2/tj-! - 2hj - hihi_1A(xi+1/2,yi-i/2) + hihi_1A(xi^l/2,yt_1/2) + 4 + fcj-i-B(gt-+1/2,yt--i/2) + fcj-i-B(gt-i/2»yi-i/2) 4fcy_i(fcy + fcy_i) C(xt-,yy) /?»• ay] + Uitj [ ( - 4 / t j _ i - 4 / i r M i - i ) x {A(gj + i / 2 ,y t - + i / 2 ) + A ( x t + 1 / 2 , y t _ 1 / 2 ) - A(gj_i/ 2 ,y t - + 1 / 2 ) - A ( x t _ 1 / 2 , y i-1/2)} 4hihi-i(hi + fot_i) (—4fcy_i — 4fcy — kjkj-i) x { j j ( x i + 1 / 2 , y i + 1 / 2 ) - g(z,-+i/2,y,--i/2) + B(a t -_ 1 / 2 ,y t - + 1 / 2 ) - -B(x^_ 1 / 2,y {_ 1 / 2)} 4 A;y /ii -y (fcy + fcy _ 1) C(x t-,yy) ft /?y *7, »,y+i -2hj-_i - 2fet- - fetfej-i^(gt-i/2> yt+1/2) + hihi-iA(xi+1f2,yi+1/2) 4hihi-i(hi + hi^i) 2 - A:yJ3(zf_i/2,y.-+i/2) + kjB{xi+1/2, y»+i/ 2 ) + 2 A/j ^ A>y -h" A*y — 1 ^ yy) ft Ty + i+ i , y - i 2 - hiA(xj+1/2,yi_1/2) 2 + fcy-i-B(xi+1/2,yt-_1/2) 4/ii(/ii + hi-i) 4kj-i(kj + fcy_i) C(x»,yy) T i «y] + 45 2 - hiA(xi+1/2,yi+i/2) ~ hiA(xi+1/2,yi_l/2) 2kj + 2fcy_i - kjkj-1B(xi+1/2,yi-1/2) + kjkj-iB(Xi+1/2,yi+1/2) Akjkj-\(kj + fcy-i) C(^t,yy) 7i /?y] + 17* +i,y+i 2 - / i t v l ( x l + 1 / 2 , y i + i / 2 ) 2 - kjB(xi+1/2,yi+1/2) 4kj(kj + /cy_i) C{xi,yj) a 7y] D{xi,yj) i = 2,..,n, j = 2,..,m. Note that the same procedure may be followed for an equation which is more general than (3.1) [G{x,y)U)xx + {F{x,y)U)yy = {A(x,y)U)x + {B{x,y)U)y + C{x,y)U + D{x,y). Now the order of the local truncation error of the nine-point formula can be found by substituting the local truncation errors from equations (3.12a), (3.12b), and from the biquadratic interpolation into equation (3.12c). - Assuming stability (which is guar-anteed to hold under standard assumptions), the error of this scheme is second-order as well, even though the direct local truncation error for (3.13) is calculated to be Th = 0{h) 46 when the Taylor expansions for U are substituted into equation (3.13) 3 . 3 N u m e r i c a l R e s u l t s The nine-point formula (3.13) was used to find approximate solutions to four two-dimensional boundary-value problems with Dirichlet boundary conditions. The exact solutions to all four problems are known. 1) Uxx + Uyy = U + exp(x + y ) , x e [ o , i ] , y e [o,i], U{0,y)=ev, U{l,y)=ev+1, U[x,0) = ex, U(x,l) = e x + \ solution U = exp(x + y ) . 2) Uxx + Uyy = -xUx -yUy + 3x2 - 3 y 2 , x E [0 ,1] , ye [0 ,1] , U(0,y)=y2, U{l,y) = l + y2, U(x,0) = x2, U[xtl) = 1 + x 2 , solution U = x2 + y 2 . 3) Uxx + Uyy = -l/xUx + xyUy - xU, x G [ l , 2 ] , y G [ l , 2 ] , 47 U{l,y)=0, C/(2,y) =y log2 , U(x,l) = logx, U(x,2) = 2 logx, solution U = ylog(x). 4) Uxx + Uyy = T T 2 ( X + y)Ux - j^jUy - Tr 2U - Tr 2 C O S ( T T X + Try) , x e [ i , 2 ] , y e [0,2], U(l,y) = ^ ^ , tf(2,y) = ^ g ± ^ , U(x,0) = s-^, C 7 ( x , 2 ) - 5 iE l 2+ x l solution U 2+x sin(x + y) x + y To create a grid of mesh points, a one-dimensional x-mesh and a one-dimensional y-mesh were created independently. The x-mesh and the y-mesh were both created in the same way as the meshes for the one-dimensional case, using the same mesh types. The number of mesh points for each of the x-meshes and y-meshes were randomly generated between 9 and 200. The maximum error and the average error of the approximate solutions are calculated from the local errors (as in the one-dimensional case). The maximum step size h is then plotted against the maximum error. The accuracy of the formula (3.13) is determined by the slope of the linear least-squares fit to the scatterplot of points. Figures 14 to 17 show the results for problems 1) to 4), respectively. The two-dimensional method given by (3.13) is obviously second-order. 48 Figure 14 Results for Two-dimensional Problem 1 slope = 2.14 MAX STEP SIZE 49 Figure 15 Results for Two-dimensional Problem 2 slope = 1.98 50 Figure 16 Results for Two-dimensional Problem 3 5] 52 4 . S I N G U L A R L Y P E R T U R B E D P R O B L E M S 4 . 1 E x p o n e n t i a l F i t t i n g S i n g u l a r l y p e r t u r b e d p rob l ems are c h a r a c t e r i z e d by the occur rence of a s m a l l pa -r a m e t e r e t h a t m u l t i p l i e s the h ighes t de r iva t ive i n the p r o b l e m . T h i s causes the p r o b -l e m t o have so lu t ions w h i c h m a y v a r y r a p i d l y i n some regions a n d s l o w l y i n o thers . T h e regions of r a p i d v a r i a t i o n are c a l l e d layers; the o ther regions are s m o o t h regions . B e c a u s e of the c o m p l i c a t e d na tu re o f t he i r d i f ferent ia l opera tors a n d so lu t ions , s i ngu -l a r l y p e r t u r b e d p r o b l e m s are t r ea t ed separa te ly f r o m the regu la r p r o b l e m s o f C h a p t e r 2. C o n s i d e r the s i n g u l a r l y p e r t u r b e d O D E ey" = {r{x)yy + q{x)y + f{x) (4.1) a < x < b w h e r e e is a s m a l l pa r ame te r , 0 < E « 1 a n d r , q , a n d f are k n o w n b o u n d e d s m o o t h func t ions as i n (1.1a). T h e s i gn o f r ( x ) over the i n t e rva l [a, b] reveals a lo t a b o u t the s o l u t i o n of the d i f ferent ia l e q u a t i o n . If r ( x ) < 0 over the ent i re i n t e r v a l , t h e n the o r d i n a r y different ia l e q u a t i o n has two 53 solution modes - one fast, one slow. These modes both exist throughout the entire interval [a, b], but near the endpoint x = a, the fast solution mode may be dominant, causing a boundary layer to occur in this region. Away from the endpoint the slow solution mode becomes dominant, and a smooth region is seen. When r(x) > 0 over the interval, the layer may possibly occur near the right endpoint and not the left. If r(x) changes sign over the interval then the fast/slow solution components cannot be separated over the entire interval. A turning point may then occur. We will not deal with the latter type of problems here; we will assume that r{x) is bounded away from zero throughout the interval [a, 6]. Numerically, the equation (4.1) can be solved by three-point schemes, like regular problems, but the schemes must be specially formulated to suit the nature of the problem. Following the methods set out in Chapter 2, we consider the conversion of (4.1) to the first-order system ey' =r(x) y + z, (4.2a) z'=q{x)y + f{x), (4.26) a < x < 6. Since (4.2b) is not singularly perturbed, and r(x) is bounded away from zero the slow solution component z is weakly coupled to the fast solution component y in the system (4.2) . Often solutions to (4.2) are found through discretization on a mesh that reflects the nature of the solution by having a high density of mesh points over the layer region, and a sparser distribution of mesh points elsewhere. This is called resolving the layer. For many common difference schemes, including those used in Chapter 2, resolving the layer is necessary because errors in the solution in one region propagate undamped elsewhere. We wish to avoid this. Constructing such meshes requires an a priori knowledge of the position of the layers. Although this information is known for the simple problem (4.1), for more complicated examples it may not be readily available. Another problem with many of the common difference schemes is that 54 they have trouble distinguishing the different solution modes of singularly perturbed problems. When one-step centered difference methods are used, fast solution modes are approximated by slow modes. This causes some difficulties. First of all, errors in the solution over the layer region will propagate into the smooth regions. Also local truncation error control is more difficult, and iterative methods for solving (1.2) will not work. Thus one-sided schemes, such as the backward Euler scheme, are often used in the solution of singularly perturbed problems like (4.2a). The one-sided schemes act to localize errors so that any errors from the layer region die out as the smooth solution becomes dominant. However, one-sided schemes cannot decouple rapidly increasing modes from rapidly decreasing modes. Since a rapidly increasing mode cannot be approximated by a rapidly deceasing mode (and vice versa), a one-sided scheme with the appropriate direction of one-sidedness must be used. This is called upwinding. The method of exponential fitting (see for example [1]) is a special one-sided scheme that does automatic upwinding for (4.1) when |r(x)| is bounded away from 0. This scheme gives accurate solutions away from the layers, even on a mesh that has not resolved the layers. Thus we will use it in the discretization of the system (4.2). To solve the system (4.2) both equations (4.2a) and (4.2b) are discretized on an arbitrary mesh 7 r n just as equations (2.1a) and (2.1b) were discretized for the regular problem. Therefore (4.2a) will be discretized at the mesh points, and (4.2b) will be discretized at the mid-points of the mesh. Since (4.2b) is not singularly perturbed, the same difference schemes that were used for the regular problem can be used in its discretization. Using the midpoint scheme on an arbitrary mesh yields the following discretization t = 2,. . ,n, where i{ is the point mid-way between X j - i / 2 and x I + 1 / 2 , and y(x t) is approximated by y t using quadratic interpolation. 55 O n a u n i f o r m m e s h th is becomes 21+1/2 I Zt~1/2 = *(*,) V, + f(x,) (4.36) i = 2 , . . , n , where h is the u n i f o r m step s ize. O r i f the t r a p e z o i d a l scheme is used o n a n a r b i t r a r y m e s h the d i s c r e t i z a t i o n of (4.2b) is zi+l/2 ~ zi-l/2 _ l/2{hi + hi_l) ~ J 9 ( ^ - 1 / 2 ) {yi + yt-i) + ^ / ( « t - i / 2 ) + 2 ? ( X i + V 2 ) ( y i+1 + y i ) + \ / ( ^ t + l / 2 ) - (4.3c) i — 2 , . . , n. O t h e r second-order difference schemes c a n a lso be a p p l i e d , j u s t as i n the case o f the r egu la r o n e - d i m e n s i o n a l p r o b l e m . E q u a t i o n (4.2a) , o n the o ther h a n d , is s i n g u l a r l y p e r t u r b e d , g i v i n g the fast s o l u t i o n c o m p o n e n t of p r o b l e m (4.1). E x p o n e n t i a l fitting is used i n s t ead o f the u s u a l cen te red difference schemes. In th is d i s c r e t i z a t i o n the func t ions r a n d z are a p p r o x i m a t e d b y piecewise cons t an t func t ions . r = r{xi+1/2) z = Zi+1/2 (4-4) o n X{ < x < X{+i, i=l,..,n. These a p p r o x i m a t i o n s are reasonable s ince r (x ) is s m o o t h a n d b o u n d e d away f r o m zero , a n d s ince z varies s lowly . O n c e these a p p r o x i m a t i o n s are m a d e , an e x p o n e n t i a l t r a n s f o r m a t i o n is a p p l i e d t o (4.2a) y = u e x p [ e _ 1 r ( x l + 1 / 2 ) {x - x,-)] g i v i n g EU = e x p [ - e _ 1 r ( x i + 1 / 2 ) (x - x t-)] zi+1/2 (4.5) 56 X{ < X < £ t - + i , = l , . . , r c . dx (4.5) c a n be i n t e g r a t e d e x a c t l y , / eu' dx= e x p [ - e " 1 r ( x i + 1 / 2 ) (x - x t ) ] zi+1/2 J Xi J Xi to g ive e [u(x t -+i) - u{xi)) = 1 + 1 / 2 [ e x p { - e - 1 r ( x t + 1 / 2 ) ( x i + 1 - x ; )} - 1]. (4.6) e r l Z i + i / 2 j R e p l a c i n g x t - + i — Xj- by i ts equiva len t hi , a n d u(xj-) by yi, a n d u ( x i + i ) b y y i + i e x p [ - e : _ 1 r ( x t + 1 / 2 ) /i;] i n (4.6) y ie lds e yi+i e x p [ - e _ 1 r ( x l + 1 / 2 ) hi] - ey{ = E Z t + i / 5 | | e x p j _ g ^ ( x - 1 / 2) /i^] - l } . (4.7) E q u a t i o n (4.7) c a n be s i m p l i f i e d by us ing the B e r n o u l l i f u n c t i o n B(x) x ex - 1 T h e d i s c r e t i z a t i o n for e q u a t i o n (4.2a) u s ing e x p o n e n t i a l f i t t i n g w i l l be Z l + 1 / 2 = h~ N o w t o get a th ree -po in t f o r m u l a B ^r{xi+1/2)hi\ ^ B / - r ( x i + 1 / 2 ) h, (4.8) « i yi-i + bi yi + Ci y l + i = gi, i = 2,.., n, (4.8) is s u b s t i t u t e d i n to the d i s c r e t i z a t i o n for (4 .2b) . If the m i d p o i n t scheme o n a u n i f o r m m e s h (4.3b) is used , the cons tan ts ai, bi, Ci, gi are g i v e n b y eB (-r(x-j»/»>h) _ £ £ f-r(xj-1/2)h^ ^ ^r(xi+1/2)h^ bi = ~ q{xi), Ci = h2 £ B ^ r ( x t + , / 2 ) h ^ h2 h2 9i = f{xi). (4.9a) 57 O n a n a r b i t r a r y m e s h th i s becomes 2sB (~r(z-^2)fe-1) - 2eB b l hi-^hi + fct-_!) ki{hi + fct-_i) ^ A ' r( I| ,+ I/2)fc.' 1 /it(/j-t + hi-i) ft = f[xi). (4.96) N o t e t h a t (4.9b) reduces t o (4.9b) o n a u n i f o r m m e s h . If the t r a p e z o i d a l scheme (4.3c) is u sed o n an a r b i t r a r y m e s h « f l ( - r ' - - y ) t ' - ) = ~ h~, i-°(*.-l/2). _ e J 5 ^ - r ( I . - i/2 ) f e ; - i ^ Eg ^r(xi+1/2)hj j 6,- = hi—i . . hi . — ^ — 9 ( ^ 1 - 1 / 2 ) - - j g ( z i + i / 2 ) , I T J /ii c* = 7^ -jQ{xi+1/2), Qi = ^ f i / ( x l _ i / 2 ) + ^ / ( x l + 1 / 2 ) . (4.9c) 58 4 . 2 D i s c r e t i z a t i o n E r r o r s Consider first the scheme (4.9a) derived on a uniform mesh. The discretization error in this scheme can be seen by looking at the local truncation error, Th, for the equation (4.8) which is found by substituting Taylor expansions for yi+i, yi, and yi-\ about the point x I + 1 / 2 into (4.8). Th =B B r{xi+i/2)h\ e e J h h y{xi+i/2) + ^ y'(xi+1/2) + o(h2) y{xi+i/2) - ^ y'{xi+1/2) + o{h2) [ey'(xi+1/2) - r(xi+1/2) y{xi+1/2)} (4.10) Obviously the order of magnitude of the local truncation error given by (4.10) depends on both the comparative sizes of h and e, and on the accuracy of the approximation (4.4). Now consider the case where e << h . If r(x) > 0 over [a,b] then 'T{xi+\/2)hs B « 1. Hence (4.10) becomes rh » - B -r{xi+i/2)h\ E y{xi+i/2) - -y'{xi+1/2) + o{h2) £y'{xi+1/2) + r{xi+1/2)y(xi+1/2) which is equivalent to / h r{xj+i/2) V e x p f - r ^ , , ) E h h , , y{xl+i/2) - -v {xi+1/2) + o{h ) Since exp ^y'{xi+1/2) + r ( x i + 1 / 2 ) y ( x i + 1 / 2 ) . -r{xi+i/2) « 1, Th = -r{xi+1/2) y{xi+1/2)-r h r{xi+1/2) , y'{xi+i/2)+r(xi+1/2) y{xi+i/2) + 0(h ) + O{E), 59 and so the local truncation error is given by rh = 0{h). (4.11) Now if r(x) < 0, then and a similar argument leads to the following simplification for (4.10) and so (4.11) still holds. The discretization scheme, then, is first-order uniformly in e. Note that the same procedures may be followed to find the local truncation error for the schemes derived on nonuniform meshes ((4.9b) and (4.9c) will also be first-order uniformly in E). Thus the scheme we have developed retains the "uniform" order of accuracy on a nonuniform mesh. The accuracy, then, is first-order regardless of the uniformity of the mesh. Note that it is the initial approximation (4.4) of r and z by piecewise constants that restricts the discretization error of the formulae (4.9) to first-order accuracy. If r(x) and z(x) are constants, then the approximation (4.4) is exact and so the discretization scheme (4.8) is also exact. The discretization error in the formulae (4.9) must be determined by the second-order local truncation error in equation (4.3). Hence, if r and z are constants, the discretization error in (4.9) is second order, regardless of the uniformity of the mesh. Also note that if hi << e , then the estimate (4.11) can be improved. If hi << E then and (4.8) is essentially a centered second-order scheme. Since (4.3) is also second order, and the discretizations were performed on a staggered mesh, the three-point formulae (4.9) will have second-order accuracy on an arbitrary mesh. 60 4.3 N u m e r i c a l R e s u l t s f o r E x p o n e n t i a l F i t t i n g T h e m e t h o d of e x p o n e n t i a l f i t t i n g was tes ted o n four s i n g u l a r l y p e r t u r b e d p r o b l e m s w i t h D i r i c h l e t b o u n d a r y c o n d i t i o n s a n d k n o w n so lu t ions . 1) e y" = —y' + y — (1 + ETV2) cos(7rx) — 7rsin(7rx) + e x p ( — x / e ) , x < E [ 0 , l ] , y (0)=0, y ( l ) = -i-exp{=±), solution cos (7rx) — e x p ( — X / E ) . 2) £ y" = [2 + c o s ( 7 r x ) ] y ' + y — (1 + sir2) cos(7rx) — 7r[2 + c o s ( 7 r x ) ] s i n ( 7 r x ) - | -N u m e r i c a l resul ts are c a l c u l a t e d u s i n g the th ree -po in t m i d p o i n t f o r m u l a (4.9b). L o c a l e r rors are f o u n d at the m e s h p o i n t s , t h e n the m a x i m u m a n d average errors are c a l c u -l a t ed . T h e m a x i m u m er ror is p l o t t e d aga ins t the m a x i m u m step s ize , a n d the s lope o f the l inea r least-squares fit is d e t e r m i n e d . T h i s s lope gives the a c c u r a c y of the m e t h o d . 61 e y" = { e - l ) y ' + y, x G [ 0 , l ] , y(0)=2, y{l) = e + exp{=^), solution y — e x p ( x ) + e x p ( — X / E ) . 4) £y" = y' + 6, x G [ 0 , l ] , y(0) = l , y ( l ) = exp{±)-6, solution y = exp(x/£r) — 6x. N o t e t h a t for p r o b l e m 4) r (x ) a n d z(x) are cons tan t s : r ( x ) = - 1 , z(x) = 6. Problems 1) and 2) were first run on uniform meshes. The results are given in Figure 18. The lower set of points are the results from problem 1), and the upper set are from problem 2). From Section 4.1, the error is expected to be first-order. Figure 18 shows that the expected results are achieved. Next each problem was run using the same meshes that were used for the one-dimensional regular problems. These non-uniform meshes do not resolve the rapid layer solution near x = 0. Figure 19 shows that the error is first-order for problems 1), 2), and 3) . The results for problem 4) are given in Figure 20. Since both r(x) and z(x) are constants for this problem, the difference scheme (4.8) is exact, and so second order errors are achieved. Lastly, all four problems were run on non-uniform meshes where the maximum step size h is much smaller than the parameter e. The graph in Figure 21 confirms the expected second-order accuracy. 62 cr o cn LUo X cr Figure 18 Results for Singular Perturbation Problems 1,2 on Uniform Meshes slope = 1.22 o o o o O o o o o o o o O <J> o° o i i i i i—i—i—r~i 1 1 1 1 — i — i — i — i — i 10-? 10-' 10° MAX STEP SIZE 63 Figure 19 Results for Singular Perturbation Problems 1,2,3 on Nonuniform Meshes 64 Figure 20 Results for Singular Perturbation Problem 4 on Nonuniform Meshes slope = 1.94 o io- J io-' io° M A X S T E P S I Z E 65 Figure 21 Results for Singular Perturbation Problems 1,2,3,4 with K « £ slope = 1.92 - i — i — i — r - T io-2 10-M R X S T E P S I Z E -i 1 — i — r i i "I 10° 66 5 . O T H E R D I F F E R E N T I A L E Q U A T I O N S T h e m e t h o d s c o n s i d e r e d i n th i s thes is for the s o l u t i o n of the l inear d i f ferent ia l equa t ions w i t h D i r i c h l e t b o u n d a r y c o n d i t i o n s (1.1) , (3.1) , a n d (4.1) m a y be ex t ended to o the r d i f ferent ia l equa t ions . C o n s i d e r f irst a quas i - l inear b o u n d a r y - v a l u e p r o b l e m y" = f{y',y,x), a < x < 6, (5.1) y(a) = 6, y(b) = #. W h e n e q u a t i o n (5.1) is l i n e a r i z e d by any o f the c o m m o n l i n e a r i z a t i o n m e t h o d s , second-o rde r a p p r o x i m a t i o n s t o i so la t ed so lu t ions c a n be f o u n d by the m e t h o d s d e s c r i b e d i n S e c t i o n 2 .1 . I so la ted so lu t ions are so lu t ions w h i c h are l o c a l l y un ique . H i g h e r - o r d e r equa t ions c a n also be s o l v e d u s ing these m e t h o d s . A n mth o rder d i f fe ren t ia l e q u a t i o n w i l l be expressed as m f i rs t -order di f ferent ia l equa t ions for t he p u r p o s e of d e r i v i n g the n u m e r i c a l scheme. Difference fo rmulas are a p p l i e d to each d i f fe ren t ia l e q u a t i o n . A s taggered m e s h Tin is used; t ha t is , each d i f ferent ia l e q u a t i o n is d i s c r e t i z e d o n different sub- in te rva ls of t he i n t e rva l [a,b] i n such a w a y tha t ( m — 1) consecu t ive s u b s t i t u t i o n s c a n be p e r f o r m e d to p r o d u c e a n ( m + l ) - b a n d e d s y s t e m A HIT = g- F o r e x a m p l e , a l inea r fou r th -o rde r d i f ferent ia l equa t i on y ( 4 ) = \ci(x) y ] ' " + [ c a ( s ) y]"+[cs{x) y}' + c4(x) y + c5(x) (5.2) c a n be r e w r i t t e n as the s y s t e m y ' =ci(x) y + p, (5.3a) 67 p' =c2{x) y + q, (5.36) q' = c 3 ( x ) y + r, (5.3c) r' =c4(x) y + c5(x). (5.3d) A m i d p o i n t scheme o n a u n i f o r m b u t s taggered m e s h c a n be a p p l i e d to the four f irs t-o rder equa t ions o n the f o l l o w i n g in te rva l s : equa t ions (5.3a) a n d (5.3c) are d i s c r e t i zed over [ x , - , X i + i ] , i = 1, ..,n a n d equa t ions (5.3b) a n d (5.3d) are d i s c r e t i z e d over [xi-i/2,xi+1/2], »' = 2 , . . , n . So t h e n V l + \ . V l = \ c i ( x i + i / 2 ) [Vi+i + Vi] + P 1 + 1 / 2 . ( 5 - 4 A ) i d , , f - = ~ e2(Xi) y i + q{, (5.46) q % = ^ c 3 ( x i + 1 / 2 ) [y,-+i + yi] + r I + 1 / 2 , (5.4c) ^ + 1 / 2 - ^ - 1 / 2 1 f \ , ( \ tz. AJ\ A n ( m + l ) - p o i n t scheme i n the u n k n o w n s yt- is d e t e r m i n e d b y ( m — 1) successive s u b s t i t u t i o n s . F o r the s y s t e m (5.4), e q u a t i o n (5.4a) is s u b s t i t u t e d i n to (5.4b) to e l i m -ina te the v a r i a b l e p. T h e n the r e s u l t i n g e q u a t i o n is s u b s t i t u t e d i n to (5.4c) a n d q is e l i m i n a t e d . F i n a l l y the new e q u a t i o n is s u b s t i t u t e d i n to e q u a t i o n (5.4d) to e l i m i n a t e r . T h i s y i e ld s a five-point f o r m u l a w i t h o n l y the m e s h f u n c t i o n yn for u n k n o w n s . W h e n a difference scheme is a p p l i e d to (5.3) o n a n o n u n i f o r m m e s h , the d i s c r e t i za t i ons b e c o m e m o r e c o m p l i c a t e d . T h e m i d p o i n t scheme o n a n o n u n i f o r m m e s h gives the f o l l o w i n g d i s c r e t i z a t i o n o f t he s y s t e m (5.3) : V l + 1h . V l = \ c i ( x i + i / 2 ) b i + i + yi] + Pi+1/2 (5.5a) 68 w, , , r - = C 2 ( x ^ y i + q i ( 5 M > 2 (ft, + / l i - l ) _ q~ = cs(xi) yi + fi (5.5c) X j ^ = c4(xi) y i + c5(xi) (5.5d) X t + l X i where x ; is the m i d p o i n t of [xi-i/2,xi+i/2}, x ; is the m i d p o i n t of [ x ; + i , x , - ] , X ; is the m i d p o i n t o f [ x l + i , X i ] . T h e m e t h o d s used i n th i s thesis c a n also be a p p l i e d t o b o u n d a r y - v a l u e p rob l ems w i t h b o u n d a r y c o n d i t i o n s t h a t are m o r e c o m p l i c a t e d t h a n the D i r i c h l e t b o u n d a r y c o n d i t i o n s c o n s i d e r e d i n the p r e v i o u s chapters . S e p a r a t e d l inea r b o u n d a r y c o n d i t i o n s of the f o r m ki y'(a) + k2 y{a) = 6X £iy(b)+£2y'(b)=62 c a n v e r y s i m p l y be i n c l u d e d i n the s y s t e m A y„ = g (1.2) w i t h o u t affect ing the t r i d i a g o n a l n a t u r e of the s y s t e m . H o w e v e r , i n c o r p o r a t i n g b o u n d -a ry c o n d i t i o n s t h a t are no t sepa ra t ed in to the s y s t e m (1.2) w i l l cause the s y s t e m to lose i ts t r i d i a g o n a l f o r m . S t i l l , a s i m i l a r s t a b i l i t y a n d cons i s tency ana lys i s c a n be s h o w n to h o l d here as for (1.2), a n d second-order accura te so lu t ions are o b t a i n e d o n a n o n u n i f o r m m e s h w h e n cen te red second-order difference schemes are used . 69 6 . C O N C L U S I O N W e have d e r i v e d a class of difference schemes to find second-order a p p r o x i m a t e so lu t ions to b o u n d a r y - v a l u e p rob lems b y d i s c r e t i z i n g o n a n a r b i t r a r y m e s h nn. T h e d e r i v a t i o n of these schemes involves , first, r e w r i t i n g the b o u n d a r y - v a l u e p r o b l e m as a s y s t e m of first-order equa t ions . E a c h e q u a t i o n is d i s c r e t i z e d b y a one-step cen te red difference scheme over s taggered sub- in te rva ls o n the m e s h nn. T h e first e q u a t i o n i n the s y s t e m , i n v o l v i n g the de r iva t ive of the m e s h f u n c t i o n (the p r i n c i p a l u n k n o w n ) , is d i s c r e t i zed at the m e s h po in t s ; the nex t e q u a t i o n , w h i c h has the de r iva t ive o f a n a u x i l i a r y u n k n o w n , is d i s c r e t i zed at the m i d p o i n t s o f the m e s h . A n y fu r the r equa t ions are d i s c r e t i z e d at po in t s m i d - w a y be tween the d i s c r e t i z a t i o n p o i n t s of the p r e c e d i n g e q u a t i o n . T h i s use of s taggered m e s h sub- in te rva l s a l lows the s u b s t i t u t i o n o f one e q u a t i o n in to the n e x t so tha t , u l t i m a t e l y , one e q u a t i o n i n v o l v i n g o n l y the p r i n c i p a l m e s h f u n c t i o n resu l t s . F o r an m t / l - o r d e r b o u n d a r y - v a l u e p r o b l e m , a n ( m + l ) - p o i n t difference f o r m u l a is p r o d u c e d i n th is way. S i m i l a r l y , an (rn + l ) 2 - p o i n t f o r m u l a c a n be o b t a i n e d for b o u n d a r y - v a l u e p rob lems posed i n t w o d i m e n s i o n s . T h e difference fo rmulas c a n also be a p p l i e d to l i nea r i zed fo rms of n o n l i n e a r equa t ions . S i n g u l a r l y p e r t u r b e d p rob l ems w i t h o u t t u r n i n g po in t s c a n a lso be so lved u s i n g s i m i l a r p rocedures , b u t some res t r i c t ions are i m p o s e d . F i r s t o f a l l w h e n r e w r i t i n g the 70 singularly perturbed boundary-value problem as a system of first-order equations, the transformation must be done in such a way that the fast and slow solution modes are weakly coupled or completely decoupled. The equation for the slow solution can be discretized using the usual difference schemes, for example,the midpoint or trapezoidal schemes. For the fast solution a fast-solution discretization method must be used, such as exponential fitting. There are two advantages to using this method over other singular perturbation methods. First, for a mesh where the maximum step size is much larger than the parameter e, the mesh need not Tesolve the solution layer to produce accurate results elsewhere; and second, explicit upwinding need not be applied because the method of exponential fitting does automatic upwinding. However, this method is onlt applicable to a limited class of problems. The three-point formulae derived in Chapter 2 for second-order linear boundary-value problems with Dirichlet boundary conditions are shown to be theoretically second-order and stable. The theory is supported by the numerical results included in this thesis. The second-order results are also shown to hold for linear boundary-value problems posed in two (or more) dimensions. For the three-point formulae derived for linear second-order singularly perturbed problems on an arbitrary mesh we have shown that the order results are obtained regardless of the uniformity of the mesh. Numerical results support all of these findings. Thus the methods presented here will be efficient and efficacious for solving differential equations that arise in mathematical modelling. 71 R e f e r e n c e s 1. U.Ascher, R. Matheij, &; R. Russel, Numerical Solution of Boundary Value Problems for Ordinary Differential Equations, Prentice Hall, New Jersey, 1988. 2. R. Burden, & J . Faires, Numerical Analysis, Prindle Weber and Schmidt, Boston, 1985. 3. T . H . Chong, " A Variable Mesh Finite Difference Method for Solving a Class of Parabolic Differential Equations in One Space Variable," S I A M J . Numer. Anal., Vol. 15, 1978, pp. 835-857. 4. J .D. Hoffman, "Relationship Between the Truncation Errors of Centered Finite-Difference Approximations on Uniform and Nonuniform Meshes," J . Comp. Phys., Vol. 46, 1982, pp. 469-474. 5. H.B. Keller, "Accurate Difference Methods for Linear Ordinary Differential Systems Subject to Linear Constraints," S I A M J . Numer. Anal. , Vol. 6, 1969, pp. 8-30. 6. H .B . Keller, &: A . B . White Jr., "Difference Methods for Boundary Value Problems in Ordinary Differential Equations," S I A M J . Numer. Anal., Vol . 12, 1975, pp. 791-802. 7. H . Kreiss, T. Manteuffel, B . Swartz, B . Wendroff, & A . B . White Jr., "Supra-convergent Schemes on Irregular Grids," Math. Comp., Vol. 47, 1986, pp. 537-554. 8. M . Lentini, &: V . Pereyra, "A Variable Order Finite Difference Method for Nonlinear Multipoint Boundary-Value Problems," Math. Comp., Vol. 28, 1974, pp. 981-1003. 9. T. Manteuffel, & A . B . White Jr., "The Numerical Solution of Second Order Bound-ary Value Problems on Nonuniform Meshes," Math. Comp., Vol. 47, 1986, pp. 511-535. 10. T. Manteuffel, &: A . B . White Jr., "On the Efficient Numerical Solution of systems of Second Order Boundary Value Problems," S I A M J. Numer. Anal . , Vol. 23, 1986, pp. 996-1006. 11. R. O'Malley, Introduction to Singular Perturbations, Academic Press, New York, 1974. 12. C .E . Pearson, "On a Differential Equation of Boundary Layer Type," J . Math. Phys., Vol. 47, 1968, pp. 134-154. 72 13. A . N . Tikhonov, &; A . A . Samarskii, "Homogeneous Difference Schemes on Nonuni-form Nets," Zh. Vychisl. Mat. i Mat. Fiz. (English Translation), Vol. 1, 1962, pp. 5-63. 14. E . Wiktor, Matched Asymptotic Expansions and Singular Perturbations, North Holland Pub. Co., New York, 1973. 15. D. Young, & R. Gregory, A Survey of Numerical Mathematics, Vol 1,2, Addison Wesley Pub., London Ont., 1984.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- A supra-convergent scheme for the solution of differential...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
A supra-convergent scheme for the solution of differential equations on an arbitrary mesh Scott, Karen 1988
pdf
Page Metadata
Item Metadata
Title | A supra-convergent scheme for the solution of differential equations on an arbitrary mesh |
Creator |
Scott, Karen |
Publisher | University of British Columbia |
Date Issued | 1988 |
Description | This thesis describes a compact numerical scheme to solve second-order linear differential equations with Dirichlet boundary conditions on arbitrary meshes. The scheme uses a staggered grid to achieve second-order accuracy on a nonuniform mesh. Regular and singularly perturbed problems in one and higher dimensions are considered. Numerical experiments are presented which support the theoretical results. The methods presented are also applicable to equations with other types of boundary conditions, and to nonlinear and higher-order differential equations. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-09-10 |
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.0079533 |
URI | http://hdl.handle.net/2429/28394 |
Degree |
Master of Science - MSc |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-UBC_1988_A6_7 S39.pdf [ 3.18MB ]
- Metadata
- JSON: 831-1.0079533.json
- JSON-LD: 831-1.0079533-ld.json
- RDF/XML (Pretty): 831-1.0079533-rdf.xml
- RDF/JSON: 831-1.0079533-rdf.json
- Turtle: 831-1.0079533-turtle.txt
- N-Triples: 831-1.0079533-rdf-ntriples.txt
- Original Record: 831-1.0079533-source.json
- Full Text
- 831-1.0079533-fulltext.txt
- Citation
- 831-1.0079533.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:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0079533/manifest