A SUPRA-CONVERGENT SCHEME FOR T H E SOLUTION OF DIFFERENTIAL EQUATIONS ON AN ARBITRARY MESH by K A R E N SCOTT B.Sc, The University of Calgary, 1986 A THESIS SUBMITTED IN PARTIAL FULLFILMENT OF T H E REQUIREMENTS FOR T H E D E G R E E 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 O F BRITISH COLUMBIA July 1988 ©Karen Scott, 1988 In presenting degree at the freely available copying this of department publication of in University of for reference this or thesis thesis by this his for or partial British Columbia, and study. scholarly her of the Department of The University of British Columbia 1956 Main Mall Vancouver, Canada V6T 1Y3 requirements I agree for purposes may representatives. be It is an advanced that the Library shall make it I further agree that permission thesis for financial gain shall not permission. DE-6G/81) fulfilment granted by the understood that for extensive head of my copying or be allowed without my written ABSTRACT T h i s thesis describes a c o m p a c t n u m e r i c a l scheme to solve second-order linear differential equations 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 on a r b i t r a r y meshes. The s c h e m e uses a s t a g g e r e d g r i d t o a c h i e v e s e c o n d - o r d e 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 . 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 o n e a n d h i g h e r d i m e n s i o n s are c o n s i d e r e d . N u m e r i c a l e x p e r i m e n t s are p r e s e n t e d w h i c h s u p p o r t t h e t h e o r e t i c a l r e s u l t s . T h e m e t h o d s p r e s e n t e d are a l s o a p p l i c a b l e t o e q u a t i o n s w i t h o t h e r t y p e s o f b o u n d a r y conditions, a n d to nonlinear a n d higher-order differential equations. 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 List of Figures 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 « 66 iv e Acknowledgement Thanks to Uri for all his patience. V 1. I N T R O D U C T I O N One of the major problems encountered in the numerical solution of boundaryvalue 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 because 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 equation 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) y + q{x)y + f{x) (1.1a) y 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 ir with n sub-intervals (or n+1 mesh points n including the end points of the interval) so that x\ = a, x TTn = {xi,X ,X , 2 3 n + i =b ...X ,X } n n+1 where X{ i = X{ + hi, + hi is called the step size of the i t h hi > 0, i=l,..,n. 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, y +i} T n 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 . 1 c ) , are used to determine the values of the mesh function. The error, e = y,-— y(xi), «' = l,..,n + 1, in these calculated values has two t 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 e- is the discretization error, that is, the t 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 hasfirst-orderaccuracy if T = O(h), and second-order accuracy if Th = 0(h ). The 2 n 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 i = 1,.., n + 1. 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 onlyfirst-orderapproximations on nonuniform meshes. 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 -_j + bi yi + c - y t t t +1 = #, i = 2,.., n, t where the a^'s, 6 's, c;'s, g^s are constants determined by the discretization; from the t three-point formula an n x n system of equations can be set up to solve for the mesh function y w Ay, = g. (1.2) The tridiagonal matrix A has the 6- 's forming the diagonal, the a; 's the lower diagonal, t and the Ci 's the upper diagonal. The right hand side g is the vector g = {0, 92, g ,-,g , s n #} 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 Previous a n d Related 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 secondorder 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 problem, 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 + i } T n b i ={6, K u K ,...,K , 2 5 n #} T and the matrix Bi has the form (B a Si Ri S Bi 2 i?2 S n R n BJ V b 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 B a and B represent the boundary conditions. On an arbitrary nonuniform mesh b 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 tridiagonal. 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 y = g. They assumed stability to show w 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 problems 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 problems. 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 t h e s e s c h e m e s w e r e d e v e l o p e d as v e r t e x - c e n t e r e d s c h e m e s , m e a n i n g t h a t u n k n o w n s a r e e v a l u a t e d o n l y a t 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 s c u s s c e l l c e n t e r e d s c h e m e s for (1.1a), w h e r e u n k n o w n s a r e o n l y e v a l u a t e d a t t h e p o i n t s m i d - w a y b e t w e e n t h e m e s h p o i n t s . T h e y s h o w t h a t , e v e n t h o u g h t h e i r s c h e m e is i n c o n s i s t e n t ( t h a t i s , t h e a p p r o x i m a t i o n t o y" is i n c o n s i s t e n t ) , i t s t i l l y i e l d s a t r i d i a g o n a l s y s t e m w i t h s e c o n d - o r d e 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 n e a 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 separated linear b o u n d a r y conditions. A g a i n they assume b u t d o n o t show stability. T h i s t h e s i s g i v e s a p r o c e d u r e t o p r o d u c e s t a b l e c o m p a c t - a s - p o s s i b l e difference schemes w i t h second-order accuracy o n n o n u n i f o r m meshes w i t h o u t t h e l i m i t a t i o n s present i n some o f t h e earlier w o r k . T h e schemes are neither vertex-centered o r cellcentered, b u t a combination of both. to be converted to a first-order T h e d i f f e r e n t i a l e q u a t i o n (1.1a) d o e s n o t n e e d s y s t e m i n order for these schemes t o b e a p p l i e d , even t h o u g h s u c h a c o n v e r s i o n is u s e d f o r d e r i v i n g t h e s c h e m e s . 7 2. 2.1 A S E C O N D - O R D E R T H R E E - P O I N T D e r i v i n g a S e c o n d - O r d e r Three-Point S C H E M E 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 o n e d i m e n s i o n is t h e e q u a t i o n (1.1) y" =(r{x)y)' + q{x)y + f{x), y(o) = y{b) = 0 , 9, (1.1) o n t h e i n t e r v a 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 t h e f i r s t - o r d e r s y s t e m o f ODE'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 t o t h e s y s t e m (2.1) h a s s o m e a d v a n t a g e s o v e r t h e m o r e c o m m o n transformation ')'•(• y'J F o r t h e s y s t e m ( 2 . 1 ) , r'(x) MCW )- 0 \q(x) + r'(x) 0 r(x)J \y'J \f(x)J <"> does n o t need t o be e v a l u a t e d , whereas t o solve t h e s y s t e m (2.2) i t is n e c e s s a r y t o e v a l u a t e r'(x). T h i s f u n c t i o n is n o t a l w a y s e x p l i c i t l y a v a i l a b l e . 8 Also for some applications, (r(x) not separate r{x) y)' may be smooth while y' is not. Since (2.1) does 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 n. n 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 = b n (1.3) x 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+ 9 1, but the values of z are approximated at the points half-way between the mesh points , ~ {Xi+l/2) Zi+l/2 Xi+1/2 z = %i + y i= l,..,n (See Figure 2 below). Figure 2 Evaluation Points on a Staggered Mesh Zi-l/2 z t'+l/2 3/t-i yt+i Xi+i ^ i - i On a staggered mesh, the discretization of the system of first-order ordinary differential equations (2.1) by one-step centered difference schemes can be represented by the banded matrix system Bv n (2.3) = b where V , r ={2/1, ^1 + 1/2, 2/2, 2 2 + 1 / 2 , — , n + l / 2 , 2 b={6, 0, / , 0, / , . . . , / , 0, # } 2 3 n and the matrix B has the form 10 T Vn+l} T (1 -1 Si *1 1 ~h "2 2 -1 B = t 2 -1 U3 _1_ V3 ^3 h 3 V W 7- n 1 n tn V 1 The s 's and t\s are determined by the equation t * i Vi + ti y i + i - z / i+1 =0, 2 i = 1, . . , n , the U i ' s , v,-'s, and w^s by Ui yi_! and hi = \(hi 1 + ^Zi_i/2 + v hi y,- + ^ - 2 hi { +/i;_i). , 1 l + 1 / 2 + Wi y i + 1 i = = 2,..,n Because the discretization was performed on a staggered mesh, the system (2.3) can be reduced to a tridiagonal system A y The i (1.2) = g. n row of the tridiagonal system, i = 2,.., n, is created by adding from (2.3) J-row[2t] + row[2i - 1] - J - r o w [ 2 ^ - 2 ] . hi hi 2 (2.4) 2 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 st and n + 1 st rows of the system (1.2). The tridiagonal system is a matrix representation of a three-point formula a; Vi-i + bi y l + c t y l + i 11 = gi i = 2,..,n. 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{, x i] i+ , i = l,..,n is . y%+1 V% = h where Zt+i/2 = Xi + a s y{xi+i/ ) + z{x / ) i i+i/2) r x a D o v e - (2.5a) i+1 2 2 Since r(x) is a given function, r(xi i/ ) + c a n D e 2 evaluated exactly. The value for y ( £ + i / 2 ) will be approximated by |(j/i+i + yi) , and t the value for will z(xi / ) +1 2 y t + \ . D e approximated by Zi i/ + V = l \ r So (2.5a) becomes 2 3 ; Equation (2.1b) can also be discretized using the midpoint scheme. stead of being discretized on i+i/2], \ i-i/2> x x i= (2.6a) ( t + i / 2 ) {Vi+i + yi) + Zi+i/2- But in- equation (2.1b) is discretized on the interval [x,-,Xi+i], 2,..,ra . This interval has length Xi = Xi + - (hi - and midpoint, called x -, hi, t 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 for y(x ) is obtained by quadratic interpolation on yi, y;+i. (For details on this t calculation, see, for example [1, 2, 14].) This yields y = a.i y t _ ! { where +& y + q { I Pi = i 1 IQhi^^hi y , { (3/i + / i _ ) ( / i _ Oil i 1 i+1 -hi) + /i;_i) (3hj + hj-i)(3hj-i IShi-yhi + hj) (3hi_i + hi)(hi — hi_i) 16hi(h + K^i) { 12 t 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)- hi (- ) 2 6b We can use the difference formulas (2.6a) and (2.6b) directly to develop a threepoint 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 i+l/2 z Zi+1/2 = V l \ . + V Substituting this definition in for 2 ( *'+i/2) (JM-i + Vi)r ~ l Zi / , +1 x Zi-\/ 2 (2.7) equation (2.6b) to eliminate all m 2 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 + r{xi-i/2) | hi~i)hihi^i ~ r(xj )(y (yi + Vi-i) +1/2 + y -) _ t i+1 hi + hi-1 (cti y -_! q(ii) t + p\ y + { n y ) f(x), + i+1 (2.8) i = 2,.., n. Equation (2.8) can be rearranged to give the three-point formula ai a . y _i + bi y t { + c, y + -2 i ~h~h— " i+1/2 im — ( > q Xi hi + hi-1 2 l ^ -r(x ) 1/2 tiihi-x C _ hi + h i - i r(xj_ ) + i = 2,.., n, u 1/2 hi_i(hi + hi-i) bi = g - r(xj_ ) 2 = i + 1 r(z htihi: + h . ) i 1 + i + 1 / 2 ht+hi^ ) q { X l ) 1 " (2.9) 9i = f{xi). For a uniform mesh, the interval [xi_i/ ,2^+1/2] has length h and midpoint X{. So 2 then discretization of (2.1b) by the midpoint scheme gives z i+l/2 _ z 7 h i-l/2 { \ , ft \ — = q{Xi) yi + f ( X i ) . 13 (o ir\\ (2.10) S u b s t i t u t i n g (2.7) i n t o (2.10) y i e l d s t h e " u s u a l " t h r e e - p o i n t f o r m u l a c» y,- i = gi, Q>i Vi-i + bi yi + i = 2 , . . , n, + where _ a 1 r(g -_ t ~h* i , -2 6i = ^ _ J L C * /* 1 / 2 ) 2h + ' , . , r(g -_ t -«(*0 + _ R 1 / 2 Y ) -r(z t + 1 / 2 , h ) (*•••+1/2) 2/* 2 ' N o t e t h a t e q u a t i o n (2.11) is t h e s a m e as e q u a t i o n (2.9) w i t h on = 7; = 0 , a n d Bi = 1. The t h r e e - p o i n t f o r m u l a s (2.9) a n d (2.11) w e r e d e r i v e d u s i n g t h e m i d p o i n t s c h e m e . O t h e r difference schemes c a n also b e used. If t h e t r a p e z o i d a l scheme were used o n ( 2 . 1 b ) , t h e n i n s t e a d o f e q u a t i o n ( 2 . 5 b ) , t h e d i s c r e t i z a t o n o f (2.1b) w o u l d b e z i+l/2 ~ i-l/2 _ z 1/2(^ + ^ - 1 ) ~ - 9(^-1/2) (yi + yi-i) + g / ( i - i / 2 ) + x ^ 9(^+1/2) {Vi+i +Vi) + \ f{x / ). (2.12) i+1 2 S u b s t i t u t i n g (2.7) i n t o (2.12) t o e l i m i n a t e z f r o m t h e e q u a t i o n leads t o t h e t h r e e - p o i n t formula where 4 a; ^ -_i(/i,: + hi-i) t _ - - t + hi-i) 2 hi + h^! -4 hi-i(hi " 2r(x -_!/ ) 1 + T — T — - 0 9(^i-i/2), 2 4 2r(a -_ t hi(hi + hi-i) 1 / 2 ) - 2r(x t + 1 / 2 ) /i,-+/i -_i t ^ 9(^-1/2) - ^ 9(^+1/2), 4 hi(hi + hi-i) 2r{x ) 1 hi + hi-i 2 i+1/2 9i = /(xt+1/2) + / ( ^ i - 1 / 2 ) 9(^i+l/2j, (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 — i-l/2 — z ^hi-i q{xi_ ) hii q{x ) y -_i) + ^/i _i /(x -_i/ ) + {yi+ l/2 t (y- + yi+i) + i+l/2 t t t 2 f(x ). (2.15) i+1/2 This yields the three-point formula where a 1 r( i-i/2) 2 1 -1 f hi-i g 1 "~fc,--i r(xj_ ) 1/2 ht + 9(^-1/2) 7 4 - r(x / ) i+1 2 hj-i g{xj_ ) + h 1/2 { g{x / ) i+1 2 4 2 - _ L _L ( « + l/2) _ fet g ( ^ t - l / 2 ) r C i ~fci ffi Z 2 =^(/(^'+l/2) + 4 /(^i-1/2)), (2.16) i = 2,..,n. 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 the O r d e r o fthe E r r o r How is the order of the discretization error of these three-point formulae determined? First the local truncation error over the i sub-interval, r^., is found for th equation (2.7). This is done using the Taylor expansions for j/(x +i) and y(x ) about t the point X j + 1 (These expansions hold for all i — 1, ..,n ). / . 2 y{x ) =y{x ) h + -^y'{x y{ ) =y{x ) h- jy'{x i+1 i+1/2 Xi t i+1/2 h ) + -£-y"{x ) h + -£y'"{x ) 2 i+1/2 3 i+l/2 h + ^-y"{x 2 i+1/2 ) + 0{h ), 4 { i+1/2 h ) - ^y"'(x 3 i+1/2 ) + 0{h ). 4 i + 1 / 2 { Then it is easily determined that h + ^Ly»'( 3 y{x ) - y( ) =ht y'{x i+l Xi i+1/2 ) + 0{hi ), ) 4 Xi+1/2 h + -^y"{x (2.17) 2 y{x ) + y{xi) =2 y{x ) i+1 ) + 0{h ), (2.18) 4 i+l/2 i+1/2 { 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). h r hi =- z + y'{x i + 1 / 2 i+1/2 ) 3 + ^y"'{x i+1/2 1 h - r{x i+1/2 ) [2y{x ) ) 0(h )4 2 + ~^y"{x i+1/2 + i+1/2 ) + 0(/i - )), 4 t (2.19) % = 1,.., n. But y' — r y + z according to equation (2.1a), or y'{x ) = r(x ) y{x i+1/2 By replacing y ' ( x i th r ht t + l i+1/2 i+1/2 ) +z i + 1 / 2 . (2.20) / ) with its equivalent (2.20), the local truncation error over the 2 sub-interval for equation (2.7) can be simplified to = -£y"(z ) i+1/2 r(x ) l+y2 + ^y'''(x f i+1/2 ) + 0(h ).y''(x 4 i ) l+1/2 r(x ) + 0(/ - ). 4 i + 1 / 2 il (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 r h = max \r . \ l < i < n 0(h ). 2 h The local truncation error for (2.6b) can be found using the Taylor expansions for z i+\/2-> As before, let hi = \(hi + /i -_i), the length of the Z i - i / 2 about the point ii. t interval and let 1; = z(xi). (Note that h = max hi = 0(h) [xi^ /2,Xi / ], 1 +1 2 .) 2 < i < n z(x l + 1 / 2 ) = ~z ±z>+ h t+ 2 (£) 1 1 2 \ 2 6 1 2 <*-V2) = « - J % %' + 0{hi), % + \(%) 2 (2.22) /- \ 3 + \ (j) + ~ "-l[j) z (- ) 2 23 Now substitute the exact expansions into (2.6b) and include the 0(h ) error for 2 quadratic interpolation to get r =~z[+ ±- h ~z'l' + 0(hi ) 2 hi - i q(ii) & + 0(h )) 4 2 A t - 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 r = ^ h z'l' + 0(hi ) 2 hi + q(ii) 0(h ). 4 (2.25) 2 { 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/ -+i, y^, y -_i, and y^ about the point x t t 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 y = w (2.3) b. Since each of the difference equations represented in (2.3) has second-order local truncation 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{x ,x ) 2 I x \ -Y(x ,x ) 3 I 2 (2.26) M -y(x ,x _i) n BJ V where the Y(XJ,X{) / n b are fundamental solutions of the system (2.1) satisfying F ( x , x ) = t t J. Thus these are 2x2 matrices that solve the homogeneous form of (2.1) fr(x) y ' ( x , *,•)=• The matrices B a 1\ \q(x) OJ )Y{x,x ). i and B represent the Dirichlet boundary conditions D B = a 1 0\ 0 0/ B =( b /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 to z(x / ), i+1 2 i + 1 / ) = z(x{) + 0(h). 2 z i+i / , an approximation 2 can therefore be considered a first-order approximation to z( i), X 18 and 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 n+3/2 z ~ n + l/2 z \{h i to the system. z / n+3 y i n+ n+ + f(x i) n+ n i Ugi w 2 q{x i) + h) n+ y e first-order approximation to z(b). If we also move a 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 f \ 1 t 1 2 -w t I 1 3 3 V lj 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, Z ,.., 2 y, n n, Vn+l, n z z + l} ', T with b = {o, i , o, / ,-, i , o, f{x ), 2 3 n 19 n+1 e, i?} , T and -1 U T- 2 ^2 «2 ^ fco B = *2 -1 -1 71+ Vo <l[ n+l) x 1 'n+1 0 o 0 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 that f. 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 (S R x \ 1 S R2 2 B = Rn \B BJ a b B can be equilibrated by multiplying by the matrix D = diag(R-\ R \..,R-\ 2 I) 2 to form the matrix (R-1S1 R B = 2 S 2 DB 2 R S 1 n \ B n I BJ a b We now show that the block entries R^~ Si approximate the fundamental solution 1 matrices Y(x;+i> i)x As discussed earlier, when the scheme represented by (2.3) 20 is a p p l i e d t o t h e f u n d a m e n t a l s o l u t i o n s , it c a n b e c o n s i d e r e d a first-order difference s c h e m e a n d so y i e l d s Si Y{xi,Xi) = 0{hi) + Ri Y{x xi) i+U (2.28) i — 1, ..n. Since H-R^Hoo = 0(h{), m u l t i p l y i n g (2.28) b y i 2 t h e m u l t i p l i c a t i o n o f (2.27) b y - 1 t D , which corresponds to 2 , gives R^Si Y{x ,xi) + = 0{h ). 2 i+1 (2.29) S e e n i n t h i s l i g h t i t b e c o m e s o b v i o u s t h a t t h e m a t r i x B is r o w - w i s e a s e c o n d - o r d e r a p p r o x i m a t i o n t o t h e 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 t h e e r r o r s i n (2.29) c a n b e represented in block diagonal form by E = diag[E , E ,..., 1 E, 2 n 0) where the last e n t r y corresponds to the errors i n the a p p r o x i m a t i o n of the Dirichlet boundary conditions. Note that \Ei\\ T h u s w e w r i t e B = M + E, 0{h ). = 2 or [M + E)\ = /J £>ib. w (2.30) 2 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 ) , w e find a s u i t a b l e b o u n d o n B B" ). 1 1 (a b o u n d o n F i r s t r e c a l l t h a t t h e i n v e r s e o f M is g i v e n b y / G{x ,x ) 2 ... G(x x ) G( ,X ) ... G( ,Xn+l) \G(x ,x ) ... 1 X2 M 2 u Y(x ,x )Q~ n+l 1 Y[x ,x )Q~ 2 X2 \ l 1 1 1 = n+1 2 G(x ,x ) n+1 Y(x ,x )Q~ l n+1 n+l l J w h e r e Q is a n o n s i n g u l a r m a t r i x r e p r e s e n t i n g t h e f u n d a m e n t a l s o l u t i o n s at t h e b o u n d ary conditions Q = B Y(x ,x ) a 1 1 + 21 B Y(x ,xi) b n+1 and G(xi,Xj) is the Green's function of the problem (see for example [l]). Hence a bound on the norm of M~ is given by H M l - 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 HM^Dalloo < K + 0(h) and | | M . f ? | | o o < c/cfo, -1 for some constant c. If c/c/i < 1 then Wil + M^E)- ^ < l + 0{h). 1 Now B- 1 = (I + M~ E)~ M~ D . 1 1 1 2 So for any mesh with h sufficiently small, l l ^ ! ! = K + o{h) - 1 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 Results Numerical results were calculated for six test boundary-value problems with Dirichlet boundary conditions. The exact solution is known for each boundary-value problem. 1) v" = \ y ' - ^ - y - l y(0) = 2, y(2) = e - solution : y = e~ 2) 4 + 4e~ * + 1 ( i + * ) - * , 31 2 e + e- , 2 + e~ . x x y" = -xy' + 4y, 2 y{l) = 2, y(2)=4.25, 2 solution : y = x 3) 2x e" 1 H—-. x x y" = xy' + 3y - 9 x - 8x - 3, 2 2 y(l) = 7, y(5)=211, solution : y = x + 3x + 2x + 1. z 4) 2 y" = A y + (A + 4TT ) cos (7rx) - 2TT , 2 2 2 2 2 y(0) = 0, y(l) = 0, . . exp(A x — A ) + exp(—A x) solution : y = ; - r : 1 + exp(-A^) 2 °) y" — — 7 r y(0) = 0, 2 2 COS 2 , * (7TX). tan(7rx) y', 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] : i x = • / n \Wi—rT\ sm \2{n t + \ J' l)J t = 0,..,n + l 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 ( — ] ( _ ) V2(n + l ) y K 0 a h « = 0,... n + 1 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 10 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 threepoint 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 Figure 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 M i d p o i n t 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 g o a l h e r e is t o d e v e l o p a n i n e - p o i n t f o r m u l a w i t h s e c o n d - o r d e r a c c u r a c y t o s o l v e a differential equation posed i n two dimensions, on an arbitrarily spaced rectangular m e s h b y u s i n g t h e s a m e p r o c e d u r e s as w e r e u s e d for t h e 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. Although a five-point f o r m u l a w o u l d have been preferrable f r o m some p e r s p e c t i v e s , t h e d e r i v a t i o n o f s u c h a f o r m u l a does n o t s e e m p o s s i b l e u s i n g t h e m e t h o d s d e s c r i b e d here. A l l of the procedures used to derive a second-order three-point f o r m u l a i n one d i m e n s i o n c a n be e x t e n d e d t o the d e r i v a t i o n of a second-order n i n e - p o i n t f o r m u l a in two dimensions. 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 b y s e c o n d - o r d e r a c c u r a c y . A d e f i n i t i o n a n a l o g o u s t o t h e o n e - d i m e n s i o n a l case is u s e d . S e c o n d - o r d e r a c c u r a c y m e a n s t h a t t h e e r r o r h a s o r d e r 0(h ), 2 w h e r e h is t h e m a x i m u m s t e p s i z e o n t h e g r i d 7r n j m in either the x - d i r e c t i o n or the y-direction. h = max(Hi, H2), w h e r e Hi is t h e m a x i m u m s t e p s i z e i n t h e x - d i r e c t i o n , a n d H 2 is t h e m a x i m u m s t e p size i n the y - d i r e c t i o n N o w , consider the two-dimensional elliptic differential equation U xx +U yy = {A{x,y)U) x + (B{x,y)U) h < x < Q y + C(x,y)U 0 2 < y < © 2 u 36 + D{x,y) (3.1) 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 TT ,m n where 7Tn,m = {{xi,yj),i = l , . . , n + l , i = l , . . , m + 1} x i+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 Q = U -B{x,y) y (3.2a) U. (3.26) D[x,y). (3.2c) Then, equation (3.1) can be rewritten as R + Q = C(x,y)U+ x y 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 secondorder 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), (x -+i,yj+i)} refer to t the rectangular region on the grid with the four corners (xi,yy), ( x ; , y y i ) , ( x i , y y ) , + t + ( x i , y y + i ) . In two-dimensions we create a staggered mesh by discretizing (3.2a) and t+ (3.2b) over box[(x»,yy), ( x i , y y + i)], and (3.2c) over i + 37 box[(xi_ / ) 1 2 2/?-i/2)> { i+i/2i Vj+i/v)]- The resulting formula will then be expected x 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, y y + 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 yy+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 ( ) /«(z.-+l/2,yy+l/2)+ 1 h{ ky y y fxy{xi+i/ ,yj+1/2) fk\ 2 + ^ (y ) 2 f y{x i/ ,yy+1/2) y i+ + 0{h ) (3.3) + 0{h ) (3.4) 3 2 h f{xi,y i) j+ y h y f[xi,yj) - y fx(x , = f{x i/ ,y / ) i+ 2 j+1 2 /y(«T+i/2»yj+1/2) + 2 k y /x (x S / l + 1 / 2 ,yy + 1 / 2 yy+i/2) + i+1/2 ( y ^ fxx{x /2,yj+l/2) 2 - 1 /k\ ) + - ( y ) f yy{xi+i/2,yj+1/2) = /(z;+i/2,yi+i/2) i+1 - y /i(a:i+i/2>yy+i/2) 38 3 y / y ( Z t i / 2 , 3 / j + l/2) + g /ii A; •f -f ("^) + fxy{Xi+i/2,yj f{zi+i,yj) 1 + l/2) + 2 X (k \ 2 ) / y y f o + i / 2 , yy+1/2) + 0 ( / i = / ( z + i / 2 , yy+1/2) + y / /i i \ Iy I . t /i- A; Y y yy+1/2) + 2 (3.5) / y y ( * t + 1 / 2 , yy+1/2) + C * ( / i ) (3.6) x , /z*(zt+i/2,y;+i/2)- 2 1 /*y(*i+i/2, ) 3 fx{ i+i/2,yj+i/2) t . l y / y ( z + i / 2 , yy+1/2) + kj ^ ^ ( i + l / 2 ' y ; + l/2) + (k\ 2 ( ) 3 N o w s u b t r a c t i n g (3.4) f r o m (3.6) g i v e s hi / x ( z i + i / 2 , y y + i / 2 ) - ^ / y ( x f{x yj) i + i/ ,yy i/ ) 2 + - f{xi,y ) i+u = 2 + 0{h ) (3.7) 3 j+1 a n d s u b t r a c t i n g (3.5) f r o m (3.3) g i v e s hi fx{xi+i/2,yj+i/2) f{x + kj f {x /2,yj i/ ) y i + u y i + 1 ) i+1 + - f{x ,y ) l j = 2 + 0{h ) (3.8) 3 A d d i n g e q u a t i o n s (3.7) a n d (3.8) w i l l g i v e a n e x p r e s s i o n f o r / ( £ { + i / 2 , y y + 1 / 2 ) x fx{Xi+i/2,yj + l/2) = - j - [ / ( z i , y y ) + f(x i + 2hi ,yj ) i+1 - +1 f(xi,y ) j+1 - f(xi, )} + 0(h ) S u b t r a c t i n g e q u a t i o n s (3.7) a n d (3.8) y i e l d s a n e x p r e s s i o n f o r f (xi / , y fy{xi+i/2, +1 (3.9) 2 yj 2 yy+1/2) yy+1/2) = ^-[/(x i + 1 ,yy ) + 1 f{x i,yj) + f{xi,y ) i+ j+1 - /(zi,yy)] E x p r e s s i o n s (3.9) a n d (3.10) c a n n o w b e u s e d as a s e c o n d - o r d e r two dimensions. 39 + 0(/i ) 2 (3.10) m i d p o i n t scheme i n 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 discretized 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) = U (x 2,yj i/ ) x i+1/ + = Uy(x /2,yj Q{Xi+l/2,yj+l/2) x i+l/2 + l/2) i+1 (3.11a) ~ M i+i/2,yj+i/2)U{x ,yj i/2) 2 i=:l,..,n - B(x /2,yj i+1 j = + + i/2)U(x /2,yj (3.116) l,..,m Equation (3.2c) is to be evaluated over b o x [ ( x _ / ) yy-1/2)) { i+i/2i x t + l/2) i+1 1 2 yy+i/2)]- (See Fig- ure 13). In the resulting equation Rx{ i,yj) x + Q {xi,yj) y = i = l,..,n C(x;,yy) C/(x;, yy) -+- D(x;,yy) j = (3.11c) l,..,m C(x,-, yy) and -D(x , yy) can be evaluated exactly because C and D are known functions. t For a uniform mesh (xi,yy) will be the point (x -,yy), however for an arbitrary mesh t i i,yj) x = {xi + hi/4 - hi-i/4 40 , yy + fcy/4 - fcy-i/4). F i g u r e 13 The Two-Dimensional Mesh 1 ! i zi-i x Vi-i i Xi+l hi. - i U(xi,yj) hi is a p p r o x i m a t e d b y Ui,j, found by biquadratic interpolation. e x a m p l e , for d e t a i l s o f t h e c a l c u l a t i o n s . ) + Pi Uij-i U j =ctj [ai Ui-ij-i it (3j [ai Ui- j + 7i Ui+ij-i] + Pi Uij + n U ij] lt i+ lj [cti Ui-ij+i + Pi Uij+i + 7t C/i+ij+i] where (Shj + hi-i){hj_i 16hi-i(hi — hi) + hi-i) (Shj + /t -_ )(3/t -_i + /ij) 1 1 t 16/ii-j/ii li + (3/lj_! + /lj)(/lj - /^t-l) I6hi(hi + /ij_i) 41 + (See [14], for _ a ~ j 0. (3kj + kj-i)(kj-i 16kj-i[kj (3fej + = t (xt r + 1 kj) fcy_i) ' kj) + Wkj-ikj _ j + fcj-i)(3fc,-_i 3 l - ' (3/cy_i + kj)(kj ~~ - kj-i) 16kj{kj + ) ' / , yy+1/2) is a p p r o x i m a t e d u s i n g t h e c o m m o n f o u r - c o r n e r f o r m u l a 2 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 u s e d i n p l a c e of t h e a r b i t r a r y f u n c t i o n / i n e x p r e s s i o n ( 3 . 9 ) , a f o r m u l a for U x 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 Ri+1/2,3 . _Ui+i,j + 1/2 - + Ui+i,j+i ^ - Uj,j+i ^M i+i/2,yj+i/2){Ui,j + Ui+ij x i = l,..,n ~ Ui,j j' = + Uij+i + Ui i + i j + i) (3.12a) 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) b y u s i n g U i n p l a c e o f / i n e x p r e s s i o n (3.10) ^ _Uj,j+i Vt+i/2,y+i/2 - + UJ+IJ+I - UJ+IJ ^ ^B{x ,y 2){Ui j i+1/2 j+1/ Ujj + U j'+Ui j > i+1> i = l,..,n j = + U ! +1 i + l t j + i) (3.126) l,..,m. T o d i s c r e t i z e ( 3 . 1 1 c ) , R is u s e d as t h e a r b i t r a r y f u n c t i o n i n (3.9) a n d Q i n ( 3 . 1 0 ) . E q u a t i o n (3.11c) is t o b e d i s c r e t i z e d o v e r b o x [ ( x _ / , y y - i / 2 ) > ( ^ + 1 / 2 , J/y+i/2)]? b u t I 1 2 e x p r e s s i o n s (3.9) a n d (3.10) w e r e d e r i v e d f o r b o x [ ( x ; , y y ) , ( x + i , y y + i ) ] - T h e r e f o r e the t e x p r e s s i o n s (3.9) a n d (3.10) m u s t b e s h i f t e d t o t h i s n e w b o x t o b e u s e d i n t h e d i s c r e t i z a t i o n o f (3.11c) f (xi,yj)= (3.9a) x [ / ( S j + l / 2 , y j - l / 2 ) + / ( S t + l / 2 , i / j + l / 2 ) ~ / ( S t - 1 / 2 . yy + 1/2) ~ / ( ^ - l / 2 , y y - l / 2 ) ] 2{l/2h i + l/2h - ) i 42 1 Q , F E [ 2^ ' (3.10a) f {xi,yj)= y [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) , + 0{h ). 2 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 ~hfcy— i C(£ -,yy) ^ y + D ^ y y ) (3.12c) t t' = 2,..,n, j ' = 2,..,m. Now (3.12a) can be used as a definition for -R;+i/2,y+i/2; for Qi+i/ y+i/22) a n d (3.12b) as a definition 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) 2 +hi - i A ( x _ 2 , y i _ i / ) i -i.y-1 1/ 2 +/c _iB(x _ 2,yi_i/ ) 4fcy_ (fcy + /cy-i) 2 J 4hi-i(h{ + /ii_i) t 1/ 2 1 C(x;,yy) cti a 3 + "4 + hi-iA(xi^ / yi-i/2) + 4/ii_i(/i -+/ij-_i) hi-iA(x _ ,yi i/ ) 1 2i i 1/2 + 2 t 2fcy + 2fcy_i - fcyfcy_i/J(x _ /2,y -i/ ) t 1 l 4kjkj_i(kj 2 +feyfcy_17J(x _ t + fcy_x) C(x;,yy) a /?y] t + 2 +/i _ ^(x _ i 1 i 1/2 ,y i+1 / ) 2 2- k B{x _ ,y /2) j i 1/2 4kj(kj + A;y_i) 4h{-i(hi + /j-i-i) C(xi,yy) a; 7y] + 44 i+1 1 / 2 ,y l+ X /) 2 -2/tj-! - 2hj - h h _ A(x / ,yi-i/ ) i Ui,j-1 4 + i 1 i+1 2 + 2 h h _ A(x ^ ,y _ ) i i 1 i l/2 t 1/2 + fcj-i-B(g - / ,y --i/ ) + fcj-i-B(gt-i/2»yi-i/2) 4fcy_i(fcy + fcy_i) t +1 2 t 2 C(x -,yy) /?»• ay] t + [(-4/tj_i-4/irMi-i) x Uj it {A(gj i/ ,y - i/ ) + A ( x + 2 t + 2 t + 1 / ,y _ 2 t 1 / 2 ) - A(gj_i/ ,y - / ) - A ( x _ 2 4hihi-i(hi (—4fcy_i — 4fcy — kjkj-i) {jj(x i + 1 / 2 ,y i + 1 / 2 t +1 2 t 1 / 2 , yi-1/2)} + fot_i) x ) - g(z,-+i/ ,y,--i/ ) + B ( a - _ 2 2 t 1/2 ,y t +1/2 ) - -B(x^_ ,y _ )} 1/2 { 1/2 4 A;y /ii -y (fcy +fcy_ 1) C(x -,yy) ft /?y t *7,»,y+i -2hj-_i - 2fe- - fetfej-i^(gt-i/2> yt+1/2) + 4hihi-i(hi + hi^i) hihi-iA(x f ,y / ) t 2 - A:yJ3(zf_i/ ,y.-+i/ ) + kjB{x / , 2 2 i+1 2 y»+i/ ) i+1 2 2 yy) ft Ty 2 A/j ^ A>y -h" A*y — 1 ^ + 2 - hiA(xj ,yi_ ) +1/2 i+i,y-i 1/2 4/ii(/ii + hi-i) 2+ fcy-i-B(x ,y -_ ) i+1/2 4kj-i(kj C(x»,yy) T i «y] + 45 + fcy_i) t 1/2 i+1 2 + 2 - hiA(x /2,y i/2) ~ i+ i+1 h A(xi / ,y _ / ) i 2kj + 2fcy_i - k kj- B(x / ,y - ) j 1 i+1 2 i +1 2 i l 2 + 1/2 Akjkj-\(kj kjkj-iB( / ,y / ) Xi+1 2 i+1 2 + fcy-i) C(^t,yy) 7i /?y] + 2 - /i vl(x t 17* +i,y+i l + 1 /2,yi i/2) 2- + kjB(x / ,y / ) i+1 2 i+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 guaranteed 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 T h = 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 Results The nine-point formula (3.13) was used to find approximate solutions to four twodimensional boundary-value problems with Dirichlet boundary conditions. The exact solutions to all four problems are known. 1) U +U xx yy = U + exp(x + y ) , xe[o,i], y e [o,i], U{0,y)=ev, U{l,y)=ev , +1 U[x,0) = e , x U(x,l) =e solution 2) U x + \ U = exp(x + y ) . + Uyy = -xU xx x E [0,1], -yU x ye y + 3x - 3 y , 2 2 [0,1], U(0,y)=y , 2 U{l,y) = l + y, 2 U(x,0) = x , 2 U[x l) t = 1 +x , solution 3) U xx 2 U = x +y . 2 + Uyy = -l/xU xG[l,2], 2 x + xyUy - xU, yG[l,2], 47 U{l,y)=0, C/(2,y) = y l o g 2 , U(x,l) = logx, U(x,2) = 2 logx, solution 4) U xx U = ylog(x). +U yy = xe[i,2], U(l,y) = TT (X + 2 y)U x - j^jUy - Tr U 2 - Tr C O S ( T T X + 2 Try), y e [0,2], ^ ^ , tf(2,y) = ^ g ± ^ , U(x,0) = -^, s C7(x,2)-5iEl2+xl 2+x solution U sin(x + y) x+ y To create a grid of mesh points, a one-dimensional x-mesh and a one-dimensional ymesh 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 = MAX STEP 49 SIZE 2.14 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 4.1 E x p o n e n t i a l P E R T U R B E D P R O B L E M S Fitting 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 are c h a r a c t e r i z e d b y t h e o c c u r r e n c e o f a s m a l l p a r a m e t e r e t h a t m u l t i p l i e s the highest derivative 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 solutions 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 slowly i n others. T h e r e g i o n s o f r a p i d v a r i a t i o n are c a l l e d layers; t h e o t h e r r e g i o n s are s m o o t h r e g i o n s . Because of the c o m p l i c a t e d nature of their differential operators a n d solutions, singul a r l y p e r t u r b e d p r o b l e m s are t r e a t e d s e p a r a t e l y f r o m t h e r e g u l a r p r o b l e m s o f C h a p t e r 2. Consider the singularly perturbed O D E ey" = {r{x) y y + q{x)y + f{x) (4.1) a < x < b w h e r e e is a s m a l l parameter, 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 f u n c t i o n s as i n (1.1a). T h e s i g n o f r ( x ) o v e r t h e i n t e r v a l [a, b] reveals a l o t a b o u t t h e s o l u t i o n of t h e d i f f e r e n t i a l e q u a t i o n . If r ( x ) < 0 o v e r t h e e n t i r e i n t e r v a l , t h e n t h e o r d i n a r y d i f f e r e n t i a l e q u a t i o n h a s t w o 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, z'=q{x)y + f{x), (4.2a) (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 just as equations (2.1a) and (2.1b) were discretized for the regular n 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 Xj-i/2 by y using quadratic interpolation. t 55 and x I + 1 / , and y(x ) is approximated 2 t O n a u n i f o r m mesh this becomes I~ = 21+1/2 Zt 1/2 *(*,) V, + f(x,) (4.36) i = 2,..,n, w h e r e h is t h e u n i f o r m s t e p s i z e . O r i f t h e t r a p e z o i d a l s c h e m e is u s e d o n a n a r b i t r a r y m e s h t h e d i s c r e t i z a t i o n of (4.2b) is z i+l/2 l/2{h ~ i-l/2 _ + h_) ~ z i i l J 9 ( ^ - 1 / 2 ) {yi + yt-i) + ^ /(«t-i/2)+ 2 ? ( i + V 2 ) (yi+1 + yi) + \ /(^t+l/2)X (4.3c) i — 2 , . . , n. O t h e r s e c o n d - o r d e r difference s c h e m e s c a n a l s o b e a p p l i e d , j u s t as i n t h e case o f t h e regular one-dimensional problem. E q u a t i o n ( 4 . 2 a ) , o n t h e o t h e r 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 t h e fast s o l u t i o n c o m p o n e n t o f p r o b l e m (4.1). E x p o n e n t i a l fitting is u s e d i n s t e a d o f t h e u s u a l c e n t e r e d difference s c h e m e s . I n t h i s d i s c r e t i z a t i o n t h e f u n c t i o n s r a n d z are a p p r o x i m a t e d b y piecewise constant functions. r = z on r{x ) i+1/2 = Zi+1/2 (-) 4 X{ < x < X{+i, 4 i=l,..,n. T h e s e a p p r o x i m a t i o n s are r e a s o n a b l e s i n c e r ( x ) is s m o o t h a n d b o u n d e d a w a y from z e r o , a n d s i n c e z varies s l o w l y . O n c e t h e s e a p p r o x i m a t i o n s are m a d e , a n 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 exp[e _ 1 r(x l + 1 / ) 2 {x - x,-)] giving EU = e x p [ - e _ 1 r(x i + 1 / 2 56 ) (x - x -)] z / t i+1 2 (4.5) = l,..,rc. X{ < X < £ - + i , t (4.5) c a n b e i n t e g r a t e d e x a c t l y , / eu' dx= exp[-e" r(x 1 J Xi J i + 1 / 2 z / dx ) (x - x ) ] t i+1 2 Xi to give e [u(x -+i) - u{xi)) = 1 + 1 / 2 [exp{-e t e t y i exp[-e: _ 1 t + 1 / 2 ) t + 1 / 2 ) hi , a n d u(xj-) + r(x r(x - 1 (x i + 1 - x ; ) } - 1]. (4.6) lZi+i/2j x - i — Xj- b y i t s e q u i v a l e n t Replacing i + r b y yi, a n d u ( x i i ) + by /i;] i n (4.6) y i e l d s ey i exp[-e i+ _ 1 EZt+i/5|| r(x e x p / ) l + 1 j_ 2 hi] - ey = { ^(x- g 1 / 2 ) /i^] - l } . (4.7) E q u a t i o n (4.7) c a n b e s i m p l i f i e d b y u s i n g t h e B e r n o u l l i f u n c t i o n x B(x) e - 1 x T h e d i s c r e t i z a t i o n f o r e q u a t i o n (4.2a) u s i n g e x p o n e n t i a l f i t t i n g w i l l b e Z l + 1 / 2 = ^r{x B h~ N o w t o get a three-point )hi\ B /-r(x i + 1 / 2 ) h, (4.8) formula « i yi-i + bi yi + Ci y (4.8) is s u b s t i t u t e d ^ i+1/2 l + i = gi, i = 2,.., n, into t h e discretization for (4.2b). If t h e 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 u s e d , t h e c o n s t a n t s ai, bi, Ci, gi a r e g i v e n b y eB (- ( -j»/»> ) r _ £ x h £ f-r(xj- )h^ bi = h £ B ^r(x Ci = ^ 1/2 t + , h / 2 2 ~ ^r(x )h^ i+1/2 h 2 q{xi), )h^ 2 (4.9a) 9i = f{xi). 57 O n a n a r b i t r a r y mesh this becomes 2sB (~ -^ - ) r(z 2)fe 1 - 2eB hi-^hi b l r 1 + ki{hi + fc -_!) t fc -_i) ^ t A ' ( | + I/2)fc.' I , /it(/j-t + hi-i) ft = f[xi). (4.96) N o t e t h a t (4.9b) r e d u c e s t o (4.9b) o n a u n i f o r m m e s h . I f t h e t r a p e z o i d a l s c h e m e (4.3c) is u s e d o n a n a r b i t r a r y m e s h « f l ( = ' - - y r ~ _ ) t ' - ) h~, J5 ^- ( .-i/2)fe;-i^ r e i-°(*.-l/2). I E g ^r(x i+1/2 )hj j 6,- = hi—i . . hi . —^—9(^1-1/2) - I c* = J T 7^ -jg(zi+i/2), /ii -jQ{x / ), i+1 2 Qi = ^ f / ( x _ i / ) + ^ / ( x i l 2 58 l + 1 / 2 ). (4.9c) 4.2 D i s c r e t i z a t i o n Errors 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 / into (4.8). 2 r{x i/ )h\ e y{x i/2) + ^y'(x / ) + o(h ) e J h i+ Th =B 2 2 i+ B i+1 y{x i/2) - ^y'{x i+ h ) - r(x / ) y{x [ey'(x i+1/2 i+1 2 i+1/2 2 ) + o{h ) 2 i+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 B 'T{xi+\/ )h 2 s « 1. Hence (4.10) becomes r » - B -r{x i/ )h\ E i+ y{xi+i/ ) - -y'{x 2 2 h ) + o{h ) 2 i+1/2 £y'{x / ) + r{xi )y(x ) i+1 2 +1/2 i+1/2 which is equivalent to / h r{xj+i/ ) E 2 e V x p h f - r ^ , , ) h,, y{x i/2) - -v {xi+1/2) + o{h ) l+ ^y'{xi ) +1/2 +r(x i + 1 / 2 ) y(x i + 1 / 2 ). Since -r{xi+i/ ) 2 exp Th = -r{x / ) i+1 2 y{x / )-r i+1 2 h r{x i+1/2 « 1, ) , y'{xi i/2)+r(x 2) y{xi i/ ) + 0(h ) + O{E), + 59 i+1/ + 2 and so the local truncation error is given by r = 0{h). 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 Results f o rE x p o n e n t i a l Fitting 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 tested 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 Dirichlet b o u n d a r y conditions a n d k n o w n solutions. 1) e y" = —y' + y x<E[0,l], (1 + ETV ) y(0)=0, 2 cos(7rx) — 7rsin(7rx) + e x p ( — x / e ) , y(l) = -i-exp{=±), cos(7rx) — e x p ( — X / E ) . solution 2) — £ y" = [2 + c o s ( 7 r x ) ] y ' + y — (1 + sir ) 2 e y" = { e - l ) y ' xG[0,l], solution cos(7rx) — 7r[2 + c o s ( 7 r x ) ] s i n ( 7 r x ) - | - + y, y(0)=2, y — exp(x) + y{l) = e + exp{=^), exp(—X/E). 4) £y" = y' + 6, xG[0,l], solution y(0) = l , y(l) = exp{±)-6, y = exp(x/£r) — 6x. N o t e t h a t f o r p r o b l e m 4) r ( x ) a n d z(x) a r e c o n s t a n t s : r(x) = - 1 , z(x) = 6. N u m e r i c a l r e s u l t s a r e c a l c u l a t e d u s i n g t h e 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 (4.9b). L o c a l e r r o r s a r e f o u n d a t t h e m e s h p o i n t s , t h e n t h e m a x i m u m a n d average e r r o r s a r e c a l c u l a t e d . T h e m a x i m u m e r r o r is p l o t t e d a g a i n s t t h e m a x i m u m s t e p s i z e , a n d t h e s l o p e o f t h e l i n e a r l e a s t - s q u a r e s f i t is d e t e r m i n e d . T h i s s l o p e gives t h e a c c u r a c y o f t h e m e t h o d . 61 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 Figure 18 Results for Singular Perturbation Problems 1,2 on Uniform Meshes slope = 1.22 o o o o O cr o cn o o o o LUo o X cr O o <J> o° i 10- o o i i i i—i—i—r~i 10-' ? MAX STEP 63 1 SIZE 1 1 1—i—i—i—i—i 10° 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- io-' J M A X S T E P 65 io° S I Z E Figure 21 Results for Singular Perturbation Problems 1,2,3,4 with K « £ slope = 1.92 -i -i—i—i—r-T io- 2 10- M R X S T E P 66 1—i—r i i "I 10° S I Z E 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 considered i n this thesis for the solution of the linear differential e q u a t i o n 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 ( 1 . 1 ) , ( 3 . 1 ) , a n d (4.1) m a y b e e x t e n d e d t o 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 . C o n s i d e r first a q u a s i - l i n e a r b o u n d a r y - v a l u e p r o b l e m a < x < 6, y" = f{y',y,x), y(a) = 6, (5.1) 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 b y a n y o f t h e 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 , s e c o n d order approximations t o isolated solutions c a n be found by the methods described i n S e c t i o n 2 . 1 . I s o l a t e d s o l u t i o n s are s o l u t i o n s w h i c h are l o c a l l y u n i q u e . H i g h e r - o r d e r e q u a t i o n s c a n also b e s o l v e d u s i n g t h e s e m e t h o d s . An m th order d i f f e r e n t i a l e q u a t i o n w i l l b e e x p r e s s e d as m f i r s t - o r d e r d i f f e r e n t i a l e q u a t i o n s f o r t h e purpose of deriving the numerical scheme. Difference formulas are a p p l i e d t o each d i f f e r e n t i a l e q u a t i o n . A s t a g g e r e d m e s h Ti is u s e d ; t h a t i s , e a c h d i f f e r e n t i a l e q u a t i o n n is d i s c r e t i z e d o n different s u b - i n t e r v a l s o f t h e i n t e r v a l [a,b] i n s u c h a w a y t h a t ( m — 1) consecutive substitutions c a n be performed to produce a n ( m + l ) - b a n d e d system A HIT = g - F o r e x a m p l e , a l i n e a r f o u r t h - o r d e r d i f f e r e n t i a l e q u a t i o n y ( 4 ) = \ci(x) y ] ' " + [ c ( s ) y]"+[cs{x) y}' + c (x) 4 a y + c (x) 5 (5.2) c a n b e r e w r i t t e n as t h e s y s t e m y ' =ci(x) y + p, 67 (5.3a) p' =c {x) y + q, (5.36) q' = c ( x ) y + r, (5.3c) r' =c (x) y + c (x). (5.3d) 2 3 4 5 A m i d p o i n t scheme o n a u n i f o r m b u t staggered m e s h c a n be a p p l i e d to the four firsto r d e r e q u a t i o n s o n t h e f o l l o w i n g i n t e r v a l s : e q u a t i o n s (5.3a) a n d (5.3c) are d i s c r e t i z e d over [x,-,Xi i], i = 1, + ..,n a n d e q u a t i o n s (5.3b) a n d (5.3d) a r e d i s c r e t i z e d o v e r [xi-i/ ,x / ], 2 So »' = 2 , . . , n . i+1 2 then V l + \. id V l c , , x f - = ~ e (Xi) 2 = ^ c (x q % 3 ^+1/2-^-1/2 An [Vi+i + Vi] = \ i( i+i/2) 1 i + 1 / 2 y i +P1+1/2. \ 5 + q, , 4 A ) (5.46) { ) [y,-+i + yi] + r f ( - ( I + 1 / 2 , \ (5.4c) tz. AJ\ ( m + l ) - p o i n t s c h e m e i n t h e u n k n o w n s y- is d e t e r m i n e d b y ( m — 1) t successive s u b s t i t u t i o n s . F o r t h e 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 t o (5.4b) t o e l i m i n a t e t h e v a r i a b l e p. T h e n t h e 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 t o (5.4c) a n d q is e l i m i n a t e d . F i n a l l y t h e n e w e q u a t i o n is s u b s t i t u t e d i n t o e q u a t i o n (5.4d) t o e l i m i n a t e r . This yields 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 y n for u n k n o w n s . W h e n a difference s c h e m e is a p p l i e d t o (5.3) o n a n o n u n i f o r m m e s h , t h e d i s c r e t i z a t i o n s 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 following d i s c r e t i z a t i o n o f t h e s y s t e m (5.3) : V l + 1 h . V l = \ i ( i + i / 2 ) b i + i + yi] + Pi+1/2 c x 68 (5.5a) w, , , 2 (ft, + r- = C 2 /li-l) ( ^ x y i + q ( i q Xj ^ = c (xi) 4 y i x 5 The l + x ; is t h e m i d p o i n t o f [ x ; i , x , - ] , x 2 [x (5.5d) + c (xi) Xi w h e r e x ; is t h e m i d p o i n t o f [ i-i/ , i+i/ }, the midpoint of > (5.5c) _ ~ = cs(xi) yi + fi Xt+l 5 M 2 + X ; is i,Xi]. m e t h o d s u s e d i n t h i s t h e s i s c a n also b e a p p l i e d t o 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 boundary conditions that are more complicated than the Dirichlet boundary conditions considered i n t h e previous chapters. Separated linear b o u n d a r y conditions of t h e f o r m ki y'(a) + k y{a) = 6 2 X £iy(b)+£ y'(b)=6 2 2 can very simply be included i n the system A y„ = g (1.2) w i t h o u t affecting t h e 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 . However, i n c o r p o r a t i n g b o u n d a r y c o n d i t i o n s t h a t a r e n o t s e p a r a t e d i n t o t h e s y s t e m (1.2) w i l l c a u s e t h e s y s t e m t o lose i t s t r i d i a g o n a l f o r m . Still, a similar stability a n d consistency analysis c a n b e s h o w n t o h o l d h e r e as f o r (1.2), a n d s e c o n d - o r d e r a c c u r a t e s o l u t i o n s a r e 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 c e n t e r e d s e c o n d - o r d e r difference s c h e m e s a r e u s e d . 69 6. C O N C L U S I O N W e h a v e d e r i v e d a class o f difference s c h e m e s t o find s e c o n d - o r d e r approximate solutions to b o u n d a r y - v a l u e problems by discretizing o n an a r b i t r a r y mesh n . n The d e r i v a t i o n o f t h e s e s c h e m e s i n v o l v e s , first, r e w r i t i n g t h e b o u n d a r y - v a l u e p r o b l e m as a system of first-order equations. 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 o n e - s t e p c e n t e r e d difference s c h e m e o v e r s t a g g e r e d s u b - i n t e r v a l s o n t h e m e s h n . n T h e first e q u a t i o n i n t h e s y s t e m , i n v o l v i n g t h e d e r i v a t i v e o f t h e 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 z e d at t h e m e s h p o i n t s ; t h e n e x t e q u a t i o n , w h i c h h a s t h e d e r i v a t i v e 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 z e d at t h e m i d p o i n t s o f t h e m e s h . A n y f u r t h e r equations are d i s c r e t i z e d at p o i n t s m i d - w a y b e t w e e n t h e d i s c r e t i z a t i o n p o i n t s o f t h e p r e c e d i n g equation. T h i s use o f s t a g g e r e d m e s h s u b - i n t e r v a l s a l l o w s t h e s u b s t i t u t i o n o f o n e e q u a t i o n i n t o t h e n e x t so t h a t , u l t i m a t e l y , o n e e q u a t i o n i n v o l v i n g o n l y t h e p r i n c i p a l mesh function results. For an m t / l -order boundary-value problem, an ( m + l)-point difference f o r m u l a is p r o d u c e d i n t h i s w a y . S i m i l a r l y , a n (rn + l ) - p o i n t f o r m u l a c a n 2 b e o b t a i n e d for b o u n d a r y - v a l u e p r o b l e m s p o s e d i n t w o d i m e n s i o n s . The difference f o r m u l a s c a n also b e a p p l i e d t o l i n e a r i z e d f o r m s o f n o n l i n e a r e q u a t i o n s . S i n g u l a r l y p e r t u r b e d problems w i t h o u t t u r n i n g points c a n also be solved using s i m i l a r p r o c e d u r e s , b u t s o m e r e s t r i c t i o n s 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 t h e 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 boundaryvalue 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 References 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 FiniteDifference 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. 791802. 7. H . Kreiss, T. Manteuffel, B . Swartz, B . Wendroff, & A . B . White Jr., "Supraconvergent 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," M a t h . Comp., Vol. 28, 1974, pp. 981-1003. 9. T. Manteuffel, & A . B . White Jr., "The Numerical Solution of Second Order Boundary Value Problems on Nonuniform Meshes," M a t h . Comp., Vol. 47, 1986, pp. 511535. 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 Nonuniform 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, Wesley Pub., London Ont., 1984. Vol 1,2, Addison
- 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 |
Aggregated Source Repository | 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