CHEBYSHEV SPECTRAL METHODSFOR POTENTIAL FIELD COMPUTATIONByShengkai ZhaoB. Sc. (Electrical Engineering) University of Science and Technology of China, 1970M. Sc. (Geophysics) University of Science and Technology of China, 1981A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFDOCTOR OF PHILOSOPHYinTHE FACULTY OF GRADUATE STUDIESDEPARTMENT OF GEOPHYSICS AND ASTRONOMYWe accept this thesis as conformingto the required standardTHE UNIVERSITY OF BRITISH COLUMBIAAugust 1993©Shengkai Zhao, 1993In presenting this thesis in partial fulfilment of the requirements for an advanceddegree at the University of British Columbia, I agree that the Library shall make itfreely available for reference and study. I further agree that permission for extensivecopying of this thesis for scholarly purposes may be granted by the head of mydepartment or by his or her representatives. It is understood that copying orpublication of this thesis for financial gain shall not be allowed without my writtenpermission.(Signature)‘ic 5 404A s-60 ,n0 nljDepartment of 4'et2The University of British ColumbiaVancouver, CanadaDate Avty. 3 DE-6 (2/88)ABSTRACTWe investigate Chebyshev spectral methods for solving Poisson's equation and thegeneralized Poisson's equations and apply them to 3-D gravity and DC resistivity modeling.When we approximate a function by a finite sum of Chebyshev polynomials, thereare two kinds of errors: truncation and aliasing. We present a cell-average discretizationscheme to reduce the aliasing error. Both theoretical analysis and numerical examples showthat when there are discontinuities in a function, the cell-average discretization is betterthan the point injection discretization.We use both the r and collocation methods to solve Poisson's equation V 2u = f. Wesolve the discrete systems by a matrix-diagonalization method. The speed and accuracy ofcollocation methods are better than those of the r method. For the generalized Poisson'sequation V • (oVu) = f, we present a new iterative method. We rewrite the equation in theform of a Poisson's equation, V2u = (f —Va•Vu)la. At each iteration we compute theright hand-side term from the current value of u first. Then we solve the resultant Poisson'sequation by the collocation method. Numerical results show that the convergence rate ofthe new method is much faster than that of the spectral multigrid method.When there are discontinuities in the source function f and/or the conductivity o, thesingle-domain Chebyshev spectral method does not converge exponentially. In this casewe use the multi-domain Chebyshev spectral method to solve the problem. We divide thewhole domain into a number of subdomains so that in each subdomain the function isinfinitely differentiable. Then we approximate the function by a separate Chebyshev seriesin each subdomain. We determine the relations between the Chebyshev polynomials inadjacent subdomains by the interface conditions. In this way we can achieve exponentialconvergence.In 3-D gravity modeling, we present a method to realize 2-D and 3-D point-injectionand cell-average discretizations, use a multipole expansion to compute the approximateiiboundary conditions and solve the equations by both single-domain and multi-domainChebyshev spectral methods. We also extend Okabe's analytic formulation for the gravi-tational field of a homogeneous, polyhedral body to the potential. Numerical results showthat the accuracy of the cell-average discretization is better than that of the point injectiondiscretization. The multi-domain solution is the best.We apply the multi-domain Chebyshev spectral method to 3-D DC resistivity modeling.We discuss the singularity removal and present a modified finite-difference method forcomparison. For a two-layered model, the multi-domain Chebyshev spectral solution is farmore accurate than the finite-difference solution. For piecewise-constant and piecewise-smooth models the solutions obtained by both methods agree with each other quite well.However, the Chebyshev spectral method is more efficient than the finite-difference method.iiiTABLE OF CONTENTSAbstract ^iiTable of Contents ^ivList of Tables ^viiList of Figures ^ixAcknowledgement ^xiChapter 1 Introduction ^11.1 Numerical Methods Review ^11.2 Motivation ^41.3 Outline of This Thesis ^5Chapter 2 Some Fundamentals of Chebyshev Expansions ^82.1 Approximation of a Function by Chebyshev Polynomials ^82.2 The Discrete Chebyshev Expansions ^132.3 The Error of Chebyshev Expansions ^172.4 Gibbs Phenomenon ^202.4.1 Cell-Average Discretization ^222.4.2 Approximation of a Function by Separate Chebyshev Seriesin Different Subdomains ^262.5 Discussion ^27Chapter 3 Solving Poisson's Equation, Part I: Single-Domain Methods^293.1 The r Method ^303.1.1 The Relations Between the Chebyshev Coefficients of a Functionand its derivatives ^303.1.2 The 1-D Problem ^313.1.3 The 2-D Problem ^343.1.4 The 3-D Problem ^36iv3.2 The Collocation Method ^383.2.1 The First and Second Order Derivative Matrices ^383.2.2 The 1-D Problem ^403.2.3 The 2-D Problem ^423.2.4 The 3-D Problem ^433.3 Discussion ^443.3.1 Accuracy ^443.3.2 Comparison of the r and Collocation Methods ^44Chapter 4 Solving Poisson's Equation, Part II: Multi-Domain Method^464.1 The Influence Matrix Method ^474.2 The 1-D Problem ^504.3 The 2-D Problem ^574.4 The 3-D Problem ^624.5 Discussion ^68Chapter 5 Solving the Generalized Poisson's Equation, Part I:Single-Domain Method ^ 705.1 Spectral Radius and Convergence Rate ^715.2 A New Iterative Method ^725.2.1 The 1-D Problem ^735.2.2 The 2-D Problem ^785.2.3 The 3-D Problem ^825.3 Discussion ^855.3.1 Convergence Rate ^855.3.2 Comparison with Chebyshev Spectral Multigrid Method ^865.3.3 Computational Operations ^86Chapter 6 Solving the Generalized Poisson's Equation, Part II:Multi-Domain Method ^ 886.1 The 1-D Problem ^886.2 The 2-D Problem ^946.3 The 3-D Problem ^986.4 Discussion ^ 104Chapter 7 3-D Gravity Modeling ^ 1057.1 The Boundary Value Problem 1067.2 Point-Injection and Cell-Average Discretization ^1077.3 The Analytic Expression of the Potential Produced by a HomogeneousPolyhedral Body ^ 1127.4 Approximate Boundary Condition ^ 1157.5 Numerical Examples ^ 1167.6 Discussion ^ 1277.6.1 The r and Collocation Methods ^1277.6.2 Cell-Average and Point-Injection Discretization ^1287.6.3 Single-Domain and Multi-Domain Solutions ^ 1287.6.4 Approximate and Accurate Boundary Conditions 129Chapter 8 3-D DC Resistivity Modeling ^ 1308.1 The Boundary Value Problem ^ 1318.2 Removal of the Source Singularity ^1338.3 The Multi-Domain Chebyshev Spectral Method ^ 1358.4 The Finite-Difference Method ^ 1388.5 Numerical Examples ^ 1458.6 Discussion ^ 159Chapter 9 Summary 162References ^ 166viLIST OF TABLESThe error norms of example 4.2 (MDCS method) ^56The error norms of example 4.2 (SDCS method) ^56The error norms of example 4.3 (MDCS method) ^61The error norms of example 4.4 (MDCS method) ^67The error norms of example 5.1 after 9 iterations on an N = 32 1-D grid^77The error norms of example 5.2 after 9 iterations on an N = 32 2-D grid^81The error norms of example 5.3 after 9 iterations on an N = 32 3-D grid^83The relative error norms of example 6.1 (MDCS method) ^93The relative error norms of example 6.2 (MDCS method) ^97The relative error norms of example 6.3 (MDCS method) ^104The Lc,. relative error norms of example 7.1 (r method) ^ 118The L2 relative error norms of example 7.1 (r method) 118The Lc*, relative error norms of example 7.1 (collocation method)^119The L2 relative error norms of example 7.1 (collocation method) .^119The relative error norms of example 7.1 (MDCS method, accurateboundary conditions ) ^ 1207.6^The Lco relative error norms of example 7.2 (r method) ^ 1227.7^The L2 relative error norms of example 7.2 (r method) 1227.8^The L03 relative error norms of example 7.2 (collocation method)•^1227.9^The L2 relative error norms of example 7.2 (collocation method) .^• 1237.10 The coordinates of the vertices of example 7.3 ^1247.11 The Lc° relative error norms of example 7.3 (r method) ^ 1267.12 The L2 relative error norms of example 7.3 (r method) 1267.13 The Lc() relative error norms of example 7.3 (collocation method) • . ^ 1267.14 The L2 relative error norms of example 7.3 (collocation method) . • • • 1274.14.24.34.45.15.25.36.16.26.37.17.27.37.47.5vii8.1 The error norms of example 8.1 (case 1, MDCS method) ^ 1488.2 The error norms of example 8.1 (case 1, FD method) 1498.3 The error norms of example 8.1 (case 2, MDCS method) ^ 1508.4 The error norms of example 8.1 (case 3, MDCS method) 1528.5 The results of example 8.1 (case 4, MDCS method) ^ 1548.6 The results of example 8.1 (case 4, FD method) 154viiiLIST OF FIGURES2.1 The Chebyshev coefficients of cos(rx), ex and 1/(1.1 — x) ^ 112.2 The function g(z) = 1/(a — x) on the interval —1 < x < 1 ^ 122.3 A 1-D Chebyshev grid (N = 4) ^ 142.4 A 2-D Chebyshev grid (N = 4) 152.5 A 3-D Chebyshev grid (N = 4) ^ 162.6 The continuous Chebyshev expansion approximation of sgn(x) obtainedby truncated series ^ 212.7 The Chebyshev expansion approximation of sgn(x) obtained by pointinjection discretization ^ 222.8 The Chebyshev expansion approximation of sgn(x) obtained by cell-average discretization ^ 254.1 1-D subdomains 504.2 The functions fi (x) and u i (x) ^ 554.3 2-D subdomains ^ 584.4 3-D subdomains 625.1 The L. norm, ems , of the differences of u between two successiveiterations of example 5.1 ^ 765.2 The spectral radius of the iteration matrix for example 5.1 ^ 775.3 The L. norm, emax , of the differences of u between two successiveiterations of example 5.2 ^ 805.4 The spectral radius of example 5.2 ^ 815.5 The L. norm, em., of the differences of u between two successiveiterations of example 5.4 ^ 846.1 The functions f3(x), Q3(x) and u3(z) ^ 927.1 Normal vector and the direction of a closed curve ^ 108ix7.2 Coordinate system rotation ^ 1137.3 The model of example 7.1 1177.4 Subdomains of example 7.1 ^ 1207.5 The model of example 7.2 1217.6 The model of example 7.3 ^ 1247.7 The relative error of example 7.3 on z = 0 plane (N = 32) ^ 1258.1 DC resistivity method configuration ^ 1328.2 A 3-D finite-difference grid ^ 1388.3 The rectangular parallelepiped Vijk for an interior point ^ 1408.4 A integration domain Vijk which contains discontinuity ^ 1428.5 The rectangular parallelepiped Viii for a point on the ground surface 1448.6 A two layered model ^ 1468.7 A dipole—dipole configuration on a two-layered model ^ 1538.8 A two-prism model ^ 1558.9 The secondary potential on the ground surface of example 8.2 ^ 1568.10 The secondary potential on the ground surface of example 8.3 ^ 158ACKNOWLEDGEMENTDedication: To my parentsI express my deepest thanks to my supervisor Professor David W. Strangway, thepresident of UBC, for providing financial support and many encouragements for most ofmy Ph.D. program. Without his help I could not have finished my studying. I also expressmy deepest thanks to my supervisor Professor Matthew J. Yedlin for his invaluable advicein every aspect of this work and financial support.I thank Drs. Yaogou Li, Dave Aldridge, David Dalton, Robert Ellis, and David Lumleyfor their help and suggestions. I thank John Amor for his help and expertise with the com-puting facilities. I'm grateful to Drs. R. Don Russell, Douglas Oldenburg, and Uri Ascherfor acting as my committee members. I thank all my friends in the Department of Geo-physics and Astronomy who have made my study so enjoyable. This work was supportedby Natural Science and Engineering Research Council of Canada (NSERC) operating grant5-80642 and Professor David W. Strangway's UBC research grant.xiCHAPTER 1INTRODUCTION1.1 Numerical Methods ReviewIn geophysics, gravity and direct current (DC) resistivity are important prospectingmethods. They reveal underground structural information from a variety of aspects. In agravity survey we measure the gravitational anomalies on the ground, air or sea and try todevelop a geological model which produces anomalies close to the observed data. In a DCresistivity survey we inject source current into the ground and measure the potential on theground surface. Then we try to find a model which produces a potential distribution closeto the observed data. In both the gravity and DC resistivity surveys, data interpretationrequires the solution of a forward problem. The gravitational potential obeys Poisson'sequation V2u = f, while the electrical potential obeys the generalized Poisson's equationV • (aVu) = f. In the two equations, f is the source function, u the potential, and a theconductivity. The forward problem requires solving these equations with the appropriateboundary conditions.For some simple models, such as the polyhedral body for gravity (Okabe 1979a) or thelayered model (Keller and Frischknecht 1966) for DC resistivity, there are analytic formulaeavailable. For more general models, numerical methods are necessary. For gravitationalmodeling the methods existing in the literature are numerical integration (Talwa,ni andEwing 1960, Talwani 1973), the Fourier transform (Bhattacharyya and Navolio 1975, 1976;Pedersen 1978; Hansen and Wang, 1988), and approximate methods (Mufti 1975). InDC resistivity modeling, the numerical methods reported in the literature are the finite-difference method (e.g. Dey and Morrison 1979a, 1979b, McGillivray 1992), the finite1CHAPTER 1: INTRODUCTION^ 2element method (e.g. Coggon 1971, Bibby 1978), and the integral equation method (e.g.Alfano 1959, Lee 1975).Since the 1970s, spectral methods have been used more and more for the solution ofpartial differential equations. Spectral methods use globally smooth functions (usuallya set of orthogonal polynomials) to approximate the unknown functions. Both Fourierand Chebyshev series are globally smooth functions. When used for solving differentialequations the Fourier series requires periodic boundary conditions. The Chebyshev seriesdoes not have any restrictions on the boundary conditions. The Chebyshev series is morecomplicated but more flexible than the Fourier series.Many authors have described the application of Chebyshev polynomials for solvingdifferential equations. Lanczos (1938) presented the r method, which expands a knownfunction as a finite sum of Chebyshev polynomials and correction terms. Clenshaw (1957)presented a direct method to solve 1-D homogeneous differential equations with Chebyshevexpansions. Gottlieb and Orszag (1977) gave a systematic description of the theory andapplications of the spectral methods. With the r approximation, Haidvogel and Zang (1979)solved the 2-D Poisson's equation on a square with homogeneous boundary conditions andproposed both alternating-direction implicit iteration and matrix-diagonalization methodsto solve the coefficient equation system. Horner (1982) solved the 2-D elliptic equationof the Laplace type with Dirichlet boundary conditions. He established an equation inthe form of a five-point finite-difference method and used an iterative method to solve thecoefficient equation system. Haldenwang et al. (1984) extended Haidvogel and Zang's 2-Dmethod to 3-D problems with general non-homogeneous boundary conditions.The multi-domain (also called domain decomposition, subdomain, substructure)Chebyshev spectral (MDCS) method and its related solution techniques have also beenextensively used to solve partial differential equations in the past decade. In the multi-domain Chebyshev spectral method, the whole domain is divided into a number of sub-domains. The unknown function is approximated by a separate Chebyshev expansion inCHAPTER 1: INTRODUCTION 3each subdomain. The connection of the Chebyshev expansions between two contiguoussubdomains is determined by the interface conditions. There are several reasons to use themulti-domain Chebyshev spectral method:1. We can convert a large problem to a number of small problems and solve them inparallel.2. When a boundary value problem is defined on an irregular domain, such as, a 2-DL shaped domain, we cannot simply expand the unknown function by a single set ofChebyshev polynomials. By dividing the whole domain into several subdomains, wecan approximate the unknown function in each subdomain by a separate Chebyshevseries provided the domain is not too irregular.3. When there are discontinuities in the source and/or the conductivity functions, therewill be discontinuities in the derivatives of the potential. In this situation the ac-curacy of the single-domain Chebyshev spectral method is low. With the multi-domainChebyshev spectral method we can obtain much higher accuracy.The multi-domain Chebyshev spectral method has been used to solve many kinds of equa-tions, such as elliptic and Navier-Stokes equations. There are also many varieties of thismethod. Delves and Hall (1979) introduced the global element method. They dividedthe whole domain into a number of subdomains. Then they mapped each curved subdo-main into a standard 2 x 2 2-D Chebyshev domains and used Chebyshev polynomials asthe finite-element shape functions. Orszag (1980) described a patching method to dealwith the interface conditions. Morchoisne (1984) developed a method based on overlap-ping multiple domains. Patera (1984) presented a spectral element method to solve 1-Dand 2-D Navier-Stokes equations, in which he divided the whole domain into rectangularsubdomains and used the Chebyshev polynomials as the finite-element shape functions.CHAPTER 1: INTRODUCTION^ 41.2 MotivationIn the finite difference, finite element, and integral equation methods, we use a low or-der local function (usually linear) to approximate the unknown function. The error of thosemethods is proportional to (1/N 2 ), where N is the number of grid points along one coordi-nate direction. To achieve high accuracy, a large number of grid points have to be used. InChebyshev spectral methods we use Chebyshev polynomials to approximate the unknownfunction. If a function is infinitely differentiable, the error of its Chebyshev expansion ap-proximation goes to zero faster than any finite power of (1/N) as N --> oo, where N is thehighest degree of Chebyshev polynomials used for the approximation (Gottlieb and Orszag1977). This is called exponential convergence. Because of the exponential convergence,Chebyshev spectral methods can achieve much higher accuracy than the finite-differencemethod with the same number of grid points. Equivalently, to obtain the same accuracy, theChebyshev spectral methods require far fewer grid points. Therefore Chebyshev spectralmethods axe more efficient than the finite-difference method. When there are discontinu-ities in the source and/or the conductivity functions, the single-domain Chebyshev spectralmethod cannot achieve exponential convergence. We can, however, use the multi-domainChebyshev spectral method to achieve the desired exponential convergence. In practicalproblems, such as DC resistivity modeling, we frequently need the function values at non-grid points. In the finite-difference and finite-element methods, this is accomplished bylocal linear interpolation. The accuracy is low. In the Chebyshev spectral methods, we cancompute the function values at any point by global Chebyshev interpolation, which is moreaccurate than the local linear interpolation. In this thesis we investigate how to solve thePoisson's and the generalized Poisson's equations by Chebyshev spectral methods. Thisincludes the r method, collocation method and multi-domain method. We also apply thesemethods to 3-D gravitational and DC resistivity modeling.CHAPTER 1: INTRODUCTION^ 51.3 Outline of This ThesisIn chapter 2 we outline some fundamentals of Chebyshev expansions as a basis forfurther discussions. When we approximate a function by a finite Chebyshev series, usuallythere are truncation and aliasing errors. When the degree of the Chebyshev expansionincreases, both the truncation and aliasing errors decrease. For a given function and afixed degree of the Chebyshev expansion, the truncation error is fixed. To reduce thealiasing error, we present a cell-average discretization scheme. We also present a theoreticalanalysis for the errors caused by the traditional point-injection discretization and cell-average discretization. By cell-average discretization, the aliasing error can be reduced tosome extent. But overall, the error norms are still proportional to (1/N). A much betterway is to divide the whole domain into a number of subdomains and in each subdomainexpand the function by a separate Chebyshev series. In this way, we can obtain exponentialconvergence.In chapter 3 we solve the Poisson's equation V 2u = f by the single-domain Cheby-shev spectral methods. In the r method, we build the discrete equation system in thetransformed space. In the collocation method, we build the discrete equation system inthe physical space directly. In both methods we solve the discrete equation system by thematrix-diagonalization method. There are some symmetry properties which can be ex-ploited to reduce the computational work. Numerical results indicate that when the sourcefunction is infinitely differentiable in the whole domain, the single-domain Chebyshev spec-tral method can achieve exponential convergence.In chapter 4 we solve the Poisson's equation by the multi-domain Chebyshev spectralmethod. We first outline the influence matrix method. Then we build the discrete equationsystem and solve it by the influence matrix method. The influence matrix is obtained bysolving Laplace's equation with special boundary conditions. For the Poisson's equationCHAPTER 1: INTRODUCTION 6the interface conditions require both the potential and its first-order normal derivative becontinuous. Numerical results indicate that if the whole domain can be divided into a finitenumber of subdomains so that in each subdomain the source function is infinitely differen-tiable, the multi-domain Chebyshev spectral method can achieve exponential convergence.In chapter 5 we present a new iterative Chebyshev spectral method for solving thegeneralized Poisson's equation V • (uVu) = f . If u is a smooth function, we can rewritethis equation in the form of a Poisson's equation V 2u = (f — Va. • Vu)/u and solve ititeratively. We give the first guess value arbitrarily. In each iteration, we compute the righthand side terms first from the current value of u. Then we solve the resultant Poisson'sequation by the direct collocation method discussed in chapter 3 to obtain the updatedvalue. We repeat the process until the norms of the differences between two successiveiterations satisfy an appropriate convergence criterion. For the 1-D problem, we give aformula to estimate the convergence rate of the new iterative method. The spectral radiusand the actual convergence curves of the numerical examples suggest that the convergencerate of the new method becomes better when the number of grid points increases. Thesame problems are also solved by the Chebyshev spectral multigrid method presented byZang et al. (1982). For the examples computed, the convergence rate of the new methodis faster than that of the Chebyshev spectral multigrid method.In chapter 6 we solve the generalized Poisson's equation by the multi-domain Cheby-shev spectral method. The procedure is similar to that for Poisson's equation describedin chapter 4. There are two differences. Firstly, the interface conditions for the general-ized Poisson's equation require both the potential and the flux, creu I On, be continuous.Secondly, we have to solve a generalized Poisson's equation instead of a simpler Poisson'sequation in each subdomain.In chapter 7 we apply the Chebyshev spectral methods to 3-D gravitational modeling.We extend Okabe's (1979a) analytical formula for the gravitational field of a polyhedralbody to its potential. We present a method to realize the point-injection and cell-averageCHAPTER 1: INTRODUCTION 7discretization for 2-D and 3-D problems. The multipole expansion technique is used tocompute the approximate boundary conditions. We present three examples. With ac-curate boundary conditions, the error norms of the single-domain solutions are roughlyproportional to (1/N 2 ). The multi-domain Chebyshev spectral solutions have exponentialconvergence. The maximum relative error caused by the approximate boundary conditionsis about 1%.In chapter 8 we apply the Chebyshev spectral methods to 3-D DC resistivity modeling.We discuss how to remove the singularity caused by the source current, give the formulationof the multi-domain Chebyshev spectral method, and present a modified finite-differencemethod for comparison. Numerical examples are given for two-layered, piecewise-constant,and piecewise-smooth models. All these examples are solved by both the multi-domainChebyshev spectral method and the finite-difference method. For the two-layered modelanalytic solutions exist. With the same number of grid points, the multi-domain Chebyshevspectral solutions are more accurate than the finite-difference solutions. For the realisticmodel which does not have an analytic solution, the results obtained by both the multi-domain Chebyshev spectral method and the finite-difference method agree with each otherquite well. However, the CPU time of the multi-domain Chebyshev spectral method isabout only one half of the finite-difference method.Chapter 9 is a summary of the results obtained in this thesis.CHAPTER 2SOME FUNDAMENTALS OF CHEBYSHEV EXPANSIONSThe Chebyshev polynomials were named after the Russian mathematician P. L. Cheby-shev (1821-1894) who first studied them. Lanczos (1938) proposed the r methods toapproximate a function by a finite Chebyshev series. Clenshaw (1957) proposed a recur-rence method to determine the Chebyshev coefficients of some known functions and to solveordinary linear differential equations. There have been many books and papers discussingthe Chebyshev polynomials and their applications. Snyder (1966), Fox and Parker (1968)and Rivlin (1974) are but a few of them.In this chapter, we outline some fundamentals of Chebyshev series as a basis for fur-ther discussion. In section 2.1, we discuss the analytic form of the Chebyshev expansionapproximation and the convergence properties. In section 2.2, we give the definition forthe 1-D, 2-D, and 3-D discrete finite Chebyshev expansions used in this thesis. In section2.3, we discuss the errors of the discrete Chebyshev expansions. In section 2.4, we discussthe Gibbs phenomenon and how to reduce the errors.2.1 Approximation of a Function by Chebyshev PolynomialsThe Chebyshev polynomial of the first kind is defined as (Lanczos 1938)Ti(x) =^cos-1 x)^— 1 < x < 1.^(2.1)From (2.1), we obtainTi(1) = 1,^ (2.2)8CHAPTER 2: FUNDAMENTALS OF CHEBYSHEV EXPANSIONS^9^Ti(-1) =^ (2.3)^T o (x) = 1, (2.4)andTi(x) =^ (2.5)Let 9 = cos -1 x. Using the trigonometric identity2 cos i0 cos 0 = cos(i 1)9 cos(i — 1)9,^(2.6)we obtain the recursion relationTi+i(x) = 2xTi(x) — Ti_1(x)^i = 1, 2, • • • .^(2.7)From (2.4), (2.5) and (2.7) we can compute Ti(x) for any i. We notice that Ti(x) is apolynomial of x with the highest degree i.Chebyshev polynomials are orthogonal with weight function 1/07- x 2 , i.e.it i=j=0,^ri Ti(z)Ti(z) dx =^i = j 0,^ (2.8)1—x2^20^i j.A function u(x) defined in —1 < x < 1 can be expanded in a Chebyshev series00u(z) =i=0From (2.8) and (2.9) the analytic Chebyshev coefficient ui is2 / 1 u(x)Ti(x) dxCor j_ i 07.7gwith=(2.9)(2.10){2 for i = 0,Ci =1 for i > 1.(2.11)CHAPTER 2: FUNDAMENTALS OF CHEBYSHEV EXPANSIONS 10Chebyshev polynomials are complete. If u(x) has continuous first-order derivativesand piecewise continuous second-order derivatives, then (2.9) is absolutely and uniformlyconvergent. Let the partial sum uN(x) be given byuN(x) = LuiTi (x ) .i.0For any given small positive number e, we can find an integer Ar e so that(2.12)iu(x) — uN(x)I e for all N > NE (-1 < x < 1). (2.13)The convergence of Chebyshev series is determined by the nature of the function to beapproximated. If u(x) is piecewise continuous and if u(x) is of bounded total variation for—1 < x < 1, thenlim uN(x) =N-►oou(1—)u(x+) u(x—)2u(-1+)for x = 1 ,for — 1 < x < 1 ,for x = —1 .(2.14)If u(P)(x) is continuous for all Ixt < 1 and for p = 0, 1,^, n —1, and if u(n)(x) is integrable,then ilk << l/kn, as k^oo. If u(x) is infinitely differentiable for Ix' < 1, the error in theChebyshev series goes to zero more rapidly than any finite power of (1/k) as k^oo (cf:Gottlieb and Orszag 1977, p28).Figure 2.1 shows the absolute values (in logarithmic scale) of Chebyshev coefficientsdefined by (2.10) for cos(rx), ex and 1/(1.1 — x) on interval [-1, 1]. Because cos(7rx) is aneven function, all the coefficients of odd order Chebyshev polynomials are zeros, which arenot shown. Since all these functions are infinitely differentiable for lx I < 1, their Chebyshevcoefficients exhibit exponential convergence. However, the convergence rates are different.The function e' has the best convergence rate and the function 1/(1.1 — x) has the poorest.From a spectral point of view, 1/(1.1 — x) contains more short-wave-length componentsthan ez.10^15^20^25xCHAPTER 2: FUNDAMENTALS OF CHEBYSHEV EXPANSIONS^11xxxxxx+^cos(irx)x^ex1/(1.1-x)xFigure 2.1 The Chebyshev coefficients of cos(rx), ex and 1/(1.1 — x)The slow convergence rate of the Chebyshev coefficients of 1/(1.1 — x) is due to thesingularity at x = 1.1, even though the singular point is outside the interval [-1, 1]. In 3-DDC resistivity problem, the primary potential 4 is of the form1= [(x x8)2 (y y,)2 ^z8)91/2On the line=z = z.,x1 . 0CHAPTER 2: FUNDAMENTALS OF CHEBYSHEV EXPANSIONS^12108g(x)a = 1.1a = 1.5 6a = 2.042.........................................—1.0 —0.5 0 0 0.5Figure 2.2 The function g(x) = 1/(a — x) on the interval —1 < x < 1.the potential can be simplified to1Xs — XTo facilitate the discussions in latter chapters, now we consider function1g(x) = a — x (-1 < x < 1 < a).At x = a, g(x) is singular. Figure 2.2 shows the values of g(x) on the interval —1 < x < 1for different a's. Intuitively, the closer is a to 1, the more rapid is the variation of g(x),CHAPTER 2: FUNDAMENTALS OF CHEBYSHEV EXPANSIONS^13that is, the more short-wave-length components exist in g(x). From a mathematical pointof view, we can explain the problem as the follows. We rewrite g(x) as1 ^1 g(z) = a 1 — (x / a) •Because of —1 < x < 1 < a, we have lx/al < 1. Function g(x) can be expressed as a Taylorseriescog(xx^tzvia 2-4 a)n=0The maximum absolute value of x/a is 1/a, which determines the convergence rate ofthe Taylor expansion (2.15). The closer a is to one, the slower is the convergence rate.Considering that a Chebyshev expansion with the highest degree N is the Nth degreepolynomial in x, it is easy to understand why the singularity at point x = a affects theconvergence rate of the Chebyshev expansion for 1/(a — x).2.2 The Discrete Chebyshev ExpansionsGenerally speaking, the Chebyshev expansions (2.9) should contain infinitely manyterms. In practice, however, we can deal with only a finite number of terms. We alwaystake a sufficiently large N and assume ui = 0 for i > N.The most commonly used 1-D Chebyshev grid points are defined as= cos —,^I=0,•••,N.^ (2.16)Figure 2.3 shows a 1-D Chebyshev grid with N = 4. Assumingul = u(x/),^ (2.17)the discrete forms corresponding to (2.9) and (2.10) are(2.15)21ui^ --ui cos= CiN E Cvi^Ni=o (2.18)CHAPTER 2: FUNDAMENTALS OF CHEBYSHEV EXPANSIONS^14 I X4^X3)X2^ XX1^00 1Figure 2.3 A 1-D Chebyshev grid (N = 4).7.14 = E ui cos 7-,^ (2.19)i=0with{2 i = 0 or N,Oi =^ (2.20)1 i=1,•••,N — 1.For 2-D problems the Chebyshev grid points are defined as11^1= 07rxi = cos — ,—, NN (2.21)ym = cos mr m = 0, ... , N.NFigure 2.4 shows a 2-D Chebyshev grid (N = 4), in which the number lm near an intersec-tion point indicates the position of (xi, y,n ). Assuminguim = u(zi, Ym),^ (2.22)the 2-D discrete Chebyshev expansions areN N4^1iiii = - - E E Oldm^N^N—uhn cos ihr_ cos j mr-(2.23)Ci Ca N2 1=0 m=0andN N ilw jmN wuhn = E E^_N co.^ (2.24)i=0 3=0-1andCHAPTER 2: FUNDAMENTALS OF CHEBYSHEV EXPANSIONS^1540^30y20^10^0031 21111 0132 22 12 02-133 23013 0334 24 14 04-1Figure 2.4 A 2-D Chebyshev grid (N = 4). The number Im representspoint (xi, y,n ).For 3-D problems the Chebyshev grid points are defined as—1/rxi = cos 1= 0,•••,NNyn = cos Mir m = 0, • • • , N^ (2.25)nx.Zn = cos — n 0,•• • , N.Figure 2.5 shows a 3-D Chebyshev grid (N = 4), in which the number lmn near anintersection point indicates the position of (z1, y n , zn). For this case we let414243441 xUlmn = U(X1) Ym, Zn)• (2.26)41002 ^-^Y^102412^- - -312^- -- _ _r^212- - - - ' -012-122 - - - 342- -232, - - - - 132 - - - - 032-^ 242^ 142 - 402^302443444401^301 ^ 201 ^ 101^001- -^ 2f1 - - lif -^111f121^- - - 021-^- - -431^-^ 1^- -^ 231^- - -^ -^031--1314f -441409^303^ 2091 IS - 103 1113003_ - -123^- - 699 _ -- - - ^133^- ^0332411 ^ 142 -^-414 - "^-314 - - -3.9.'^: "^ 2f4 - - -^ 11- ^104 -014 - - 004204- - - - - -^ - - - - - -124^- - -^024- - - - - - - - - - - - - - - -^ - - - -- - - - 134^- - - - 1342C4 ^ ' 144" - - - -644 - --Figure 2.5 A 3-D Chebyshev grid (N = 4). The number lmn indicates the position of point (zt, yin , zn).N N N ilw jmw knrUlmn = ^ilijk cos — cos — cos N •i=0 j=0 k=0(2.28)CHAPTER 2: FUNDAMENTALS OF CHEBYSHEV EXPANSIONS^17The 3-D discrete Chebyshev expansions are8 ^N N Nn=1sjk O'^N3k E E E C, nzO a tihnn COS 27. COSs 1=0 m=0 n=0andjmir cosNknirN (2.27)We call (2.18), (2.23) and (2.27) forward Chebyshev transforms; (2.19), (2.24) and (2.28)inverse Chebyshev transforms. Equations (2.18), (2.19), (2.23), (2.24), (2.27), and (2.28)are simply Fourier cosine series, which can be computed by the fast Fourier transformtechnique (Cooley and Tukey 1965; Cooley et al. 1970; Haidvogel and Zang 1979).2.3 The Error of Chebyshev ExpansionsWhen we approximate a function u(x) by a finite discrete Chebyshev series, there aretwo kinds of errors: truncation and aliasing. The truncation error comes from the omissionof short-wave-length components. The aliasing error comes from the superposition of short-wave-length components on the long-wave-length components. The amplitude of the errorsis measured by their L2 and L00 norms. The L2 error norm is the global mean square valueof the error and the L00 error norm is the maximum error over the domain, which is theabsolute error at a particular point.The truncation error of the partial sum (2.12) is 00e(t)(x) = uN(x) — u(x) = — E tliTi(X).i=N+1(2.29)Let 11.112andand 11.11 00 represent the L2 and L00 norms. We haveJ-1 1-^ = 1Il e(t)(x) 11 2 = 1ir i ^dx2 . E=N+1 ui(2.30)Ile(t)(x)110. = max E fiiTi(x)coi=EN+1(2.31)CHAPTER 2: FUNDAMENTALS OF CHEBYSHEV EXPANSIONS^18In practice we usually compute the discrete Chebyshev coefficients ui by (2.17) and(2.18). If u(x) contains components with Chebyshev wave number higher than N, therewill be an abasing error. From (2.9) the discretized value of u(x) at xi is00^00ui = Efinz,(.0 = E fin cos nir (2.32)n=0^n=0We can express the discrete coefficients "'Pik in terms of the continuous Chebyshev coefficientsfin . Using (2.18) the discrete Chebyshev coefficients fik are given byBecause^2 N 1 oo^nix)^kiruk =Ck N 4.4 Zri coss=0^n=02 co ( N E n cos N1^nix^kircos —N— cos —N—CkN n44=0 Ldi=0 (2.33)NE C1 nix kir^N k-.---- COS — cos^=^2s=o^N^N^0we havewhen n = 2mN f k,^(2.34)otherwise,(2.35)00iik iik -F E (112rnN–k -F^(k == 0, 1, • • • ,N).m=1The second term on the right-hand side of (2.35) is the abasing error which represents theshort-wave-length components in (2.32) appearing as long-wave-length components. If weapproximate u(x) bythe error e(d)(x) is given byNUN(d) (X) = E fliTi(x),i=0(2.36)e(d) (X) = E aiTi(x)i=000— EfIkTk (X)k=0^= E^2mN E E (112mN-Fk fi2mN–k) I Tk(Z)^E iikTk(x) . (2.37)m=1 k=1 m=1^ J^k=N+1oo N [oo^ 00CHAPTER 2: FUNDAMENTALS OF CHEBYSHEV EXPANSIONS 19On the right-hand side of (2.37), the first two terms represent the aliasing error and thethird term is the truncation error, which is the same as (2.29). The L2 and L. norms ofe(d)(x) are given respectively byIle (d) (x)11 2 =m=1^k=1„ 21 fl le(d) (Z ) .1 dIT sk i 1,/17.70"=oo 2ii2mN)00^ 2012m N 1-k + t2mN —le) +[m=1 and(2.38)(ie(d) (x ) L max le(d) (x ) 100Eg2mNm=iN+ Ek=100E (112mN-Ek tl2mN—k)m=100+ E Iuil. (2.39)i=N+1 The error norms given in (2.30), (2.31), (2.38), and (2.39) are useful in theoretical ana-lysis. However, they are rarely used in numerical computation. In practice the definitionsare slightly different. Let the superscript "(a)" represent the accurate value and the super-script "(n)" represent the discrete approximations. Suppose there are (N. + 1), (Nv + 1),and (N. + 1) grid points along x-, y-, and z-direction respectively. For discrete Chebyshevexpansions we define the maximum absolute value um , the relative error e, the L.. and L2error norms 'lel'. and 11e11 2 as follows.For 1-D problemsu,n = max luCal (0 i N.),sei = —^— Ui(n)^(a)Umiiefico = niPx led (0 i 5_ N.),$(2.40)(2.41)(2.42)andNa11e112 = [ ^1 ^ei21/2 (2.43)Nz +1 i,r)CHAPTER 2: FUNDAMENTALS OF CHEBYSHEV EXPANSIONS^20 For 2-D problems= max lu(7) I (0^< i <^0 j Ny ) ,1 (n)^(aeij =^)—Um(2.44)(2.45)(2.46)(2.47)andMoo = nta:xleiil (0^Nz, 0^NO,is1/2N. N,1^211e112 = [ (Nz 1)(Ny + E e•is=0 .J=0For 3-D problemsum = mai( u (i)k I (0 < < Arz , 0 < j < Ny, 0 < k < N.)^(2.48)eijk =—^ (2.49)UraHello° = max leol (0 i < Nx , 0 < j < Ny , 0 < k < Nz ),^(2.50)i,j,kand]1121^N. N, 22= ( ./Vm 1)(Ny 1)(Nz + 1)^E= k = eijk (2.51)2.4 Gibbs PhenomenonWhen u(x) contains discontinuities, its Chebyshev expansion (2.9) does not convergeuniformly. This is the well-known Gibbs phenomenon. Let us construct a simple example.Suppose u(x) is the signum function:—1 — 1 < < 0 ,u(x) = sgn(x) =^ (2.52)1^0 < < 1 .From (2.10), the analytic Chebyshev coefficient of u(x) isui = i)k ^4 7r(2k + 1) for i = 2k -I- 1 (k = 0,1, • • •).^(2.53)0^for i = 2kFrom (2.53), the analytic Chebyshev coefficient ui is proportional to 1/i. The convergencerate is quite slow. Figure 2.6 shows the partial sum (2.12) of sgn(x) for different N's. WhenN increases, the maximum errors near x = 0 remain the same, about 8.9% of the jump ofCHAPTER 2: FUNDAMENTALS OF CHEBYSHEV EXPANSIONS^211.0 \^10.5-- N = 65- - - - N = 33^ N = 17x-1.0^-0.5 0.5^1.0-0.5-1.0vFigure 2.8 The continuous Chebyshev expansion approximation ofsgn(x) obtained by the truncated series.u(x) at x = 0, even though the width of the maximum error area become narrower. Ata fixed point x (except x = 0) the error goes to zero as O(1 /N). From (2.53), there areinfinitely many short-wave-length components. Therefore the aliasing error in the finiteChebyshev expansion can become serious. Figure 2.7 shows the finite Chebyshev seriesapproximation (2.36) to sgn(x) (2.52). The discrete Chebyshev coefficients are obtainedby (2.17) and (2.18). The curves in figure 2.7 are obtained by changing (2.19) toNu(x) = EfikT,(.).k=0CHAPTER 2: FUNDAMENTALS OF CHEBYSHEV EXPANSIONS^22I‘1.0 ^r II0.5 Ir^- - N = 65- - - - N = 33^ N = 17-1.0^-0.5 0.5^1.0-0.5-1.0Figure 2.7 The Chebyshev expansion approximation of sgn(z) obtainedby point-injection discretization.Comparing figures 2.6 and 2.7 we find that the error of the discrete Chebyshev series islarger than that of the truncated series. This is due to the aliasing error.2.4.1 Cell-Average DiscretizationOne way to reduce the aliasing error is to use the so-called cell-average discretizationscheme. Let1 i = 0,di = cos Pgr ( — 1 ^i= 1,•••,N,—1^i=N-1-1.(2.54)CHAPTER 2: FUNDAMENTALS OF CHEBYSHEV EXPANSIONS^23We then haveJ .^- 1 J^7Ccos- 1 di -1— cos^=NCi •The cell-average discretization is defined asu(c) = dx.d;+, v 1 — x2(2.55)(2.56)In contrast we call (2.17), (2.22), and (2.26) point-injection discretization. Substituting(2.9) into (2.56) and using (2.16), we obtainulc) = fLo + - sin nr— fin cos nir (2.57)nw 2N^Nn=iFrom (2.18) and (2.57), the Chebyshev coefficient obtained by cell-average discretization is...(c)^2^1 (c)^kiruk = chN NE 6,7ui cos Ni=o[2^N 1^°° 2N nr.n67. 'do + E —nir si — fin ni7r^kir= )^C kN 2N cos IV cosCOS N^s=o s^n=1[2^N 1 kir °° 2N . nr N 1n x^kir _n . (2.58)^ilo E -„,..,,i^n=1^s=0^Ncos — + E — sm _ E -, co. — cos^U= C kN i=0 ‘j N^nr 2N^Ci^N From (2.34) and (2.58),Uo = U0^ (2.59)and_Jo^Br-= 2N^kir ['ilkT 12°(—^x sin^E 1r, ( 2mN —k 112mN+k km=1^2mN - k 2mN k(k = 1, • • • ,N).^(2.60)If we approximate u(x) byNU(;) (z) =^fikc)Th(x),^ (2.61)k=0( 11 2mN —k^1712mN +k + 11e(t)I I co (2.64)2mN — k 2mN kCHAPTER 2: FUNDAMENTALS OF CHEBYSHEV EXPANSIONS^24the error, e(c)(x) ise(c)(x) =t4car) (X) — ti(X)=E{(kr sin 2N — 1 uh2N . kir^) _—2N kr °°— sin — E ,w 2N(ii2mN —k 112mN 4-k )1 Tk (x )2mN — k 2mN km=1E flkTk(Z)k=N +1The L2 and Loo norms of e(c)(x) are given respectively by^2^1^1 [e(c)(x)1 2^Ile(c)(x)112 ^dx1 ,c--NN r(2N sin kr _ ,) fik^=i^kr 2Nk=1(2.62)2 ^2sin —2N^kr °°E (_1) m_i (2-12mN—k^fi2mN —r^2N m=1^-1 — k^2mN k^Ile(i) 112 ,(2 '63)andIle (c) (x)11 00 =max le(c)(x)I2N kr—kr sin 2N — 1) iikk =12N kr °°+—r sin gr- E (-1)m-1m=1Lanczos (1957, p221-228) discussed how to reduce the Gibbs oscillations of Fourierseries by ak method. The cell-average discretization defined by (2.54)—(2.56) is an extensionof Lanczos' methods to Chebyshev series. From (2.57) we see that by the cell-averagediscretization we actually multiply the kth Chebyshev coefficient by a factorsin(kr/N) crk = kir/Nwhich is the same as the Lanczos' ak factor.(2.65)CHAPTER 2: FUNDAMENTALS OF CHEBYSHEV EXPANSIONS^25Figure 2.8 The Chebyshev expansion approximation of sgn(x) obtainedby cell-average discretization.Figure 2.8 shows the finite Chebyshev series approximation (2.61) to sgn(x) obtainedby the cell-average discretization. Comparing figures 2.6, 2.7, and 2.8 we find that theerrors of the finite Chebyshev expansion obtained by the cell-average discretization is smal-ler than that obtained by point-injection discretization, especially for points far from thediscontinuity.Comparing (2.30), (2.38) and (2.63), we see that for a given N, the truncated analyticChebyshev series give the best accuracy, i.e., give the smallest L2 error norm. If 92,i = 0 forn > N, there are no aliasing errors in the approximation. In this case the point-injectionCHAPTER 2: FUNDAMENTALS OF CHEBYSHEV EXPANSIONS^26discretization works better than the cell-average discretization. When u(x) contains discon-tinuities whether in itself or in its derivatives of any order, it will contain short-wave-lengthcomponents. Its Chebyshev coefficients will converge slowly. When we approximate u(x)by a finite Chebyshev series, there are aliasing errors. For the cell-average discretization,there is no aliasing error for uo . This makes it more accurate than the point injectiondiscretization, especially when we use spectral methods to solve differential equations.We will discuss how to realize the cell-average discretization for some 2-D and 3-Dproblems and give examples of their application in chapter 7.2.4.2 Approximation of a Function by Separate Chebyshev Series in DifferentSubdomainsWe can convert an arbitrary 1-D domain Si =^I t2 < x < i2} to the standardChebyshev domain SY = I - 1 < < by= 2x — (21 + 22) xl —(2.66)Therefore we can approximate a function defined in an arbitrary 1-D domain by Chebyshevexpansions. In the previous subsection we improved the accuracy of Chebyshev expansionsby cell-average discretization. However, the improvements are small. To achieve the desiredexponential convergence, we divide the whole domain into a number of subdomains sothat in each subdomain the function is infinitely differentiable. Then we approximate thefunction by a separate Chebyshev series in each subdomain.We consider the Chebyshev approximation of sgn(x) again. From (2.53) the analyticChebyshev coefficient fii of sgn(x) is of 0(1/4 To obtain an accuracy of 10 -6 we have touse 106 terms! However, if we divide the whole domain into two subdomains [-1, 0] and[0, 1], and approximate sgn(x) in each subdomain by a separate Chebyshev series, becausein each subdomain sgn(x) is constant, only one term is required to represent the functionwithout any error!CHAPTER 2: FUNDAMENTALS OF CHEBYSHEV EXPANSIONS^27Let us consider another example. Suppose thatZ 2 + X1 < X < 0 ,(2.67)2If we expand u(x) into a single Chebyshev series (2.9), from (2.10) the analytic Chebyshevcoefficient fh isfti =4^131-(-1)k+1 47-(2k — 1)(2k + 1)(2k + 3)0for i = 1for i= 2k+1 (k > 1)for i = 2k (k > 0).(2.68)Because u(x) and its first-order derivative are continuous in [-1, 1], the Chebyshev coeffi-cients is of 0(03 ). The convergence rate is better than that in example 2.1. However,it is still slow. In order to obtain an accuracy of 10 -6 , we need to use 100 terms. If wedivide the whole domain into two subdomains and approximate u(x) in each subdomainby a separate Chebyshev series, then because u(x) is the second degree polynomial of x,only three terms are required to represent u(x) in each subdomain.The idea to represent a function by separate Chebyshev series in different subdomainscan be extended to 2-D and 3-D functions.2.5 DiscussionThe Chebyshev polynomial Ti(x) is the ith degree polynomial of x. If u(x) is a polyno-mial with the highest degree equal to or less than N, the Chebyshev expansion with degreeN is accurate everywhere. If u(x) is infinitely differentiable, its Chebyshev expansion con-verges exponentially.When there are discontinuities in u(x) or its derivatives, its single domain Chebyshevexpansion approximation does not converge exponentially. By using cell-average discret-ization, we can reduce the error to some extent. However we cannot achieve exponentialconvergence. To obtain the desired exponential convergence, we should divide the wholeU(X) = x2 _2x0 < x < 1 .CHAPTER 2: FUNDAMENTALS OF CHEBYSHEV EXPANSIONS 28domain into a number of subdomains so that in each subdomain the function is infinitelydifferentiable. Then we approximate the function in each subdomain by a separate Cheby-shev series. In this way, we can use fewer terms to obtain higher accuracy. It is thisvery idea that leads to the multi-domain methods which will be discussed in the followingchapters.CHAPTER 3SOLVING POISSON'S EQUATION, PART I:SINGLE-DOMAIN METHODSIn this chapter we discuss how to solve the boundary value problemV2u = f^in SZ ,^ (3.1a)u = known on (3.1b)by Chebyshev spectral methods on a single domain. In section 3.1 the r method is discussed.In section 3.2 the collocation method is discussed. A comparison is given in section 3.3.In the 1-D case equation (3.1a) becomesd2udx2 = f(x) •In the 2-D case equation (3.1a) becomes02u a2uOz2 ey2 = f(x, Y)•In the 3-D case equation (3.1a) becomes02u a2u 82uOx2 ay2 az2 = f(x ,Y1z).^(3 .4)(3.2)(3.3)29CHAPTER 3: SOLVING POISSON'S EQUATION, PART I^ 303.1 The r MethodIn the r method, proposed by Lanczos (1938), a function is expanded into a finite sumof Chebyshev polynomials and correction terms. The coefficients of the Chebyshev poly-nomials and the correction terms are determined from the function to be expanded. Haid-vogel and Zang (1979) used the r method to solve the 2-D Poisson's equation on a squarewith homogeneous boundary conditions, and proposed both alternating-direction implicititeration and matrix-diagonalization method to solve the coefficient equation system. Hald-enwang et al. (1984) extended Haidvogel and Zang's 2-D method to the 3-D Helmholtzequation with general inhomogeneous boundary conditions. We will discuss the 1-D prob-lem in detail, to illustrate the basic ideas, and then generalize the method to 2-D and 3-Dproblems. In this section we assume that N is even and -1 < x, y, z < 1.3.1.1 The Relations Between the Chebyshev Coefficients of a Function and itsDerivativesWe expand the sth derivatives of u(x) into a Chebyshev seriesoodsu(x) =^ (3.5)i.0Fromwe have/Ti(x)2 I Ti(x) dx = 1 772(x)Ti+i(x) Ti-i(x)i + 1^i - 1for i = 0 ,for i = 1 ,for i > 1,(3.6)4+1)^.48) 48+1)= 2sui^ui+1 (3.7)where Cl is defined by (2.11). Using (3.7) repeatedly, we obtain the Chebyshev coefficientsof the first-order derivative of u(x) (cf. Orszag 1971),1^°3-zu• = — E 21iii^ (3.8)Ci i=i+i1-1-i oddand1+i odd1=i+1U* = 1^21a,Ci (3.10)N-, z = 1 E l(12 i2 )ut (3.11)CHAPTER 3: SOLVING POISSON'S EQUATION, PART I^ 31and the Chebyshev coefficients of the second-order derivative of u(z),u = — -x.^1v 1(12_ i2)111.^ (3.9)Ci 1=i+21+i evenThe finite discrete forms corresponding to (3.8) and (3.9) are1=i+21+i evenwhere Cl is defined by (2.20).3.1.2 The 1 -D ProblemFor the 1-D equation (3.2), the Chebyshev grid is given by (2.16). We expand f(x) anduzz(z) into the 1-D Chebyshev series (2.19). Equating the coefficients of the like Chebyshevpolynomials on both sides of (3.2), we obtain equationsfir = A^= 0,1, ... N — 2) .^(3.12)From (3.11) and (3.12) we obtain1 N i2 A = A = 0,1, ... N — 2).^(3.13)1=i+21+i evenLetu(1) = b1^ (3.14)andu(-1) = b2 .^ (3.15)Combining (2.2), (2.3), (2.19), (3.14) and (3.15), we obtainao-Fai+ii2 1-•-•+t-IN=bi^ (3.16)and1 .4 — il(bi + b2)71i8i = .-f: — 2 (bi — boeii even,i odd.CHAPTER 3: SOLVING POISSON'S EQUATION, PART I^ 32and—^+ '112 — • • +^= ba •^ (3.17)ThereforeUN-1 = 1(bi — 62) — (ui + u3 + • • + eN-3)andux = 2 (bi^) — (110 + 712 + • • • + eN_2).For notational convenience we define71i = N(N2 — i2 )(3.18)(3.19)(3.20)and6 = (N — 1) [(N — 1)2 — i2} .^ (3.21)Eliminating^and eN in (3.13) with (3.18) and (3.19), we find thatN-2E (Lae: = si^(i = 0,1, .. N — 2) ,^(3.22)b=0whereail =1--2 77i i even,1 = 0,i 1 even, 1 < 1 <- +1(12 — i2 ) i,1 even, l>—6^i, 1 odd, 1 <—6 + 1(12 — i2 ) i, 1 odd, l>0^otherwise,(3.23)From (3.23), all the entries of the matrix [aii] with i even, j odd or i odd, j even are zero.System (3.22) can be solved by many different methods such as Gauss eliminationmethod or Gauss-Seidel iterative method. For a better understanding of the 2-D and 3-DCHAPTER 3: SOLVING POISSON'S EQUATION, PART I^ 33problems to follow, we solve (3.22) by the matrix-diagonalization method. We factorize thematrix [ail] asN-2ail = E a-A .a311si (3.24)j=owhere Aj is the jth eigenvalue of matrix [ail], [aii] is the corresponding eigenvector matrix,and [aji] -1 is the inverse of [a4]. Therefore^N-2N-2^{1 = 1,E = (3.25)j=0^J=0 0 i 1.All the eigenvalues are real, non-negative, and distinct (Canuto et al. 1988, p.409). Thecondition number of the matrix [aii] computed by numerical experiments up to N = 64 is0(N4 ). Therefore the matrix-diagonalization method is not suitable for large N.Substituting (3.24) into (3.22) we obtainN-2 N-2E E1=0 j=0= 8i^(i = 0, 1, . . . N - 2) .^(3.26)We then multiply the ith equation of (3.26) by aml and sum over i to obtainN-2^N-2 N-2^N-2E a,n1 E E aijAjajllul = Ei=0^1=0 5=0^i=0(3.27)LetN-2wj = E ail ui1=0(3.28)andN-2Tm = Ei=0Because of (3.25), (3.28) and (3.29), we rewrite (3.27) asN-2E amAjw; = rm1=0Therefore(3.29)(3.30)r„,^(3.31)^CHAPTER 3: SOLVING POISSON'S EQUATION, PART I^ 34From (3.28),N-2uj = E „win .^ (3.32)m=0In summary, the system (3.22) can be solved using the following two steps:N-2Wi = —Ai^(i = 0, • • • , N — 2),^(3.33)1=0andN-2ul = Ealiwi^(/ = 0, • • • , N — 2).^(3.34)i.0With ao, ui, . .. , fi-N_2 determined by (3.33) and (3.34), we can calculate CiN_1 and 2N by(3.18) and (3.19). Then ui can be calculated by (2.19).From (3.11) and (3.23) the coefficients of even and odd order Chebyshev terms are notcoupled. As a consequence of this, half of the entries of both [aii] and [aki] — i are zeros.We can take advantage of this property to reduce the computational work.3.1.3 The 2-D ProblemFor the 2-D equation (3.3) the Chebyshev grid is given by (2.21). We expand uzz, um/and f into 2-D Chebyshev polynomials (2.24). Equating the coefficients of the same orderChebyshev polynomials on the two sides of (3.3) we obtain=^(0 < i , j < N — 2).^(3.35)From (3.11) and (3.35) we obtain equationsN^ N1 V 1(12 ;2)15^E m(m2 — j2 ) .aini =Cs^ m1=42 =j+21+1 even^m+j even(0 < i,j < N — 2).^(3.36)Letu(1, y) = h(y) ,u(-1, y) = 1(y) ,u(x,1) = b(x),andCHAPTER 3: SOLVING POISSON'S EQUATION, PART I^ 35u(x, -1) = k(x).From (3.1b) the functions h(y), 1(y), b(z) and k(x) are known. We compute the Chebyshevcoefficients hi, li , 3i and hi by (2.18). By a derivation similar to that for (3.18) and (3.19),we have1=^- ki) -^+ ui3 + • • • +^(0^N - 2 ) ,uiN = 1 + ki) - (aio + ui2 + • • + ei(N_2))^(0^N - 2) ,17(N-1); =^-^- kik; +113; + • • • + u(N-3)j)(0 < j < N),(3.37)(3.38)(3.39)anda^1 uNj =^+ hi) — (uoi + uzi + • • • + u(N-2)j)^(0 j < N).^(3.40)Substituting (3.37)-(3.40) into (3.36) to eliminate^eiN, 17(N-1)j, and aN; andmoving the known terms to the right-hand side we obtain the final systemN-2^N-2ailulj +^= 8ij^(i,i =0,1,— N - 2),^(3.41)i=0^i=owhere the matrix [aii] is given by (3.23) andiii - 1 [0; +TIAN+ (Ei + TodIii - i [0; + Tovi + 6i - reofd4 = ii.; - 1 [(I; - li.ofi + 0-, + iondIv - 2[(1.i - rti)fi + (bi - Icoedi even, j even,i even, j odd,i odd, j even,i odd, j odd,with vi and fi defined by (3.20) and (3.21). We solve the system (3.41) by the matrix-diagonalization method using the following two steps:1 N-2 N-2 —1 —1 nzwij^ail aimsi= m2d=oandN-2 N-2elm = LEi=0 .J=0(1,m = 0, • • • N - 2).^(3.43)(i, j = 0,^N - 2),^(3.42)CHAPTER 3: SOLVING POISSON'S EQUATION, PART I^ 36We can find iii(N-1), uiN, 1-1(N-1)i and 11.:Ni by (3.37)—(3.40). Then /hi can be computedby (2.24).3.1.4 The 3-D ProblemFor the 3-D equation (3.4), the Chebyshev grid is given by (2.25). We expand u",u" and f into 3-D Chebyshev polynomials (2.28). Equating the coefficients of the sameorder Chebyshev polynomials on the two sides of (3.4) we obtain equationsIL'gC +^+^•jk^(0 < j, k < N — 2).— ijk^ijk^i jk^sFrom (3.11) and (3.44) we obtain equations(3.44)E 1(12 .:2)frCif 1=421-I-i evenm(m2 - J ) 0irnkm=j-1-2m+j even1 N+- E n(n2 — k2 )1Iijn = fijk^(0 < i, j < N — 2) .^(3.45)Ck n=k+2evenLetu(1,y , z) = h(y, z),u(-1, y, z) = 1(y, z),u(x,1, z) = b(z , z),u(x, —1, z) = k(x, z),u(x,y, 1) = t(x,y),u(x, y, —1) =^, y).From (3.1b) the functions h(y, z), 1(y, z), b(x, z), k(z, z), t(x,y), and d(z, y) are known.We compute the 2-D Chebyshev coefficients he,, /kb bij , kij , tij and dij by (2.23). By a^CHAPTER 3: SOLVING POISSON'S EQUATION, PART I^ 37derivation similar to that for (3.18) and (3.19), we have1^dij)^ILO + • • • + ilij(N-3))^(0 < j < N - 2) , (3.46)/sof =^3) — (usio + uij2 + • • • + aij(N-2))^(0 i,^N - 2) ,^(3.47)ai(N-1)k = 12 (bik - kik) - (aiik fii3k + • • ' tii(N-3)k)^(0 5_ i , k < N),^(3.48)eiNk =^(17i0k fii2k • • • + ei(N-2)k)^(0 5 k N),^(3.49)1 ft.fl(N-1)jk =^— ilk) — (ui.jk U3jk + • ' • + ri(N-3)jk)(0 < i < N - 2, 0 < j < N),^(3.50)andfiN';k = 2 31 (it •h^3k) (110jk^2jk ' • • + 11(N-2)jk)(0 < < N - 2, 0 < j < N).^(3.51)Substituting (3.46)-(3.51) into (3.45) to eliminate eij(N-1), aqN,^aiNk,and eNjk, and moving the known terms to the right-hand side we obtain the final equationsystemN-2^N-2^N-2E E E akiew = 8ijk^(i,j,k = 0,1 ... N - 2),^(3.52)1=0^1=0^1=0where the matrix [aij] is given by (3.23), andfijk^Pijk +Tik) 7/i (bik +iik) 771^clij)7/ki^Khik +Ijk)ni (gik + 17ik)ni^- it/Mk]fijk - -21 [(iiik +rik)77/ (bik -^(4)7/k1fijk — 5 [viik + &Joni + toile -,Tfijk — 2 [lajk^/CO%fijk[Oijk ' 1/06 +^+ kik )771 +^- Zigkifijk iVijk^ + (Fs.; + J0;ik3ijkj, k even,j even, k odd,i, k even, j odd,i even, j, k odd,i odd, j,k even,i, k odd, j even,i odd, j,k even,k odd, j even,and^E E E aimakn.'imn_,wick Ai 4. Ai -1-^1=0 m=0 n=0N —2 N —2 N —2= • • • ,N — 2)^(3.53)CHAPTER 3: SOLVING POISSON'S EQUATION, PART I^ 38with ni and 6 defined by (3.20) and (3.21). We solve system (3.52) by the matrix-diagonalization method using the following two steps:N —2 N —2 N —2ulmn = E E E aiiamjank Wilk^(1, m, n = 0, • , N — 2).^(3.54)i=0 j=0 k=0We can find 47--(N-1)jk; UNjk, fii(N-1)k,^and uijN by (3.46)—(3.51). Then uqkcan be computed by (2.28).3.2 The Collocation MethodIn the r method, the equation system is derived for the Chebyshev coefficients. Inorder to build the equation system we have to apply forward Chebyshev transforms tof to convert the quantity from physical space to transformed space. After obtaining theChebyshev coefficients of u, we have to apply inverse Chebyshev transform to convert thequantity from transformed space to physical space. In the collocation method we requireequation (3.1a) be satisfied exactly at the interior points and we directly solve the equationsystem in physical space. No Chebyshev transforms are required.3.2.1 The First and Second Order Derivative MatricesIn the application of Chebyshev expansions, we frequently need to compute the deriv-atives of a function u(x). When the values are given on the Chebyshev grid points, thefirst-order derivatives can be computed as follows. We first compute the Chebyshev coeffi-cients iii by a forward Chebyshev transform (2.18), then compute u by (3.10). Finally wecompute the first-order derivatives by an inverse Chebyshev transform (2.19). The aboveprocess can be expressed in the form of a matrix multiplication. We denote the forwarddi) -i.; -2(1 —2N2 + 162N2 -I- 160 < i,j < N, j1 <i=j<N-1,i = j = 0,i=j=N.(3.58)CHAPTER 3: SOLVING POISSON'S EQUATION, PART I^ 39Chebyshev transform operator by matrix [cii] and the inverse transform operator by matrix[co. From (2.18) and (2.19), we obtain2cif =^ cosCiCfN N (3.55)and—1^ijlrcif•• = cos^. (3.56)Let{2jv ) = Zi; j > i and i + j odd,0^otherwise.From (3.55), (3.56), and (3.57) the first-order derivative matrix [e] is given byN Nc4; ) = E E C 1 1)1, 1,?cmi1=0 m=0N^N m mi4^1^ihr^n-NO.^,,,_E co. COS.7\ -r - E an N •=2 1=0 ' 1^m=i+il+m oddGottlieb et al. (1984) gave the following explicit expression for PO as(3.57)Let 2 _ i2V(?) =%I^/ ici a )) 3. > i.and I+ j even,0 otherwise.From (3.55), (3.56), and (3.59) the second-order derivative matrix [d (j)] is given byN N2) _ E E —1 (2)Egli —^Cii Vim Cm)1=0 m=0N^2 N 1^its-^m(m2 _12)^mill.= — E — cos — E ^ co.No. 0 N Om^N3 i=o 1^m=1+21-1-m even(3.59)4(2)^j(2)"(N -i)(N - j) - "ij •We can take advantage of (3.63) and (3.64) to reduce the computational work.3.2.2 The 1-D Problem(3.61)(3.62)(3.63)(3.64)1ui = E ) uiiti=0andIt is easy to verify thatand2(1)"(N -i)(N - j) = 611Pxx^1(2)u = E a„ u i .i=0CHAPTER 3: SOLVING POISSON'S EQUATION, PART I^ 40Peyret (1986) gave the following explicit expression for [d (in as(-1)i+i (x? zixi — 2) 1 < i < N —1, 0 < j < N, j i,Ci( 1 x?)(zi xi) 2(N2 — 1)(1 — x?) -I- 3 1 <i=j<N-1,3(1— 4)2do) = (-1)2 2 [(2N2 + 1)(1 — xi) — 6]ij^ i = 0, 1 < j < N,^(3.60)301(1 — xj) 2(-1)N+1 2 [(2N2 -I- 1)(1 + xi) — 6] i = N, 0< j < N —1,30i( 1 + zi)2N4 — 1 k = j =0,N.15With (3.58) and (3.60) the first and second-order derivatives of tt at grid points can befound byIn section 3.1, we assumed —1 < z ,y, z < 1. In practice, we frequently face arbitraryintervals. In these cases, in order to solve the problem by Chebyshev spectral methods wehave to convert the arbitrary domains into the standard Chebyshev domains by coordinatemapping. For the convenience of further application, from now on we will consider problemsin more general domains.CHAPTER 3: SOLVING POISSON'S EQUATION, PARTIn the 1-D equation (3.2), we now suppose 0 = {xi— (21 + 22)Ix2 < X < ill. We can use41(3.65)(3.66)= ^ -_xi — x2to map 0 to 0' = le I — 1 < e < 1}. Assuming that2r =_^_x1 — x2we havedu^du.dx —andd2u^2 d2u=r(3.67)(3.68)dx2^42 •Therefore equation (3.2) becomes2 d2 u (3.69)=r eNow the Chebyshev grid is21 + 22^1 — . (3.70)xi = 2^+^cosThe discrete form of (3.69) is• • , N — 1), (3.71)r2^= fi^(i = 1, •i=owhere 4; ) is defined by (3.60). From (3.1b) the values of uo and uN are known. Movingthe known terms to the right-hand side of (3.71), we obtain the system of equationsN-17.2 E gut^(i = 1,•••,N— 1),^(3.72)i=iwherece) = 4?) (1 < i , 1 < N — 1),^ (3.73)andai = A— r2 (402)U0 42pitiN) •^ (3.74)CHAPTER 3: SOLVING POISSON'S EQUATION, PART I^ 42As is the case of r method, we factorize the matrix [42)] asN-142)1 = E (3.75)i=1where vi is the jth eigenvalue of matrix [c4 12)], [j9,] is the corresponding eigenvector matrix,and [Pi] -1 is the inverse of [AA ThereforeN-1^N-1^{1 i = 1,E^= E pijlQjl = Sil =^ (3.76)5=1 1=1 0 i I.All the eigenvalues are real, negative, and distinct (Gottlieb and Lustman 1983). Thecondition number of matrix [42) ] is about 0.019N4 . Because of (3.64) we havepi(v- j) =^ (3.77)andflCALoi — (_i)N-1-i+lki 1 .^ (3.78)We can use (3.77) and (3.78) to reduce computational work. We solve the system (3.72)by the following two steps:1 N-1wi = r2vi E Pi7 131^(i = 1 , • • • ,N — 1),1=1andN-1U/ = 2 riwi^(i = 1, • • • ,N — 1).i=i3.2.3 The 2-D Problem(3.79)(3.80)In the 2-D equation (3.3), we now suppose Cl = {(z, y)I 22 < X < 21, P2 < y <Let 2r =_ _^ ,z1 — Z228 = .Y2(3.81)1r =^8 = 21 2— 2 '22yl —2t ^(3.87)CHAPTER 3: SOLVING POISSON'S EQUATION, PART I^ 43Then the 2-D Chebyshev grid is1 ii + 22 + r X1 irXi = 2 cos — ,-I — cos — .yi = Pi +92 , 1 yr2 a NAs is the discussion in subsection 3.2.2, we obtain an equation system(3.82)N-1^N-1^r2 E c4,2)„, + 82 E^= sq^(i, j = 1,. , N —1),^(3.83)1=1^1=1where_7(2^/2)^, _4^ 2)8ij = fi; - r2 kaki) UOI aiNuNI) — 2 kujo2) ^+U80 itiNtigN) • (3.84)We solve the system (3.83) by the matrix-diagonalization method using the following twosteps:andWii =1^N-1 N-1r2vi 821/4 EEi3i71/3118inzi=1 m=1 j = 1,• • • ,N —1),^(3.85)N-1 N-1Ultn = E^puomiwii^(1,m = 1, • • • , N — 1 ).^(3.86)i=1 j=13.2.4 The 3-D ProblemIn the 3-D equation (3.4) we suppose SI = {(x, y, z)I xa < x < 21, h <Y < gi, xa <z < il}. LetThe 3-D Chebyshev grid is^+ 22 1^ir1 = 2 + cos —r X '_^1 cosya — 2^— ,8 N^21 + 22 1^kr= 2 -I- t cos N— .(3.88)CHAPTER 3: SOLVING POISSON'S EQUATION, PART IAs is the discussion in subsection 3.2.2, we obtain an equation system^N-1^N-1^N-1r2 E + 82 E^+ t2 E 42,)„„, = 8ijk^1=1^1=1 1=144(i, j, k = 1,whereaiik = fijk — r2(402)uojk^duNik) — 82(C420)UiOk^d(2)u.—t2(d120)uijo^di2),TuiiN)•• • , N — 1), (3.89)(3.90)System (3.89) can be solved using the following two steps:N-1 N-1 N-11, • • • , N — 1) , (3.91)Wijk = 82y j^t2vk E E E^(i, ;7, k =1=1 m=1 n=1andN-1 N-1 N- 1Ulmn =^PliSmAnkWijk^(1, m,n = 1, • • • ,N — 1). (3.92)i=i j=1 k=13.3 Discussion3.3.1 AccuracyThe condition number is about 0.123N 4 for the matrix [aq] of the r method (section3.1.2) and 0.019N4 for the matrix [ii of the collocation method (section 3.2.2). WhensjN increases, the condition number increases rapidly. For example, when N = 64, thecondition number is about 2 x 10 6 for the r method and is about 3 x 10 5 for the collocationmethod. To ensure that the solutions have enough accuracy, we should use double precisionin the computation.3.3.2 Comparison of the r and Collocation MethodsIn the r method, the matrix-diagonalization is performed in the transformed space.Because the even and odd order Chebyshev coefficients are uncoupled, half the entries areCHAPTER 3: SOLVING POISSON'S EQUATION, PART I^ 45zero in matrices [aii] and [a0. We can use these properties to reduce the computationalwork. However, both forward and inverse Chebyshev transforms are required.In the collocation method, the discrete equation system is built in the physical spacedirectly. The matrices Ai] and [137 1] are full, but also have some symmetry properties (seesj(3.77) and (3.78)). We can use these properties to reduce the computational work. Sinceno Chebyshev transforms are required, the collocation method is faster and simpler thanthe r method.In the r method we need to know the values of f at the boundaries (which are notrequired in the collocation method) in order to obtain the Chebyshev coefficients fi in(3.12), A i in (3.35), and liik in (3.44). When the boundary conditions and the right-hand side terms are self-consistent, the r method works well. When we solve the 2-D and3-D Poisson's problems with the multi-domain Chebyshev spectral method in chapter 4,initially the interface conditions are given arbitrarily. Because of a lack of consistency, ther method does not work but the collocation method does. The collocation method is moregeneral than the r method and can be used to solve more general equations.CHAPTER 4SOLVING POISSON'S EQUATION, PART II:MULTI-DOMAIN METHODIn chapter 3 we solved the Poisson's equation V 2u = f by the single-domain Chebyshevspectral (SDCS) method. When the source function f is infinitely differentiable, the con-vergence is exponential. However, when there are discontinuities in the source function, theaccuracy deteriorates seriously. This will be evident from the example 4.2 to follow. Thereason is that when the source function f is discontinuous, the second-order derivative ofthe potential function u is discontinuous too. Near the end of chapter 2 we discussed thatwhen there are discontinuities in the function or its derivatives, the error in its Chebyshevexpansion approximation does not go to zero exponentially.In this chapter we solve Poisson's equation by the multi-domain Chebyshev spectral(MDCS) method. In this method the whole domain is divided into a number of subdomainsso that in each subdomain the source function is infinitely differentiable. Then the potentialfunction is approximated by a separate Chebyshev series in each subdomain. The relationsof the Chebyshev series in adjacent subdomains are determined by the matching conditionsat the interface points. In this way we can obtain higher accuracy and use fewer gridpoints than the single-domain Chebyshev spectral method. In section 4.1 we outline theinfluence matrix method for solving a linear equation system, which is an essential part ofthe multi-domain Chebyshev spectral method to follow. In sections 4.2-4.4 we present themulti-domain Chebyshev spectral method for 1-D, 2-D, and 3-D Poisson's equations. Wediscuss the 1-D problem in detail to illustrate the main ideas. Then we generalize these to46CHAPTER 4: SOLVING POISSON'S EQUATION, PART II^ 472-D and 3-D problems. We present numerical examples for each case. In section 4.5 wediscuss the results and present the conclusions.4.1 The Influence Matrix MethodIn this section, we outline the basic idea of the influence matrix method, which wassuggested by Buzbee et al. (1971) to solve a linear equation system arising from thenumerical solution of partial differential equations. Following Buzbee et al. we use acapital letter to represent a matrix and a bold lower case letter to represent a vector in thissection. Suppose we are given a linear equation systemAx = y,^ (4.1)where A is an n x n matrix and y is a vector with n elements. The matrix A is partitionedasAiA=A2where Ai is an ni x n matrix, A2 an n2 x n matrix, and ni + n2 = n. The vector y ispartitioned accordingly asYlY=Y2where yi and y2 are vectors with ni and n2 elements respectively. Then we replace Al byB1 to obtain another n x n matrixB = ( B1 )A2and replace yi by y1 to obtain another vector 7 with n elementsCHAPTER 4: SOLVING POISSON'S EQUATION, PART II^ 48Given an arbitrary nonsingular ni x ni matrix W, we create an n x ni matrix W by addingan matrix of n2 x ni zeros. ThusW =The influence (capacitance) matrix C is defined as the ni x ni matrixC =provided det B 0. Assume that the vector z with ni elements is the solution of equationCz = Yi —^ (4.2)Then the solution of system (4.1) is given byx = B-1(y+ Wz)^ (4.3)See Buzbee et al. (1971) for a rigorous proof. Now we present a simple example todemonstrate the procedures discussed above.Example 4.1 For a given system(1 1 1) (xi^61 —1 1^x2 = 2 ,1 1 2^x3^9(4.4)instead of solving (4.4) directly, we replace the first two equations in (4.4) byxi = 2andx2=11to obtain a simplified system(1 0 0 ) (xi^( 20 1 0^x2 = 1^.1 1 2^.z3^9CHAPTER 4: SOLVING POISSON'S EQUATION, PART II^ 49Using the notations in the above discussion, here n = 3 , n1 = 2, n2 = 1,1 0 0B= (0 1 0)1 1 2A1 = (1 11 1 )'2y= (1),9andYi = (6)2,It is easy to find that1^0 0B-1 =^0^1 0-0.5^-0.5 0.5Let1^0W^(0 1) .0^0ThenC = A 1s-1W = ( 0 : 5Yi - 11 1 13-1if = °2 )^ (4.9)Substituting (4.8) and (4.9) into (4.2) and solving the resultant equation system we obtainz =^ 11) (4.10)Substituting (4.5)-(4.7) and (4.10) into (4.3) yields the correct solution of system (4.4)(1x = B ---1 (y+Wz) = 2) .3The procedure for this example seems tedious. But when applying the influence mat-rix method to solve the equation system arising from multi-domain Chebyshev spectralmethod, we do not need to compute B -1 explicitly. We will discuss this problem later.,(4.5)(4.6)(4.7)(4.8)CHAPTER 4: SOLVING POISSON'S EQUATION, PART II^ 504.2 The 1-D ProblemThe basic ideas of the multi-domain Chebyshev spectral method can be illustrated bythe following 1-D model problem:d2udx2 = .f(z)^n,u = b(x)^on r ,where f(x) and b(x) are known functions, 1 = {x^< x <^r is the boundary,B. is an integer, and^< it (1 = 1, • • • , B.). We assume that^f(z) = fl(x) for 41_1 < z < xj (1= 1, • • • ,^(4.12)where fl(x) (1 = 1, • • • , B,,) is a smooth function. There are bounded discontinuities inf(x) at x = xi (1 = 2, • • • , B.). Motivated by (4.12) we divide the whole domain CI intoB. subdomains fli = I 2,44 < x < (1 = 1, • • • , B.) (see figure 4.1, where Bo = 3). Insubdomain (2 1 , u = u1 and f =^The boundary value problem (4.11) is equivalent tod2u1= f^in 11' (1 = 1, • • • , B.),0 3 2^1I^I374^ 373^)72dx2u = u^at x=dul-1= duidxdx at x = xi (I = 2, • — ,^,u = b(x)^on r.(4.13a)(4.136)(4.13c)(4.13d)Figure 4.1 1-D subdomains.CHAPTER 4: SOLVING POISSON'S EQUATION, PART II 51In the multi-domain Chebyshev spectral method, we approximate the potential functionu by a separate Chebyshev polynomial in each subdomain. For simplicity we use the samedegree N in all subdomains. Let221 — 211-1 .The 1-D Chebyshev grid is/^2/ -I- 2/+1^1^ix2^—= cos N— (1 = 1, • • • , B.; = 0, • • • , N).7./ = (4.14)(4.15)For the 1-D problem, we call the points that belong to r external boundary points, thosebelonging to only one subdomain interior points, and those belonging to two subdomainsinterface points. Let /4 = u(xi) and^= Az!). Similar to the discussion in subsection3.2.2, from (4.13) we haver? E (42Y = fi^(i = 1, • • • , Br=07,1_1 E^_ rt E euri = 0r=0 r=04,1-1^4 ,1^n ty-N^-0 =^= • • • /^/where 41) and 42) are defined by (3.58) and (3.60). From (4.13d) ul and uk are known.Boundary conditions (4.16c) can be satisfied by taking the same value for u li 1 and u& (1 =1, • • • , From (3.63) and (4.16), we obtain the system of equations:r?^d1r2)tir' = ft^(1 = 1, • • • , B.; i = 1, • • • ,N —1),^(4.17a)r=0E + E iori)uri = 0^(1 = 2, • • • , B.)•^(4.176)r=0^r=0Equations (4.17a) are for the interior points; (4.176) for the interface points. The unknownsin the adjacent subdomains 0 1-1 and 11 1 are coupled by (4.176).In principle the system (4.17) can be solved by the Gauss elimination method. Butfor illustrating the idea of 2-D and 3-D problems, we solve (4.17) by the influence matrixmethod outlined in section 4.1 using the following steps:= 1,• • •,N —1), (4.16a)(1= 2, • • • , (4.166)(4.16c)CHAPTER 4: SOLVING POISSON'S EQUATION, PART II^ 521. We change (4.17b) to1- 1^1= uo = 0^= 2 , • • • , B.). (4.18)In the new system formed by (4.17a) and (4.18) the unknowns in different subdomainsare no longer coupled. The new system can be separated into B. smaller systems, eachof which can be solved independently by the single-domain Chebyshev spectral methoddiscussed in subsection 3.2.2. The solution obtained in this way will not satisfy (4.17b),since the values at the interface points are not given correctly. The left-hand side of(4.170 gives the differences of the first-order derivatives at the interface points. Wedenote the difference by di, that is,N^ N_1( ^1 J^1+1oil = T 1 E ..or) tiN_r + 11+1 E uO).or ur^(1 = 1, • • • , B. — 1).^(4.19)r=0 r=02. We have to find the correct values vg for the interface points from d1. That is, we haveto build the influence matrix [cpq] and solve the resultant systemBr _1E cpqvq = 4^= 1,• • • ,B. — 1).^(4.20)q=1If the whole coefficient matrix and its inverse had to be computed to obtain the influencematrix as implied by (4.3), it would be too expensive computationally. Fortunately,the influence matrix can be obtained in a much simpler way. Both the first and thesecond-order derivative operators are linear with respect to the interface values. Therelation between the interface values and the jumps of the first-order derivatives intwo contiguous subdomains is linear too. In the influence matrix system (4.20), vgrepresents the difference between the assigned value and the correct value at the qthinterface point, while dp represents the difference of the first-order derivatives on thetwo sides of the pth interface point. The coefficient cpq represents the influence ofa unit change in vg to the right-hand side term 4. Therefore the influence matrixcan be found by solving a Laplace problem with special boundary conditions. Weassign 1 to one interface point and 0 to all others, solve the Laplace problem on theCHAPTER 4: SOLVING POISSON'S EQUATION, PART II 53two subdomains connected to the non-zero interface point, compute the jumps of thefirst-order derivatives at that point and the two neighbour interface points to obtainthe entries of the influence matrix. In principle we can repeat the above procedurefor all the interface points in turn to obtain the whole matrix. For the 1-D problemunder consideration there is yet an even easier way to build the influence matrix. In asubdomain the solution of a 1-D Laplace equation is a straight line. When the boundaryconditions are set to 1 and 0 at the two ends, the first-order derivative of the solutionat any point is the slope k = ±1/(ii —3. For the 1-D problem because a change of an interface value affects only the first-orderderivatives at that point and the two adjacent points, the matrix [cpq] is tridiagonal.We solve system (4.20) to obtain vq . Then we correct the interface values with vq .4. We solve system (4.17a) in each subdomain again to obtain the correct solution.If a function consists of piecewise-smooth (infinitely differentiable) functions, its multi-domain Chebyshev expansion approximation should display exponential convergence. Theexamples given in the rest of this chapter are designed to show these properties. Now weconsider an 1-D example, in which the analytic solution is piecewise-smooth function.Example 4.2 In boundary value problem (4.11) we assume— cos x^for — 3 < x < —1,f (z) = h (z) =u(3) = 1,—(sin x + cos x)— sin xfor — 1 < x < 2,for 2 < x < 3,(4.21)(4.22)andu(-3) = —1. (4.23)CHAPTER 4: SOLVING POISSON'S EQUATION, PART II^ 54We can find the analytic solution as follows. We integrate (4.21) twice to obtain thegeneral solutioncos Z + cix -I- c2^for — 3 < x < —1,u(x) = ui(x) = cos x + sin x + csx + c4 for — 1 < x < 2,^(4.24)sin x + c5x + c6^for 2 < < 3.Then we determine the constants in (4.24) by the interface conditions (4.13b), (4.13c) andthe external boundary conditions (4.22) and (4.23). Because the analytic expressions forthese constants are quite complicated, we present the numerical values here only. Theconstants are:• = 0.58567833289620,c2 = 1.7470274952890,• = 0.045376027028056,• = 2.0481961742288,cs = —0.86392139979763,andce = 3.4506441913330.We denote the source function f(x) in (4.21) as f1(x) and the potential function u(x)in (4.24) as ui(x). They satisfyd2ui(x) dx2 — fi(x) (4.25)The functions ui(x) and fi(x) consist of piecewise-smooth functions sin x, cos x, and linearfunctions ax + b. These functions are infinitely differentiable in each subdomain. At theinterface points x = 2 and x = —1, ui(x) and ui(x) are continuous, whereas fi(x) isdiscontinuous. We plot functions fi(x) and ui(x) in figure 4.2. The functions ul(x) andf1(x) will be used to construct 2-D and 3-D piecewise-smooth functions later.CHAPTER 4: SOLVING POISSON'S EQUATION, PART II^ 55Figure 4.2 The functions fi(x) and ul(x).CHAPTER 4: SOLVING POISSON'S EQUATION, PART II^ 56Table 4.1The error norms of example 4.2 (MDCS method)N B. Pt L2 Loo4 3 13 7.018 x 10—°3 6.690 x 10—°35 3 16 2.673 x 10 —°3 2.059 x 10 -036 3 19 1.118 x 10—°4 1.137 x 10—°47 3 22 3.697 x 10 -09 2.964 x 10—°58 3 25 9.861 x 10 -07 1.048 x 10 -069 3 28 2.941 x 10—°7 2.408 x 10—°210 3 31 5.566 x 10—°9 6.074 x 10 -09Table 4.2The error norms of example 4.2 (SDCS method)N Pt L2 Lc°8 9 8.462 x 10 -02 1.310 x 10 -0116 17 2.546 x 10 -02 4.364 x 10 -0232 33 2.111 x 10—°3 3.777 x 10—°3Tables 4.1 and 4.2 show the relative error norms of the solutions obtained by the 1-Dmulti-domain Chebyshev spectral method and the single-domain Chebyshev spectralmethod discussed in subsection 3.2.2 respectively. In tables 4.1 and 4.2, N is the degree ofChebyshev polynomial used in each subdomain, B. the number of subdomains along thex-direction, Pi the number of total grid points, and L2 and Lao are the corresponding errornorms. From table 4.1, we see that the error norms of the multi-domain Chebyshev spectralsolutions go to zero exponentially. On a grid with 31 points, the relative error norms are10-9 . From table 4.2, the error norms of the single-domain Chebyshev spectral solutionsare roughly proportional to (1/N 2 ). On a grid with 33 points, the relative error normsCHAPTER 4: SOLVING POISSON'S EQUATION, PART II^ 57of the single-domain Chebyshev spectral solutions are f•-• 10 -3 . The multi-domain Cheby-shev spectral solutions are far more accurate than the single-domain Chebyshev spectralsolutions.4.3 The 2-D ProblemNow we consider the 2-D boundary value problem82u 82u= f (2%1Y)0 z2 ay2U = b(z, y)^on r ,in fl, (4.26a)(4.26b)where Il = {(x, y) I Be+1 < Z < xl, YB,+1 Y <^r is the boundary, B. and By areintegers, functions f(x,y) and b(x, y) are known. We divide the whole domain Il into B.Bysubdomains (see figure 4.3, where B. = By = 2) with dividing lines x = xi (1 = 2, • • • , B.)and y = ym (m = 2, • • • , By ). We denote the subdomain located at the lth column and themth row by Olm. In Stlm u = ulm and f = flm. The boundary value problem (4.26) isequivalent to82 ulm 82 ulm8x2^0y2 = f^in O lm (1 = 1, • • • , B.; m = 1, • • • , By) ,^(4.27a)U(1-1)114 = ulm at x = xl (1 = 2, • • • , B.; m = 1, • • • , By ) ,^(4.27b)ul(m-1) = ulm at^=ym (1 1, • • • , B.; m = 2, • • • , Be,),^(4.27c)Ou(1-1)m = Oulmex^at x = xl (1 = 2, • • • , B.; m = 1, • • • , By ) ,^(4.27d)eul(m-1 )^Oul'n^at y = ym (1 = 1, • • • , B.; m = 2, • • • , By ),^(4.27e)ey^eyu = b(x, y)^on r. (4.271)Assuming that2Iri - _ _xi -28m, = _ _^Ynt Ym+1(4.28)CHAPTER 4: SOLVING POISSON'S EQUATION, PART II^ 580 21 0 11022012x3^X2^ X,Figure 4.3 2-D subdomains.the 2-D Chebyshev grid is{ x1 = 21 -I-:1+1 + 1 cos i.ffri^N'Yi --=^m §m + Vm+1 + 1^i .r2^-I-^cos —^s,„^N (4.29)For the 2-D problem, we call the points that belong to r external boundary points, thosebelonging to only one subdomain interior points, those belonging to two subdomains inter-face points, and those belonging to four subdomains corner points. There are B.137,(N —1) 2interior points, (B.By — B. — By )(N — 1) interface points, and (B. — 1)(Bv — 1) cornerpoints. At the corner points, if we require the first-order derivatives be continuous alongboth x- and y-directions, the system will be over determined. Therefore we require the371Y2ST3CHAPTER 4: SOLVING POISSON'S EQUATION, PART II^ 59first-order derivatives be continuous along only one direction. From (3.61), (3.62), (3.63),(3.68), and (4.27) we haveNE iir)utn = firr=0E +r=0(1 = 1, • • • , B.;Nrl—i ' 1)^rnu(1-1)Or (N-r)1r=0(1 = 2, • • • , B.;Nam-1 E dvitil(;:r1r=0(1= 1,•••, B.;Nri—i E d( u01) —1)mOr (N—r)0r=0^r=0(1 = 2, • • • , Bz ; m = 2, • • • , By ),i= 1,•••,N— 1),1), (4.30a)(4.306)(4.30c)(4.30d)m= 1, • • • , By ; i , j = 1, • • • ,N —NE cioir)u;.7 = 0r=07n, = 1, • • • , By ; j = 1,• • • ,N — 1),8771^do;.)utn = 0r=0m= 2,•••,By ;\—•lm = 0+ *1^Or Ur0where d19 ) and^are defined by (3.58) and (3.60). Equations (4.30a) are for the interiorpoints, (4.306) and (4.30c) for the interface points, and (4.30d) for the corner points. Theunknowns in different subdomains are coupled by (4.30b) and (4.30c). The corner pointsdo not appear in (4.30a)—(4.30c).Similar to the 1-D problem we solve the subsystem (4.30a)—(4.30c) by the influencematrix method using the following steps:1. We change (4.306) and (4.30c) toutril =0^(1 = 2, • • • ,Ba ; 711. = 1, • • • , By ; j = 1, • • • , N — 1)^(4.31a)andu!'on = 0^(1 = 1, • • • ,Bz ; rn. = 2, • • • , By ; = 1, • • • ,N — 1),^(4.31b)solve the new system formed by (4.30a) and (4.31) in each subdomain independentlyby the single-domain Chebyshev spectral method discussed in subsection 3.2.3, andcompute the jump of the first-order normal derivatives at all the interface points.CHAPTER 4: SOLVING POISSON'S EQUATION, PART II^ 602. We build the influence matrix by solving Laplace's equations. There are some symme-tries which can be exploited to reduce the computational work. For example, in figure2.4 (here we view it as a subdomain), the influence of point 01 on point 20 is equal tothat of point 03 on point 24. The influence of point 02 on point 10 is equal to that ofpoint 42 on point 34. For general rectangular subdomains, we need to solve 2[(N+1)/2]([x] is the integer part of x) Laplace equations in each subdomain to build the wholeinfluence matrix. If all subdomains are squares with the same size and are divided intothe same N2 cells, we need to solve the Laplace problem only [N/2] times. All theother entries of the influence matrix can be obtained by making use of the symmetryproperty.3. The influence matrix is sparse and diagonally dominant. We solve it by Gauss-Seideliterative method to obtain vk. This method works. But by no mean is it the mostefficient method to solve the problem. Then we correct the interface values in (4.31)with vk.4. We solve the system formed by (4.30a) and (4.31) again to obtain the correct solutionfor all the interior and interface points.We find the values of u at the corner points by solving much smaller systems. Substi-tuting the interface values into (4.30d) and moving the known terms to the right-hand side,we obtain (By — 1) smaller tridiagonal systems, each of which contains (B. — 1) equationsfor (B. — 1) unknowns and is easy to solve.Now we consider a 2-D example, in which both the source function f and the potentialfunction u are the products of piecewise-smooth functions.Example 4.3 In the boundary value problem (4.26) we assume the domain is= {(x, y) I - 3 <^< 3}and the analytic solution isu(x,y) = ui(z)ui(Y),^ (4.32)CHAPTER 4: SOLVING POISSON'S EQUATION, PART II^ 61where ui is given by (4.24). From (4.25), (4.32) and (4.26a) we havef (x y) = (z )ai (Y) + (z )fi (Y ).where h is defined by (4.21). The boundary conditions are given accordingly. In the multi-domain Chebyshev spectral method, we divide the whole domain into 3 x 3 subdomainswith the dividing lines i2 = 2, 23 = -1, g2 = 2, and g3 = -1. Table 4.3 shows theerror norms of the numerical solutions obtained by the multi-domain Chebyshev spectralmethod. In table 4.3, N is the order of Chebyshev polynomials used in each subdomainalong one direction, B. and By the number of subdomains along the z- and y-directionsrespectively, Bt the total number of subdomains, Pt the total number of grid points, K thenumber of iterations for solving the influence matrix system, ce the relaxation parameter,L2 and Lac, are the corresponding error norms, and the CPU time is obtained on a SPARCstation 2 computer. From table 4.3, the error norms of multi-domain Chebyshev spectralsolutions go to zero exponentially. On a grid with 961 points, the error norms are ti 10-9 .Table 4.3The error norms of example 4.3 (MDCS method)N Bx By Bt Pt K w L2 CPU(sec)4 3 3 9 169 11 1.394 1.760 x 10-03 4.879 x 10-03 0.055 3 3 9 256 17 1.458 7.176 x 10-04 1.571 x 10-°3 0.086 3 3 9 361 24 1.505 2.972 x 10 -135 8.847 x 10-°5 0.137 3 3 9 484 27 1.540 1.047 x 10 -05 2.226 x 10-°5 0.198 3 3 9 625 39 1.570 2.718 x 10-°7 8.108 x 10-07 0.299 3 3 9 784 43 1.591 8.581 x 10-08 1.825 x 10-°7 0.3410 3 3 9 961 58 1.610 1.563 x 10 -139 4.629 x 10 -49 0.51CHAPTER 4: SOLVING POISSON'S EQUATION, PART II^ 624.4 The 3-D ProblemNow we consider the 3-D boundary value problem02u 02u 82u^ f8x2 oy2 az2 — f kz, Y , z)^in 1 ,^(4.33a)u = b(x,y, x)^on t ,^ (4.33b)= {(x,Y,z)1 2B.-1-1 < x < x1, VB,+1 < Y < gl, 213,,+1 < z < 21}, F is theboundary, B., By and B. are integers, functions f(x, y, z) and b(x, y, z) are known. Wedivide the whole domain f/ into Be .13yB, subdomains (see figure 4.4, in which B. = By =B, = 2) with dividing planes x = xl (1 = 2, • • • ,B.), y = (in = 2, • • • , Be,), andz = in (n = 2, • • • ,134. We denote the subdomain located at the lth column, mth row, andnth layer by In f lmn , u = ulmn and f = fimn. The boundary value problem (4.33)is equivalent to the following problem:where 11z iz 2z•Sr4X3. • •^ca.XI3Figure 4.4 3-D subdomains.CHAPTER 4: SOLVING POISSON'S EQUATION, PART II^ 63492ulmn a2ulmn 492ulmnOZ2^0y2^az2^ = flmnin Simn (1= 1,• • • , B.; m = 1, • • • ,By ; n = 1,• • • ,B,),^(4.34a)U^= ulmnlmul(m--1)n = ulmnlm(n-1) = u lmnat x = xt (I = 2, • • • , B.; m = 1,• • • , By ; n = 1, • • • ,^(4.34b)at y = gm (I = 1,• • • , B.; m = 2, • • • , By ; n = 1, • • • ,B,),^(4.34c)at z = in (1= 1, • • • , B.; m =1,• • • By ; n = 2, • • • ,E4),^(4.34d)Ou(1-1)mn &LimnOx = Oxat x =^(1 = 2, • • • , Bz ; m = 1,• , By ; n = 1,- • • ,B,),^(4.34e)Oul("1-1)n &Limnay^ayat y =^(1 = 1,• • • , B.; m = 2, • • • By ; n = 1, • • • ,Bz ),^(4.34f)auin(n-1) &limnay^ayat z = in (1= 1, • • , B.; m = 1, • • • By ; n = 2, • • • ,^(4.34g)u = b(x, y, z)^on r.^ (twoLetThe 3-D Chebyshev grid is2ri = _ _^^7x1 — zi+i2sm - _2 _^,Yrn —* Ym+1to - ^.Zn — Zn+1it + it+i , 1^iirzi =^-i- cos —2^ri^N 7m gm+ Pm-I-1 , 1^jriY =^-t-^cos 2^sm Nn in + 4+1 , i^kirzk =^-I- to cos —2^N .(4.35)(4.36)^CHAPTER 4: SOLVING POISSON'S EQUATION, PART II^ 64For the 3-D problem we call the points that belong to 1" external boundary points, that be-longing to only one subdomain interior points, that belonging to only two subdomains in-terface points, that belonging to four subdomains edge points, and that belong to eight sub-domains corner points. There are a total of B.ByBz (N - 1)$ interior points, (3B.ByB. -B.By - BxBz - ByB.)(N -1)2 interface points, (3B.ByB, - 2B. By 2B.B.-2ByBzB. + By + B.)(N - 1) edge points, and (B. - 1)(By - 1)(B. - 1) corner points. At theedge and corner points, for the same reason as that for the 2-D problem in section 4.3, werequire the first-order derivatives be continuous along only one direction.From (4.34h) the values of u at the external points are known. The interface conditions(4.34b)-(4.34d) can be satisfied by taking the same interface value for the two contiguoussubdomains. From (3.61), (3.62), (3.63), and (4.34) we haveN^N^NE 4r2)ur,74, sm2 E 42,)utne^E ik2r) urir = fixr=0^r=0^r=0(1 = 1, • • • , B.; m = 1,• • • , By ; n = 1,• • • , Bz;^=^• , N - 1),^(4.37a)7.1_1 E d(ori,u(wi_lznk^E 41)47: = 0r=0^ r=0^, B.; m =1,• • • , By ; n = 1,• • • ,B.; j,k = 1,• • • ,N -1),^(4.37b)+ E 4„.)utr = 0r=0, B.; m = 2, • • • , By ; n = 1, • • • ,^i,k = 1,• • • ,N -1),^(4.37c)t^d(1)ul.m(n-1) + tn v■ d(1)orinn- i^Or sj(N-r)^Or Ur = or=0 r=0^(1 = 1, • • • , B.; m = 1,• • , By ; n = 2, ••• Bz; j = 1, • • • , N - 1),^(4.37d)d(1)u(1-1)mn j_ r . E iori)urv: = 0ri-i Li Or (N-r)Okr=0^ r=0(1 = 2, • • • , B.; m = 2, • • • , By ; n = 1, • • • ,^k = 1,• • • ,N -1),^(4.37e)7.1_1 E d(ori,u(wi_!):,);^E iori)ur,Ton = 0r=0^ r=0(1 = 2, • • • , B.; m = 1,• • • , By ; n = 2, • • , Bz; j = 1, • • , N - 1),^(4.37f)N 11T8m-1 E^+ 3,n E ioir)utnon = 0r=0 r=0(1 = 2, • • •8m-1 Er=0(1 = 1, • ••CHAPTER 4: SOLVING POISSON'S EQUATION, PART II^ 65^(1= 1, • • • ,B.; m = 2, • • • , By ; n = 2, • • • , Bz ; i = 1, • • • ,N — 1),^(4.37g)8m-1 E ciorl'uokinN— lr;n0 + E 41)uurnji = 0f=0^,•=0^(1 = 2, • • • , Bz ; m = • • • , Bei ; n = 2, ••• , Bz ),^(4.37h)where dW and 4) are defined by (3.58) and (3.60). Equations (4.37a) are for the interiorpoints; (4.376)—(4.37d) for the interface points; (4.37e)—(4.37g) for the edge points; and(4.37h) for the corner points. The edge and corner points do not appear in equations(4.37a)—(4.37d).Similar to the 1-D and 2-D problems, we solve the subsystem (4.37a)—(4.37d) by theinfluence matrix method using the following steps:1. We change equations (4.376)—(4.37d) tou107: = 0(1= 2, • • • , Bz ; m = 1,— • , By ; n = 1, • • , Bz;Imnuiok = 0( 1 = , • • • , B.; m = 2, • • • , By ; n = 1, • • • ,Bz ; i,k = 1, • • • ,N 1),^(4.38b)u10 =0^(1 = 1, • • • , Hz ; m = 1, • • • , By ; n = 2, • • • , Bz ; i,j =1,• • • ,N —1).^(4.38c)In the new system formed by (4.37a) and (4.38), the unknowns in different subdomainsare no longer coupled. We solve them in each subdomain independently by the 3-D single-domain Chebyshev spectral method discussed in subsection 3.2.4. Then wecompute the left-hand side terms in (4.376)—(4.37d) to obtain the differences betweenthe first-order normal derivatives in two contiguous subdomains related to a particularinterface point.2. We build the influence matrix by solving Laplace's equations with special boundaryconditions. Similar to 2-D problems, there are some symmetries on the influence of onepoint on another to be exploited to reduce the computational work. In figure 2.4 (herej, k = 1, • • • , N — 1),^(4.38a)CHAPTER 4: SOLVING POISSON'S EQUATION, PART II 66we view it as a subdomain), for example, the influence of point 110 on point 011 isequal to that of point 130 on point 031. The influence of point 210 on point 101 is equalto that of point 210 on point 301. For general rectangular parallelepiped subdomainswe need to solve the Laplace problem 3[(N + 1)/2] 2 times in each subdomain to obtainthe entries of the influence matrix. If all the subdomains are cubes of the same size,and are divided into the same N3 cells as shown in figure 2.5, we need to solve theLaplace problem only [N/2]([N/2] + 1)/2 times. All the other entries of the influencematrix can be obtained by making use of the symmetry property.3. We solve the influence equation system by Gauss—Seidel iteration to obtain the correc-tion vector vq . Then we correct the interface values in (4.38) with vq .4. We solve the system formed by (4.37a) and (4.38) again to obtain the correct solutionfor all interior and interface points.We compute the values at the edge points as follows. We substitute the known valuesinto (4.37e)—(4.37g) and move the known terms to the right-hand side. Equation (4.37e) canbe separated into (By - 1)Bz(N - 1) independent smaller systems. Equation (4.37f) can beseparated into By (B, — 1)(N — 1) independent smaller systems. Each of the small systemscontains (B.-1) equations for (B.-1) unknowns. Equation (4.37g) can be separated into(B. —1)Bz (N —1) smaller systems, each of which contains (B y — 1) equations for (By — 1)unknowns. All these smaller systems are tridiagonal and are easy to solve.We compute the values at the corner points in a similar way. Substituting the knownvalues into (4.37h) and moving the known terms to the right-hand side, we obtain (B. —1)(B, — 1) independent tridiagonal smaller systems, each of which contains (B y — 1) equa-tions for (By - 1) unknowns. Solving them, we obtain all the unknowns.Now we consider a 3-D example, in which the analytic solution is the product ofpiecewise-smooth functions.CHAPTER 4: SOLVING POISSON'S EQUATION, PART II^ 67Example 4.5 In the boundary value problem (4.33), we assume the analytic solution isu(x,y , = tii(x) ui(y)ui(x),^ (4.39)where ui is defined by (4.24). From (4.25), (4.33a) and (4.39), we havef (x, y, z) = h(z) ui (y) ui(z) + ui(x) fl (y) ^+ ui(x) ui(y) h(z).The boundary conditions are given accordingly. The domain is divided into 27 subdomainswith dividing planes 22 = 2, 23 = -1, g2 = 2, gs = —1, 22 2 , and 4 = — 1. Therelative error norms are shown in table 4.4. In table 4.4, N is the degree of the Chebyshevpolynomial used in each subdomain along each coordinate direction, B., By , and B. thenumber of subdomains along the x-, y-, and z-directions respectively, Bt the total number ofthe subdomains, Pt the total number of grid points, K the number of iterations for solvingthe influence matrix system, w the relaxation parameter, L2 and L. the correspondingerror norms, and the CPU time is obtained on a SPARC station 2 computer. From table4.4, the relative error norms of the multi-domain Chebyshev spectral solution go to zeroexponentially. On a grid of 15,625 points, the relative error norms are 10 -7 .Table 4.4The error norms of example 4.4 (MDCS method)N B. By By Bt Pt K w L2 L. CPU(sec)4 3 3 3 27 2, 197 19 1.400 1.251 x 10 -03 3.828 x 10—°3 2.375 3 3 3 27 4,096 25 1.470 4.935 x 10—" 1.295 x 10—°3 8.046 3 3 3 27 6,859 34 1.520 1.856 x 10 -05 7.415 x 10—" 22.827 3 3 3 27 10, 648 37 1.560 6.729 x 10 -06 1.801 x 10—°5 50.398 3 3 3 27 15, 625 52 1.590 1.607 x 10 -07 6.401 x 10 -07 114.10CHAPTER 4: SOLVING POISSON'S EQUATION, PART II^ 684.5 DiscussionIn this chapter we solve the 1-D, 2-D, and 3-D Poisson's equation using the multi-domain Chebyshev spectral method. We divide the whole domain into a number of sub-domains and approximate the potential function by a separate Chebyshev series in eachsubdomain. The relations of the Chebyshev series in adjacent subdomains are determinedby the matching conditions at the interface points. We solve the discrete system by theinfluence matrix method. We build the influence matrix by solving Laplace's equation withspecial boundary conditions.For 1-D problem, the influence matrix is tridiagonal and is easy to solve. For 2-D and3-D problems, we solve the influence matrix system by the Gauss-Seidel iterative method.The number of iterations depends on the particular problem and the accuracy requirements.In examples 4.3 and 4.4, the analytic solutions are piecewise-smooth functions. For a fixeddegree of Chebyshev expansion, the accuracy is limited. A convergence criterion that istoo stringent will only waste time.In the 1-D example 4.2, the source function is discontinuous. Therefore the second-order derivative of the potential is also discontinuous. The relative error norms of thesingle-domain solutions are roughly proportional to (1/N 2 ) (see table 4.2). On a grid with33 points, the relative error norms of the single-domain Chebyshev spectral solutions areti 10-3 . This is consistent with the conclusion in chapter 2. For the same problems, themulti-domain Chebyshev spectral method displays exponential convergence. On a gridwith 31 points, the relative error norms of multi-domain solution obtained are 10 -9(see table 4.1), which are about 6 orders better than those obtained by the single-domainChebyshev spectral method on a grid with 33 points. In examples 4.3, and 4.4, the analyticsolutions are piecewise-smooth functions. The multi-domain Chebyshev spectral solutionsalso display exponential convergence.CHAPTER 4: SOLVING POISSON'S EQUATION, PART II 69From the above discussion we conclude that for Poisson's equations, when the sourcefunction is discontinuous, the multi-domain Chebyshev spectral method can achieve ex-ponential convergence, while the single-domain Chebyshev spectral method cannot. Inother words the multi-domain Chebyshev spectral method can achieve higher accuracywith fewer grid points than the single-domain Chebyshev spectral method. Although theformulation of the multi-domain Chebyshev spectral method is more complicated than thesingle-domain method, the multi-domain Chebyshev spectral method is more efficient thanthe single-domain method.CHAPTER 5SOLVING THE GENERALIZED POISSON'S EQUATION, PART I:SINGLE-DOMAIN METHODThere have been many papers describing the solution of the generalized Poisson's equa-tions by Chebyshev spectral method. Zang et al. (1982, 1984) solved the 2-D generalizedPoisson's problem with Chebyshev spectral multigrid method. They used incomplete LUdecomposition method to solve the finite difference preconditioning equations and gaveseveral examples. Brandt et al. (1985) proposed several modifications to the method pre-sented by Zang et al. for periodic problems. Deville and Mund (1985, 1990) discussed theChebyshev pseudospectral algorithm for the second-order generalized Poisson's equationsusing finite element preconditioning. They computed the spectral radius of the generalizedPoisson's solvers with finite difference and finite element preconditioning and with differ-ent interpolation elements. Canuto et al. (1988, chapter 5) and Zang et al. (1989) gavea comprehensive discussion and review on the Chebyshev spectral method in solving thegeneralized Poisson's equations.In this chapter we present a new iterative Chebyshev spectral method for solving thegeneralized Poisson's equation in a single domain. In section 5.1 we outline the relationbetween the convergence rate and the spectral radius of the iteration matrix. In section 5.2we present the formulation of the new method for 1-D, 2-D, and 3-D generalized Poisson'sequations. Numerical examples are presented for each situation. A discussion is given insection 5.3.70CHAPTER 5: SOLVING GENERALIZED POISSON'S EQUATION, PART I^715.1 Spectral Radius and Convergence RateLet [gii] be an N x N matrix with eigenvalues (1 < i < N). ThenP([NA) = max NI1<i<N (5.1)is the spectral radius of the matrix [gii] (Varga 1962). Suppose we solve equationNui = fi +^girur^= 1,•••,N)^(5.2)r=1by the iteration relationNui(k) =^E giruVz-1) .r=1Subtracting (5.2) from (5.3) yieldsNu(k) — ui =^gir (uf.k-1) — ur) .r=1Equation (5.4) indicates that the error el k) = ulk) — ui satisfiesNei(k) = E gir 4") .r=1SupposeN= E cr (er )i,r=1where (g,.)i is ith element of the rth eigenvector of matrix [gii] and cr 's are constants. Afterk iterations the error is(k) (k—i)ei =^gipep^= • • • = E cr Airep=1 r=1(5.3)(5.4)(5.5)(5.6)(5.7)When k becomes large enough, elk) will be dominated by the largest eigenvalue(s). Thatis, the convergence rate is determined by the spectral radius of the iteration matrix.CHAPTER 5: SOLVING GENERALIZED POISSON'S EQUATION, PART I^725.2 A New Iterative MethodWe consider the generalized Poisson's problemV . (oVu) = f in 1, (5.8a)u = b on r, (5.8b)where a and b are known functions, a > e > 0, SZ is the domain, and F is the boundary. Inthe new iterative method, equation (5.8a) is rewritten asv2 f — Vu u (5.9)aprovided a is a smooth function and Vala is not too large. If the right-hand side of (5.9)were known, it could be solved by the direct method discussed in chapter 3. Based on thisidea we solve boundary value problem (5.8) by the following procedure:1. For a given guess value of u, compute the right-hand side of (5.9).2. Solve the resultant Poisson's equation by the collocation method discussed inchapter 3.3. Take the solution obtained as new guess and repeat the above process until the normsof the differences between two successive iterations satisfy an appropriate convergencecriterion.Since a is given, Vu can be computed analytically or numerically. However Vu has to becomputed numerically. In the Chebyshev spectral method the first-order derivatives alongany coordinate direction are quite easily obtained by formula (3.61) with high accuracy.Now we apply the above method to 1-D, 2-D, and 3-D problems. Numerical examplesare given for each situation. In each example we compute the solution by both the presentmethod and the Chebyshev spectral multigrid method presented by Zang et al. (1982). Inthe Chebyshev spectral multigrid method, a simple V cycle with one relaxation on each gridCHAPTER 5: SOLVING GENERALIZED POISSON'S EQUATION, PART I^73level is used. The correction vector is obtained by solving a finite difference preconditioningequation system.5.2.1 The 1 -D ProblemFor the 1-D problem, equation (5.8a) becomesddu\dx kffax) =' in II,^ (5.10)where Si = {x Ix2 < x < 21}. Equation (5.10) can be rewritten asd2u f — Q(x)u(x)(5.11) dx2^a.^•Let2r =_^_ .xi — X2 (5.12 )The 1-D Chebyshev grid isxi =21 + 22^1^1:71-+ —2^rcos .N (5.13)From (3.61), (3.62) and (5.11) we obtain equationsN^ (2) NE12) 11,3 _ ^Ai)— — — r ^.^= 1, • • • , N —1),^(5.14)o7^o-i=o i=owhere 4 1 ) and d1P are defined by (3.58) and (3.60). The values of uo and uN are knownfrom (5.8b). They can be combined into the first term on the left-hand side of (5.14).Equation (5.14) can be simplified toN- 1^(z) N-1r2 E ceu, = si —^E^= 1, • • • ,N — 1),^(5.15)ofj=i 1=1where(41) = 41 ) (i, j = 1, • • • , N — 1),^(5.16)4?) = d(i) (i,j= 1 ,^N — 1 ), (5.17)and(z)Si = —fi — r^(d(1)^d(1)d(1) )^2 (P )^)+ CPQs^ sO 0^iNUN r^io uo^iNuN) •Qs (5.18)CHAPTER 5: SOLVING GENERALIZED POISSON'S EQUATION, PART I 74Before the iteration process, we compute ai according to (5.18). The first guess valueis given arbitrarily. In each iteration we calculate the right-hand side of (5.15) from thecurrent ui first. Then we solve the resultant system of equations by the procedure describedby (3.79) and (3.80). From (5.15) the iteration process isN-1^N-1_.ye) = E [7.2 (g2)17 1 tip + E gi ju(ik-1) .us 1 spp=1^j=1The iteration matrix [gij1 is given byN-1[r2gi; = —P= 1From (5.20) we can compute the spectral radius of [Ni]. Because there is a factorap(r)^1 din arap r2)E 1^(z)ap 11)r = N-1 fp,' —1^(z)ffp^/1).ip^Crp^RI ip rap "i — p=1 L^.1X=M.(5.19)(5.20)the iteration matrix is model dependent, as is the spectral radius. We can build the iterationmatrices and compute the spectral radii for 2-D and 3-D problems similarly.Now we give an estimate of the convergence rate of the 1-D problem. We exchange theorder of the matrices on the right-hand side of (5.20) to obtain another matrixN-1(x)qij = E cri c) [pirlrai^ kjThe matrices [gii] and [Di] have the same eigenvalues (Wilkinson 1965, p54) and thereforethe same spectral radii. Let^N-1^—1bii =^11k) [12)1z_d wi "'^k=1^kjThenCr(x)^— ^iraj kJ.Numerical experiments show that the spectral radius of the matrix [k] is bounded by 1/7r.Since matrix [c4z)/rcri] is diagonal, the spectral radius is the maximum absolute value ofits elements. From Varga (1962, p5), we have^( z) ^In aPaqii]) P(Prai DPa^(dbiji) — max dx ) •rw(5.21)k=1CHAPTER 5: SOLVING GENERALIZED POISSON'S EQUATION, PART I^75That is, the spectral radius of the iteration matrix [gii]1^ a.c/lnPagii]) 5- ir— max I dx (5.22)The convergence rate R is defined as (Canuto et al. 1988, p138)R = —hip.The reciprocal of R measures the number of iterations required to reduce the error bya factor of e. The larger the convergence rate is, the fewer iterations that are requiredto obtain a solution to a given accuracy. For the 1-D problem under consideration, theconvergence rate is bounded byR .> ln(rx) — In (max din adx)(5.23)Example 5.1^In equation (5.10), we assume that0 = Ix I — 1 < x < 1},u(x) = 1+ ex 2 ,for the exact solutionu(x) = sin('x cos x).The boundary conditions and f(x) are given accordingly. Figure 5.1 shows the Lao norm,em., of the differences of u between two successive iterations for N = 4,8,16, and 32. Infigure 5.1, (a) is for 6 = 0.2 and (b) is for e = 1. In both (a) and (5) the solid lines onthe top show the results obtained by the Chebyshev spectral multigrid (CSMG) methodpresented by Zang et al. (1982) on an N = 32 grid. Figure 5.2 shows the spectral radiusfor N = 4 to 16. Table 5.1 shows the L2 and Lc,. error norms (comparing with the exactsolutions) and the CPU time on a SPARC station 2 computer after 9 iterations for boththe present method and the Chebyshev spectral multigrid method.155^ 10Number of iterations(a) E = 0.2tiN = 4N = 8N = 16N = 32N = 32 CSMGI CHAPTER 5: SOLVING GENERALIZED POISSON'S EQUATION, PART I^760N = 4N = 8N = 16N = 32N = 32 CSMG-5-1 0-155^ 10^15Number of iterationsFigure 5.1 The Lco norm, e., of the differences of u between twosuccessive iterations of example 5.1. (a) e = 0.2. (b) e = 1...........CHAPTER 5: SOLVING GENERALIZED POISSON'S EQUATION, PART I^770.25E = 0.20.200.15E= 1.0X0.100.054^ 8^ 12^16NFigure 5.2 The spectral radius of the iteration matrix for example 5.1.Table 5.1The error norms of example 5.1 after 9 iterations on an N = 32 1-D gridpresent method spectral multigrid methodL2 0.2 4.658 x10 -13 1.125 x 10—°5Lco 0.2 1.146 x 10 -12 2.295 x10 -05L2 1.0 4.200 x 10 -09 1.922 x10—°4Lco 1.0 4.950 x10—°9 3.595 x 10 —°4C PU (sec) 0.08 0.110.0021r =_^_ ,xi — x22=_^_^ .y1 — Y2CHAPTER 5: SOLVING GENERALIZED POISSON'S EQUATION, PART I^78For this 1-D example, the maximum absolute value of din cr/dx is f. From (5.23), theestimated convergence rate of the iteration matrix [gii] is R(e) > ln(r7r)— 0.51n e. ThereforeR(e) > 1.9494 for e = 0.2 and R(e) > 1.1447 for e = 1. From figure 5.2, the asymptotic valueof the spectral radius of the iteration matrix is 0.0325 for e = 0.2 and 0.125 for e = 1. Thereal convergence rate is R. = 3.4265 for e = 0.2 and R = 2.0794 for e = 1. The real values ofthe convergence rate are larger than the estimated bound. This can be explained as follows.The estimation (5.23) is based on an inequality condition. If dln o/dx is a constant, theequality sign is valid in (5.23). Only in this situation, is the estimate accurate.5.2.2 The 2-D ProblemFor the 2-D problem, equation (5.8a) becomesox( Ou Oucr02;) ( er .G) = i(x,Y) in fl, (5.24)where SI = {(x, y) I 22 < x < 21, f2 < y <fil}. We rewrite equation (5.24) as82u 02u f — ti(z)0.(z) — u(a) er(Y)= ^oz2 Oy2Let(5.25)(5.26)The 2-D Chebyshev grid is21 + 22 1^iir2 + r cos X ,yi = "Vi. + h2 + —1 cos 3--T .a N(5.27)Similar to the derivation for the 1-D problem in subsection 5.2.1, from (3.61), (3.62), and(5.25) we obtain a system of equationsN-1^N-1^cr(!) N-1^0.(1) N-1 11)r2^42)uti + 82^= ski _^ctiPuji — 8^al tia1=1 1=1 Crii 1=1^crii z 1=1(0 < i,j,k < N),^(5.28)CHAPTER 5: SOLVING GENERALIZED POISSON'S EQUATION, PART I^79where1sii = — [fii - ra,C.7 ) (4 10 )uoi duNi) - sal! ) (c401)ui0-1- guiN)}crii_2 (P)._^82 (p)uT VS0 'kW 'AIN "IV^k 30 su ' a3NuiN (5.29)We begin the iteration by computing^according to (5.29). The first value is givenarbitrarily. In each iteration we compute the right-hand side of (5.28) from the currentvalues of uji first. Then we solve the resultant equation system by the 2-D collocationmethod for Poisson's equation described by (3.85) and (3.86) to obtain the new values ofExample 5.2^In equation (5.24), we assume thatSi = {(x,y) I - 1 < x,y < 1},rr(x,Y) = 1 + 6(x2 + Y2 ),for the exact solutionu(x, y) = sin(r cos x) sin(r cos y).The boundary conditions and f(x, y) are given accordingly. Zang et al. (1982) solved thisproblem by the Chebyshev spectral multigrid method. Figure 5.3 shows the L oo norm,emax , of the differences of u between two successive iterations for N = 4,8,16, and 32. Infigure 5.3, (a) is for e = 0.2 and (b) is for e = 1. In both (a) and (b) the solid lines onthe top show the results obtained by the Chebyshev spectral multigrid method presentedby Zang et al. (1982) on an N = 32 grid. Figure 5.4 shows the spectral radius for N = 4to 16. Table 5.2 shows the L2 and L. error norms (comparing with the exact solutions)after 9 iterations and the CPU time on a SPARC station 2 computer for both the presentmethod and the Chebyshev spectral multigrid method. Probably due to the difference inprogramming and the selection of relaxation parameters, in table 5.2 the L2 error norm(1.708 x 10-5 ) for e = 0.2 of the Chebyshev spectral multigrid method after 9 iterations isbetter than that (3.16 x 10 -2 ) presented by Zang et al. (1982, p499).(a) E = 0.210Number of iterations(b) E = 1.0jp.0N = 4N = 8N = 16N = 32N = 32 CSMG10^ 15Number of iterationsCHAPTER 5: SOLVING GENERALIZED POISSON'S EQUATION, PART I^80Figure 5.3 The Lc° norm, en., of the differences of u between twosuccessive iterations of example 5.2. (a) 6 = 0.2. (b) 6 = 1.0.250.20Cl) 0.100.050.00CHAPTER 5: SOLVING GENERALIZED POISSON'S EQUATION, PART I^81E = 0.2E = 1.0X4^ 8^ 12^ 18NFigure 5.4 The spectral radius of example 5.2Table 5.2The error norms of example 5.2 after 9 iterations on an N = 32 2-D grid6 present method spectral multigrid methodL2 0.2 2.816 x 10 -13 1.708 x 10 -05L.. 0.2 1.525 x10 -12 8.399 x 10—°5L2 1.0 1.329 x10 -1° 1.172 x 10 -04L03 1.0 1.938 x10-1° 6.113 x10—°4CPU(sec) 0.68 2.28Let(5.32)2r=21 2— 22'CHAPTER 5: SOLVING GENERALIZED POISSON'S EQUATION, PART I^825.2.3 The 3-D ProblemFor the 3-D problem, equation (5.8a) becomesfox (er 0u) + 8 (0. 8u) + 8 (0. 8u) = f(x,y,z)^in II,^(5.30)where SI = {(x, y,^i2 < x < 21)^< y < gi, 22 < x < ill. We rewrite equation (5.30)as02u 82u 82uf — u(00.(2) — u(00.(v) — u(z)Q(z)0232 + 0y2 8x2 =^ (5.31)The 3-D Chebyshev grid is21 + 22 1^iirxi = ^2 ^-I- r cos --g, ,g 1 + v2 1^3wyi = 2 + ; cos —N— , (5.33)21 + 2.2 1^kwZk =^2 ^-I- t cos N— .Similar to the derivation for the 1-D problem in subsection 5.2.1, from (3.61), (3.62), and(5.31) we obtain an equation systemN-1^N-1^N-1r2 E gp,ul + 82 E 421)uilk + t2 E /2) .....kir us, /1=1^1=1^1=1(z) N-10. 0.C1.1) N-1sjk E p. ) . k -^e,scz _i= siik — r.)k NEli)....tNI U4^8 sik^--^ula mu.2,= E ii/i)u. — tcriik 1=1^crijk 1=1^Oink 1=1(0 < i, j, k < N),^(5.34)where1^+ J uNik) — sO)silk^ (y) (11)= — — 7.41 (41)tiojk -r csiN criik Jo 'Umcriik(z) Al)J ^JO)^1^2 td(2)^2(2)^)- "rijk kak0 Uji0 clieNthiN)1 r^uoik ceiNUNikid(2)^ 2)^id(2)^,(„2)— 2 k jo tii0k u3NUiNk) — 42 k k0 Uji0 "kNUjiN) • (5.35)CHAPTER 5: SOLVING GENERALIZED POISSON'S EQUATION, PART I^83Before the iteration process, we compute so according to (5.35). In each iteration wecompute the right hand side of (5.34) from the current values of ujik first. Then we solvethe resultant equation system by the procedure described by (3.91) and (3.92) to obtainthe new values of ujik.Example 5.3^In equation (5.30), we assume thatCl = {(x,y^-1 < x,y,z <11,z) = 1 + e(x2 + y2 + z2),for the exact solutionu(x, y, z) = sin(ir cos x) sin(ir cos y) sin(ir cos z).The boundary conditions and f(x, y, z) are given accordingly. Figure 5.5 shows the Lo„norm, em.„, of the differences of u between two successive iterations for N = 4,8,16, and32. In figure 5.5, (a) is for e = 0.2 and (b) is for c = 1. In both (a) and (b) the solid lines onthe top show the results obtained by the Chebyshev spectral multigrid method presentedby Zang et al. (1982) on an N = 32 grid. Table 5.3 shows the L2 and L, error norms(comparing with the exact solutions) and the CPU time on a SPARC station 2 computerafter 9 iterations.Table 5.3The error norms of example 5.3 after 9 iterations on an N = 32 3-D gridpresent method spectral multigrid methodL2 0.2 2.139 x 10 -13 3.049 x 10 -05Loo 0.2 2.266 x 10 -12 1.529 x 10 -04L2 1.0 3.189 x 10 -12 1.341 x10 -°4Loo 1.0 8.194 x 10 -12 7.113 x 10 -04CPU(sec) 32.38 112.28(a) e = 0.210—10 N =N = 8N = 16N = 32N = 32 CMG15—15^ `*10Number of iterations(b) e = 1.0N = 4N = 8N = 16N = 32N = 32 CSMG—1510^ 15Number of iterations— 10CHAPTER 5: SOLVING GENERALIZED POISSON'S EQUATION, PART I^8400Figure 5.5 The Loi, norm, emax , of the differences of u between twosuccessive iterations of example 5.4. (a) e = 0.2. (b) e = 1.CHAPTER 5: SOLVING GENERALIZED POISSON'S EQUATION, PART I^855.3 Discussion5.3.1 Convergence RateFrom figures 5.2 and 5.4 we see that for small N (< 8) as N increases, the spectralradius decreases. When N > 8, the spectral radius essentially reaches an asymptotic value.For the 1-D example 5.1, the asymptotic spectral radius is 0.0325 for 6 = 0.2 and 0.125for e = 1. For the 2-D example 5.2, the asymptotic spectral radius is 0.0335 for c = 0.2and 0.120 for 6 = 1. We have not computed the spectral radius for the 3-D example 5.3.From figures 5.1-5.5, however, we expect that the asymptotic spectral radius for the 3-Dexample 5.3 is smaller than 0.033 for e = 0.2 and smaller than 0.120 for e = 1. From figures5.1, 5.3, and 5.5, we see that unlike most iterative methods the convergence rate of thepresent method does not deteriorate as N increases. This is consistent with the foregoingdiscussion above the spectral radius. The convergence rate of the present method does notdeteriorate as the dimension increases either.The convergence rate depends on the behaviour of a. In examples 5.1-5.3 the constant6 determines the variations of Vu/Q. If e = 0, Val°. = 0. The problems become exactPoisson's equations. Thus, no iterations are required. When e increases, the variation ofVu/u also increases. The convergence rate decreases. This can be seen by comparing theresults for 6 = 0.2 and 6 = 1.0 in figures 5.1-5.5 and tables 5.1-5.4.In figures 5.1, 5.3, and 5.5, the convergence curves for N = 16 and N = 32 are almostidentical. This is due to the nearly identical spectral radius. When N = 4 all the eigenvaluesof the iteration matrix are real. The convergence curves for N = 4 are nearly straight lines.The wave-like convergence curves for N = 8 are presumably due to the dominant complexeigenvalues of the iteration matrix. This can be explained as follows. Because the iterationmatrix is real, all the coefficients of the characteristic equation IG — alb = 0 are real too.If there are complex eigenvalues, they must appear in conjugate pairs. The correspondingCHAPTER 5: SOLVING GENERALIZED POISSON'S EQUATION, PART I^86eigenvectors are also in conjugate pairs. In the iteration process, the imaginary parts in theerror cancel each other due to the conjugate property. When there are dominant complexeigenvalues, from (5.7) the error will be dominated by them. A complex number z = a + ibcan be written as reie, where r = (a + b)'! 2 and 8 = tan-1 (b/a). The real part of z" iscos(n0). It is the phase factor, cos(n0), that makes the convergence curves wave-like.5.3.2 Comparison with Chebyshev Spectral Multigrid MethodComparing the convergence curves of N = 32 cases of the present method with thecorresponding curves of the Chebyshev spectral multigrid method in figures 5.1, 5.3, and5.5, we see that the convergence rate of the present method is much better than that of theChebyshev spectral multigrid method. In tables 5.1-5.3, after 9 iterations the error normsobtained by the present method are about 5-8 orders better than those obtained by theChebyshev spectral multigrid method. The speed of the present method is also superior tothe Chebyshev spectral multigrid method. For the 2-D and 3-D problems, the CPU timeof the present method is less than one third of that of the Chebyshev spectral multigridmethod for the same number of iterations and the same number of grid points. Because theproblem is solved on a single grid, the programming for the present method is much simplerthan that for the Chebyshev spectral multigrid method, especially for 3-D problems.5.3.3 Computational OperationsSuppose that d represents the dimension of the problem. The solution for the Poisson'sproblem requires 2d(N —1)d+ 1 d(N —1)d operations (one addition and one multiplicationcount for one operation). The evaluation of Vu requires d(N — 1)d+1 operations. In eachiteration the present algorithm requires total 2dN(N —1) 1 operations. The factorization ofmatrix 4t ) requires 8(N-1)s operations (Golub and Van Loan 1983, Algorithm 7.6-3). Forthe 1-D problems the factorization may be too expensive. For the 2-D and 3-D problemsthe cost of the factorization process is relatively small. The total operation count for 2-DCHAPTER 5: SOLVING GENERALIZED POISSON'S EQUATION, PART I^87and 3-D problems is of 0(Na+ 1 ). In a practical application we can build a lookup table tostore the factorized matrices and the eigenvalues for commonly used Ns.CHAPTER 6SOLVING THE GENERALIZED POISSON'S EQUATION, PART II:MULTI-DOMAIN METHODIn chapter 5 we solved the generalized Poisson's equation V • (crVu) = f by the single-domain Chebyshev spectral method. When the potential u is infinitely differentiable, thesingle-domain Chebyshev spectral method can achieve exponential convergence. Whenthere are discontinuities in a, because Vo does not exist at the discontinuities, the single-domain Chebyshev spectral method does not work. To overcome this problem, we nowdiscuss how to solve the generalized Poisson's problem by multi-domain Chebyshev spectralmethod. The ideas are similar to those discussed in chapter 4 for Poisson's equations. Insections 6.1-6.3, we discuss 1-D, 2-D, and 3-D problems respectively. Numerical examplesare given in each situation. We discuss the results in section 6.4.6.1 The 1-D ProblemWe consider the 1-D boundary value problemd I duNVia ; ) = f(x)u = b(x)^on r,in 0, (6.1a)(6.1b)where a(x), f(x), and b(x) are known functions, a > e > 0, 12 = {z I iBr+i < x < 21},r is the boundary, B. is an integer, and ii+i < i i (1 = 1, • • • , Be ). We divide 12 into B.subdomains: Ill = {x I ii+i < x < ii} (1 = 1, • • • , B., see figure 4.1, where B. = 3).In subdomain 0 1 , u = ul(x), a = a'(x), and f = fl(x). We assume that al(x) and88CHAPTER 6: SOLVING GENERALIZED POISSON'S EQUATION, PART II^89fl(x) are smooth functions and that there are some discontinuities in a(x) and/or f(x) atx = xi (I = 2, • • • , Be ). The boundary value problem (6.1) is equivalent to the followingproblem.( duldx eri dx) = -fi^in Il l (I = 1, • • • , BO,^(6.2a)ui—i = ui^at x = xi (1 = 2, • • • , BO , (6.2b)i_ i dul-1^i dula^dx— = a dx^at x = xi (1 = 2, • • • , BO ,^(6.2c)u = b on r. (6.2d)Equation (6.2a) is the governing equation (6.1a) in each subdomain. Equation (6.2d) isthe same as (6.16). Equations (6.26) and (6.2c) are the interface conditions, that is, thecontinuity of the potential u and the flux a du/dx.For simplicity, we use the same Nth degree Chebyshev polynomials for all subdomains.Letting2ri =_^ 7X1 - X1+1the 1-D Chebyshev grid is given byxs^2= 21 + 1+1 -I- —1 cos N— (1 = 1,^, B.; i = 0, • • • , N).^(6.3)We assumeand(Gq)(z) = da •dxIn numerical computation, equations (6.26) can be satisfied by taking the same value forulp.7 1 and tito (1 = 2, • • • , Be ). By a derivation similar to that in subsection 5.2.1, from (3.61),(3.62), (3.67), (3.68) and (6.2), we obtain the system of equations:A21 1^(04^^N^r? E^-^ri )(z)s = 1, • • • ,Bz ; i = 1, • • • ,N —1),^(6.4a)^j=0 s j=0CHAPTER 6: SOLVING GENERALIZED POISSON'S EQUATION, PART H^90N^Nri_ia1N^O_1 E dop)Up1 H) + rjao E dOpo)Up —_ 0^(1 = 2, • • • , B.).^(6.4b)p=0 p=0We solve system (6.4) by the influence matrix method outlined in section 4.1 using thefollowing steps:1. Assign 0 to all the interface points. Then separate equations (6.4a) into B. independ-ent smaller systems, and solve them by the single-domain Chebyshev spectral methoddiscussed in subsection 5.2.1. Compute the jump of the flux at the interface points,that is, the left-hand side of (6.4b) and denote them by a vector di.2. Build the influence matrix [cif] by solving the generalized Poisson's equation with spe-cial interface conditions. We compute the jth column of the influence matrix [cii] asfollows. We assign 1 to the jth interface point and zero to all others, solve the boundaryvalue problem, compute the jumps of the flux, and denote them by another vectorThe difference ti — di is the jth column of the influence matrix [cii].3. For the 1-D problem, the influence matrix [cii] is tridiagonal. We solve the systemBr -1E =j=1to obtain v. Then we correct the interface values with v.4. Solve the Dirichlet problem again to obtain the correct solution.Now we consider a 1-D example. In this example, f(x) is a smooth function on the wholedomain, a(x) and u(x) are piecewise-smooth functions. We find the analytic solution first.Then we solve the problem by the multi-domain Chebyshev spectral method and comparethe numerical solution with the analytic solution.Example 6.1 In the boundary value problem (6.1), we assumefi={x1-3<z<3},f(x)= h(z)= sin(z),^ (6.5)CHAPTER 6: SOLVING GENERALIZED POISSON'S EQUATION, PART II^91anda(z) = 72(x) =u(3) = 1,u(-3) = —1.10 1 + 0.8 cos x11 + 0.8 cos x51 -I- 0.8 cos x2 < < 31 < x < 23 < < —1,(6.6)We first solve this problem analytically. In each subdomain we integrate both sides of(6.2a), then divide both sides by o, and finally integrate both sides again to obtain a generalsolution. We determine the constants in the general solution by the interface conditions(6.2b), (6.2c) and the external boundary conditions (6.2d). The analytic solution iswith1r1)ct(z) + ciu(x) = u2(x) = a(x) + c2a(x)i N + C302 < x < 3—1 < < 2—3 < x < —1(6.7)a(x) = co(x + 0.8 sin x) — (sin x + 0.2 sin 2x + 0.4x),co = 1.0421556809681,ci = 0.80411146510375,= —0.35191875057152,andc3 = —0.60822293020750.In (6.5)—(6.7), we denote the functions by special names f2(x), cr2(x), and u2(x). Todistinguish them from those defined in chapter 4, here we use the index 2. These functionssatisfyi17; V727 ) = n(z).d^du2\^ (6.8)CHAPTER 6: SOLVING GENERALIZED POISSON'S EQUATION, PART II^92f2(x)Figure 6.1 The functions 12(x), cr2(x) and u2(x).CHAPTER 6: SOLVING GENERALIZED POISSON'S EQUATION, PART II 93We plot their values in figure 6.1. The function f2(x) is smooth on the whole domain. Thefunction a2(x) is discontinuous at points x = 2 and x = —1. From (6.2c), the first-orderderivative du/dx is also discontinuous at these two points. We will use f2(x), 0r2(x), andu2(x) to construct 2-D and 3-D piecewise-smooth functions later.In the multi-domain Chebyshev spectral method, we divide the whole domain CZ into3 subdomains: = 2 < x < 3}, n2 - - {xI - < < 2}, and f/3 = I — 3 < <—1}. In each subdomain we approximate u(x) by an Nth degree Chebyshev polynomials.Table 6.1 shows the relative error norms of the numerical solutions obtained by the multi-domain Chebyshev spectral method. In table 6.1, N represents the degree of Chebyshevpolynomials used in each subdomain, B. the number of subdomains along the x-direction,Pt the total number of grid points, L2 and Lo. the corresponding error norms, and CPUtime is obtained on a SPARC station 2 computer. From table 6.1, the multi-domainChebyshev spectral solutions exhibit exponential convergence. On a grid with 37 points,the error norms are ti 10-8 .Table 8.1The relative error norms of example 6.1 (MDCS method)N B. Pt L2 Loo4 3 13 1.839 x 10 -02 2.952 x 10 -025 '3 16 2.867 x 10-03 3.557 x 10 -836 3 19 1.213 x 10-83 1.856 x 10-037 3 22 1.882 x 10-04 2.316 x 10-048 3 25 4.651 x 10 -88 6.723 x 10 -859 3 28 7.059 x 10 -88 8.468 x 10 -8810 3 31 1.133 x 10 -88 1.585 x 10 -0811 3 34 1.613 x 10 -07 1.901 x 10-0712 3 37 1.897 x 10 -88 2.604 x 10-88a — = 0^at y =^(1 = 1, • • • , B.; m = 2, • • • , By), (6.10e)(6.10f)In each subdomain, we approximate the potential function by the same Nth degree Cheby-shev polynomials. Assumingwhere f (x, y), b(x, y), and cr(x, y) are known functions, o(x, y) e > 0, SZ = {(x, y) 2,13.41< x < 21, < Y < gi}, r is the boundary, and B. and By are integers. We divide thedomain into B.By subdomains with the subdomain boundaries x = it (1 = 1, • • • , B. + 1)and y = gni (m = 1, • • • , By + 1) (see figure 4.4, in which B. = By = 2). We denote thesubdomain located at the lth column and mth row as Slim. In film, u = utm, o. = all", andf = tm. The boundary value problem (6.9) is equivalent to the following problem.0 ( im Oulm) 49^hn &Lim^thnex VY Ox ) ay e7 ey =in Sl im (1 = 1, • • • , B.; m = 1, • • • , By )^(6.10a)u(1-1)m = ulm^at x = xj (1 = 2, • • • , B.; m = 1, • • • , By ),^(6.10b)km-1)^ntu^= uIo0,(1-1)m u(i-i)mOx0.l(m-1) Oul(m-1)ay^ayu = b(x, y)^on r.at y =^(1 = 1, • • • , B.;crIm^ =0^at x =^(1 = 2, • • • , B.; m = 1, • • • , By ), (6.10d)lm aulmm = 2, • • • , By ),(6.10c)CHAPTER 6: SOLVING GENERALIZED POISSON'S EQUATION, PART II^948.2 The 2-D ProblemThe basic ideas for the 1-D problem described in the previous section can be extendedto the 2-D and 3-D problems. Certainly there are special problems to deal with in eachcase. We consider the 2-D generalized Poisson's boundary value problemf Ou\ 0 f Ou\) = f (;Y)u = b(x, y)^on r,in CZ, (6.9a)(6.9b)2rj = _^zi -28m = _ ^,Ym - Ym-1-1CHAPTER 6: SOLVING GENERALIZED POISSON'S EQUATION, PART II^95the 2-D Chebyshev grid is given byIx! = 21 +:1+1 + 1 cos ilrri^N 'm gm + gm+1 , 1^j's-Yj =^-t- — cos —^2^am N .utri = u (.!, yr) ,o•qr = cr (4, yr) ,^( 0.!,,,)(x)^I= 8asj ,°° KIT)((Y)^8atrtrt) — ____8y (4 ,Yr) From (6.10f) the values of u at the external boundary points are known. Equations (6.105)and (6.10c) can be satisfied by taking the same interface value for the two contiguoussubdomains. By a derivation similar to that in subsection 5.2.2, from (3.61), (3.62), (3.63),(3.67), (3.68) and (6.10), we obtain the system of equations:N^N^ (lm^(2)aim) N. .2 E do) lm 2 E d(2) lma,.. __ 4i, ,ri^ip up' + 3 TA^" Sp^i m ri elm E dVup17. •p=0 p=0 13^ ii^p=0I^(Y)(Or) N3711 elm^ E c6p1)4";-ii^p=0(1= 1, • • • , B.; m = 1, • • • , By ; i, j = 1, • • • ,N —1),^(6.11a)N Nri_laWm E d(oip) u(i;ii)„7, + riot! E ioip) uip7 = 0p=0^ p=0(1 = 2, • • • , B.; m = 1, • • • , By; j = 1, • • • , N —1),^(6.11b)N N8m_iail( -1) E cpuirN:p1 + „motEdo)4,77 = 0p=0^ p=0(1= 1, • • • , B., m = 2, • • • , By ; i = 1, • • • ,N — 1), (6.11c)N N8m_ia ign-1) E d(1)u1471,),) + 8„,,rg E d(10„) uionpi = 0p=0^ p=0(1= 2, • • • , B.; m = 2, • • • , By ),^ (6.11d)Letandup =0^= 2, • • • ,Bx; m= 1, • • • , By ; j = 1, • • • ,N — 1)u1S =0^(1 = 1, • • ,Bm ; m = 2, • • • , By ; i = 1, • • • , N — 1).andCHAPTER 6: SOLVING GENERALIZED POISSON'S EQUATION, PART II^96where 41) is defined by (3.58) and d (j) is defined by (3.60). Equations (6.11a) are forthe interior points, (6.11b) and (6.11c) for the interface points, and (6.11d) for the cornerpoints. The corner points do not appear in equations (6.11a)—(6.11c).We solve the sub-system consisting of (6.11a)—(6.11c) by the influence matrix methodfirst using the following steps:1. Change (6.11b) and (6.11c) toSolve the new system formed by (6.11a) and (6.12) in each subdomain independentlyby the single-domain Chebyshev spectral method discussed in subsection 5.2.2. Thencompute the jump of the flux at all the interface points and denote them by a vector2. Build the influence matrix [cii] by solving the generalized Poisson's equation with spe-cial interface conditions. We compute the jth column of the influence matrix [cif] asfollows. We assign 1 to the jth interface point and zero to all others, solve the boundaryvalue problem, compute the jumps of the flux, and denote them by another vectorThe difference ti - di is the jth column of the influence matrix [45].3. For the 2-D problems, the influence matrix [cii] is sparse and diagonally dominant. Wesolve the influence matrix system by Gauss-Seidel iterative method to obtain v. Thenwe correct the interface values in (6.12) with qi.4. Solve the system formed by (6.11a) and (6.12) again to obtain the correct solution forall the interior and interface points.The values at the corner points can be found by substituting the known values of theinterface points into (6.11d) and solving the resultant smaller systems.CHAPTER 6: SOLVING GENERALIZED POISSON'S EQUATION, PART II^97Now we consider a 2-D example. In this example, the analytic solution is a piecewise-smooth function.Example 6.2 In the boundary value problem (6.9), we assume thatSZ = {(z, y) I — 3 < x,Y < 3},cr(x,y) = 02(x)472(Y),^ (6.13)and the analytic solution isu(x,y) u2(x) u2(Y).^ (6.14)where u2 and o.2 are defined by (6.6) and (6.7). From (6.8) and (6.9a) we havef(x,y) = f2(x)u2(y)cr2(y) + u2(x) h(Y) cr2(x),where 12 is defined by (6.5). The boundary conditions are given accordingly. Table 6.2shows the relative error norms of the multi-domain Chebyshev spectral solutions. In tableTable 6.2The relative error norms of example 6.2 (MDCS method)N Bx By Bt Pt K w L2 L. CPU(sec)4 3 3 9 169 9 1.150 7.304 x 10—°3 3.451 x 10 -02 0.225 3 3 9 256 11 1.240 1.190 x 10—°3 4.073 x 10 -03 0.506 3 3 9 361 16 1.290 4.425 x 10 -04 2.086 x 10 -03 1.097 3 3 9 484 20 1.340 6.365 x 10—°3 2.489 x 10 -04 2.228 3 3 9 625 26 1.390 1.674 x 10—°3 7.454 x 10 -05 3.969 3 3 9 784 28 1.400 2.225 x 10 -06 9.159 x 10-06 5.7010 3 3 9 961 33 1.430 4.122 x 10—°7 1.749 x 10 -06 9.5411 3 3 9 1,156 39 1.450 4.996 x 10 -08 2.051 x 10—°7 15.0412 3 3 9 1,369 41 1.470 7.033 x 10—°9 2.898 x 10 -08 20.82CHAPTER 6: SOLVING GENERALIZED POISSON'S EQUATION, PART II^986.2, N is the degree of Chebyshev polynomials used in each subdomain along one direction,B. and By the number of subdomains along the x- and y-directions respectively, Bt thetotal number of subdomains, Pt the total number of grid points, K the number of iterationsfor solving the influence matrix system, w the relaxation parameter, L2 and Loo are thecorresponding error norms, and CPU time is obtained on a SPARC station 2 computer.From table 6.2, we see that the error norms of the 2-D multi-domain Chebyshev spectralsolutions go to zero exponentially. On a grid with 1,369 points, the error norms are ti 10-8 .6.3 The 3-D ProblemWe consider the 3-D generalized Poisson's boundary value problemOu\ 0 Ou\ 0 Ou\ax VTG.)^1 cr $3-i) = f(z ' Y ' z)u = bi(x, y, z)^on r (except on z = xi plane),o-OnOu = 62(x, y, z)^on z = xi planein ft,^(6.15a)(6.15b)(6.15c)where SZ = {(x, y, z) 2B.+1 < x < 21, gB,+i < Y <^2B:+1 < z < z1}, r is theboundary, B., By and B. are integers, f , y, z), cr(x,y, z),^z), and b2(x, y, z) areknown functions, and cr(x,y, z) > e > 0. For the convenience of further application ofthis method to 3-D DC resistivity modeling, we specify the boundary conditions on the topplane by the Neumann type. We divide the domain into B.ByB, subdomains with dividingplanes x = ti (1 = 2, • • • , y = ym (m = 2, • • • , By), and z = in (n = 2, • • , B.). Thesubdomain located at the lth column, mth row, and nth layer is denoted by Oh'. In Oh',u uhnn , almn, and f flmn . Boundary value problem (6.15) is equivalent to thefollowing problem.Ox^ax^ay^ay ) az k^az )a( ahnnauhnn) + a (ahnnatilmn) + 8 ( crimnOulmn ) = flmn^in lmn (1 = 1, • • • , B.; m =^, B; n = 1,• • • , B4,^(6.16a)U(1-1)ffin = 241mnCHAPTER 6: SOLVING GENERALIZED POISSON'S EQUATION, PART II 99at x =^(1 = 2, • • • , B.; m = 1,u l(m— 1)n = ulmnat y =^(1 = 1, • • • , B.; m = 2,ulm(n-1) = ulmnat z =^(1 = 1,• • • , B.; m = 1,• • • ,By ; n = 1,•• • , By ; n = 1,•• • , By ; n = 2,•••,•• • , Bz ),•• • , Bz ),(6.16b)(6.16c)(6.16d)(1- 1 )ffin 8u(1-1)nin int?, Oulmn 00.=• • • , By ; n = 1, • • • , Bz ), (6.16e)Oxat x =^(1= 2, • • • , .13z ; m = 1,00.1(m-1)n 8likm-1)n^almn°Ulnin =• • • , By ; n = 1,•• • , By ; n = 2,• • • ,• • • ,(6.16f)(6.16g)Oy Oyat y = §,,, (1 = 1, • • • , B.; m = 2,Oulm(n-1)^°Ulmn 0cri,n(n-1) =ezat z = in (1 = 1,• • • , B.; m = 1,°uhniahni^=^Y) z)^(1=^...,Bs; m = 2 , • • • , By), (6.16h)b2(z, 1,u = bi(x, y, z)^on r (except on the ground surface). (6.16i)In each subdomain, we approximate the potential function by the same Nth degree Cheby-shev polynomials. Assuming^2^= ^_^21 — 21+1Ym Ym+1^t o =^=22-the 3-D Chebyshev grid is1 xl = 21 + 2 21+1 + 71 cos1^ilrN 'Yi =^1-^m gm + gm+1 ,^cos —1^jr2^8m^N 'n 2n + 4+1 . 1^kirZic =^1- COS — •2^tn^N 2n 2n+1(4yr,z;)CHAPTER 6: SOLVING GENERALIZED POISSON'S EQUATION, PART II^100LetutIgn = u (4, yr,^,az7gn = (x!, yr, .41) ,imn (z)VT° =(Y)^acr(orkn) =ayand)(z)(criikn^=49Z (4,yr,q)From (6.16i) the values of u at the external boundary points (except for those points onz = x1 plane) are known. Equations (6.16b)—(6.16d) can be satisfied by taking the sameinterface value for the two related subdomains. By a derivation similar to that in subsection5.2.3, from (3.61), (3.62), (3.67), (3.68) and (6.16), we obtain the system of equations:N^N^NE 42 )upiw sm2 E (42p)47 tn2 E di2p)utnpnp=0^p=0^p=0fz\ N1^lmn ftir^(44rnk) ‘ E dTti'p'N —cro^ p=0N(a!rkn) ‘ E ,gyp uipkp=0— tn (41,n) — ciklp)u!rpnl( . 1 Np=0(1= 1, • • • , B.; m = 1,. • • , B; n = 1,• • , Bz;N^ N(1-1)mn V-N j1) ..0-1)mn^njmn V• 11)^= 0r1- 1 47Nik^u'Op^• 1" Ojk L. Op "pjlep=0 p=0(1 = 2, — • , B.; m = 1, • • • , By ; n = 1, • • • , Bz; j,Nkm-1)n VN /1) 1(m-1)n^Irnn [11 P)4,1Inn 08m—laiNk a^0—p)k Sm (7i0k La 'bop "'iP14p=0^ p=0(1 = 1, • • • ,^m = 2, • • • , By; n = 1 , • • • , Bz;N Nlm(n-1)^4(1)..1m(n-1)^,ltrin^d(1)ulynn1n-1 eriiN^U'Op NAN —p) i'n‘'ijO ^Op s3Pp=0 p=0(1= 1, • • • ,13.; m = 1,• • • ,B; n = 2, • • • 7 Bz;j, k = 1, • • • , N — 1), (6.17a)k = 1, • • • ,N — 1), (6.17b)k = 1, • • • , N — 1), (6. 17c)j = 1, • • • ,N —1), (6.17d)CHAPTER 6: SOLVING GENERALIZED POISSON'S EQUATION, PART II^101N44V E 101p)utp; =^yr ,p=0(1= 1, • • • , B.; m = I,. • • , By ; i,j = 1,• • • ,N — I),NdiCi.p)u(ovi i)pm)onk riot:^Al)E= 0 op 7dpOkp=0• • , B.; m = 2, • • , By ; n = 1, • , Bz ; k = I, • ,N — I),NriakTon E utrze 0p=0• • , By ; n = 2, • • • , Bz ; j =NBmotionE 419,40zn 0p=0• • , By ; n = 2, • • • , B.; i = 1, • • • ,N — I),N°moron dWp) u icon = 0p=0• • , By ; n = 2, • • • , /3,),(1-1)mn^(l-1)mnri- laNJO^"'Op U(N -p)j0p=0(1 = 1, • • • , B.; m = I,.Nl(m-1)n V-■ p)ul(m-l)nsm_icrilVO^Op i(N-p)0p=0(1 = 1, • • • , B.; m = 2, •N8m-100N01(7-1)n^p) ul(m-l)n,14 Op 0(N-p)0p=0( 1 = 2, • • • , Bz ; m = 2, •,N — I),(6.17e)(6.171)(6.I7g)(6.17h)(6.17i)where dIP is defined by (3.58) and dli ) is defined by (3.60). Equations (6.17a) are for theinterior points, (6.17b)—(6.17d) for the interface points, (6.17e) for the points on z =plane, (6.17f)—(6.17h) for the edge points, and (6.17i) for the corner points. The edgepoints and corner points do not appear in equations (6.17a)—(6.17e).We solve the sub-system consists of (6.17a)—(6.17e) by the influence matrix methodfirst using the following steps:1. Change equations (6.17b)—(6.17e) tou 107: = 0(1 = 2, • •lmn•, B.; m = 1, • • •, By ; n = 1, • • • , Bz ; j, k = 1, • • • , N — I), (6.18a)Ui0k = 0(1 = 1, • • • , Bz ; m = 2, • • • , By; n = 1, • • • , Bz;^k = 1, • • • , N — 1), (6.180Utrzon = 0(1 = 1, • • • ,^m= 1,• •• ,By ; n = 2, • • • , Bz; i,j = 1,• •,N —1). (6.18c)CHAPTER 6: SOLVING GENERALIZED POISSON'S EQUATION, PART II^102u re 0(1= 1,• • • ,B. ; m = 1,• • • ,By; j = 1,• • • ,N — 1). (6.18d)In the new system formed by (6.17a) and (6.18), the unknowns in different subdomainsare no longer coupled. We solve them in each subdomain independently by the single-domain Chebyshev spectral method discussed in subsection 5.2.3. Then we compute theleft-hand side terms in (6.17b)—(6.17e) to obtain the jump of the flux in two contiguoussubdomains related to a particular interface point and denote them by a vector di.2. Build the influence matrix [cii] by solving an generalized Poisson's equation with specialinterface conditions. We compute the jth column of the influence matrix [cif] as follows.We assign 1 to the jth interface point and zero to all others, solve the boundary valueproblem, compute the jumps of the flux, and denote them by another vector ti. Thedifference ti — di is the jth column of the influence matrix [cii].3. For the 3-D problem, the influence matrix is sparse and diagonally dominant. Solvethe influence matrix system by Gauss—Seidel iteration to obtain the correction vectorvi. Then we correct the interface values in (6.18) with vi.4. Solve the system formed by (6.17a) and (6.18) again to obtain the correct solution forall interior and interface points.We compute the values at the edge points as follows. We substitute the known valuesinto (6.17e)—(6.17g) and move the known terms to the right-hand side. Equation (6.17e) canbe separated into (By — 1)B.(N— 1) independent smaller systems. Equation (6.171) can beseparated into By(Bz —1)(N-1) independent smaller systems. Each of the smaller systemscontains (B. —1) equations for (B.-1) unknowns. Equation (6.17g) can be separated into(B, —1)Bz (N — 1) smaller systems, each of which contains (By — 1) equations for (By — 1)unknowns. All these smaller systems are tridiagonal and are easy to solve.We compute the values at the corner points in a similar way. Substituting the knownedge values into (6.17h) and moving the known terms to the right-hand side, we obtainCHAPTER 6: SOLVING GENERALIZED POISSON'S EQUATION, PART II^103(B. — 1)(13, — 1) independent tridiagonal smaller systems, each of which contains (B y — 1)equations for (By — 1) unknowns. Solving them, we obtain all the unknowns.Now we consider a 3-D example. In this example, the analytic solution is a piecewise-smooth function.Example 6.3 In the boundary value problem (6.15), we assume that= {(x z) — 3 < xl Y/ z < 3 }7the analytic solution isu(x , y, = u2(x) u2(y) u2(x),^ (6.19)andcr(x , y, = o2(x) o2(y) 02(x) ,^ (6.20)where u2 and 02 are defined by (6.7) and (6.6). From (6.15a) and (6.8) we havef (x, y, z) =f2(x)u2(y)u2(z)(72(y) 472(z) + u2(m) h(y)u2(z) (72(x) 02(z)u2(x)u2(Y) .f2(z) 02(x) 02(y),where f2 is defined by (6.5). The boundary conditions are given accordingly. The domain isdivided into 3 x 3 x 3 subdomains with dividing planes 22 = 2, 23 = —1, in = 2, = —1,z2 = 2, and 23 = —1. Table 6.3 shows the relative error norms of the multi-domainChebyshev spectral solutions. In table 6.3, N is the degree of Chebyshev polynomialused in each subdomain along each coordinate direction, B., By , and B. the number ofsubdomains along the y-, and x-directions respectively, Bt the total number of thesubdomains, Pt the total number of grid points, K the number of iterations for solvingthe influence matrix system, w the relaxation parameter, L2 and L. the error norms, andCPU time is obtained on a SPARC station 2 computer.CHAPTER 6: SOLVING GENERALIZED POISSON'S EQUATION, PART II^104Table 6.3The relative error norms of example 6.3 (MDCS method)N Bx By Bz Bt Pt K w L2 Lc° CPU(sec)4 3 3 3 27 2, 197 14 1.390 3.340 x 10 -02 6.055 x 10-03 3.615 3 3 3 27 4,096 21 1.460 5.535 x 10 -03 1.003 x 10-03 17.056 3 3 3 27 6,859 25 1.510 2.038 x 10 -03 3.741 x 10-04 52.577 3 3 3 27 10, 648 33 1.530 2.766 x 10-04 4.947 x 10-°8 177.278 3 3 3 27 15,625 46 1.570 7.448 x 10-08 1.443 x 10 -08 535.006.4 DiscussionIn this chapter, we solve the 1-D, 2-D, and 3-D generalized Poisson's equation by multi-domain Chebyshev spectral method. We divide the whole domain into many subdomainsand approximate the unknown function by a separate Chebyshev expansions in each sub-domain. We solve the discrete system by the influence matrix method. We build theinfluence matrix by solving the generalized Poisson's equations with special boundary con-ditions. Unlike the procedure for Poisson equation in chapter 4, the entries of the influencematrix are the difference between the flux jump vector ti and the reference vector di. Wesolve the influence matrix system by Gauss-Seidel iteration.In examples 6.1, 6.2, and 6.3, the analytic solution is infinitely differentiable in eachsubdomain. From tables 6.1-6.3 we see that when N increases, the error norms go to zeroexponentially. In the 1-D example 6.1, the relative error norms on a grid of 37 points(N = 12) are ti 10-8 (see table 6.1). In the 2-D example 6.2, the relative error norms ona grid of 1,369 points (N = 12) are e•-• 10-8 (see table 6.2). In the 3-D example 6.3, therelative error norms on a grid of 15,625 points (N = 8) are •-•-• 10 -5 (see table 6.3). From theabove discussion we conclude that when there are discontinuities in u, the multi-domainChebyshev spectral method can achieve exponential convergence for the canonical modelof a cube.CHAPTER 73D GRAVITY MODELINGIn the geophysical literature, forward gravity modeling falls into one of three classes.The first class consists of analytic methods, in which solutions are obtained by assuminga constant density distribution and integrating over certain geometrically simple spatialdomains. Formulae have been published for simple geometrical bodies, such as the circulardisk, straight rod, solid ellipsoid, sphere, spheroid, rectangular parallelepiped, horizontalcylinder, vertical line element, vertical cylinder, vertical fault, vertical sheet, and the rightcircular cone (Nettleton 1942; MacMillan 1958; and Talwani 1973). Nagy (1966) calculatedthe gravitational attraction due to a rectangular body. Barnett (1976) modeled a polyhedralbody composed of triangle facets. Okabe modeled an arbitrary polyhedral body (1979a)and a 3D solid of revolution (1982). Cady (1980) computed the vertical gravity field of abody with constant polygonal cross-section and finite strike-length. Gaze and Lahmeyer(1988) computed the potential of a triangular prism and suggest that a polyhedral bodycan be represented as a combination of many triangular prisms.In the second class of gravity modeling, solutions are obtained using approximate meth-ods, in which an arbitrary body is subdivided into geometrically simple bodies such ascubes, rectangular parallelepipeds and prisms. Talwani and Ewing (1960) cut an arbitrary3D body into horizontal slices, approximated each slice by a polygon and then computedthe gravity field by integration along the edges of the polygon. Botezatu et al. (1971) usedmany cubes to construct an arbitrarily shaped body. Mufti (1975) developed an approxi-mate formula for the gravitational field of a cube and used it to compute the gravitationalfield of an arbitrarily shaped body.105CHAPTER 7: 3D GRAVITY MODELING^ 106The Fourier transform technique is used in the third class of gravity modeling. Odegardand Berg (1965) developed methods for analyzing gravity data in the frequency domain.Bhattacharyya and Navolio (1975, 1976) used a digital convolution method to compute thepotentials and fields due to arbitrarily shaped bodies with homogeneous density distribu-tions. Pedersen (1978) derived the spatial Fourier transform of Barnett's results, pointingout that the required calculations are considerably simpler in the frequency domain andare therefore correspondingly faster. Hansen and Wang (1988) recast Pedersen's Fouriertransform expression in a simpler, coordinate-invariant form.In this chapter we apply the Chebyshev spectral methods discussed in chapters 3 and4 to 3D forward gravity modeling. In section 7.1 we outline the boundary value problem.In section 7.2 we present the methods to realize 3D point-injection and cell-average dis-cretizations. In section 7.3 we extend Okabe's (1979a) analytic expression of the gravityfield of a homogeneous polyhedron to the potential field. In section 7.4 we discuss how tocompute the boundary conditions approximately by multipole expansions. In section 7.5we present three examples of 3D gravity modeling. In section 7.6 we discuss the resultsand give the conclusions.7.1 The Boundary Value ProblemThe gravitational potential u satisfies Poisson's equation (Telford et al. 1976). Thefundamental problem of forward gravity modeling is to solve boundary value problemV2u = —47rGp^in 0 .^ (7.1a)u = known^on r, (7.1b)where G is the universal gravitational constant and p is the density. The gravitational field,g, can be found byg = Vu.CHAPTER 7: 3D GRAVITY MODELING^ 107On the interface common to two media the potential and the first-order normal derivativesare continuous. That is,ui =112,^(7.2)andOui^0112On^On ' (7.3)where n is the normal vector of the interface.7.2 Point-Injection and Cell-Average DiscretizationIn order to solve boundary value problem (7.1) by numerical methods, we have todiscretize the density function p. For 3D problems the point-injection discretization isdefined as(P)^lPiik = "Xi, Y j t Zk) (0 < i, j, k < N).^(7.4)For simplicity in the mathematical treatment, the 3D cell-average discretization is definedas follows. Let i = 0,di =N ( — il.)11COSy^i=1,•••,N, (7.5)—1 i=N+1,andhi = di — di+1^(i = 0, • • • , N).At a point (xi, yi, 4), we take a small rectangular parallelepiped cell confined by di+i <x < di, (11+1 < y < d1, and dk+i < z < 4. The volume of this cell isVijk = hihihk .The discrete value is1p13 = T . . ^p(z , y, z)dV .rv i.7 k lilts(7.6)CHAPTER 7: 3D GRAVITY MODELING 108The models to be considered in this chapter have a piecewise-constant or piecewise-smooth density distribution. Each piece will be confined by a polyhedron. What followsare some fundamentals of analytic geometry. We briefly outline them here as a basis forfurther discussion.A plane has two sides. We can arbitrarily define one side as positive corresponding toan outward-normal and the other negative. The direction of a closed curve in the plane isdetermined by the right-hand rule. In figure 7.1, n represents the normal vector of planeP2P3. The arrows on polygon PiP2P3 indicate the direction of the curve. They obey theright-hand rule. The vertices of a polygon are always numbered in the positive directionin this thesis.P2Figure 7.1 Normal vector and the direction of a closed curveFor piecewise-constant models, the 2D point-injection discretization can be simplifiedto determine if a point is internal or external to a polygon. On a plane, two points P1(zi, yi)and P2 (z2, y2) define a lineax + by = c,whereCHAPTER 7: 3D GRAVITY MODELING^ 109a = yz — Y1b=zi—x2,andc = aziThe positive direction is from P1 to P2. The distance from a point P(z, y) to the line is3 =c — a(z — xi) — b(y — yi) (az + bz)1/z^• (7.7)If s > 0, P is on the right-hand side ( we face to the positive direction). If a = 0, P is onthe line. If a < 0, P is on the left-hand side. If P(x, y) is an interior point of a polygon,the distances (defined by (7.7)) from P to all the edges of the polygon will be less thanor equal to zero or vice versa. This property can be used to realize the 2D point-injectiondiscretization.The 2D cell-average discretization can be simplified to determine the common area oftwo polygons. Two linear equations consist of a system1 aiz + bly = ci(7.8)azr + b2y = czIf aib2 — a2bi 0, system (7.8) has a unique solution, which is the intersection point of thetwo lines. The common area of two polygons can be determined by the following steps:1. If all the vertices of one polygon are inside the other, the common part is the smallerone.2. If not all the vertices of one polygon are internal points of the another, keep all theinner (if any) points. Then find the intersection points of each edge of one polygonwith all the edges of the another by solving (7.8). If an intersection point is internal toboth polygons, keep it.3. Remove the repeatedly-counted points. If the number of the remaining points is greaterthan 2, there is a common polygon. The vertices are the remaining points.CHAPTER 7: 3D GRAVITY MODELING^ 1104. Rearrange the remaining points in proper order so that they satisfy the right-hand rule.5. Suppose K is the number of the remaining points and (xi, yi) are the coordinates ofpoint i. The area of the common polygon isa = —2 E(zi — zi-Fi)(Yi Yi+i1i=1(7.9)When i = K, we should replace (K + 1) by 1 in (7.9).The 3D point-injection and cell-average discretizations are similar to the 2D problemsbut they are little more complicated. The 3D point-injection discretization can be simplifiedto determine if a point is internal to a given polyhedron. The normal vectors of thesurfaces of a polyhedron are defined outwards. Three non-collinear points Pi(xi, yi, zi),P2 (X21 Y2, Z2) and P3 (x3, y3, z3) determine a planeax + by -I- cz = d ,wherea = (y2 — Yi)(z3 — zi) — (y3 — Yi)(z2 — zi),b = (z2 — zi)(x3 — x1) — (x3 — zi)(x2 — xi),c = (x2 — xi)(Y3 — yi) — (xs xl)(Y2 Y1),andd = axi byi + czi.The normal vector isa^n = (n., ny , nz) = ( (a2 62 + 0)1/2 (a2 b2 c2)1/2 (a2 b2 c2)1/2)^(7.10)The distance from a point P(x,y,z) to the plane is^d — a(x — xi) — b(y — yi) — c(z — zi)s ^ (7.11)(a2 b2 c2 ) 1 /2If s > 0, P is on the positive side. If a = 0, P is on the plane. If s < 0, P is on the negativeside. If P(x, y, z) is an interior point of a polyhedron, all the distances (defined by (7.11))CHAPTER 7: 3D GRAVITY MODELING^ 111from P to the surfaces are less than or equal to zero or vice versa. This property can beused to realize the 3D point-injection discretization.The 3D cell-average discretization can be simplified to determine the common area oftwo polyhedrons. The three plane equations consist of a systemaix + bly + ciz = d1a2x + bzy + caz = d2^ (7.12)a3Z b3y C3Z = d3.Ifal b1 cla2a3b2b3 C30,the system (7.12) has a unique solution, which is the intersection point of the three planes.The intersecting volume of the two polyhedra can be determined as follows.1. If all the vertices of one polyhedron are the interior points of the other, the intersectionis the smaller one.2. If not all vertices of one polyhedron are the interior points of the other, keep all theinner (if any) points. Then we find the intersection points of each two planes of onepolyhedron and one plane of the other (by solving (7.12)). Keep the common interiorpoints.3. Remove the repeatedly-counted points. If the number of the remaining points is greaterthan 3, there is a common polyhedron. The remaining points are the vertices of thecommon polyhedron.4. Find the coplanar vertices and number them according to the right-hand rule.5. Suppose N. is number of surfaces, Ki is the number of vertices of the ith surface,and (x .(ii), zn is the jth vertex of the ith surface. The volume of the commonpolyhedron isIC;1V =^E(x`ii)-J;21)(y(i)+1,(i) )i=i j=1^•^Z^•31=1CHAPTER 7: 3D GRAVITY MODELING^ 112For 2D problems if all the internal boundaries are parallel to the coordinate axes, the2D discretization can be resolved into two 1D discretizations. For 3D problems if all theinterfaces are parallel to the coordinate planes, the 3D discretization can be resolved intothree 1D discretizations or one 1D and one 2D discretizations as in the examples to followin section 7.5.7.3 The Analytic Expression of the Potential Produced bya Homogeneous Polyhedral BodyThe gravitational potential of a homogeneous polyhedral body can be computed analyt-ically. Okabe (1979a) gave an analytic formulation for the first and second-order derivativesof the potential due to a homogeneous polyhedral body. Here we extend his method tothe potential. We consider the potential at the origin. Any other point (xo, yo, zo) can beeasily converted to the origin by the coordinate translationx = x — XoYi = y — zoz' = z — zo.We assume the polyhedral body has M facets and the ith facet is a polygon with Ni edges.Let r = (x2 + y2 + z2)1/2 and p be the density of the polyhedron. The potential at theorigin iswhere^u = G —pdV = 1 Gp V • IdV = —1 Gp — • ds = 1^,vr^2^r^2^s r^2(7.13)= r ! • ds .^ (7.14)s. rWe rotate system (x, y, z) to a new position (X, Y, Z) so that the Z-axis points along theCHAPTER 7: 3D GRAVITY MODELING^ 113Figure 7.2 Coordinate system rotationnormal vector of facet Si. The rotation is achieved by a two-step method (see figure 7.2)as follows:1. Rotate the x- and y-axis around the z-axis until the rotated x-direction is coincidentwith the projected direction of the normal vector of facet Si onto the x-y plane. Bythis counterclockwise rotation through the angle 0 (0 < 0 < 2r) in the x-y plane, wehave the transient system (x', Y, z). If the z-axis is perpendicular to the facet, 0 = 0.2. Rotate the z- and x'-axes around the Y-axis until the rotated z-direction is coincidentwith the direction of the outward normal. By this counterclockwise rotation throughthe angle 0 (0 < 4 < r), we obtain the target coordinate system (X, Y, Z).and8i =^Yi-F1 — — Xj) 2 (yi+i — 13)91/2. (7.18)CHAPTER 7: 3D GRAVITY MODELING^ 114The coordinate transformation can be written ascos 0 0 — sin 0 cos 9 sin 9 0 x[Xy0[1 0 — sin 9[cos 0 0 y[Z sin 0 0 cos 0 0 0 1cos 0 cos 0 cos 0 sin 9 — sin 0^x— sin 0^cos 0^0^y •sin 0 cos 0 sin 0 sin 9 cos 0(7.15)In practical computation all the vertices are given. The normal vector (n.,n y ,n.) canbe determined for each facet by (7.10). It is easy to find thatcos 0 = n, ,sin 0 = (1 — n!) 1 /2 ,1cos 0 =^n.(n2 + n2 )1 /2and0sin 0 = nv(4 + n:) 1 /2In the new system (X, Y, Z) the equation of plane Si is Z = Zi = constant. Theintegration in (7.14) becomesZidXdY Ii = (X2 + y2 + 4)1/2.9;I .X;+1^=Zi E^log [y + (x2 + y2 + zi2)1/2] dx^(7.16)j=i XiSupposing that the jth edge of the polygon Si has vertices j and j 1 (for the last edge,vertex j 1 corresponds to vertex 1). Let^C'== X j+i - Xi(7.17)[(Xj+i — X1)2 (1'j+1 — Y)91/2if In.1= 1,if In.1 < 1,if inz I = 1 ,if In.l < 1.—2Zi tan-1 Zici(xi,n)xci (Xi+i,n+i)(1 + aj) [y (K2 y2 + 4)1/21(7.19)CHAPTER 7: 3D GRAVITY MODELING^ 115Following Okabe (1979a), we haveIi =Zi Ef(Xsi — Yci) log [Xci Ysi + (X2 + y2 + 4)1/21j=1All the (X, Y, Z)' a can be found by (7.15).7.4 Approximate Boundary ConditionFor most practical gravity problems, the boundary conditions are not known. In thosecases we compute the potential values at the boundary approximately by multipole expan-sions. Let r = (x, y, z) be the exterior observation point, = , y' , z') the source point,and V' the source volume. Thento,^=^_P r, I dV(,.., G It 1$ pi 0^—1 i i di. 02 1r + 2.8xir + 6. 3 8xi8xi ri=1^s=1 j=1where x 1 x, x2 = y, x3 = z. On the right-hand side of (7.20) the quantity(7.20)m = p dVrepresents the total mass. The first term represents the potential of the total mass locatedat the origin. The quantitypi =^xipdV^(i = 1,2,3)JVrepresents the mass dipole moment vector. The second term represent the potentialproduced by the mass dipole. The quantitydii =^(3x1ixi — ri2) p dV^(i, j = 1,2, 3)CHAPTER 7: 3D GRAVITY MODELING^ 116represents the mass quadrupole moment tensor. The third term represents the potentialproduced by the quadrupole moment. We shift the origin to the center of mass, so thatpi = 0. The quadratic terms can be computed by either analytic or numerical integrationmethods. The expression for u as given by (7.20) is then used to compute the approximateboundary conditions on the exterior bounding surface r.'7.5 Numerical ExamplesNow we present three numerical examples. In example 7.1, the source body is a uniformcube. In example 7.2, the source body is a uniform tetrahedron. Example 7.3 simulates ageological normal fault. The source body consist of four blocks and with a spatial densityvariation. In each example we compute the potential by the analytical formulation pre-sented in section 7.3 as a standard to compare with the numerical solutions. We solveeach problem by both the r and the collocation methods on a single-domain and with boththe accurate and the approximate boundary conditions. We also present the multi-domainsolution for example 7.1.Example 7.1 For this example in equation (7.1) we take the source body to be a cubewith a unit density (figure 7.3), that is,P(x,Y,z) = {1 for — 0.5 < x, y, z < 0.5 ,0 otherwise.We first solve this problem by the single-domain method discussed in chapter 3. Becauseall the interfaces are parallel to the coordinate planes, both point-injection and cell-averagediscretizations can be expressed as a production of three 1D discretizations. For the point-injection discretization we define1 if —0.5 < xi < 0.5 ,ti =0 otherwise.di — di+iwhere di is defined by (7.5). The discrete values of p for this discretization areti — di+i if di+i < 0.5 < didi di + 0.5 if di+i < —0.5 < di .CHAPTER 7: 3D GRAVITY MODELING^ 117Figure 7.3 The model of example 7.1. The big cube is the computationaldomain. The small cube (thick lines) is the source body.The discrete values of p for point-injection are()^{ 4piik = tiqqc •For the cell-average discretization we define0^if di+i > 0.5 or di < —0.5 ,1^if di+i > —0.5 and di < 0.5 ,(C)^4 4 sPo = tlit,j6kCHAPTER 7: 3D GRAVITY MODELING^ 118The accurate solutions and boundary conditions are computed by the methods discussed insection 7.3. The approximate boundary conditions are computed by the multipole expan-sion method discussed in section 7.4. We solve for potential u by both the r and collocationmethods. The relative error norms are shown in tables 7.1-7.4, in which N is the degreeof Chebyshev polynomials and Pt the number of total grid points.Table 7.1The Lc,„ relative error norms of example 7.1 (r method)N PtAccurate BCscell-average^point-injectionApproximate BCscell-average^point-injection8 729 2.269 x 10 -02 1.812 x 10 -01 2.170 x 10-02 1.825 x 10 -0116 4,913 6.046 x 10 -03 5.085 x 10 -02 5.213 x 10 -03 4.967 x 10 -0232 35,937 1.284 x 10 -03 2.852 x 10 -02 5.213 x 10 -43 2.970 x 10 -0264 274,625 3.180 x 10-" 1.305 x 10 -02 5.213 x 10 -03 1.179 x 10 -02Table 7.2The L2 relative error norms of example 7.1 (r method)N PtAccurate BCscell-average^point-injectionApproximate BCscell-average^point-injection8 729 4.114 x 10 -03 2.912 x 10 -02 4.203 x 10 -03 2.950 x 10 -0216 4,913 6.680 x 10-" 1.180 x 10 -02 1.153 x 10 -03 1.142 x 10 -0232 35, 937 1.762 x 10 -04 7.226 x 10-43 9.288 x 10 -04 7.684 x 10 -0364 274, 625 3.711 x 10-" 3.424 x 10 -03 9.162 x 10 -04 3.103 x 10-133CHAPTER 7: 3D GRAVITY MODELING^ 119Table 7.3The Lo. relative error norms of example 7.1 (collocation method)N PtAccurate BCscell-average^point-injectionApproximate BCscell-average^point-injection8 729 2.032 x 10 -02 1.383 x 10 -01 1.918 x 10 -02 1.394 x 10 -0116 4, 913 4.193 x 10-°3 4.722 x 10 -02 5.213 x 10-03 4.603 x 10 -0232 35, 937 9.733 x 10 -04 2.809 x 10 -02 5.213 x 10 -1)3 2.935 x 10-0264 274, 625 2.929 x 10-" 1.300 x 10 -02 5.213 x 10 -°3 1.174 x 10-02Table 7.4The L2 relative error norms of example 7.1 (collocation method)N PtAccurate BCscell-average^point-injectionApproximate BCscell-average^point-injection8 729 2.838 x 10 -03 2.864 x 10 -02 2.882 x 10 -03 2.901 x 10-0216 4,913 4.979 x 10-" 1.177 x 10 -02 1.077 x 10-°3 1.139 x 10 -0232 35, 937 1.661 x 10 -04 7.224 x 10-03 9.259 x 10-" 7.682 x 10 -0364 274, 625 3.627 x 10 -1)5 3.424 x 10 -513 9.162 x 10-" 3.103 x 10 -113Now we solve the problem by the multi-domain Chebyshev spectral method describedin chapter 4. The whole domain is divided into 3 x 3 x 3 subdomains (figure 7.4). Thedividing planes are x = y = z = ±0.5. Each subdomain is divide into N3 cells. Theaccurate solution and boundary conditions are computed by (7.13)-(7.19). Table 7.5 showsthe relative error norms obtained by the multi-domain Chebyshev spectral method withaccurate boundary conditions. In table 7.5, B., By , and B. represent the number ofsubdomains along the x-, y-, and z-directions respectively, Bt the number of subdomains,L2 and L, the corresponding error norms, and the CPU time is obtained on a SPARCstation 2 computer.CHAPTER 7: 3D GRAVITY MODELING^ 120. • . - • ...... • 'fr^.^;. ^::.--'":....,.......i. ....... ,,,„i...L.^-,..^:^: ..00".^I::: : .. .* : ...^ .... ....Im...-^--, --..0 ^_ .,r^....,- IIII.••--....•---........—I^...-— r.o...^.••••• .....• . ^• •..^.^—...,. .^•..••^. .,-• .^..-^...-- .,•••• .-^....- ^--Figure 7.4 Subdomains of example 7.1Table 7.5The relative error norms of example 7.1(MDCS method, accurate boundary conditions )N B. By Bz Bt Pt L2 Lo. CPU(sec)2 3 3 3 27 343 8.014 x 10 -43 6.557 x 10 -02 0.253 3 3 3 27 1,000 4.378 x 10—°3 1.898 x 10 -02 0.254 3 3 3 27 2,197 3.168 x 10 -04 1.683 x 10 -133 1.435 3 3 3 27 4,096 1.826 x 10 -134 9.745 x 10—" 4.506 3 3 3 27 6,859 4.467 x 10 -45 2.864 x 10 —" 12.727 3 3 3 27 10,648 2.778 x 10 -1)5 2.019 x 10—" 28.968 3 3 3 27 15,625 1.094 x 10 -135 9.242 x 10 -05 62.56CHAPTER 7: 3D GRAVITY MODELING 121Example 7.2 For this example in equation (7.1), we take the source body to be a tetra-hedron defined by four points Pi(-0.7, 0.7, —0.7), P2(0.7, —0.7, —0.7), P3(0.7, 0.7, 0.7) andP4(-0.7, —0.7, 0.7) and with a unit density (figure 7.5). We discretize the density functionby both point-injection and cell-average schemes and compute the potential with both the rand the collocation methods with both the accurate and approximate boundary conditions.The relative error norms are shown in tables 7.6-7.9.Figure 7.5 The model of example 7.2. The thick lines indicate the sourcebody. The thin lines indicate the computational domain.CHAPTER 7: 3D GRAVITY MODELING^ 122Table '7.6The Lc° relative error norms of example 7.2 (r method)Accurate BCs^Approximate BCsN Pt cell-average point-injection cell-average point-injection8 729 2.999 x 10 -02 1.341 x 10 -01 2.996 x 10 -02 1.324 x 10-0116 4,913 7.285 x 10 -03 4.377 x 10-02 9.865 x 10-°3 4.307 x 10 -0232 35,937 2.273 x 10 -03 1.258 x 10 -02 9.865 x 10-°3 1.451 x 10 -02Table 7.7The L2 relative error norms of example 7.2 (r method)Accurate BCs^Approximate BCsN Pt cell-average point-injection cell-average point-injection8 729 3.355 x 10-°3 2.601 x 10 -02 5.687 x 10 -03 2.601 x 10 -0216 4,913 9.713 x 10 -04 1.093 x 10 -02 4.640 x 10 -°3 1.152 x 10 -0232 35,937 2.747 x 10-" 1.600 x 10-°3 4.514 x 10-°3 5.124 x 10-03Table 7.8The Lor, relative error norms of example 7.2 (collocation method)Accurate BCs^Approximate BCsN Pt cell-average point-injection cell-average point-injection8 729 2.345 x 10 -02 1.301 x 10 -01 2.190 x 10-02 1.284 x 10 -0116 4,913 6.868 x 10 -°3 4.295 x 10 -02 9.865 x 10-03 4.225 x 10 -0232 35,937 2.242 x 10-°3 1.266 x 10 -02 9.865 x 10 -03 1.458 x 10 -02CHAPTER 7: 3D GRAVITY MODELING^123Table 7.9The L2 relative error norms of example 7.2 (collocation method)Accurate BCs^Approximate BCsN Pt cell-average point-injection cell-average point-injection8 729 2.904 x 10 -03 2.600 x 10 -02 5.452 x 10 -133 2.600 x 10-0216 4,913 9.539 x 10-04 x 10-02 4.637 x 10 -03 1.152 x 10 -0232 35,937 2.738 x 10 -04 1.599 x 10-°3 4.514 x 10-°3 5.124 x 10-03Example 7.3 This model consists of four blocks simulating a geological normal fault(figure 7.6). The coordinates of the vertices are shown in table 7.10. This model has aspatial density variation. The densities of blocks I and II are p = 2 -I- cos iry, and thoseof blocks III and IV are p = 1. In order to employ the techniques discussed above, wecompute the center of mass (zo,yo, zo) first and then shift the origin to (xo, yo,For this model there are no known analytic formulae to compute the potential field. Wecompute the "accurate" potential as the follows. For the blocks III and IV, we computethe potential with the analytic formulation in section 7.3. For the blocks I and II, wedivide each blocks into many thin vertical slices along y =constant planes. In each slicewe approximate the density by a constant. Then we compute the potential by the analyticformulation presented in section 7.3. The maximum relative error of the "accurate" solutionis within 10 -5 .We solve this problem by both the r and collocation methods. Both point-injection andcell-average discretizations can be resolved into one 1D discretization (along the x-direction)and one 2D discretization (on y-z plane). All the other computations are the same as before.The relative error norms are shown in tables 7.11-7.14. Figure 7.7 shows the relative errormaps in the z = 0 plane (after coordinate shifting) for N = 32. Comparing figure 7.7(a)with figure 7.7(b), and figure 7.7(c) with figure 7.7(d), we see that for both accurate andapproximate boundary conditions, the relative error norms obtained by cell-averageX YFigure 7.6^The model of example 7.3. The stars indicate the locationCHAPTER 7: 3D GRAVITY MODELING^ 124of the z = 0 plane (after coordinate shifting) which is shown in figure 7.7.Table 7.10The coordinates of the vertices of example 7.3No. ( 0, y, z) No. ( z, y, z) No. ( z, y, z)1 ( 0.5,-0.6, 0.4) 9 ( 0.5,-0.6,-0.2) 17 ( 0.5, 0.1,-0.1)2 ( 0.5,-0.4, 0.4) 10 ( 0.5, 0.2,-0.2) 18 ( 0.5, 0.6,-0.1)3 (-0.5,-0.4, 0.4) 11 (-0.5, 0.2,-0.2) 19 (-0.5, 0.6,-0.1)4 (-0.5,-0.6, 0.4) 12 (-0.5,-0.6,-0.2) 20 (-0.5, 0.1,-0.1)5 ( 0.5,-0.6, 0.1) 13 ( 0.5,-0.2, 0.2) 21 ( 0.5, 0.4,-0.4)6 ( 0.5,-0.1, 0.1) 14 ( 0.5, 0.6, 0.2) 22 ( 0.5, 0.6,-0.4)7 (-0.5,-0.1, 0.1) 15 (-0.5, 0.6, 0.2) 23 (-0.5, 0.6,-0.4)8 (-0.5,-0.6, 0.1) 16 (-0.5,-0.2, 0.2) 24 (-0.5, 0.4,-0.4)1^1^1^1L 1^ 0.5-0.5 10x............^1^..............................^2^--..^•0•.^• •• • •›, 0• • • ..... ......^- • -„ ...................... - 1 -0.••........^3 .......^2^............... , .........^1^.......^........... ^......^..... . .............^............. . ... ,.. .*^-^......."-1-2..^..."„••"'........ ..^ ....1 ............. ^3 \^1..................4,1lf?0), 0IL0CHAPTER 7: 3D GRAVITY MODELING^ 125(a)(c)(b)x(d)L 1 -0.5^0xFigure 7.7I '"0.5^1^L1^-0.5^0^0.5^1xThe relative error of example 7.3 on z = 0 plane (N =tr?00032). The contour levels are 0.5% of the number on the line. The solidrectangles are the cross sections of the source body at the z = 0 plane(after coordinate shifting) which is indicated by the stars in figure 7.5. Thepanels are obtained by: (a) accurate BCs and cell-average discretization.(b) accurate BCs and point-injection discretization. (c) approximate BCsand cell-average discretization. (d) approximate BCs and point-injectiondiscretization.CHAPTER 7: 3D GRAVITY MODELING^ 126Table 7.11The Lo, relative error norms of example 7.3 (r method)Accurate BCs^Approximate BCsN Pt cell-average point-injection cell-average point-injection8 729 4.580 x 10-02 1.906 x 10 -01 4.785 x 10 -02 1.927 x 10-0116 4,913 2.909 x 10 -02 1.399 x 10 -01 2.846 x 10 -02 1.387 x 10 -0132 35,937 2.030 x 10 -02 2.625 x 10 -02 2.376 x 10-02 2.848 x 10 -02Table 7.12The L2 relative error norms of example 7.3 (r method)N PtAccurate BCscell-average^point-injectionApproximate BCscell-average^point-injection8 729 5.388 x 10-°3 2.411 x 10 -02 6.076 x 10-03 2.500 x 10 -0216 4,913 2.002 x 10 -03 1.533 x 10 -02 3.456 x 10-°3 1.493 x 10 -0232 35,937 1.789 x 10-03 2.736 x 10 -°3 3.370 x 10-03 4.063 x 10-°3Table 7.13The Loc, relative error norms of example 7.3 (collocation method)Accurate BCs^Approximate BCsN^Pt^cell-average point-injection^cell-average^point-injection8 729 2.393 x 10 -02 1.921 x 10 -01 2.515 x 10 -02 1.941 x 10 -0116 4,913 2.675 x 10 -02 1.398 x 10-01 2.612 x 10-02 1.385 x 10 -0132 35,937 2.003 x 10-02 2.607 x 10-02 2.348 x 10 -02 2.783 x 10 -02CHAPTER 7: 3D GRAVITY MODELING^127Table 7.14The L2 relative error norms of example 7.3 (collocation method)Accurate BCs^Approximate BCsN Pt cell-average point-injection cell-average point-injection8 729 3.272 x 10—°3 2.409 x 10 -02 4.285 x 10—°3 2.498 x 10 -0216 4,913 1.942 x 10—°3 1.532 x 10 -02 3.424 x 10—°3 1.492 x 10'232 35, 937 1.787 x 10 -03 2.735 x 10 -03 3.369 x 10—°3 4.062 x 10 -03discretization are less than those obtained by point-injection discretization. Comparingfigure 7.7(a) with figure 7.7(c), and figure 7.7(b) with figure 7.7(d), we see that the relativeerror norms for the case of the approximate boundary conditions are only slightly largerthan those for accurate boundary conditions. In figure 7.7(c) and figure 7.7(d), due to theapproximate boundary conditions, the contour lines intersect the boundaries. The relativeerror caused by the approximate boundary conditions is about 1%.7.6 DiscussionIn this chapter we solve the 3D forward gravitational problem by Chebyshev spectralmethods. We present a method to realize the point-injection and cell-average discretizationfor 2D and 3D problems. We also extend Okabe's (1979a) analytic expression of thegravity field of a homogeneous polyhedron to the potential field. We compute the boundaryconditions approximately by multipole expansions.7.6.1 The r and Collocation MethodsComparing the results of the r method with the corresponding results of the collocationmethod shown in tables 7.1-7.4, 7.6-7.9, and 7.11-7.14, we find that the accuracy of thecollocation method is slightly better than the r method. The CPU time is 4.06 secondsfor the r method and 2.53 seconds for the collocation method on a 33 x 33 x 33 grid. TheCHAPTER 7: 3D GRAVITY MODELING 128speed of the collocation method is also faster than the r method. More comparisons aboutthe r and the collocation methods can be found in subsection 3.3.2. The results of the rmethod have been published in our paper (1991) and we include them in this chapter.7.6.2 Cell-Average and Point-Injection DiscretizationIn practical gravitational modeling, we are interested in models with piecewise-constantor piecewise-smooth density distribution, such as the examples presented in this chapter.There are usually discontinuities in the density function. In the single-domain Chebyshevspectral method, these discontinuities cause many short-wave-length components. Fromthe analysis in chapter 2, in this case there are both truncation and aliasing errors in thefinite Chebyshev expansion approximation. For a given function and a given degree Nof the Chebyshev expansion, the truncation error is fixed. We can use the cell-averagediscretization scheme to reduce the aliasing error. From tables 7.1-7.4, 7.6-7.9, 7.11-7.14,and figure 7.7 we see that for a given model with the same grid and boundary conditionsthe relative error norms obtained by cell-average discretization is 4-80 times smaller thanthat obtained by point-injection discretization. However, overall the error is still roughlyproportional to (1/N2 ).7.6.3 Single-Domain and Multi-Domain SolutionsIn the multi-domain Chebyshev spectral method, the functions are smooth in eachsubdomain and we approximate the function in each subdomain by a separate Chebyshevseries. Therefore there is no Gibbs phenomenon. The point-injection discretization isaccurate enough. It is unnecessary to use the cell-average discretization scheme. The multi-domain Chebyshev spectral method is more accurate than the single-domain Chebyshevspectral method. From table 7.5, the error norms of the multi-domain Chebyshev spectralsolutions go to zero exponentially. On a grid with 15,625 points, the error norms are•-•-, 10-5 . Comparing table 7.5 with tables 7.3-7.4 the relative error norms of the multi-domain solution on a grid with total 10,648 points are better than those of the single-domainCHAPTER 7: 3D GRAVITY MODELING 129collocation method on a grid with total 274,625 points. The CPU time is 28.96 secondsfor the multi-domain Chebyshev spectral method and 37.32 seconds for the single-domaincollocation method. To obtain the same accuracy, the multi-domain Chebyshev spectralmethod is more efficient than the single-domain Chebyshev spectral method.7.6.4 Approximate and Accurate Boundary ConditionsThe error caused by the approximate boundary conditions is about 1% for the threemodels computed. In tables 7.1 and 7.3, the L„, relative error norms for both the r andcollocation methods and for both N = 16 and N = 32 are 5.213 x 10-3 . In tables 7.6and 7.8, the L co relative error norms for both the r and collocation methods and for bothN = 16 and N = 32 are 9.865 x 10-3 . These values are the maximum errors caused bythe approximate boundary conditions. In numerical experiments with the approximateboundary conditions, we have found that if the ratio of the maximum dimension of thedensity distribution to the dimension of the computational volume is less than 70%, thenthe maximum relative error of the potential is about 1%. As this dimension ratio decreases,because the mass is concentrated in a relatively smaller volume, the error of the multipoleexpansion approximation at the exterior boundaries also becomes smaller.CHAPTER 83-D DC RESISTIVITY MODELINGIn DC resistivity modeling, analytic formulas exist for some simple configurations, suchas a layered model, and a conducting sphere buried in a uniform half space (Grant and West,1965, chapter 14). For more general 2-D and 3-D models numerical methods are necessary.The numerical methods reported in the literature are the finite-difference method, finite-element method, and the integral equation method. With the finite-difference method,Mufti (1976) computed 2-D models. Dey and Morrison (1979a and 1979b) computed 2-D and 3-D models. James (1985) combined the Polozhii decomposition with the finite-difference method to solve the 3-D DC resistivity problem. McGillivray (1992) computed 2-D pole-pole models. With the finite-element methods, Bibby (1978) computed the apparentresistivity of a body with vertical axis of symmetry embedded in a uniform half space. Fox etal. (1980) computed 2-D topographic effects for DC resistivity and IP surveys. Pridmore etal. (1981) computed 3-D DC resistivity models. Holcombe and Jiracek (1984) computed the3-D terrain correction by replacing the near-surface layers of regular hexahedron elementswith contiguous distorted isoparametric elements. With integral equation formulations,Alfano (1959) computed 2-D models with horizontal and vertical interfaces. Dieter etal. (1969) computed IP and DC resistivity type curves for 3-D bodies. Pratt (1972)computed 3-D DC resistivity models. Lee (1975) computed 2-D and 3-D models for IPand DC resistivity. Daniels (1977) computed 3-D IP and DC resistivity models for buriedelectrodes. Spiegel et al. (1980) calculated 2-D and 3-D geological anomalies and terraineffects. Okabe (1979b) suggested a method to remove the singularity arising in the integralequation formulation. Xu et al. (1984) computed 3-D models with piecewise-constant130CHAPTER 8: 3-D DC RESISTIVITY MODELING 131conductivity and a horizontal ground surface. Poirmeur and Vasseur (1988) computed 3-Dborehole models for electrical method. Xu et al. (1988) computed 3-D terrain effects forDC resistivity problem. Eskola et al. (1989) computed DC resistivity models containinga thin conductor. All the finite-difference , finite-element, and integral equation methodsuse polynomial-based interpolation to approximate the unknown function. If we want toobtain high accuracy, a large number of grid points has to be used. This is costly, both inboth memory and time consumption.In this chapter we solve 3-D forward DC resistivity problems by the multi-domainChebyshev spectral method discussed in chapter 6. This method can achieve much higheraccuracy with fewer grid points than the finite-difference method. In section 8.1 we outlinethe boundary value problem. In section 8.2 we discuss how to remove the singularitycaused by the source current. In section 8.3 we present the formulation of the multi-domainChebyshev spectral method for 3-D DC resistivity modeling. In section 8.4 we present amodified finite-difference method. In section 8.5 we present three numerical examples,which include the comparison of the multi-domain Chebyshev spectral method with thefinite-difference method. In section 8.6 we discuss the results and present the conclusions.8.1 The Boundary Value ProblemIn DC resistivity survey, we inject a source current I into the ground and measurethe potential difference Au between electrodes near the source current (figure 8.1). Thenwe try to find a geoelectrical structure that produces a potential distribution close to theobserved data. Suppose the conductivity structure is known. The fundamental task forforward DC resistivity modeling is to solve the boundary value problemV • (aVu) = —I 8(x — z a , y — y,, z — z,)On—On _ 0 on the ground surfaceu = b(x, y, z)^on other boundaries.in IICHAPTER 8: 3-D DC RESISTIVITY MODELING^ 132Figure 8.1 DC resistivity method configurationwhere the conductivity or(x, y, x) and the boundary condition b(x, y, x) are known functions,u the potential, (x, y, x) the observation point and (x,, y,, x a ) the source point. The homo-geneous Neumann boundary condition is given on the air-earth surface because no currentcrosses it. The external boundary conditions (except for the air-earth surface) are givenby the Dirichlet type. On the interface between two media with different conductivities er1and 02, the matching boundary conditions areui = ti2^ (8.2)andOui Ou2cri^= —472— •Oni^On2 (8.3)CHAPTER 8: 3-D DC RESISTIVITY MODELING 133The negative sign in (8.3) is due to the fact that the normal vectors on two sides are inopposite directions. If the same normal vector n is used for both sides of the interface,equation (8.3) can be rewritten as&al^0214On — 472 8n = u (8.4)8.2 Removal of the Source SingularityAt the source point the potential becomes infinite. This causes large errors in thefinite-difference method (Lowry et al. 1989, McGillivray 1992). The same problem existsfor Chebyshev spectral methods. Let us consider a simple case. Suppose a current I isinjected into the surface of a uniform half-space with conductivity a. The potential isIu= 27rar'wherer= Rz z8)2 +(3l — Ya)2 + (z z8)2 1 1/2In (8.5) and (8.6) if I = 22-a, x, = 1.1, and y, = z, = 0, the potential along x-axis isu = 1/(1.1 — x). From figure 2.1 the Chebyshev coefficients of 1/(1.1 — x) on the interval[-1, 1] converge much slower than those of cos Tx and ez. This is due to the singularityof 1/(1.1 — x) at x = 1.1. Even though the singular point is outside of the computationaldomain, it still causes large errors.Lowry et al. (1989) suggested the removal of the singularity by splitting the potentialinto the primary potential 4 caused by the source current and the secondary potential &caused by the inhomogeneities. Thus(8.5)(8.6)u=4+0^ (8.7)CHAPTER 8: 3-D DC RESISTIVITY MODELING^ 134The idea is correct and classical (see, for example, Coggon 1971, Charbeneau and Street1979). However, Lowry et al. incorrectly gave the primary potential asI= 2war(8.8)where o was defined as the average conductivity of the whole model. With the aboveformula the singularity cannot be removed completely in the neighbourhood of the sourcepoint where the singularity effects are the greatest. The correct formula should be0 = 2iror (8.9)if (z., yg , z3) is on the ground surface, or0 = 47rarI(8.10)if (z5 , ya , z8 ) is underground, where a is the conductivity of the media at the source point.In this chapter we always assume that the source point is on the ground surface and take(8.9) as the primary potential field.Substituting (8.7) into (8.1) yieldsV •(aVik) = — [M(x — zg, Y — Ys, z — z8) + V • (aV0)]On = 0^on ground surface= known on other boundaries.in SZBecause u and are known functions, the right-hand side of (8.11a) is known too. Substi-tuting (8.7) into (8.2) and (8.4) yields the interface conditions for the secondary potentials—^= 0 (8.12)001^42^00Cr2 Tr/ + (Cri (72)^= °^(8.13)CHAPTER 8: 3-D DC RESISTIVITY MODELING^ 1358.3 The Multi-Domain Chebyshev Spectral MethodThe boundary value problem (8.11) can be solved by the multi-domain Chebyshevspectral method. We expand the right-hand side terms of (8.11a) to obtainV • (crVii)) = —Va. • VO^in O.^(8.14)Then we solve the boundary value problem consisting of (8.14), (8.11b), and (8.11c). Theprocedure is very similar to that for (6.15). Here the matching boundary conditions (8.12)and (8.13) have to be used at the interface points. Let2ri =_ -^^,xi — xi+128m = _ -^^1Ym — Ym+12to =_^_^.zn — Zn4.1The 3-D Chebyshev grid is given by/ = 2/ + 2/4.1 +r cos ir2^ri^N 'Vnt + girs+1^1^ixYr = 2^+ COS —^am^N '2n + 4-1-1^1^kr4 =^ + — cos —2^to^N .Let= [(xi x42 + (yr ys)2^zet^1/2^— — Z.)2b 7 = (4 ,yr ,^,_Iranujik(8.15)(8.16)(8.17)CHAPTER 8: 3-D DC RESISTIVITY MODELINGQs'^ (z) =j3k ax 1(4yr,.1)(m) 00 I^i(16tpktt) = _ A Zi' Z s1Ox - 270" 4./.71n) 3 '(4rq)^uk )hnn (0 80ijk^= —^—0Y (z! ,Yr ,zr)^27ro• (../ nn \ 3 'vsjk )and136(8.18)(8.19)By the same subdomain partition and derivation as that presented in section 6.3, from(3.58), (3.60), (8.11), (8.13) and (8.14) we obtain the system of equations:N 4^E 4 .4.7 d)^Er? E 4 47 + 4p)orp=0 p=0(.1 NE d^Nvopbx + am (0.1r.nn) (19 E d(1) hn^lmn)(z)sjk^jp ?kip: + to (0% •k^E d( 1 ),Altnns3 kp Y'i.7PP=0 p=0^ p=0(1 = 1, • • • , Bx ; m = 1, • • • , By ; n = 1, • • • , B.; 1, j, k = 1, • • • , N — 1),^(8.21a)N^ N(l-l)mn V-• 4(1)./.(l-l)mn +^ hnn > 40.)„hintnri_ierNjk^Zd'Op '(N-p)jk -r- ricrOjk Z., 'Op irpjkp=0 p=0p=0(m)^Inn (z)^hnn (Y)^1^(11)^inns (z)^lni (z)=^[ (Tr) (oi k)(cri ik^ oat)^(criik^(oipin)(z) 00^I^— z, (4:7r) = =(z!,yr,q)^27r0 (ri nn)ijk(8.20)1 ) (z) — 0(ri_iff(pli.711. )mn — ricrZn) (47:^—(1 = 2, • • • , B.; m = 1, • • • , By ; n = 1, • • • , Bz ; j,k = 1, • • • ,N —1),^(8.21b)l(tn -1)n3m CriNkjo)./.1(rn_i)n^a^Op 1)n E d(i) 1 rill“Op i(N -P)k smalp=0^ p=0(ans — 1 ark-1)n Sm (rtin) (Olnkn) (V) = 0^(1 = 1, • • •^m = 2, • • , By ; n =^1, • • • , Bz ; i,k = 1, • • • ,N — 1 ),N Nn-1 sjN E ioip)ov;:;; + tnair E 41p)tkr^P=0^ p=0(8.21c)(tn -1 ers1•71-1) tn air ) (41)(z) = 0CHAPTER 8: 3-D DC RESISTIVITY MODELING^ 137^(1 = 1, • • • , B.; m = 1, • • • , By; n = 2, • • • , B.; i , j = 1, • • • , N — 1),^(8.21d)E d(oip).47;),. 0 (1 = . . . , B.; m = 1, • • • , By ; 1, j = 1, • • • , N — 1),^(8.21e)p=0N^ N(1-1)mn^JO) (1-1)mn^lmnr1—lallok d(Lawop (N—p)ok + 00 E 1)01mnOp pokp=0 p=0(7.1-1a(Pli011s )mn ria0017) (010NI) (2) = 0(1 = 1, • • • , B.; m = 2, • • • , By; n = 1, • • • , B.; k = 1, • • • , N — 1),^(8.211)Nri_iogr—joi)mn E^+ notion E iolp)oplmonp=0 p=0— TjObio^wojk(1-1)mn^lmn) (.1Imn) (2)^(1 = 1, • • • , B.; m = 1, • , By ; n = 2, • • • , /4; j = 1, • • • , N — 1),^(8.21g)N^ N1(m-1)n^d(1)./.1.(m-1)n,8171- 1 47iN0^Z..d Op Vs(N—p)0 "W.100^1/4Clp Yspop=0 p=0(8m-1C15"1-1)m^Linn (Almn (Y) 0sNO^8mai00 ) ViOk )^u^(1 = 1, • • • , B.; m= 2, • • • , By ; n = 2, • • • ,Bz ; i = 1, • • • , N — 1),^(8.21h)N^ N8m-10.01(70-1)n^14,001((raNzlpro arnalgon E 10p1)0trAnp=0 p=0(sni_10.10Z-1)st smaimn) (ohnn)(Y) . 0^(1 = 2, • • • , Bz ; m = 2, • • • , By; n = 2, • • • ,B,),^(8.21i)where dli ) and dIP are defined by (3.58) and (3.60) respectively. Equations (8.21a) are forthe interior points, (8.21b)—(8.21d) for the interface points, (8.21e) for the ground surfacepoints, (8.21f)—(8.21h) for the edge points, and (8.21i) for the corner points. The system(8.21) can be solved by the same procedure as that for system (6.17). When computing thejump of the flux, o014/On, the contribution due to the primary field should be computedby the analytic formula (8.18)—(8.20) to reduce errors.xy6zYCHAPTER 8: 3-D DC RESISTIVITY MODELING^ 1388.4 The Finite-Difference MethodFor comparison, we present a modified integrated finite-difference method for solvingthe boundary value problem (8.11) in this section. Although the basic idea is similar to thatpresented by Dey and Morrison (1979b) and Lowry et al. (1989), there are two significantdifferences. One is the singularity removal, which has been discussed in section 8.2. Theother is the computation of the right-hand side terms of the finite-difference equations.Figure 8.2 A 3-D finite-difference gridCHAPTER 8: 3-D DC RESISTIVITY MODELING^ 139McGillivray (1992) suggested the analytic computation of the right-hand side terms of thefinite-difference equations for the 2-D problem. Here, we present an analytic method tocompute the right-hand side terms for the 3-D problem. In the finite-difference method,we divide the whole domain into N.NyN, rectangular cells. The dividing planes arex = Xi^(i = 1 , • • • ZN.1-1))Y^Yi^11 • • • ,YN,+1),andz = zis^(k = 1, • • • ,zivz+i),where xi > zi-Fi, > Yi+i, and zk > zk+i (see figure 8.2, in which N. = Ny = NN = 5).There are total of (N. + 1)(Ny + 1)(N. + 1) grid points. Following Dey and Morrison(1979b), the conductivity in each cell is approximated by a constant which is chosen asthe conductivity at the cell center. At each interior grid point (xi, yj , zk), we define arectangular parallelepiped Vijk = {(x, y, z) I (xi + xi+i)/2 (xi_l xi)/2, (y, +Y5+1)12^y^+ WW2 , (zk + zk+i)/2^z 5_ (zk_ 1 zk)/2 } (see figure 8.3).Integrating both sides of (8.11a) we obtainV • (crVO) dv =^[I8(x — za,Y Ya, z — z.) + V • (crV^dv.^(8.22)Using Gauss's theorem, we obtain6. (crVO) • ds = — r [1.8(z — zoY — Y8, Z zg) + V • (crV4)] dv.^(8.23)roThe volumetric integration in (8.23) is performed in Vijk and the surface integration isperformed on its closure. LetAzi = xi — zi+i,=^Y1+1,andAzk = zk — zk+1.CHAPTER 8: 3-D DC RESISTIVITY MODELING^ 140ijk....••^ .•^ij.....•i(j+ 1)k^-4,^....•••.•‘.•••••.•(i+1)jkii(k+1)Figure 8.3 The rectangular parallelepiped Vijk for an interior pointOn the surfaces x = (xj_i xi)/2, the normal derivative of the secondary potential bi is.•'..., ....•.•.•••••^ .••"...,^..^if.• ,..•approximated by8x^bai_i^•The integration on this surface is given approximately by1 (00- 1)/k '00)^[6(i-1)(/-1)(k-1)41-1Azk-1 '(i- 1)/(k- 1 )A Acr(i-i)jkAYibizig] •The integrations on all other surfaces can be computed in a similar way. Adding all theseresults together we obtain an equation+ 2^3^4 _Lcijkyq^-ri-1)jk Cijk ,V(i+1.)jk CijkiVi(j-1)k^Ciippi(j+i)k^5CHAPTER 8: 3-D DC RESISTIVITY MODELING^ 141+c?jk ibij(k+1) cijklkips = fijk•^ (8.24)The coefficients are:1 1^= 4AXi-1 hi--1)(j-1)(k-1)AW-lAzk-1 0'(i—iV(k—i)AYjAzk-1cro_ ivkAyjAzki ,2^rCi^1ik -4Azi lui(j-1)(k-1)■6%-lAzk-1 iTij(k_i)AyjAZk_icri(j-1)kAYj-10zk 47iik4j2Zki 74Ac^1^r• •k =^" yi-1(70-1)(j-1)kAzi-10ZIg orri(i_i)kAziAzki1^r—^" — 4AyioijkLxiAzki ,1 [(1(i-1)(i-1)(k-1)Azi—ibqj-1 cri(j-1)(k—i)Azi4j—icrii(k_i)AziAyil1 ^reVk =4Azk +criikAxiAYjiand7^ 6= - (Clik + 4k Ci3ps ci4ik + 5 + Co).The fijk in (8.24) is obtained by computing the volumetric integration on the right-handside of (8.23). At the source point, the first term and the second term cancel each other.At any other point, the first term is zero. Therefore we need to consider the second termonly. This integration can be separated in two parts:iv V • (crVO) dv =^krV20 Vo • Vq] dv.^(8.25)JvExcept at the source point, the first term on the right-hand side of (8.25) is zero. Thereforeiv V • (crVO) dv = jv Vcr • VO dv.^(8.26)CHAPTER 8: 3-D DC RESISTIVITY MODELING^ 142If the conductivity is constant within Vijk, the result of the integration is zero. If theconductivity is a smooth function within Viik, the integration can be computed analyticallyor numerically. At a point where the conductivity is discontinuous, such as at an interfacepoint, Ver does not exist. The integration on the left-hand side of (8.26) is, however,still computable. Let us consider a simple case, where Vijk contains one interface on an= constant plane (figure 8.4). Across the interface the conductivity jumps from a l tocr2. We divide Vij k into three parts Vi, V2 , and Vg so that V2 contains the interface andhas a thickness 2e. In Vi and Vg, the conductivity is smooth function and the integrationis easily done. In Vg we convert the volume integration into a surface integration80• (crVO) dv = I cr—On • ds.In the limiting case of e –+ 0 we obtain00lim^V • (crVO) dv = (cri – cr2) Tda-.0 v2where S is the area of the interface within the small volume. Therefore00fijk = (01 – az) —d8 + r Vcr • V0 dv.Ox^Vijk (8.27)z Figure 8.4 A integration domain VJk which contains discontinuityCHAPTER 8: 3-D DC RESISTIVITY MODELING 143Since both a and 4 are known functions, the integrations in (8.27) are computable. Innumerical computation, a is assumed to be constant within each cell. Therefore Va iszero. Equation (8.27) can be rewritten as0 ,!kik = (al — az) —aa.s Ox(8.28)The integration in (8.28) can be computed analytically. We consider the x = constantinterface shown in figure 8.4, which is confined by al < z < a2 and b1 < y < b2. Assumingthat a =1 and I = 2ir in (8.9), we have1[(x — x,)2 (ii Ya) 2^Z )9112and^00^X — X,Oz [(x x8)2^— Y8)2 (z — x8 )93/2.It is easy to verify thatdz (1 + z2)3/2 — (1 + z2)1/2anddy^1 a2)1/2J (a2 y2)02 + y2)1/2^ a(b2 + y2)1/2Using (8.29) and (8.30), we have^= a(b2 — a2)1/2 tan a(b2y(b2 — fa2 dz dyvi^Oz^/62 fa2^ 1(x xa) Jai [(z — z,, )2 + (y y8)2 (z z. )213/2^dz dy(8.29)(8.30)b2^1 Z — Z.(° — x, )^(x x,)2 (y y^i )2 [(x x8)2 (y yg )2 (z zs)211/2= — (z — z8 ) fb2^1 az — Z8(x^8)2 (y — y4 2 [E(0 — 08)2 (y y42 (a2 _ za)9 1/2al — x,[(x — 8)2 + (y — 42 + (al — z8)91/21= tan-1 ^(b2 — Y8)(a2 — z8) (x — x8)^— x8)2 + (1)2 — y8 )2 (a2 — z8 )2 1 1/ 2a2dyalCHAPTER 8: 3-D DC RESISTIVITY MODELING^ 144— tan-1^— y8)(a2 — z, )(X x8) [(x x,)2 + (bi y8)2 (a2 — z8)91/2— tan-1 ^(62 — Y8)(ai — z8) (x — x8 ) [(x — z8 )2 (62 — y8 )2 + (al — z8)2[1/2+ tan-1^— ya)(ai — z8) (8.31)(x xs) [(z z8)2 (14 y42 + (al _ z8 )2[ 1 /2.The integrations on y = constant and z = constant planes can be computed in the sameway. This approach is different from that presented by Lowry et al. (1989).For a ground surface point (xi, y, , 0), we define a rectangular parallelepiped 145 1 ={(x,Y,z) I (xi + xi+i)/2 5_ x (xi-i + xi)/2, + y,i+i)/ 2 5_ Y 5 + yi)/2,22 /2 < z < 0 } (figure 8.5). Because the normal derivative of the potential on the groundsurface is zero, the integration calculated for this surface is also zero. The resultant equationisjld. + 3 , _4Y(i-1)11 tijivi(j+1)1 =0-1b11j2Figure 8.5 The rectangular parallelepiped 3/451 for a point on the groundsurface1Eiji2Cull3Ciiic116CiiiandCHAPTER 8: 3-D DC RESISTIVITY MODELING^ 145The coefficients are1cro_1v10yi421] ,1^r= 4Azi 10'4 bai]+ aijiAY 1 ^r44i-1 1.(i-1)(1-1)1Axi—lbal + cri(i--1)1AziAgi]1 ,^4Ayi^+14Azi r=^1.47(1-1)(j-1)1Axi—lAW-1 ++ Gro-iviAzi-lAy; + aiiiAziAxdcal 1 2 3= — k cip + + c ii1 + 4 +c 6 pBecause CZik = C(i+i )jk , Ctik = 4i+i )k , and ctlik = 4.(k44), the coefficient matrix issymmetric and sparse. We need to compute and store only the upper half of the coeffi-cient matrix. The resultant equation system is solved by the conjugate gradient method(Kershaw 1978, Dey and Morrison 1979b, Meijerink et al. 1981) with a simple incompleteLDLT factorization preconditioner.8.5 Numerical ExamplesWe now present three numerical examples. We solve each problem by both the multi-domain Chebyshev spectral method and the finite-difference method. In the followingexamples, the x- and y-axes are horizontal and the z-axis upwards with the origin on theground surface.zxFigure 8.6 A two layered modelCHAPTER 8: 3-D DC RESISTIVITY MODELING^ 146Example 8.1 A two-layered modelThe model is shown in figure 8.6. The upper layer has a conductivity al and a thicknessa. The lower layer is a half space with a conductivity al. A current I = 2w A is injectedat the origin (0, 0, 0).For the two-layered model, an analytic solution is available. LetR2 =(x— w8)2 + (y — Y3)2andk al — =01 + 02(8.32)le l w^kal a2/ICHAPTER 8: 3-D DC RESISTIVITY MODELING^ 147At the point (x, y, z) the analytic solution is (Keller and Frischknecht 1966, pp104-112)00. ....._^I {^kn ^00 kn u1^(8.33)27ra1 n.1 [R2 + (2na — z)2]1/2 +nE.0 [R2 + (2na -I- z) 2 ] 1 /2if —a < z < 0 andir(cri 02) nE=0 [R2 (2na — z)91/2if z < —a. The analytic secondary potential b' can be obtained by= Ua — 4.,ita =ka(8.34)(8.35)where 4 is given by (8.9) with o = Q1. We use (8.33)—(8.35) to compute the accuratesecondary potential at grid points for comparison and also for the boundary conditionsin numerical modeling. In computing the analytic solution, when the absolute value of aterm (for a given n ) is less than 10 -10 , we stop. The accuracy of the truncated "analytic"solution is expected to be better than 10 —s. This can be proved by the following derivation.Suppose we truncate the series at a number M and z < —a. From (8.34), the absoluteerror of the truncated solution is00 kpP=M+1 [R2 + (2pa — zoo)2p/2ikim<11- (01 + 0'2) [R2 + (2Ma — z)91/2^'kr= ^I^ikim ^iki 2(01 01) [R2 + (2Ma — z)9 1 /2 1 — ikiFor this example, lki < 0.99 and IkI/(1 — lkl) < 100. SinceIikim ^< 10 -10 ,ir(cri 01) [R2 + (2Ma — .091/2 —the absolute error is given bylel < 10-2 .We can do a similar analysis for the case of —a < z < 0 and obtain the same results.CHAPTER 8: 3-D DC RESISTIVITY MODELING 148Now we report the results for four different cases. The differences between these casesare the size of the computational domain Si, the number of subdomains, as well as thevalues of the conductivities. In cases 1-3, we investigate the effects of the geometricalparameters. In case 4, we compare the results of our methods with those presented by Deyand Morrison (1979b) and Lowry et al. (1989).Case 1. We take CI = {(x,y, z) I -1 m < x < 1 m, -1 m < y < 1 m, -4 m < z < 0 m},a = 2 m, 0i = 1 S/m, and 02 = 0.01 S/m.In the multi-domain Chebyshev spectral method, we divide the whole domain into twosubdomains with the dividing plane z = -2 m. Both subdomains have a standard 3-DChebyshev domain 2 x 2 x 2 m 3 . We use (8.33)-(8.35) to compute the accurate boundaryconditions. Table 8.1 shows the error norms of the solutions obtained by the multi-domainChebyshev spectral method. In table 8.1, N is the degree of Chebyshev polynomials usedin each subdomain along one direction; N., N y , and N. the number of grid points alongthe x-, y-, and z-directions respectively; Pi the total number of grid points; L2 and Lc,.Table 8.1The error norms of example 8.1 (case 1, MDCS method)N N. Ny Pt L2 Loo CPU (sec)4 5 5 9 225 1.377 x 10-" 1.203 x 10 -°3 0.045 6 6 11 396 2.039 x 10-°3 1.128 x 10-" 0.096 7 7 13 637 2.979 x 10 -06 3.490 x 10 -4)3 0.207 8 8 15 960 4.792 x 10-°7 3.333 x 10 -06 0.368 9 9 17 1, 377 8.755 x 10-08 1.153 x 10-" 0.839 10 10 19 1, 900 2.182 x 10-08 1.541 x 10-" 1.3810 11 11 21 2,541 1.132 x 10 -08 7.919 x 10-08 2.46CHAPTER 8: 3-D DC RESISTIVITY MODELING^149Table 8.2The error norms of example 8.1 (case 1, FD method)h Nx Ny Nz Pt L2 CPU (sec)0.5000 5 5 9 225 6.978 x 10 -84 4.835 x 10 -83 0.180.2500 9 9 17 1, 377 1.990 x 10-04 1.439 x 0.330.1250 17 17 33 9,537 5.243 x 10 -85 3.776 x 10-84 4.080.0625 33 33 65 70,785 1.444 x 10 -05 9.566 x 10-85 52.25the corresponding error norms of the potential; and the CPU time is obtained on a SPARCstation 2 computer. From table 8.1 the error norms of the multi-domain Chebyshev spectralsolutions go to zero exponentially. On a grid with 2,541 points, the error norms are ti 10-8 .In the finite-difference method, we also use (8.33)—(8.35) to compute the accurateboundary conditions. We use a uniform grid. Table 8.2 shows the error norms of thefinite-difference solutions. In table 8.2, h is the grid size. All the other symbols have thesame meaning as those in table 8.1. From table 8.2 the error norms of the finite-differencesolutions are roughly proportional to h 2 . On a grid with 70,785 points, the error norms are10-5 .Comparing table 8.1 with table 8.2, the error norms of the multi-domain Chebyshevspectral solution on a grid with 637 points are better than those of the finite-differencesolution on a grid with 70,785 points. The CPU time is 0.20 seconds for the multi-domainChebyshev spectral method and 52.25 seconds for the finite-difference method. The multi-domain Chebyshev spectral method is far more efficient than the finite-difference method.Case 2. We take = {(z, z)i —1m<o<lm, —1 m<y<lm, —16 m<z<0 m}, a = 2 m, vi = 1 Sim, and al = 0.01 Sim. We divide the whole domain into twosubdomains with the dividing plane z = —2 m. The lower subdomain has a height 7 timesas large as that in case 1.CHAPTER 8: 3-D DC RESISTIVITY MODELING^ 150Table 8.3The error norms of example 8.1 (case 2, MDCS method)N Ary N. Pt L2 L00 CPU (sec)4 5 5 9 225 2.225 x 10—" 1.603 x 10 -03 0.045 6 6 11 396 8.631 x 10 -05 4.868 x 10—" 0.086 7 7 13 637 3.796 x 10 -05 2.467 x 10—" 0.197 8 8 15 960 2.046 x 10 -05 1.409 x 10—" 0.358 9 9 17 1, 377 1.292 x 10 -05 9.399 x 10 -135 0.839 10 10 19 1,900 7.686 x 10 -06 6.391 x 10 -05 1.4110 11 11 21 2,541 4.102 x 10 -06 4.377 x 10 -05 2.46Table 8.3 shows the error norms of the solutions obtained by the multi-domain Cheby-shev spectral method. All the symbols in table 8.3 have the same meaning as those in table8.1. Comparing table 8.3 with table 8.1, the error norms of the multi-domain Chebyshevspectral solutions in this case go to zero much slower than that in case 1. The poorer ac-curacy of the multi-domain Chebyshev spectral solutions can be explained as follows. Theonly difference between this case and case 1 is the vertical dimension of the lower subdo-main, therefore the degraded accuracy in this case is caused by the elongated shape of thelower subdomain. An analogous situation occurs in the finite-element method when ele-ments with large aspect ratios are employed. From (8.34) and (8.35) the analytic secondarypotential in the lower layer is/^oo^kn^ /IP = N^ E740'1 + g2) n=0 [R2 + (2na _ z)2]1/2 27.01 (R2 + z2)1/2*The term for Ti = 0 is{ ^I 1 ^Ik^1 01,) 7r(ai + 01)^21r0•1] (R2 + z2 ) 1 /2^27ra1 (R2 + z2)1/2 •(8.36)2z — (21 + 22) =•ii — 22(8.40)CHAPTER 8: 3-D DC RESISTIVITY MODELING^ 151All the other terms with n > 1 are smoother than 14. The dominant contribution to theerror in the Chebyshev expansion is principally determined by the approximation of then = 0 term in the analytic solution, to On the linez = zslY = Yi,R = 0 and (8.36) can be further simplified asIk 1114 = — 2r0.1 —z (x < —2).^ (8.37)The constant Ik/2iroi in (8.37) has little effect on the accuracy of the Chebyshev expan-sions. For further insight we consider a simple but more general functiong(z) = t 1^(8.38)wheret > li > z > i2.^ (8.39)When approximating g(z) by a Chebyshev expansion, we actually map the interval 22 <z < il to —1 < E < 1 byWith (8.40), point t is mapped toThe function g(z) becomes2t — (21 + 22) 6 =^_ _^•z1 — .z2 (8.41)9(0 —221 — 22 6 1— •The constant 2/(21 — 22) does not affect the convergence rate. At the point 6, g(0 issingular. This singularity affects the accuracy of the Chebyshev expansions. From thediscussion in section 2.1, the closer is 6 to 1, the more prominent of the effects of thesingularity on g(). In case 1, t = 0, 21 = —2, 22 = —4, and 6 = 3. In this case, t = 0,21 = —2, z2 = —16, and 6 = 9/7. In the mapped domain the singular point 6 of this caseCHAPTER 8: 3-D DC RESISTIVITY MODELING^ 152is closer to 1 than that of case 1. Therefore the accuracy in this case is poorer than thatin case 1 for a given N.Case 3. There are two ways to overcome the degraded accuracy in case 2. One way is toincrease the number of grid points in each subdomain and keep the number of subdomainsfixed, such as the results shown in table 8.3 for case 2. From table 8.3, as N increases,the error norms go to zero exponentially. The other is to divide the whole domain intomore subdomains. In this case, the model parameters are the same as those in case 2.Here we divide the whole domain into 4 subdomains with dividing planes z = —2, —4, and—8 m. Table 8.4 shows the error norms of the multi-domain Chebyshev spectral solutions.Comparing table 8.4 with table 8.1, we see that the error norms in case 1 and this case arealmost the same for a given N. Because the whole domain in this case is larger than thatin case 1, we have to use more grid points to achieve the same accuracy.Table 8.4The error norms of example 8.1 (case 3, MDCS method)N N. Ny NN Pt L2 L.. CPU (sec)4 5 5 17 425 2.658 x 10 -4)5 3.769 x 10 -04 0.075 6 6 21 756 5.106 x 10-06 5.809 x 10—°5 0.166 7 7 25 1225 7.968 x 10—°7 1.230 x 10 -05 0.427 8 8 29 1856 1.434 x 10—°7 2.505 x 10 -06 0.798 9 9 33 2673 1.945 x 10 -08 3.729 x 10—°7 1.799 10 10 37 3700 7.701 x 10 -09 8.478 x 10—°8 2.9410 11 11 41 4961 7.910 x 10 -09 5.536 x 10 -08 5.22Case 4. In this case, we assume f2 = {(x,y,z)10 m < x < 16 m, —16 m < y <16 m, —16 m < z < 0 m}, a = 2 m, = 0.01 Sim, and er2 = 0.1 Sim. Dey and MorrisonCHAPTER 8: 3-D DC RESISTIVITY MODELING^ 153iaa, = 0.01U2 = 0.1Figure 8.7 A dipole—dipole configuration on a two-layered model(1979b) and Lowry et al. (1989) used this model to test the accuracy of their methods.They computed the apparent resistivity for a dipole—dipole configuration (figure 8.7). Nowwe solve this problem by both the finite-difference and the multi-domain Chebyshev spectralmethods described in this chapter. For this model, we compute the apparent resistivityfor a dipole—dipole configuration along the x-axis as follows. Suppose we know the totalpotential (analytically or numerically) ui at xi = ia, where a is the length of the dipole andi = 1, 2, 3, • Because I = 2w A, the apparent resistivity pi at a separation ia is given byPi^ •Li xi —21zi+i +11zi+2In the multi-domain Chebyshev spectral method, we divide the whole domain into4 x 8 x 4 = 128 subdomains. The dividing planes are z = 2,4, 8 m; y = —8, —4, —2,0, 2, 4, 8 m; and z = —8, —4, —2 m. We use the same N th degree Chebyshev polynomialsin all subdomains. After solving the problem, we compute the secondary potential th onthe x-axis by 1-D Chebyshev interpolation. The total potential is given by (8.7). Theresults are shown in table 8.5, in which MARE is the maximum absolute relative error ofthe apparent resistivity and all the other symbols have the same meaning as those in tableUi — 214+1 ui-F2CHAPTER 8: 3-D DC RESISTIVITY MODELING 1548.1. In the finite-difference method, we use a uniform grid. The results are shown in table8.6, in which h is the grid size, all the other symbols have the same meaning as those intable 8.5.Comparing table 8.5 with 8.6, we find that on a grid with 9,537 points the finite-difference results are better than the multi-domain Chebyshev spectral results. However,when we increase the number of grid points, the maximum absolute relative error of theapparent resistivity obtained by the multi-domain Chebyshev spectral method goes to zeroexponentially as N increases. On an N = 6 grid with 30,625 points, it is 8.459 x 10-5 . TheCPU time on a SPARC station 2 computer is 178.21 seconds. In contrast, the maximumabsolute relative error of the apparent resistivity obtained by the finite-difference methodgoes to zero slowly, roughly proportional to h 2 . On a grid with 545,025 points, it is5.323 x 10-4 . The CPU time is 904.55 seconds. When we want to obtain higher accuracy,such as 10-5 , the multi-domain Chebyshev spectral method becomes more efficient thanthe finite-difference method.Table 8.5The results of example 8.1 (case 4, MDCS method)N N. Ny N.^Pt MARE CPU (sec)45617212533414917^9,53721^18, 08125^30,6251.871 x 10 -022.300 x 10 -038.459 x 10 -0520.1966.54178.21Table 8.6The results of example 8.1 (case 4, FD method)h N. Ny N.^Pt MARE CPU (sec)1.00 17 33 17^9,537 8.066 x 10 -03 4.920.50 33 65 33^70,785 2.105 x 10 -03 65.620.25 65 129 65^545,025 5.323 x 10 -04 904.55CHAPTER 8: 3-D DC RESISTIVITY MODELING^ 155Example 8.2 A two-prism modelFigure 8.8 shows a model composed of two prisms, Pi and P2, embedded in an otherwiseuniform half space. The conductivity of the half space is 0.01 Sim. The dimension andconductivity of the two prisms are shown in table 8.7. In the multi-domain Chebyshevspectral method, the whole domain is divided into 5 x 5 x 3 = 75 subdomains. Thedimensions of the subdomains are 70, 20, 20, 20, 70 m along the x- and y-directions and20, 30,50 m along the z-direction. The Chebyshev grid is given by (8.15) and (8.16).For this model no analytic formulae are available for the secondary potential. We set thexFigure 8.8 A two-prism modelTable 8.7The dimension and conductivity of the two-prism modelx(m)^y(m)^x(m)^a(S/m)Pi 30 — 10 30 —^10 0 — — 20 0.001P2 - 10 - - - 30 30 — — 30 0 — — 20 0.100CHAPTER 8: 3-D DC RESISTIVITY MODELING^ 156xFigure 8.9 The secondary potential on the ground surface of example 8.2(a) The multi-domain Chebyshev spectral solution on a grid with 31 x 31 x 19 =18,259 points. (b) The finite-difference solution on a grid with 81 x 81 x 41 =269,001 points. The contour intervals are 5 mV. The solid lines represent thepositive values and the dash lines represent the negative values. The star in thecenter indicates the location of the source current. The two rectangles near thecenter indicate the location of the the inhomogeneous bodies.CHAPTER 8: 3-D DC RESISTIVITY MODELING^ 157secondary potentials at the external boundaries to zero. Figure 8.9(a) shows the secondarypotentials on the ground surface obtained by the multi-domain Chebyshev spectral methodon an N = 6 grid with 31 x 31 x 19 = 18, 259 points. In the finite-difference method,we use a uniform grid and the same boundary conditions as those in the multi-domainChebyshev spectral method. Figure 8.9(b) shows the secondary potential on the groundsurface obtained by the finite-difference method on a grid with 81 x 81 x 41 = 269, 001points. In figure 8.9, the solid lines are for the positive values and the dash lines are for thenegative values. The contour intervals are 5 mV. The results of the two different methodsare in good agreement. The maximum difference is less than 2 mV.Example 8.3 A two-prism model with piecewise-smooth conductivity variation.Geometrically this model is the same as that in example 8.2 except that in the twoprisms the conductivity is smooth spatial function. The conductivity in P1 is0 1 = 0.001 [1 + 9(y 40+ 10)1 Sint.When y varies from —10 to 30, al varies from 0.001 Sim to 0.01 Sim. The conductivityin P2 is—0'2 = 0.1 [1 + 0.9( y60 30)]When y varies from —30 to 30, o2 varies from 0.01 Sim to 0.1 Sim. Figure 8.10 (a) showsthe secondary potential on the ground surface obtained by the multi-domain Chebyshevspectral method on an N = 6 grid with 31 x 31 x 19 = 18, 259 points. Figure 8.10(b) showsthe results obtained by the finite-difference method on a grid with 81 x 81 x 41 = 269, 001points. The results shown in the two panels are also in good agreement. The maximumdifference is less than 2 mV.CHAPTER 8: 3-D DC RESISTIVITY MODELING^ 158x0^50^100xFigure 8.10 The secondary potential on the ground surface of example 8.3Legend is the same as that of figure 8.9CHAPTER 8: 3-D DC RESISTIVITY MODELING^ 1598.6 DiscussionIn this chapter, we applied the multi-domain Chebyshev spectral methods to 3-D DCresistivity modeling. To compare the results, we modified the finite-difference methodpresented by Dey and Morrison (1979b) and Lowry et al. (1989). We presented moreaccurate formulae for the removal of the source singularity and for the computation of theright-hand side terms of the finite-difference equations.For the two-layered model in example 8.1, there is an analytic solution for comparison.From tables 8.1, and 8.3-8.5, the error norms of the multi-domain Chebyshev spectralsolutions have exponential convergence. In case 1, with a total 2,541 grid points, the errornorms are ti 10-8 . When the geometrical aspects of the domain varies, the accuracy doesvary, but for all the four cases considered, the convergence is exponential. From tables 8.2and 8.6, the error norms of the finite-difference solutions are roughly proportional to h 2 ,where h is the grid size.For the same two-layered model considered in example 8.1 case 4, Dey and Morrison(1979b, p763) reported that the maximum absolute relative error of the apparent resistivitywas better than 5 x 10 -2 on a grid with 57 x 17 x 12 = 11, 628 points. In their results,the error came from the effect of the singularity and the mixed boundary conditions. Withsingularity removal, Lowry et al. (1989) reported that the accuracy had been improved toabout 1.5 x 10-2 . In their results, the error came from the incorrect formula for the removalof the source singularity and the mixed boundary conditions. In addition, they computedthe right-hand side terms by numerical differentiation of the primary field, which also causederrors. From table 8.6, the maximum absolute relative error of the apparent resistivity ofthe present finite-difference solutions is 8.006 x 10 -3 on a grid with 17 x 33 x 17 = 9,357points. It is better than the result reported by Lowry et al. (1989).CHAPTER 8: 3-D DC RESISTIVITY MODELING 160In examples 8.2 and 8.3, because there are no analytic solutions available, we approx-imately set the secondary potential on the boundary to zero. In the two examples, theinhomogeneities are located within a small area. The secondary potential is mainly a di-pole field and is proportional to the inverse square of the distance. The maximum absolutevalue of the secondary potential is less than 52 mV. The minimum distance from the in-homogeneous body to the boundary is 70 m. Therefore the error caused by the approximateboundary conditions is less than 0.02 mV.Comparison of figure 8.9(a) with figure 8.9(b) and figure 8.10(a) with figure 8.10(b)illustrates that the solutions obtained by the multi-domain Chebyshev spectral methodand finite-difference method agree with each other quite well. The difference between themis less than 2 mV.For all the models presented in this chapter, the multi-domain Chebyshev spectralmethod is more efficient than the finite-difference method. In example 8.1 case 1, theerror norms of the multi-domain Chebyshev spectral solution on a grid with 637 points arebetter than those of the finite-difference solution on a grid with 70,785 points. On a SPARCStation 2 computer, the CPU time is 0.20 seconds for the multi-domain Chebyshev spectralmethod and 52.25 seconds for the finite-difference method. The multi-domain Chebyshevspectral method is about 260 times more efficient than the finite-difference method. Inexample 8.1 case 4, the maximum absolute relative error of the apparent resistivity of themulti-domain Chebyshev spectral solution on a grid with 30,625 points is 8.549 x 10 -5 .The maximum absolute relative error of the apparent resistivity of the finite-differencesolution on a grid with 545,025 points is 5.323 x 10 -4 . The CPU time is 178.21 seconds forthe multi-domain Chebyshev spectral method and 904.55 seconds for the finite-differencemethod. The multi-domain Chebyshev spectral method is about 5 times more efficient thanthe finite-difference method. In examples 8.2 and 8.3, the CPU time is 110.98 seconds forthe multi-domain Chebyshev spectral method and 225.86 seconds for the finite-differenceCHAPTER 8: 3-D DC RESISTIVITY MODELING^ 161method. The multi-domain Chebyshev spectral method is about 2 times more efficientthan the finite-difference method.The issue of efficiency still needs more qualification, however, since different linear-system solvers were used in the finite-difference and multi-domain Chebyshev spectral meth-ods. The most effective preconditioned conjugate gradient method (Saad and Shultz 1986)was used to solve the finite-difference equations. On the other hand, a simple Gauss-Seideliteration was used to solve the influence matrix system in the multi-domain Chebyshevspectral method. This may put the multi-domain Chebyshev spectral method at a disad-vantage. The generalized conjugate residual method (Saad and Shultz 1986) may be moreefficient than the simple Gauss-Seidel iteration in solving the influence matrix system. Thepreconditioned biconjugate gradient method (Press et al. 1992, p.77) is another alternativefor solving the discrete equation system arising from the multi-domain Chebyshev spectralmethod. This will be the subject of further research.CHAPTER 9SUMMARYThe purpose of the research presented in this thesis was to develop and apply Chebyshevspectral methods to solve the partial differential equations arising in potential field prob-lems. In the existing finite difference, finite element, and integral equation methods, a locallow order function is used to approximate the unknown function. To obtain high accuracy,a large number of grid points have to be used. In the Chebyshev spectral methods, we usea set of globally smooth functions, Chebyshev polynomials, to approximate the unknownfunction. When f and a are piecewise-smooth functions with interface discontinuities par-allel to the coordinate planes (that is, a box-like structure), the Chebyshev spectral methodcan achieve higher accuracy with fewer grid points than the finite-difference method. Inthese situations, the Chebyshev spectral methods can be more efficient than the finite-difference method. But when a and f are generally discontinuous or the domain is farfrom being a box-like structure, this is no longer true. Then local methods such as thefinite difference or finite element methods are more applicable. However, in geophysicalapplications, when a smooth node-based model parameterization rather than a generallydiscontinuous cell-based model parameterization is employed, then the Chebyshev spectralmethods may again be used. The choice of separate model blocks in the multi-domainapproach must be such that regions with large conductivity gradients coincide with thedensely-spaced Chebyshev interpolation points in the interior of the block.In chapter 2, we outline the fundamentals of Chebyshev expansions. A function can beapproximated by a Chebyshev series. When there are discontinuities in the function, itsChebyshev expansion will contain infinitely many modes. If we approximate the function162CHAPTER 9: SUMMARY 163by a finite Chebyshev series, there will be truncation and aliasing errors. For a fixed numberof terms, the truncation error is fixed. We can use a cell-average discretization scheme toreduce the aliasing error. However, the overall accuracy is still of 0(1/N). A better schemeconsists of dividing the domain into a number of subdomains so that in each subdomain thefunction is infinitely differentiable and then approximating the function in each subdomainby a separate Chebyshev series.In chapters 3 we use both the r method and the collocation method to solve the Pois-son's equation in a single domain. In both cases the discrete equation system can be solvedefficiently by the matrix-diagonalization method. The condition number of the coefficientmatrix of the equation system is proportional to N 4 . To obtain a satisfactory accuracy, weshould use double precision in the computation. The accuracy of the collocation methodis slightly better than that of the r method. In the r method we have to use forwardand inverse Chebyshev transforms to convert data between physical space and transformedspace. In the collocation method, no Chebyshev transforms are required. Therefore thecollocation method is also faster than the r method.In chapter 4, we solve the Poisson's equation by the multi-domain Chebyshev spec-tral method. When there are discontinuities in the source function, the single-domainChebyshev spectral method cannot achieve exponential convergence. However, we can usethe multi-domain Chebyshev spectral method to achieve this convergence rate. In themulti-domain Chebyshev spectral method, the whole domain is divided into a number ofsubdomains so that in each subdomain the function is infinitely differentiable. Then thepotential function is approximated by a separate Chebyshev series in each subdomain. Therelations of the Chebyshev series in adjacent subdomains are determined by the matchingconditions at the interface: the continuity of the potential and the first order normal de-rivatives. In this way we can achieve the desired exponential convergence. We solve theequation system by the influence matrix method. This method works. But there is noguarantee that it is the best method to solve the problem.CHAPTER 9: SUMMARY 164In chapter 5, we presented a new iterative Chebyshev spectral method for solving thegeneralized Poisson's equation V • (uVu) = f in a single domain. If u is smooth, thegeneralized Poisson's equation can be rewritten in the form of a Poisson's equation V 2u(f — Va • Vu)/a. The initial values are given arbitrarily. In each iteration, the right handside terms are computed first from the current values of the potential. Then the resultantPoisson's equation is solved by the direct collocation method to obtain the updated values.The process is repeated until the norms of the differences between two successive iterationssatisfy an appropriate convergence criterion. For the examples computed, the convergencerate of the new method is better than that of the Chebyshev spectral multigrid methodpresented by Zang et al. (1982). For 1-D problem, the convergence rate is bounded by> ln(r7r) — In [max (dln a/dx)] . When a contains discontinuities or when R. is too smallto be acceptable, we can divide the whole domain into a number of subdomains so thatin each subdomain a is continuous and R. is large enough, and solve the problem by themulti-domain Chebyshev spectral method.In chapter 6, we solve the generalized Poisson's equation by the multi-domain Cheby-shev spectral method. The procedure is similar to that for Poisson's equation in chapter4. There are two differences. Firstly, for the generalized Poisson's equation the interfaceconditions require that both the potential and the flux, aOu/On, be continuous. Secondly,we have to solve a generalized Poisson's equation instead of a simpler Poisson's equation ineach subdomain. We solve the discrete equation system by the influence matrix method.Like that for Poisson's equation, this method works. But there is no guarantee that it isthe best method to solve the problem.In chapter 7, we apply the Chebyshev spectral methods to 3-D gravity modeling. Weextend Okabe's (1979a) analytic formula for the gravitational field of a polyhedral body tothe potential. We present a method to realize the point-injection and cell-average discret-ization for 2-D and 3-D problems. The multipole expansion technique is used to computethe approximate boundary conditions. Numerical results suggest that the error norms ofCHAPTER 9: SUMMARY 165the single-domain solution go to zero as order (1/N 2 ) while the multi-domain Chebyshevspectral method retains exponential convergence. The error caused by the approximateboundary conditions is about 1%.In chapter 8, we apply the multi-domain Chebyshev spectral method to 3-D DC resis-tivity modeling. We discuss the singularity caused by the source current. For comparisonwe also present a modified finite difference method which includes a more accurate formulafor the removal of the source singularity and an analytic formula for the computation ofthe right-hand terms of finite-difference equations. For a two-layered model, there is ananalytic solution available for comparison. The accuracy of the modified finite differencemethod is better than that reported by Dey and Morrison (1979b) and Lowry et al. (1989).The multi-domain Chebyshev spectral solutions are more accurate than the finite differencesolutions. For the piecewise-constant and piecewise-smooth models, there are no analyticsolutions available for comparison. The results obtained by both the multi-domain Cheby-shev spectral method and the finite-difference method agree with each other quite well forthese cases. For the numerical examples presented, the multi-domain Chebyshev spectralmethod is faster than the finite-difference method.The Chebyshev multi-domain spectral methods implemented in the thesis may be usedas the forward modelling engine when considering that portion of model space exploredvia a smooth model inversion. When there are a significant number of piecewise-smoothconductivity blocks, the computations can be done in parallel, with the calculations foreach block assigned to a computational node. Inter-node communication is necessary onlywhen computing the interface boundary error. These foregoing topics will be studied infuture research efforts.REFERENCESAlfano, L., 1959, Introduction to the interpretation of resistivity measurements for compli-cated structural conditions. Geophysical Prospecting, 7, 311-366.Brandt, A., Fulton, S. R., and Taylor, G. D., 1985, Improved spectral multigrid methodsfor periodic elliptic problems. Journal of Computational Physics, 58, 96-112.Barnett, C. T., 1976, Theoretical modelling of the magnetic and gravitational fields of anarbitrarily shaped three-dimensional body. Geophysics, 41, 1353-1364.Bhattacharyya, B. K. and Navolio, M. E., 1975, Digital convolution for computing gravityand magnetic anomalies due to arbitrary bodies. Geophysics, 40, 981-992.Bhattacharyya, B. K. and Navolio, M. E., 1976, A fast Fourier transform method for rapidcomputation of gravity and magnetic anomalies due to arbitrary bodies. GeophysicalProspecting, 24, 633-649.Bibby, H. M., 1978, Direct current resistivity modeling for axially symmetric bodies usingthe finite-element method. Geophysics, 43, 550-562.Botezatu, R., Visarion, M., Scurtu, F., and Cucu, G., 1971, Approximation of the gravita-tional attraction of geological bodies. Geophysical Prospecting, 19, 218-227.Buzbee, B. L., Dorr, F. W., George, J. A., and Golub, G. H., 1971, The direct solution of thediscrete Poisson equation on irregular regions. SIAM Journal on Numerical Analysis,8, 722-736.Cady, J. W., 1980, Calculation of gravity and magnetic anomalies of finite-length rightpolygonal prisms. Geophysics, 45, 1507-1512.Canuto, C., Hussaini, M. Y., Quarteroni, A., and Zang, T. A., 1988, Spectral Methods inFluid Dynamics. Springer-Verlag, New York.166REFERENCES 167Charbeneau, R. J., and Street, R. L., 1979, Modeling groundwater flow fields containingpoint singularities: a technique for singularity removal. Water Resources Research, 15,583-594.Clenshaw, C. W., 1957, The numerical solution of linear differential equations in Chebyshevseries. Proceedings of Cambridge Philosophical Society [A], 53, 134-149.Coggon, J. H., 1971, Electromagnetic and electrical modeling by the finite element method.Geophysics, 36, 132-155.Cooley, J. W. and 'Palmy, J. W., 1965, An algorithm for the machine calculation of complexFourier series. Mathematics of Computation, 19, 297-301.Cooley, J. W., Lewis, P. A. W., and Welch, P. D., 1970, The fast Fourier transform al-gorithm: programming considerations in the calculation of sine, cosine and Laplacetransforms. Journal of Sound and Vibration, 12, 315-337.Daniels, J. J., 1977, Three-dimensional resistivity and induced polarization modeling usingburied electrodes. Geophysics, 42, 1006-1019.Delves, L. M. and Hall, C. A., 1979, An implicit matching principle for global elementcalculations. Journal of the Institute of Mathematics and Its Applications [A pub], 23,223-234.Deville, M. and Mund, E., 1985, Chebyshev pseudospectral solution of second-order ellipticequations with finite element preconditioning. Journal of Computational Physics, 60,517-533.Deville, M. 0. and Mund, E. H., 1990, Finite-element preconditioning for pseudospectralsolutions of elliptic problems. SIAM Journal on Scientific and Statistical Computing,11, 311-342.Dey, A. and Morrison, H. F., 1979a, Resistivity modeling for arbitrarily shaped two-dimensional structures. Geophysical Prospecting, 27, 106-136.Dey, A. and Morrison, H. F., 1979b, Resistivity modeling for arbitrarily shaped three-dimensional structures. Geophysics, 44, 753-780.REFERENCES^ 168Dieter, K., Paterson, N. R., and Grant, F. S., 1969, IP and resistivity type curves forthree-dimensional bodies. Geophysics, 34, 615-632.Eskola, L., Soininen, H., and Oksama, M., 1989, Modelling of resistivity and IP anomaliesof a thin conductor with an integral equation. Geoexploration, 26, 95-104.Fox, L. and Parker, I. B., 1968, Chebyshev Polynomials in Numerical Analysis. OxfordUniversity Press, London.Fox, R. C., Hohmann, G. W., Killpack, T. J., and Rijo, L., 1980, Topographic effects inresistivity and induced-polarization surveys. Geophysics, 45, 75-93.Golub, G. H. and Van Loan, C. F., 1983, Matrix Computations. John Hopkins UniversityPress, Baltimore.Gottlieb, D. and Lustman, L., 1983, The spectrum of the Chebyshev collocation operatorfor the heat equation. SIAM Journal on Numerical Analysis, 20, 909-921.Gottlieb, D. and Orszag, S. A., 1977, Numerical Analysis of Spectral Methods: Theory andApplications. SIAM-CBMS, Philadelphia.Gottlieb, D., Hussaini, M. Y., and Orszag, S. A., 1984, Theory and applications of spectralmethods. In Spectral Methods for Partial Differential Equations, ed. by Voigt, R. G.,Gottlieb, D., and Hussaini, M. Y., SIAM-CBMS, Philadelphia, 1-54.GOtze, H. J. and Lahmeyer, B., 1988, Application of three-dimensional interactive modellingin gravity and magnetics. Geophysics, 53, 1096-1108.Grant, F. S. and West, G. F., 1965, Interpretation theory in applied geophysics. McGraw-Hill, New York.Haidvogel, D. B. and Zang, T., 1979, The accurate solution of Poisson's equation by ex-pansion in Chebyshev polynomials. Journal of Computational Physics, 30, 167-180.Haldenwang, P., Labrosse, G., Abboudi, S., and Devine, M., 1984, Chebyshev 3-D spectraland 2-D pseudospectral solvers for the Helmholtz equation. Journal of ComputationalPhysics, 55, 115-128.REFERENCES^ 169Hansen, R. 0. and Wang, X., 1988, Simplified frequency-domain expressions for potentialfields of arbitrary three-dimensional bodies. Geophysics, 53, 365-374.Holcombe, H. T. and Jiracek, G. R., 1984, Three-dimensional terrain corrections in resis-tivity surveys. Geophysics, 49, 439-452.Horner, T. S., 1982, A double Chebyshev series method for elliptic partial differentialequations. In Numerical Solutions of Partial Differential Equations, Noye, J. (editor),North-Holland Publishing Company.James, B. A., 1985, Efficient microcomputer-based finite-difference resistivity modeling viaPolozhii decomposition. Geophysics, 50, 443-465.Keller, G. V. and Frischknecht, F. C., 1966, Electrical methods in geophysical prospecting.Pergamon Press, New York.Kershaw, D. S., 1978, The incomplete Cholesky-conjugate gradient method for the iterativesolution of systems of linear equations. Journal of Computational Physics, 26, 43-65.Lanczos, C., 1938, Trigonometric interpolation of empirical and analytical functions.Journal of Mathematics and Physics, 17, 123-199.Lanczos, C., 1957 , Applied Analysis. Sir Isaac Pitman & Sons, Ltd., London.Lee, T., 1975, An integral equation and its solution for some two- and three-dimensionalproblems in resistivity and induced polarization. Geophysical Journal of the Royal As-tronomical Society, 42, 81-95.Lowry, T., Allen, M. B., and Shive, P. N., 1989, Singularity removal: A refinement ofresistivity modeling techniques. Geophysics, 54, 766-774.MacMillan, W. D., 1958, The Theory of the Potential. Dover Publications, Inc., New York.McGillivray, P. R., 1992, Forward modeling and inversion of DC resistivity and MMRData, Ph. D. thesis, Department of Geophysics and Astronomy, University of BritishColumbia.REFERENCES^ 170Meijerink, J. A. and Van der Vorst, H. A., 1981, Guidelines for the usage of incompletedecompositions in solving sets of linear equations as they occur in practical problems.Journal of Computational Physics, 44, 134-155.Morchoisne, Y., 1984, Inhomogeneous flow calculations by spectral methods. In SpectralMethods for Partial Differential Equations, ed. by Voigt, R. G., Gottlieb, D., and Hus-saini, M. Y., SIAM-CBMS, Philadelphia, 181-208.Mufti, I. R., 1975, Iterative gravity modelling by using cubical blocks. Geophysical Pro-specting, 23, 163-198.Mufti, I. R., 1976, Finite-difference resistivity modeling for arbitrarily shaped two-dimen-sional structures. Geophysics, 41, 62-78.Nagy, D., 1966, The gravitational attraction of a right rectangular prism. Geophysics, 31,362-371.Nettleton, L. L., 1942, Gravity and magnetic calculations. Geophysics, 7, 293-310.Odegard, M. E. and Berg, J. W., 1965, Gravity interpretation using the Fourier integral.Geophysics, 30, 424-438.Okabe, M., 1979a, Analytical expressions for gravity anomalies due to homogeneous poly-hedral bodies and translations into magnetic anomalies. Geophysics, 44, 730-741.Okabe, M., 1979b, Treatment of singularities in the integral equation approach. Geophysics,44, 2004-2006.Okabe, M., 1982, Analytical expressions for gravity due to homogeneous revolutional com-partments in the Gaussian divergence approach. Geophysical Prospecting, 30, 166-187.Orszag, S. A., 1971, Accurate solution of the Orr-Sommerfeld stability equation. Journalof Fluid Mechanics, 50, part 4.689-703.Orszag, S. A., 1980, Spectral methods for problems in complex geometries. Journal ofComputational Physics, 37, 70-92.Patera, A. T., 1984, A spectral element method for fluid dynamics: Laminar flow in achannel expansion. Journal of Computational Physics, 54, 468-488.REFERENCES^ 171Pedersen, L. B., 1978, Wavenumber domain expressions for potential fields from arbitrary2-, 21-, and 3-dimensional bodies. Geophysics, 43, 626-630.Peyret, R., 1986, Introduction to spectral methods, von Karman Institute Lecture Series,1986-04 Rhode-Saint Genese, Belgium.Poirmeur, C., and Vasseur, G., 1988, Three-dimensional modeling of a hole-to-hole electricalmethod: Application to the interpretation of a field survey. Geophysics, 53, 402-414.Pratt, D. A., 1972, The surface integral approach to the solution of the 3-D resistivityproblem. Bulletin, Australian Society Exploration Geophysicists, 3, No. 4,33-50.Press, W. IL, Teukolsky, S. A., Vetterling, W. T., and Flannery, B. P., 1992, NumericalRecipes in FORTRAN, the Art of Scientific Computing, second edition. CambridgeUniversity press.Pridmore, D. F., Hohmann, G. W., Ward, S. H., and Sill, W. R., 1981, An investigationof finite-element modeling for electrical and electromagnetic data in three dimensions.Geophysics, 46, 1009-1024.Rivlin, T. J., 1974, The Chebyshev polynomials, John Wiley Sc Sons, New York.Snyder, M. A., 1966, Chebyshev Methods in Numerical Approximation, Prentice-Hall, Inc.Saad, Y. and Schultz, M. H., 1986, GEMRES: a generalized minimal residual algorithmfor solving nönsymmetric linear system. SIAM Journal on Scientific and StatisticalComputing, 7, 856-869.Spiegel, R. J., Sturdivant, V. R., and Owen, T. E., 1980, Modeling resistivity anomaliesfrom localized voids under irregular terrain. Geophysics, 45, 1164-1183.Talwani, M., 1973, Computer usage in the computation of gravity anomalies. Methods inComputational Physics, 13, 343-389.Talwani, M. and Ewing, M., 1960, Rapid computation of gravitational attraction of three-dimensional bodies of arbitrary shape. Geophysics, 25, 203-225.Telford, W. M., Geldart, L. P., Sheriff, R. E., and Keys, D. A., 1976, Applied Geophysics.Cambridge University Press, Cambridge.REFERENCES^ 172Varga, R. S., 1962, Matrix iterative analysis. Prentice-Hall Inc., New Jersey.Wilkinson, J. H., 1965, The Algebraic Eigenvalue Problem. Clarendon Press, Oxford.Xu, S. Z., Gao, Z. C., and Zhao, S. K., 1988, An integral formulation for three-dimensionalterrain modeling for resistivity surveys. Geophysics, 53, 546-552.Xu, S. Z., Zhao, S. K., and Lou, Y. J., 1984, The boundary element method for solving thethree dimensional electric field under the conditions of horizontal ground. ComputingTechniques for Geophysical and Geochemical Exploration, 6, No. 3,53-60 (in Chinese).Zang, T. A., Wong, Y. S., and Hussaini, M. Y., 1982, Spectral multigrid methods for ellipticequations. Journal of Computational Physics, 48, 485-501.Zang, T. A., Wong, Y. S., and Hussaini, M. Y., 1984, Spectral multigrid methods for ellipticequations II. Journal of Computational Physics, 54, 489-507.Zang, T. A., Streett, C. L., and Hussaini, M. Y., 1989, Spectral Methods for CFD. ICASEreport, No. 89-13.Zhao, S. K. and Yedlin, M. J., 1991, Chebyshev expansions for the solution of the forwardgravity problem. Geophysical Prospecting, 39, 783-802.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Chebyshev spectral methods for potential field computation
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Chebyshev spectral methods for potential field computation Zhao, Shengkai 1993
pdf
Page Metadata
Item Metadata
Title | Chebyshev spectral methods for potential field computation |
Creator |
Zhao, Shengkai |
Date Issued | 1993 |
Description | We investigate Chebyshev spectral methods for solving Poisson's equation and the generalized Poisson's equations and apply them to 3-D gravity and DC resistivity modeling. When we approximate a function by a finite sum of Chebyshev polynomials, there are two kinds of errors: truncation and aliasing. We present a cell-average discretization scheme to reduce the aliasing error. Both theoretical analysis and numerical examples show that when there are discontinuities in a function, the cell-average discretization is better than the point injection discretization. We use both the r and collocation methods to solve Poisson's equation V2u = f. We solve the discrete systems by a matrix-diagonalization method. The speed and accuracy of collocation methods are better than those of the r method. For the generalized Poisson's equation V • (oVu) = f, we present a new iterative method. We rewrite the equation in the form of a Poisson's equation, V2u = (f —Va•Vu)la. At each iteration we compute the right hand-side term from the current value of u first. Then we solve the resultant Poisson’s equation by the collocation method. Numerical results show that the convergence rate of the new method is much faster than that of the spectral multi grid method. When there are discontinuities in the source function f and/or the conductivity o, the single-domain Chebyshev spectral method does not converge exponentially. In this case we use the multi-domain Chebyshev spectral method to solve the problem. We divide the whole domain into a number of subdomains so that in each subdomain the function is infinitely differentiable. Then we approximate the function by a separate Chebyshev series in each subdomain. We determine the relations between the Chebyshev polynomials in adjacent subdomains by the interface conditions. In this way we can achieve exponential convergence. In 3-D gravity modeling, we present a method to realize 2-D and 3-D point-injection and cell-average discretizations, use a multi pole expansion to compute the approximate boundary conditions and solve the equations by both single-domain and multi-domain Chebyshev spectral methods. We also extend Okabe's analytic formulation for the gravitational field of a homogeneous, polyhedral body to the potential. Numerical results show that the accuracy of the cell-average discretization is better than that of the point injection discretization. The multi-domain solution is the best. We apply the multi-domain Chebyshev spectral method to 3-D DC resistivity modeling. We discuss the singularity removal and present a modified finite-difference method for comparison. For a two-layered model, the multi-domain Chebyshev spectral solution is far more accurate than the finite-difference solution. For piecewise-constant and piecewise-smooth models the solutions obtained by both methods agree with each other quite well. However, the Chebyshev spectral method is more efficient than the finite-difference method. |
Extent | 6959227 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2008-09-17 |
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.0085632 |
URI | http://hdl.handle.net/2429/2110 |
Degree |
Doctor of Philosophy - PhD |
Program |
Astronomy |
Affiliation |
Science, Faculty of Physics and Astronomy, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 1993-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_1993_fall_phd_zhao_shengkai.pdf [ 6.64MB ]
- Metadata
- JSON: 831-1.0085632.json
- JSON-LD: 831-1.0085632-ld.json
- RDF/XML (Pretty): 831-1.0085632-rdf.xml
- RDF/JSON: 831-1.0085632-rdf.json
- Turtle: 831-1.0085632-turtle.txt
- N-Triples: 831-1.0085632-rdf-ntriples.txt
- Original Record: 831-1.0085632-source.json
- Full Text
- 831-1.0085632-fulltext.txt
- Citation
- 831-1.0085632.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0085632/manifest