Orthogonal collocation revisited

Accepted Manuscript Orthogonal collocation revisited Larry C. Young

PII: DOI: Reference:

S0045-7825(18)30522-X https://doi.org/10.1016/j.cma.2018.10.019 CMA 12124

To appear in:

Comput. Methods Appl. Mech. Engrg.

Received date : 15 July 2017 Revised date : 16 August 2018 Accepted date : 12 October 2018 Please cite this article as: L.C. Young, Orthogonal collocation revisited, Comput. Methods Appl. Mech. Engrg. (2018), https://doi.org/10.1016/j.cma.2018.10.019 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

*Manuscript Click here to download Manuscript: OCRevisit.pdf

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Click here to view linked Referenc

Orthogonal Collocation Revisited Larry C. Young Dedicated to Professor Bruce A. Finlayson Table of Contents 1. Abstract ...............................................................................................................................1 2. Introduction..........................................................................................................................2 3. Historical Perspective ..........................................................................................................6 3.1 Galerkin Method.............................................................................................................7 3.2 Method of Moments .......................................................................................................7 3.3 Least Squares ................................................................................................................8 3.4 Collocation .....................................................................................................................8 3.5 Other Weight Functions ...............................................................................................11 3.6 Tau Method ..................................................................................................................11 3.7 Boundary Conditions ....................................................................................................12 3.8 Finite Elements. ...........................................................................................................13 4. Boundary Value Problem - Catalyst Pellet Problem ..........................................................14 4.1 Orthogonal Collocation Method....................................................................................15 4.1.1 Linear Source, Dirichlet B.C. .....................................................................................16 4.1.2 Linear Source, Third Kind B.C. .................................................................................19 4.2 Method of Moments .....................................................................................................21 4.3 Galerkin Method...........................................................................................................23 4.3.1 Linear Source, Third Kind B.C., Galerkin/Moments ..................................................26 4.4 Mass Conservation and Fluxes ....................................................................................27 5. Parabolic Problem - Falling Liquid Film .............................................................................29 5.1 Continuous Solutions in z ............................................................................................31 5.2 Summary Remarks - Global Methods and Fluxes ........................................................37 6. Multidimensional Applications ...........................................................................................38 7. Orthogonal Collocation - Finite Elements ..........................................................................38 7.1 Quadrature in Finite Element Methods ........................................................................40 7.2 C1 Collocation at Gauss Points ....................................................................................41 7.3 C0 Collocation at Lobatto Points ..................................................................................42 7.4 Convergence Rate and Efficiency ................................................................................45 Summary ...............................................................................................................................48 Notation .................................................................................................................................49 References ............................................................................................................................50 Appendix A – Galerkin Method with Flux Boundary Conditions.............................................57 Appendix B – Links to Available Software .............................................................................58

1. Abstract [1]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

2018 marks the 80th anniversary of Lanczos (1938) and last year was the 50th anniversary of Villadsen and Stewart (1967), both seminal works in development of the Orthogonal Collocation Method, also known as the Pseudospectral Method and Differential Quadrature Method. Due to these milestones, a retrospective on development of the method seems appropriate. The method’s development as a Method of Weighted Residuals (MWR), with numerical integration, is reviewed. It is placed into historical perspective and developments in disparate disciplines are reconciled. We show the equivalence of methods that were thought by many to be different. Guidance in navigating the minefield of terminology is also provided. In addition to reviewing the development of the method, we show some refinements which are not known to most practitioners. Chief amongst these are the correct treatment of boundary conditions involving fluxes. We demonstrate that a proper treatment of flux boundary conditions can sometimes reduce errors by several orders of magnitude.

2. Introduction Orthogonal Collocation (OC), synonymous with the Pseudospectral Method (PS), also goes by the name of Differential Quadrature Method (DQ). It is a method used for the approximate solution of differential equations. This method differs from finite difference methods because the solution is represented as a continuous or piecewise continuous function. For example, if one is solving for a variable y, the solution is approximated by: 𝑛

𝑦 ≅ 𝑦̃ = ∑ 𝑎𝑘 𝜓𝑘 (𝑥 )

(1)

𝑘=1

Where x represents one or more spatial coordinates, the a are adjustable parameters which could be time dependent, and the ψ are called trial functions or basis functions. The trial functions usually, but not always, satisfy the boundary conditions. The collocation procedure satisfies the differential equation exactly at a set of collocation points, xk. This criterion produces values for the adjustable parameters. By its various names (OC, PS and DQ), the Orthogonal Collocation method has been applied extensively in such diverse fields as: solid and structural mechanics, fluid mechanics, elasticity, computational chemistry and chemical physics, quantum mechanics, biomedical modeling, optimal control, acoustics, geophysics and seismic modeling, heat transfer, mass transfer, water resource modeling, petroleum reservoir simulation, chemical reaction engineering, meteorology and atmospheric modeling, and electromagnetism. The number of books and articles number more 100,000. For this reason, comprehensive coverage is not feasible. This paper’s primary focus is on development and formulation of the method rather than applications. It is further restricted to problems where a trial function representation is applied to one or more spatial coordinates, such as in a boundary value problem. The application of collocation methods to initial value problems is another large area of study [Hairer, et al.

[2]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

(1993), Hairer and Wanner (1996)]. Only nonperiodic problems in finite domains are considered, but these other cases are addressed in the many references provided. We distinguish between several different types of trial functions. In this article, all trial functions are polynomials. In early applications, the trial functions were often simple monomials, xk. In later applications the trial functions were orthogonal polynomials, usually Legendre or Chebyshev polynomials. We call this a modal approximation since the adjustable parameters are analogous to the modes in a Fourier series. The method is more intuitive and finitedifference-like if the adjustable parameters represent the approximate solution at the collocation points rather than polynomial coefficients, i.e. 𝑎𝑘 = 𝑦̃(𝑥𝑘 ). In this case it is called a nodal approximation and the trial solution is: 𝑛+1

𝑦̃ = ∑ 𝑦̃(𝑥𝑖 )ℓ𝑖 (𝑥)

(2)

𝑖=0

where the trial functions, ℓi(x), are Lagrange interpolating polynomials: 𝑛+1

ℓ𝑖 ( 𝑥 ) = ∏ 𝑗=0 𝑗≠𝑖

(𝑥 − 𝑥𝑗 ) (𝑥𝑖 − 𝑥𝑗 )

The choice of modal or nodal representations is largely a matter of convenience. In the absence of computer rounding errors all choices produce absolutely identical approximations. This important fact is obscured in many texts that emphasize the importance of orthogonal polynomial trial functions. The easiest way to treat multidimensional problems is by using combinations of the onedimensional trial functions. For example, Eq. (2) can be extended by the addition of a second, z, coordinate by: 𝑛𝑥 +1 𝑛𝑧+1

𝑦̃ = ∑ ∑ 𝑦̃(𝑥𝑖 , 𝑧𝑗 )ℓ𝑖𝑥 (𝑥) ℓ𝑗𝑧 (𝑧)

(3)

𝑖=0 𝑗=0

Trial functions composed of one dimensional trial functions like this are called tensor product trial functions. Finally, we distinguish between global trial functions where a single continuous polynomial approximates the solution and piecewise continuous trial functions called finite element or spectral element approximations. These trial functions are also called shape functions. In finite elements, polynomial trial functions are continuous within subdomains called elements but have limited continuity at the interface between elements. When the elements are distorted using a polynomial to describe their boundaries, they are called isoparametric elements. Lanczos (1938,1956) was an early adopter of modal approximations. Villadsen and Stewart (1967) were the first to use nodal approximations in this context. Neither one considered finite

[3]

element trial functions. In this article, we will concentrate on global approximations, but will show how the method extends to finite elements in a natural way. The distinguishing feature of orthogonal collocation is that the collocation points are chosen to be roots of orthogonal polynomials. The collocation points are usually the roots of Chebyshev polynomials of the 1st or 2nd kind, or the base points of Gaussian, Radau or Lobatto quadrature. Gaussian quadrature base points are the roots of Legendre polynomials. All these choices are the roots of Jacobi polynomials, which possess the orthogonality: 1

(𝛼,𝛽) (𝛼,𝛽) ∫ (1 − 𝑥 )𝛼 (1 + 𝑥)𝛽 𝑃𝑛 (𝑥)𝑃𝑚 (𝑥)𝑑𝑥 = 0 𝑓𝑜𝑟 𝑛 ≠ 𝑚

(4)

−1

where α = β = -½, 0, +½, +1 for the Chebyshev 1st kind, Legendre (Gauss), Chebyshev 2nd kind polynomials, and for those whose roots are Lobatto quadrature base points, respectively. For increasing n, these polynomials are alternating symmetric and antisymmetric about the centerline. Radau points are asymmetric with α = 0, β = 1 or α = 1, β = 0. The right half of the symmetric polynomials along with their extrema are compared in Fig. 1 for n = 8. What should be noted in Fig. 1 is that as α increases the roots are shifted away from the boundary and there is a greater variation in the extrema. The endpoint for the Jacobi (1,1) polynomial is off scale at 9. In all cases the points are more tightly clustered near the boundaries, with spacing 4 proportional to 1/n2. The Chebyshev 1st Chebyshev polynomials of the Legendre 3 Chebyshev 2nd 1st kind are the favorite ones Jacobi (1,1) Chebyshev 1st, extrema for interpolation, since in Legendre, extrema 2 Chebyshev 2nd, extrema general they minimize the Jacobi (1,1), extrema maximum error [Lanczos 1 (1956) p. 245, Hildebrand 0 (1987) p. 469]. Note that in Fig. 1 the Chebyshev -1 polynomials of the 2nd kind are normalized like other Jacobi -2 0 0.2 0.4 0.6 0.8 polynomials rather than the x Fig. 1 Polynomials and extrema for, n = 8,  =  = -½, 0, +½, 1 traditional way. Pn

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Orthogonal polynomials are closely tied to the theory of accurate numerical integration methods. For a general reference to orthogonal polynomials and quadrature, the reader is directed to Hildebrand (1987) and Krylov (1962). The selection of these points is related to numerical integration of the form: 1

𝑛

∫ 𝑓(𝑥 )𝑑𝑥 ≅ 𝑊0 𝑓(0) + ∑ 𝑊𝑖 𝑓(𝑥𝑖 ) + 𝑊𝑛+1 𝑓(1) 0

(5)

𝑖=1

Fig. 2 illustrates the point locations and quadrature weights for n = 6. All but Radau points are symmetric about the centerline. For Gaussian quadrature, endpoint weights are not used, W0 = Wn+1 = 0, while the Radau-Left quadrature shown has W0 > 0 and Wn+1 = 0. There is a mirror [4]

image Radau-Right formula where the weights and points are switched around. For Lobatto and Chebyshev (of the 2nd kind) points W0 = Wn+1 > 0. Some claim the approximation of boundary conditions is improved by nonzero weights at the boundaries; however, we will show that the exact opposite is sometimes true.

0.3 0.25 0.2

w

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

0.15 0.1

Lobatto Gauss Chebyshev Radau Left

0.05

The various quadrature formulas can calculate exact integrals when the 0 0 0.2 0.4 0.8 1 x 0.6 integrand of Eq. (5) is a polynomial of Fig. 2 Quadrature points and weights degree 2n+1, 2n, 2n-1 and n+1 for Lobatto, Radau, Gauss and Chebyshev points, respectively. For Chebyshev points, the Clenshaw-Curtis (1960) quadrature is formally less accurate, which is certainly not obvious from examination of Fig. 2. Some claim it is just as accurate in many applications [Trefethen (2008)]. Clenshaw-Curtis quadrature is distinct from Chebyshev-Gauss quadrature which has accuracy O(2n) when the radical 1/√(1 − 𝑥 )2 is included in the integrand of Eq. (5). With this introduction, we can proceed with the main part of the paper. We begin with Section 3, a historical review of the development of collocation and related methods called Methods of Weighted Residuals (MWR). Then two linear one dimensional example problems are considered with global trial functions: a boundary value problem in section 4, the catalyst pellet problem, and a parabolic problem in section 5, the falling liquid film problem. These problems are intentionally simple to enable a clear description of the method and some potential pitfalls if it is not correctly applied. Chief amongst the pitfalls is the correct treatment of flux boundary conditions. The problems are solved with global collocation, Galerkin and moments (or tau) methods. Although the examples are simple, the methods presented carry over to nonlinear multidimensional problems, section 6. In section 7 we show how the global methods extend naturally to finite elements using the catalyst pellet problem as an example. One goal of this article is to present a unified and consistent description of methods, and to show relationships between methods, including the equivalence of methods which are normally considered to be different. Also, this subject area can be a virtual minefield of conflicting terminology due to different naming conventions and so forth. To aid the reader, we have endeavored to include alternate names for various quantities and italicize their definition when first presented. Designation of the roots or collocation points is one of the more confusing terminology issues one encounters. For examples, the names for the roots considered are:

[5]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

− Chebyshev polynomials of the 1st kind also called: Chebyshev-Gauss or abbreviated CG points or simply Chebyshev points or a roots grid − Chebyshev polynomials of the 2nd kind also called: Chebyshev Extrema or Chebyshev-Gauss-Lobatto or abbreviated CGL points or simply Gauss-Lobatto or Chebyshev points − Legendre polynomials also called: Gaussian quadrature base points or abscissa or Legendre-Gauss or abbreviated LG points or simply Legendre or Gauss points − Jacobi Polynomials (α = 0, β = 1 or α = 1, β = 0) also called: Radau points, Radau quadrature base points or abscissa, Legendre-Gauss-Radau or abbreviated LGR points − Jacobi Polynomials (α = β = 1) also called: Legendre Extrema points or Jacobi points or Lobatto quadrature base points or abscissa or Legendre-Gauss-Lobatto abbreviated LGL points or simply Gauss-Lobatto or Lobatto points Beware of the terms Chebyshev and Gauss-Lobatto since both are thrown around loosely for two different sets of points. In this article, we will not consider collocation at the roots of Chebyshev polynomials of the 1st kind, so Chebyshev with no clarification designates Chebyshev polynomials of the 2nd kind. The other points in Fig. 2 will be called by the names found in most standard numerical analysis texts, i.e. Gauss, Lobatto or Radau (left or right) points [Hildebrand (1987)].

3. Historical Perspective The collocation method is one of a family of methods collectively called Methods of Weighted Residuals (MWR) [Finlayson (1972)]. These methods approximate the solution with a linear combination of trial functions like Eq. (1) or (2). For example, suppose we seek a solution of the equation: 𝑦𝑥𝑥 + 𝑔(𝑥, 𝑦, 𝑦𝑥 ) = 0 (6) for 0 < x < 1, with y(0) = y0 and y(1) = y1. The residual is formed by substituting the approximate solution into the differential equation: 𝑛

𝑅(𝑥, 𝒂) = ∑ 𝑎𝑘 𝜓𝑘𝑥𝑥 + 𝑔(𝑥, 𝑦̃, 𝑦̃𝑥 )

(7)

𝑘=1

If the residual is zero for all x, an exact solution has been achieved. In general, this will not be possible, so the parameters, a, are selected so the residual will be small in some sense. If the residual is small for all x, Eq. (1) will be a good approximation to the true solution. MWR achieves this goal by forcing the function to zero in a weighted average sense: 1

∫ 𝑤𝑘 (𝑥)𝑅(𝑥, 𝒂)𝑑𝑥 = 0 for 𝑘 = 1, … , 𝑛

(8)

0

The wk are called weight or test functions. Different choices for the weight functions lead to fundamentally different forms of MWR. Finlayson and Scriven (1966) review the early development of these basic procedures. The MWR were created in the early 20th century to generalize the Ritz (1908) variational procedure, which has limited applicability. The name [6]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

MWR came about from comments by Courant (1924) at the First International Congress of Applied Mechanics. Crandall (1956, p. 147) was first to illustrated the methods in a unified way and following Courant stated the basic principal is that “weighted averages of the residuals” should vanish. These comments led to the name Methods of Weighted Residuals. They were called error distribution principals by Collatz (1961) and Ames (1965). The early literature sometimes called them direct methods. They are also called projection methods, orthogonalization methods or orthogonal projection methods. The popular name spectral methods if of unknown origin but seems to have originated about the time of Orzag (1972) and was initially applied only to global trial functions. Different choices for the n weight functions produce different methods, we will consider the following methods: 1. 𝑤𝑘 = 𝜓𝑘 (𝑥), Galerkin method 2. 𝑤𝑘 = 𝑥 (𝑘−1) or 𝑃𝑘 (𝑥 ), method of moments 𝜕𝑅

3. 𝑤𝑘 = 𝜕𝑎 , least squares 𝑘

4. 𝑤𝑘 = 𝛿(𝑥 − 𝑥𝑘 ), collocation method There are a few other methods which are not listed above. Several texts include the tau method of Lanczos (1938, 1956) as a MWR, which is discussed in section 3.6. 3.1 Galerkin Method. The Galerkin method weights the residual by the trial functions. For problems which can be treated with a variational procedure, the Galerkin method is equivalent. It was described by Galerkin (1915) in a study of elastic equilibrium and the stability of rods and plates. An earlier paper by Bubnov (1913) had many similarities, but Galerkin was the first to describe the method without any reliance on a variational principal. Due to these relationships it is called the Bubnov-Galerkin method in Russian literature and often the Raleigh-Ritz-Galerkin method in western literature. The Galerkin method is the most heavily used MWR and can be considered the gold standard for these methods. 3.2 Method of Moments. The basic ideas behind the method of moments were first described by N.M. Krylov at a session of the All Ukraine Academy of Sciences in 1926 [Krylov (1926,1961), Lucka and Lucka (1992)]. The ideas behind the original method are far more general than weighting by simple monomials. Originally, given some linear differential operator, L, the weight functions were chosen as 𝑤𝑘 = 𝐿 𝜓𝑘 . The method was further developed by Kravchuk (1926,1932). Shuleshko (1959) provided additional elaboration. Kravchuk1 spent much of his short and tragic career studying the various conditions on the linear operator and

1

Mykhailo Kravchuk (1892-1942) was the foremost Ukrainian mathematician of the 20th century. He entered the University at Kyiv in 1910 studying mathematics and physics, teaching after 1917. He endured many hardships following the Bolshevik revolution and the subsequent loss of Ukranian independence in 1922, but despite the obstacles, rose to full professor in 1925. In 1938, after years of public service as an educator and academician, he fell victim to Stalin’s purges and was sent to a Gulag death camp in Siberia. He succumbed to the harsh treatment in 1942.

[7]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

trial functions and resulting convergence rates. The generality afforded by this choice of weight functions includes special cases of the Galerkin and least squares methods, while it also admits eigenfunctions of other linear operators, e.g. Legendre polynomials. Of course, if the residual is orthogonal to Legendre polynomials it is also orthogonal to all the monomials, x(k-1), up to the same degree. The idea of using monomials is apparently due to Yamada (1947) as described by Fujita (1953) in his study of diffusion in gels. Weighting by monomials can lead to ill conditioned approximations and some have rejected it for that reason [Boyd (2000) p. 62]. In the western literature today, the method of moments is generally considered to be weighting by monomials, but in the current paper we include equivalent methods, such as weighting by Legendre polynomials. A basic idea behind the Galerkin and moments methods are that if the weight functions are members of a complete set, then as n ⟶ ∞ the residual must converge to zero. Of course, only a finite number of trial functions are used in practice, but the property of completeness guarantees that the sequence will converge. 3.3 Least Squares. Least squares is a very old idea for approximating functions. Here, it is used to force the residual to approximate zero. In the context of solving differential equations, Finlayson and Scriven (1966) attributed its first use to Picone (1928). However, we note the method was discussed in an earlier paper by Krylov (1926). The least squares method works well for simple differential equations. However for equations with nonlinearity and many terms it is more difficult to implement. For this reason, its popularity has waned as more complex problems have been tackled. 3.4 Collocation. The collocation method dates from the 1930’s. Early work on the method is due to Slater (1934), Kantorovich (1934), Frazer, et al. (1937) and Lanczos (1938). Kantorovich called it the interpolation method, Lanczos called it the method of selected points, while Frazer, et al. applied the name collocation, citing its definition in the Oxford dictionary. Using all three names we can say the method interpolates the residual to zero at selected collocation points. In the context of MWR, the weighting functions can be viewed as Dirac delta functions. Since integration is not required, it is simpler to apply. Simplicity is its most attractive feature. Slater applied the method to study energy bands in metal. A translation of the Kantorovich article is not available, but it is discussed in other articles. Frazer, et al. compared equally spaced collocation to the Galerkin and least squares methods on several problems. Frazer, et al. argued that the method would converge whenever the Galerkin or least squares methods converge. However, Karpilovskaya (1963) points out flaws in their argument. Since collocation [8]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

is basically interpolation of the residual, the Runge (1901) phenomenon can occur. A good example of problems associated with equally spaced collocation can be found in Bert and Malik (1996). Frazer, et al. were apparently aware of the Runge phenomenon, but did not observe it for their problems. Lanczos used Chebyshev polynomial trial functions and collocated at the roots of Chebyshev polynomials of the 2nd kind. Clenshaw and Norton (1963) collocated at the same points and discuss several implementation details for nonlinear problems. Wright (1964) collocated at the roots of Chebyshev polynomials of the 1st kind, arguing that this choice would minimize the maximum deviations of the residual from zero. However, a uniform distribution of the residual does not produce a uniform error in the solution. The use of Chebyshev roots was an important advance, since the Runge phenomenon does not occur. Villadsen and Stewart (1967) made a further improvement by selecting the roots of Jacobi polynomials that are the base points of accurate numerical integration schemes, e.g. Gaussian or Lobatto quadrature. They chose the name Orthogonal Collocation for their variation on the method. They explained that collocation at Lobatto points was equivalent to numerical integration of the Galerkin method. An accurate (sometimes exact) approximation of the Galerkin method is achieved. They also improved the method by formulating it with nodal values rather than polynomial coefficients. This modification produced a more intuitive, finite difference like method and facilitated automation of the method for general applications. Despite the title of the article, they demonstrated the method for: a boundary value problem, the parabolic transient catalyst pellet problem, eigenvalues for the Graetz problem and the elliptic 2-dimensional Poisson equation. All but hyperbolic equations were covered. In the 1970’s, the development of collocation methods occurred in three threads: (1) Orthogonal Collocation (OC), (2) Pseudospectral (PS) and (3) Differential Quadrature (DQ). The OC thread started with the Villadsen and Stewart article, with further development reported in Villadsen (1970), Finlayson (1972) and Villadsen and Michelsen (1978). These references describe collocation at Gauss, Lobatto and Radau points. Collocation at Gauss points was shown to be equivalent to numerical integration of the moments method. In addition to standard cartesian coordinates, they describe methods for axisymmetric problems in planar, cylindrical and spherical coordinates. Nodal differentiation matrices are used exclusively. These developments led to a flurry of activity applying the method to different problems [see references in Michelsen and Villadsen (1972,1981), Finlayson (1974)]. Many of these early papers demonstrated very favorable comparisons of the method to finite differences. A later text describing the method is Finlayson (2003). The Pseudospectral thread began with the work of Orzag (1972), with further developments by Gottlieb and Orzag (1977). The name pseudospectral is normally synonymous with collocation, but occasionally it is used to describe other approximations to exact integration in MWR. Orzag applied collocation to a periodic problem using trigonometric functions. He also [9]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

considered a nonperiodic linear first order hyperbolic problem with collocation at the roots of Chebyshev polynomials of the 2nd kind. An important contribution of this work was the use of fast Fourier transforms (FFT) to perform the calculations. For the hyperbolic example, he demonstrated that collocation approximated the Galerkin method accurately. For the nonperiodic problem Chebyshev trial functions were used. Nodal approximations were not considered. It is not known whether Orzag was aware of the earlier work of Villadsen and Stewart, but it was not referenced. FFT was the spark that ignited work in the PS thread. Subsequent work investigated various methods of solution using FFT. A strong mathematical foundation for the method was laid and the solution of periodic problems was highly developed. Most of the early work on nonperiodic problems employed collocation at Chebyshev points with a modal approach using Chebyshev trial functions. Nodal trial functions and collocation at Lobatto points were adopted later on, but Villadsen and Stewart are rarely credited for these innovations. Some excellent texts describe the method, solution algorithms and applications [Canuto, et al. (1988, 2006, 2007), Bernardi and Maday (1997), Trefethen (2000), Boyd (2000) , Peyret (2002), Shen, et al. (2011)]. The Differential Quadrature Method originated with Bellman and Casti (1971), Bellman, et al. (1972), and Bellman (1973, p. 244). The first was a short paper advancing the idea of using a nodal differentiation matrix in the solution of differential equations. The second paper developed the method further and demonstrated it for both first and second order equations in one dimension. However, the paper provided few details regarding boundary condition treatment. Bellman, et al. (1972) developed the differentiation matrices for collocation at Gauss points. The use of a nodal differentiation matrix was not a new idea. Nielsen (1956) presents formulas for them with arbitrary nodal locations. Nodal collocation at Gauss points was already well established in the OC literature. Nevertheless, the method was adopted for many engineering applications. Bert and Malik (1996) and Bellomo (1997) give excellent reviews of developments in DQ. The method was initially developed independently of work in the OC and PS communities. Apparently, it was not recognized as a form of collocation until the late 1980’s [Quan and Chang, 1989)]. A reference text on the method is Shu (2000). Ideas have been exchanged between the different threads of development; however, their source are not always referenced. Chebyshev points have historically been popular in PS applications and are used in DQ as well. Villadsen and Stewart were the first to advocate collocation at Gauss and Lobatto points. Lobatto points, also called Legendre-Gauss-Lobatto or LGL points, have become very popular in the PS literature. Gauss points, also called Legendre-Gauss or LG points, are popular in OC applications, but not as popular in the PS and DQ communities. The nodal formulation, pioneered by Villadsen and Stewart, is popular with “the engineers” due to its similarity to finite differences. Nodal approximations are used in OC, PS and DQ applications. Symmetric problems in planar, cylindrical and spherical coordinates are common in OC applications, while problems with periodic solutions are covered almost exclusively in the PS literature. [10]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

3.5 Other Weight Functions. Some authors take a more general view and refer to virtually all MWR weightings as Galerkin or variational methods. Alternative weighting schemes are sometimes called Petrov Galerkin methods. For Chebyshev methods the weights in Eq. (8) usually include the radical 1/√1 − 𝑥 2. This term is necessary to exploit the orthogonality of the polynomials and to achieve accurate integration with Gauss-Chebyshev quadrature. In naming methods, most authors do not distinguish between methods which include the radical. It would be more appropriate to add Chebyshev to the name, e.g. Chebyshev-Galerkin method. The naming of methods with alternative weight functions is yet another area where terminology can be confusing. In the current article, methods which include the radical will not be used. 3.6 Tau Method. The tau method was originated by Lanczos (1938, 1956) as a method for approximating functions which can be described by differential equations. The basic idea behind the method is that rather than find an approximate solution to the differential equation, why not find an exact solution to an approximate differential equation. For simple problems, he represented the approximate differential equation as the original one with the addition of a small perturbation term of the form τPm(x), where τ is a constant and Pm(x) is a high order orthogonal polynomial, i.e. m ≈ n. Depending on the complexity of the problem, several terms may be required. These perturbation terms are the residual as defined by Eq. (7), so in general Lanczos’ idea was to set: 𝑅(𝑥, 𝒂) = 𝜏0 𝑃𝑚 (𝑥 ) + 𝜏1 𝑃𝑚+1 (𝑥) + ⋯

(9)

where the number of terms required depends on the nature of the differential equation. Lanczos describes a recursive procedure for determining the values of the τ parameters and the coefficients in the approximation, i.e. a in Eq. (1). The tau method was further developed in a series of articles by Ortiz and coworkers [e.g. Ortiz (1975)]. Several solution methods are described which do not involve integration or orthogonalization as required by most MWR. Later papers [El-Daou, et al. (1993), El-Daou and Ortiz (1997)] explore the equivalence between the tau method and other methods. The relationship between the tau method and the integrated MWR is obvious. Since most MWR make the residual orthogonal to the first few members of a complete set, the residual will be proportional to the first neglected term and possibly higher terms. If tau methods are defined to include all methods with a residual term like Eq. (9), then virtually all MWR would be tau methods. This idea seems silly. Instead, we view the tau method as an alternate solution procedure for MWR which does not involve integration. The lack of a need for integration is the method’s distinguishing feature, but it also causes its applicability to be somewhat limited. The tau method as presented in the spectral literature [e.g. Gottlieb and Orzag (1977, p. 11), Canuto, et al. (1988, p. 10, 2006, p. 21), Boyd (2000) p. 473] bears little resemblance to the method described by Lanczos (1956). Unlike the algorithms of Lanczos and Ortiz and coworkers, the spectral literature describes a method which requires integration. The Legendre-tau method makes the residual orthogonal to the first few Legendre polynomials, so [11]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

it is equivalent to the method of moments. The Chebyshev-tau method makes the residual orthogonal to Chebyshev polynomials with the additional weight factor of 1/√1 − 𝑥 2 [Johnson (1996)], i.e. a Chebyshev-moments method. In the spectral literature, the boundary conditions are usually satisfied using side conditions, creating a messy matrix structure. Whether the boundary conditions are satisfied by the trial functions or whether they are imposed by side conditions makes no difference in the final result. Section 4.2 describes a nodal form of the method with a symmetric matrix structure. The tau method offers an interesting view of the residual which we will revisit when discussing the examples. However, as practiced in the spectral literature it is identical to the method of moments. In this article, we will use the traditional name method of moments, but the reader should be aware that the tau method is equivalent. The residual for all the MWR with polynomial trial functions can be represented by Eq. (9), so the distribution of residual errors can be shifted by the choice of weighting function, cf. the polynomials in Fig. 1. 3.7 Boundary Conditions. Descriptions of MWR distinguish between three different approaches: (1) interior methods, (2) boundary methods and (3) mixed methods. An interior method employs trial functions which obey the boundary conditions, so the parameters are chosen to approximate the differential equation. In boundary methods, the trial functions obey the differential equation and the parameters are chosen to approximate the boundary conditions. In mixed methods, the parameters are adjusted to approximate both the differential equation and boundary conditions. The foregoing sections have primarily addressed approximations to the differential equation, which is usually the biggest challenge. When the boundary conditions are approximated, the residual of the boundary conditions may be required to meet their own independent weighted residual criteria or they may be enforced with criteria based on a weighted combination of interior and boundary residuals. For the moments (or tau) method, all boundary conditions must be satisfied exactly. They may be satisfied through the selection of the trial functions or by using side conditions [Boyd (2000), pp 111-115]. For the Galerkin method, the boundary conditions can be divided into two types – essential boundary conditions and natural boundary conditions [Finlayson (1972)]. For second order equations, Dirichlet conditions or other conditions on the value of the solution variable are essential and must be satisfied exactly. Conditions of the 2nd or 3rd type (Neuman or Robin) or others involving derivatives or fluxes are natural boundary conditions. Natural conditions can be handled in one of two ways. They can be satisfied exactly, or they can be approximated using a combined weighting of interior and boundary residuals. The natural treatment is also called a weak formulation of the boundary conditions. Exact satisfaction of the boundary condition is called a strong treatment or boundary collocation since, like collocation in the interior, the residual of the boundary condition is zero. [12]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Although not a fundamental requirement, collocation implementations normally use a strong treatment of boundary conditions [e.g. Finlayson (1972), Villadsen and Michelsen (1978), Bert and Malik (1996), Bellomo (1997), Trefethen (2000), Boyd (2000), Peyret (2002)]. Note that these references include the OC, PS and DQ threads of development. There have been a few references that describe an alternative weak formulation of flux boundary conditions [Canuto, et al. (1988, 2006), Funaro (1992), Shen, et al. (2011)], but these texts also cover the strong treatment, and none claim a significant benefit to a weak formulation. Consequently, most applications use a strong treatment or boundary collocation for all boundary conditions. Although not widely publicized, early applications of OC uncovered a difficulty with flux boundary conditions [Ferguson (1971), Finlayson (1972), Elnashaie and Cresswell (1973)]. Ferguson (1971) described an integral procedure which circumvented the problem, but its applicability is limited. It was also suggested that, when the method approximates a Galerkin method, flux boundary conditions should be formulated as natural/weak boundary conditions [Young (1977)]. The problem with flux conditions was absent when quadrature weights were zero at the boundary, i.e. Gauss points, which flies in the face of the commonly held notion that nonzero boundary weights somehow aid in the approximation of boundary conditions. The examples in sections 4 and 5, provide a close look at this issue. 3.8 Finite Elements. In the 1970’s there was an explosion of interest in finite element methods [Zienkiewicz (1971), Strang and Fix (1973), Hughes (1987)]. Its success in structural mechanics lead to an interest in other application areas. Lanczos and Villadsen and Stewart considered only global approximations. However, Villadsen and Stewart’s idea of collocating at quadrature points was soon used to extend OC for use with finite elements. First and most popular is collocation at Gauss points with early articles by DeBoor and Swartz (1973), Douglas and Dupont (1973), and Cary and Finlayson (1975). Finite element collocation at Gauss points is a method with continuous derivatives (C1 continuity) at element interfaces. Several texts are available which describe the method [Davis (1984), Lapidus and Pinder (1999), Finlayson (2003)]. Collocation at Lobatto points extends to a finite element method with simple continuity, C0 continuity. The earliest developments are due to Diaz (1975), Dunn and Wheeler (1976), and Wheeler (1977) who describe a method called the Hybrid-Collocation-Galerkin method. Independently, Gray (1977), Young (1977), Hennart (1982), and later Leyk (1986,1997) developed a method Young called the Lobatto-Galerkin method. Heretofore, these two methods were considered distinctly different, but in section 7 we show they are equivalent. This method was later popularized as the hp Spectral Element method [Maday and Patera (1989), Vosse and Minev (1996), Canuto, et al. (2007), Karniadakis and Sherwood (2013)].

[13]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Since the discussion of finite element collocation methods gets into the large area of reduced quadrature (i.e. less than exact integration) in finite element methods, further discussion is deferred until section 7, after the background information is presented in sections 4 and 5.

4. Boundary Value Problem - Catalyst Pellet Problem Consider reaction and diffusion in a porous catalytic slab. This differential equation is one of the most ubiquitous ones in engineering and physics. It is analogous to heat conduction in a slab with a source dependent on temperature and position. If the source, r, is constant, the equations describe laminar flow between parallel plates. It also describes the position of an elastic string, where r characterizes the load. The governing equations are: 𝑑2𝑦 (10) + 𝑟(𝑥, 𝑦) = 0 𝑑𝑥 2 where 𝑟(𝑥, 𝑦) = 4𝜑2 𝑟̂ (𝑥, 𝑦) and y is the fractional conversion, 0 for no reaction and 1 for complete reaction. φ is called the Thiele modulus. The rate dependence on x permits a reactivity that depends on position. φ is defined such that the average value of r̂ (x,0) is 1. The boundary conditions are either first kind, Dirichlet, boundary conditions: 𝑦(0) = 𝑦(1) = 0 (10a) or third kind, Robin, boundary conditions: 𝑑𝑦 𝑑𝑦 | = 2𝐵𝑖0 𝑦(0) and − | = 2𝐵𝑖1 𝑦(1) 𝑑𝑥 𝑥=0 𝑑𝑥 𝑥=1

(10b)

In Eq. (10b), the Biot numbers, Bi, account for an external transfer resistance. As Bi tends to infinity, Eq. (10b) reduces to Eq. (10a). This problem is often symmetric about the centerline, but here we allow for a nonsymmetric profile. Symmetric problems in various geometries can be efficiently treated using polynomials in x2. The Thiele modulus and Biot numbers are defined using half the thickness of the slab, so the factors of 4 and 2 are required when the full slab is considered. For a simple kth order reaction with the reactivity independent of position, the source term is: 𝑟̂ (𝑥, 𝑦) = (1 − 𝑦)𝑘

(11)

For a first order reaction, k = 1, and a single value for Bi, the analytical solution is: cosh[𝜑(2𝑥 − 1)] 𝑦 =1− cosh(𝜑) + (𝜑/𝐵𝑖)sinh(𝜑)

(12)

The quantity of interest from the solution is called the effectiveness factor, η, which gives the overall rate of reaction relative to that with no diffusional resistance: 1

𝜂=

∫0 𝑟(𝑥, 𝑦)𝑑𝑥 1

∫0 𝑟(𝑥, 0)𝑑𝑥

1

= ∫ 𝑟̂ (𝑥, 𝑦)𝑑𝑥

(13)

0

The effectiveness factor is basically a normalized boundary flux, since the divergence theorem in one dimension gives: [14]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

−

1 𝑑𝑦 1 | = 4𝜑2 ∫ 𝑟̂ (𝑥, 𝑦)𝑑𝑥 = 4𝜑2 𝜂 𝑑𝑥 0 0

(14)

Eq. (14) is given the generic name average energy equation, since in a heat transfer setting it is an overall energy balance and the right-hand side is the average energy generated. Using the analytical solution, Eq. (12), this normalized flux is: 1 𝜂= (15) φ[coth(𝜑) + 𝜑/𝐵𝑖] η is approximately one for φ < ½ and becomes asymptotic when φ > 2. The asymptotic state corresponds to a condition where all the reactants are consumed near the boundary and the conversion, y, approaches one at the center. 4.1 Orthogonal Collocation Method Early developments of the Orthogonal Collocation method state that the trial functions are orthogonal polynomials, but then use monomials and through transformations formulate problems in terms of nodal values. Although the approximation turns out the same, we prefer to dispense with the transformations and develop the method directly. To solve the problem with a nodal formulation, we start with the trial solution Eq. (2). The interpolation points in Eq. (2) include the endpoints, x0 = 0 and xn+1 = 1, while the interior points are the roots of an orthogonal polynomial. Most authors seem almost paranoid about these boundary points when used with Gauss points in the interior, because the associated quadrature weights are zero, see Eq. (5) and Fig. 2. They state that nonzero weights are needed on the boundary to help meet the boundary conditions. There is never any justification for these statements and, in fact, we will show that in some instances nonzero boundary weights are a detriment rather than an asset for meeting the boundary conditions. The residual, Eq. (7), is formed by substitution of the approximate solution into the equation: 𝑛+1

𝑛+1

𝑑 2 ℓ𝑖 ∑ 𝑦(𝑥𝑖 ) 2 + 𝑟 (𝑥, ∑ ℓ𝑖 (𝑥)𝑦(𝑥𝑖 )) = 𝑅(𝑥, 𝒚) 𝑑𝑥 𝑖=0

(16)

𝑖=0

Where y is the vector of nodal values y(xi). With the collocation method, the residual is set to zero at the interior collocation/interpolation points, xj, j = 1,…,n. Since ℓi(xj) = δij (the dirac delta function), the resulting equation simplifies to: 𝑛+1

∑ 𝑦 (𝑥 𝑖 ) 𝑖=0

𝑑 2 ℓ𝑖 | + 𝑟 (𝑥𝑗 , 𝑦(𝑥𝑗 )) = 0 𝑑𝑥 2 𝑥

(17)

𝑗

By defining B as indicated below and letting 𝑦𝑖 = 𝑦(𝑥𝑖 ) the equation simplifies to: 𝑛+1

∑ 𝐵𝑗𝑖 𝑦𝑖 + 𝑟(𝑥𝑗 , 𝑦𝑗 ) = 0

(18)

𝑖=0

The boundary conditions provide two additional conditions, either first kind, Dirichlet, boundary conditions: [15]

𝑦0 = 𝑦𝑛+1 = 0 or third kind, Robin, boundary conditions: 𝑛+1

𝑛+1

∑ 𝐴0,𝑖 𝑦𝑖 = 2𝐵𝑖0 𝑦0 and − ∑ 𝐴𝑛+1,𝑖 𝑦𝑖 = 2𝐵𝑖1 𝑦𝑛+1 𝑖=0

(19)

𝑖=0

where: 𝑛+1

𝑑ℓ𝑖 𝑑 2 ℓ𝑖 | = ∑ 𝐴𝑗𝑘 𝐴𝑘𝑖 𝐴𝑗𝑖 = | and 𝐵𝑗𝑖 = 𝑑𝑥 𝑥𝑗 𝑑𝑥 2 𝑥 𝑗

(20)

𝑘=0

For Dirichlet conditions, Eq. (10a), the boundary values (whether zero or finite) are simply substituted into Eq. (18). For the third kind conditions, Eq. (10b), most texts and articles recommend that the boundary conditions be satisfied exactly as in Eq. (19), i.e. boundary collocation or a strong treatment [Finlayson (1972), p. 101; Villadsen and Michelsen (1978), p. 137; Bert and Malik (1996); Bellomo (1997); Trefethen (2000), p. 137; Boyd (2000), p. 111, Peyret (2002), p. 59]. Below, boundary derivatives or fluxes are discussed further and a method superior to Eq. (19) is demonstrated. To solve the problem, we need only the collocation points, i.e. the roots of the orthogonal polynomial, x, the derivatives of the Lagrange interpolating polynomials, A and B, and the quadrature weights, W, for approximating integrals, e.g. Eq. (13). A and B are called differentiation matrices. The original methods for computing these quantities was crude, but improvements were quickly developed [Michelsen and Villadsen (1972)]. The roots can be calculated iteratively, while the weights and differentiation matrices can be calculated from the polynomial derivatives evaluated at the roots. Alternatively, the roots and quadrature weights can be calculated from the eigenvalues and eigenvectors of a symmetric tridiagonal matrix formed from the polynomial recurrence relationships [Golub and Welch (1969)]. Computation of the differentiation matrices with arbitrary point locations is described by Nielsen (1956). Many of the texts cited discuss procedures for performing these calculations [e.g. Villadsen and Michelsen (1978)]. For software availability in a variety of computer languages see Appendix B. 4.1.1 Linear Source, Dirichlet B.C. Fig. 3 shows solutions of Eqs. (10) and (10a) for a first order reaction and φ = 5, i.e. Eq. (11) with k = 1. Approximate solutions are shown with n = 4 for collocation at Lobatto, Gauss and Chebyshev points. With this relatively high reaction rate most of the reaction occurs near the boundary. The fifth order polynomial can only approximate the sharp profile by oscillating about the exact

1

0.8

0.6 n=4, Lobatto n=4, Gauss n=4, Chebyshev Exact

y

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

0.4

0.2

0

0

0.2

0.4

x

0.6

0.8

1

Fig. 3 Solutions of Eq. (9), Dirichlet, n = 4,  = 5, k = 1

[16]

solution. Since this problem is symmetric about x = 0.5, the coefficient of the fifth order term is zero.

4 6 8

12 Lobatto n = 6 Chebyshev n = 6 Gauss n = 6

10 8

Residual

n

Table 1 Coefficients of Residuals Gauss Chebyshev Lobatto (0.0) (0.5,0.5) (1,1) 𝑃𝑛 (𝑥 ) 𝑃𝑛 𝑃𝑛 33.6927 16.4722 9.0058 5.6737 2.8395 1.5242 0.5711 0.2892 0.1544

6 4 2

Fig. 4 shows the residual for the same 0 problem and n = 6. The nature of the -2 residual was discussed briefly in section 3.6. -4 Due to the symmetry of the problem, each 0 0.2 0.4 0.6 0.8 1 x th Fig. 4 Residual for Eq. (9), n = 6, Dirichlet  = 5, k = 1 of the residuals is an n order Jacobi polynomial with α = β = 0, +½, +1 for Gauss, Chebyshev and Lobatto points, respectively. The residuals follow Eq. (9), but each choice of points makes the residuals orthogonal to a different Jacobi polynomial, see Fig. 1, so the residual errors are distributed differently. The values of the proportionality constant, τ, given in Table 1 decrease at a rapid rate as n is increased. 100

Fig. 5 shows the L2 error norms versus n for three choices of points. A plot of the L1 error norm looks very similar. All methods show the typical exponential convergence with Lobatto points giving slightly better results. The error with Chebyshev and Gauss points averages 1.3 and 1.9 times that with Lobatto points, which is relatively small since the error drops by almost an order of magnitude with each increment of n. The largest residual values occur at the boundary with Lobatto points, but, nevertheless, this choice produces the smallest errors in the solution.

10-2

10

L2 Error

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

-4

10

-6

10

-8

10

-10

10

-12

Lobatto Gauss Chebyshev

2

4

6

8

n

10

12

14

Fig. 5 L2 Error Norm, Dirichlet,  = 5, k = 1

Using Eq. (15), the normalized boundary flux, η = 0.199983. Numerical quadrature can be used to calculate this value from the numerical solutions using Eq. (13): 𝑛+1

𝜂 = ∑ 𝑊𝑖  𝑟̂ (𝑥𝑖 , 𝑦𝑖 )

(21)

𝑖=0

Using Eq. (21) for the cases in Fig. 3, the calculated fluxes are in error by 0.8, -3.1 and -4.3 percent for Lobatto, Chebyshev and Gauss points, respectively. The same quantity can be calculated using Eq. (14), with the derivatives approximated by: [17]

𝑛+1

𝑑𝑦 | = ∑ 𝐴0𝑖  𝑦𝑖 𝑑𝑥 𝑥=0

𝑛+1

and

𝑖=0

𝑑𝑦 | = ∑ 𝐴𝑛+1,𝑖  𝑦𝑖 𝑑𝑥 𝑥=1

(22)

𝑖=0

Using Eq. (22), the derivatives of the solution at the boundaries are in error by -14.3, -10.2 and -4.3 percent for Lobatto, Chebyshev and Gauss points. These errors are much larger and the relative accuracy of the methods is a complete reversal of that found with Eq. (21). Only Gauss points give the same result with either method of flux calculation. Shortly, we will explain why this is so. For the other two methods, integration gives a far more accurate result. Fig. 6 shows the flux errors for increasing n. Relative to Fig. 5, this graph shows a much greater difference between the calculations. For n < 4, the error with all three methods is relatively large, while for n > 14, some of the results are affected by rounding errors. Engineering accuracy is obtained with 4 to 8 points depending on the method and accuracy required. In this range, the errors with Gauss and Chebyshev points are similar, while Lobatto points give somewhat greater accuracy.

10

Flux Error

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

0

10

-2

10

-4

10

-6

10

-8

Lobatto Gauss Chebyshev Lobatto Der. Chebyshev Der.

10-10 10

-12

10

-14

2

4

6

8

n

10

12

14

Fig. 6 Flux Error, Dirichlet  = 5, k = 1

The understanding of convergence properties of the method has improved greatly since it was first developed. Ferguson and Finlayson (1972) presented error bounds for Eqs. (10) and (10a), while Michelsen and Villadsen (1981) give an error expression for the linear problem, k = 1. Examples of other early work on convergence analysis are in Gottlieb and Orzag (1977) and Canuto and Quarteroni (1981). More recent work is summarized in many of the reference books cited in section 3.4 [e.g. Canuto, et al. (2006)]. Since Figs. 5 and 6 are so different let us take a closer look at these measures of the error. Given the exact solution, y*, the Lp error norm is a measure of the error in the internal profile: 1

1 𝑝

(23)

𝜖𝑝 = [∫ |𝑦 − 𝑦 ∗ |𝑝 𝑑𝑥 ] 0

The error in the normalized flux for this linear source function is: 1

𝜖𝜂 = |∫ (𝑦 − 𝑦 ∗ )𝑑𝑥 |

(24)

0

These two error measures look similar, especially when p = 1. However, upon closer examination, it is clearly possible to achieve an exact flux with an imperfect solution which oscillates about the exact solution but in a way that gives the correct solution on average. If one plots y – y* it is obvious that the error at the collocation points is smaller than the average

[18]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Lp error. The combination of greater accuracy at the points and accurate quadrature makes for the accurate flux calculations in Fig. 6. Villadsen and Michelsen (1978) (p. 85) also note that other methods may produce similar profile errors, such as shown in Fig. 5, but the moments and Galerkin methods do an especially good job of balancing the error, so that the flux calculation is very accurate. As we shall see, collocation at Gauss and Lobatto points have equivalence to moments and Galerkin methods, respectively. The only study which seems to address this better than expected performance of certain global approximations is the one by Lanczos (1973). It is likely the more rapid convergence rate for fluxes with Lobatto and Gauss points is related to the long known superconvergence property of finite element methods [see Křížek and Neittaanmäki (1998)]. The superconvergence phenomenon is one where the solution and/or its derivative converge faster at certain points than the solution on average. Superconvergence for orthogonal collocation finite element procedures is discussed in Section 7. Since the global methods considered here are but a single element of a finite element procedure, it is likely the exceptional accuracy observed here in Fig. 6 is a related phenomenon. The error curves in Figs. 5 and 6 converge at a supergeometric rate, i.e. with (n)log(n), but are plotted versus n as is customary. Lobatto points give the best convergence rate. Gauss points converge at the same rate, but the errors are about 10 to 15 times larger, equivalent to about one increment of n. The convergence rate with Chebyshev points is roughly half that with Gauss or Lobatto points. If derivatives are used to calculate the flux, Eq. (22), the results are poor with Lobatto or Chebyshev points. Gauss points produce the same result regardless of calculation method. Since all the methods give exponential convergence and the L1 and L2 error norms are similar, one could easily be misled by an incomplete comparison. Significant differences show up only when fluxes are compared. One of the advantages of orthogonal collocation or pseudospectral methods is that virtually exact solutions are feasible. The method with Lobatto points is superior for achieving high accuracy for this example, but only if fluxes are accurately calculated. 4.1.2 Linear Source, Third Kind B.C. Now, consider the same problem as above, but with the third kind boundary conditions, Eq. (10b). It is clear from Eqs. (12) and (14) that φ/Bi gives the relative importance of internal and external resistance. The shape of the profile is not changed, but the boundary value is scaled up to account for the external resistance. We first present calculations with φ = 5 as above and Bi = 10. For these conditions both are important, but the external resistance is less important than the internal resistance. Later, other values of Bi are considered to determine its effect on the results. Most texts recommend that boundary collocation, Eq. (19), or an equivalent strong method be used to approximate the boundary condition. Given that Fig. 6 shows poor accuracy of [19]

derivatives (Lobatto and Chebyshev) for calculating fluxes, one might question the suitability of this approach for all but Gauss points. 1

y

Fig. 7 compares the profiles for solutions calculated with n = 4. Note 0.8 that the boundary value of y is about 0.3, indicating roughly 30% of the resistance is external. The relative 0.6 accuracy of the methods appear similar Lobatto, n = 4 to that shown in Fig. 3 with Dirichlet Gauss, n = 4 Chebyshev, n = 4 conditions. Fig. 8 shows the L2 error Exact 0.4 norms are almost the same for Chebyshev and Gauss points and averages about 50 percent greater for 0.2 0 0.2 0.4 0.6 0.8 Lobatto points. If only the norms are x Fig. 7 Solutions Eq. (6), Bi = 10,  = 5, k = 1 compared, one might conclude that all methods give similar accuracy. However, Fig. 9 shows a significant disparity in fluxes calculated with Eq. (21). Since Eq. (12) shows that the shape of the profile is the same, the greater error must be due to a problem in the formulation. Although, the differences are not large for n < 6, the differences in convergence rate are almost a factor of 2. As discussed below, similar flux errors persist even when Bi is so large that a Dirichlet condition is approached. 100

10

-2

10

-4

10

0

10-2

Flux Error

10-4

L2 Error

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

10-6

10

10

Gauss Lobatto Coll. Chebyshev Coll.

-8

4

6

8

n

10

-6

10

-8

10

-10

2

10

12

10-12

14

Fig. 8 L2 Error Norm, Bi = 10  = 5, k = 1

Gauss Lobatto Coll. Cheb. Coll.

-10

2

4

6

8

n

10

12

Fig. 9 Flux Error, Bi = 10  = 5, k = 1

Although, significant errors occur only in the fluxes for this problem, other equations show larger errors in the internal profiles as well [e.g Young (2018)]. This problem with flux boundary conditions was pointed out in the early 1970’s [Ferguson (1971), Elnashaie and Cresswell (1973)]. The problem is clearly evident for two different problems in Finlayson (1972) (see Table 5.7 and Fig. 5.7). The problem is only briefly [20]

14

1

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

discussed by Villadsen and Michelsen (1978, p. 248). Other articles attributed the problem to the strong treatment of flux boundary conditions (see Appendix A) and recommended collocation at Gauss points for these problems [Young and Finlayson (1976), Michelsen and Villadsen (1981)]. Unfortunately, many are apparently unaware of the problem, since later texts continued to recommend Eq. (19) or its equivalent for Chebyshev and Lobatto points. This problem caused collocation at Lobatto points to fall from favor in the OC community. A superior alternative to the strong treatment, Eq. (19), can be found by examining other Methods of Weighted Residuals (MWR) and the foundation of orthogonal collocation. 4.2 Method of Moments As discussed in section 3.6, the method of moments is equivalent to the tau method as practiced in the spectral literature [Boyd (2000), Canuto, et al. (2006)]. To solve Eq. (10) by the method of moments, the residual is weighted by monomials xk for k = 0,…,n - 1; however, as discussed in section 3.2, weighting by any linearly independent set of n polynomials through degree n - 1 will give identical results. For example, we could use the first n Legendre polynomials and the results would be identical except possibly for rounding errors. Another suitable set of linearly independent polynomials are the Lagrange interpolating polynomials for only the n interior points. These polynomials are related to those in Eq. (2) by: 𝑥 𝑖 (1 − 𝑥 𝑖 ) (25) ℓ∗𝑖 (𝑥) = ℓ𝑖 (𝑥) 𝑥(1 − 𝑥) where the asterisk indicates the reduced polynomial. Weighting the residual function, Eq.(16), by these polynomials and integrating numerically gives the following problem: 𝑚

𝑛+1

𝑛+1

𝑑 2 ℓ𝑖 | ) + 𝑟 (𝑥𝑘 , ∑ ℓ𝑖 (𝑥𝑘 )𝑦𝑖 )] 𝑊𝑘 ℓ𝑗∗ (𝑥𝑘 ) = 0 ∑ [(∑ 𝑦𝑖 𝑑𝑥 2 𝑥

𝑘=1

𝑖=0

𝑘

(26)

𝑖=0

for j = 1,…,n, where xk and Wk designate quadrature base points and weights, respectively, which are yet to be specified. The method of moments requires that all boundary conditions be satisfied exactly. For Dirichlet conditions, Eq. (10a), the boundary values are substituted. For third kind conditions, Eq. (10b), Eq. (19) provides the two extra equations. For m > n the quadrature base points are naturally different from the nodal interpolation points used to define the trial functions, Eq. (2). However, if m = n, and the interpolation points coincide with the quadrature points, some wonderful simplifications occur. For this case, Eq. (26) simplifies to: 𝑛+1

∑ 𝑊𝑗 𝐵𝑗𝑖 𝑦𝑖 + 𝑊𝑗 𝑟(𝑥𝑗 , 𝑦𝑗 ) = 0

(27)

𝑖=0

Eq. (27) is equal to Eq.(18) when each row is multiplied by the quadrature weight, Wj, so the equations are equivalent. To make this approach work, we need a quadrature formula that will give an exact or accurate approximate integration of Eq. (26) with m = n and 0 < xi < 1. The most accurate quadrature of this type is Gaussian quadrature. [21]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Let us now examine the accuracy of Gaussian quadrature for integration of Eq. (26). The trial functions, ℓi(x), are polynomials of degree n + 1, so the second derivative is of degree n – 1. The weight functions, ℓj*(x) are of degree n – 1. Combining the two terms, the diffusion term is of degree 2n – 2. The first order reaction term is of degree 2n. Since Gaussian quadrature is exact for polynomials through degree 2 n – 1, integration of the diffusion terms is exact, while the source term misses exact integration by one degree. So, Eq. (27) evaluated with Gaussian quadrature is both an accurate approximation of the moments method and equivalent to collocation at Gauss points, Eq.(18). Lobatto and Radau quadrature with m = n have sufficient accuracy to integrate all the terms in Eq. (26) exactly for a first order reaction. However, neither one reduces to a collocation method because end point terms would appear in Eq. (27). These terms are zero with Gauss points, because the quadrature weights are zero at both boundaries. This is a case where quadrature weights on the boundary are not an asset as many claim. Orthogonal collocation at Lobatto or Radau quadrature base points bear no direct relationship to the moments method. With Chebyshev points, the associated Clenshaw-Curtis quadrature also has nonzero quadrature weights at the endpoints. In addition, Clenshaw-Curtis quadrature is not accurate enough to produce a good approximation to the moments method. For these reasons, it also bears no direct relationship to the moments method. To achieve an exact representation of the moments method with a first order reaction, 𝑟̂ (𝑥, 𝑦) = 1 − 𝑦, the following integration must be performed more accurately: 1

𝐷𝑗𝑖 = ∫ ℓ𝑗∗ (𝑥 )ℓ𝑖 (𝑥 ) 𝑑𝑥 ≈ 𝛿𝑗𝑖 𝑊𝑗

(28)

0

We call D the mass matrix and define D0,i = Dn+1,i = 0. With Gauss points it is approximated by the diagonal matrix of quadrature weights shown at the far right. In finite element methods, the reduction of the mass matrix to a diagonal one is called lumping. Lumping is achieved automatically here due to the approximate integration. D is also called a capacity matrix, since it is often associated with time derivative terms. To list the complete set of equations, it is convenient to combine the boundary and interior equations by defining: 𝐶𝑗𝑖 = 𝛿𝑗,𝑛+1 𝐴𝑛+1,𝑖 − 𝛿𝑗,0 𝐴0,𝑖 − 𝑊𝑗 𝐵𝑗𝑖 (29) With this definition, the complete set of equations for the first order reaction with constant coefficients and third kind boundary conditions is: 𝑛+1

𝑛+1 2

𝛿𝑗,0 2𝐵𝑖0 𝑦0 + 𝛿𝑗,𝑛+1 2𝐵𝑖1 𝑦𝑛+1 + ∑(𝐶𝑗𝑖 + 4𝜑 𝐷𝑗𝑖 )𝑦𝑖 = 4𝜑 ∑ 𝐷𝑗𝑖 = 4𝜑2 𝑊𝑗 𝑖=0

2

(30)

𝑖=0

C is a symmetric matrix which we call the stiffness matrix. For the full moments method, the complete matrix problem is not symmetric because of D. However, with the diagonal [22]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

approximation of D (far right of Eq. (28)), the equations are symmetric and positive definite and equivalent to collocation. One deficiency of orthogonal collocation, pseudospectral or differential quadrature methods is that self adjoint operators do not lead to symmetric matrix problems. This development shows that by a simple rescaling of the equations this desirable feature is achieved with Gauss points. A symmetric matrix problem cuts the calculations for solution almost in half [Faddeeva, 1959]. Symmetric matrices also have theoretical and computational benefits for other problems, e.g. the eigenvalue problem in section 5. 4.3 Galerkin Method To solve the problem with the Galerkin method, the residual is weighted by the trial functions ℓj(x). Since the Robin boundary conditions reduce to Dirichlet conditions for large Bi, we will consider only the more general conditions, Eq. (10b). As discussed in section 3.7, the Galerkin method permits flux boundary conditions to be treated strongly or weakly [Finlayson (1972)]. The weak or natural treatment is most commonly used. However, the normal procedure with the orthogonal collocation, pseudospectral or differential quadrature method is to satisfy the boundary conditions exactly, i.e. boundary collocation, Eq. (19). The short development in Appendix A shows the difficulties with this approach, especially for nonzero quadrature weights at the boundary. For this reason, we will use the natural boundary condition treatment here. The natural treatment allows the polynomial greater flexibility to adapt to the solution, while the boundary condition is satisfied approximately along with the rest of the differential equation. To solve the problem with the Galerkin method, the residual, Eq. (2), is weighted by the trial functions ℓj(x), for j = 0,…,n + 1. The equations are converted to the weak formulation by integrating the second derivative term by parts: 𝑛+1

𝑛+1

𝑖=0

𝑖=0

1 𝑑ℓ𝑗 𝑑ℓ𝑖 𝑑ℓ𝑖 1 ∫ (∑ ∑ ℓ𝑗 | 𝑦𝑖 − 𝑦 − ℓ𝑗 𝑟(𝑥, 𝑦)) 𝑑𝑥 = 0 𝑑𝑥 0 𝑑𝑥 𝑑𝑥 𝑖 0

(31)

The first term contains the two boundary derivatives, so we substitute the boundary conditions and integrate the other terms using a suitable quadrature formula: 𝑚

𝑛+1

𝛿𝑗,0 2𝐵𝑖0 𝑦0 + 𝛿𝑗,𝑛+1 2𝐵𝑖1 𝑦𝑛+1 + ∑ 𝑊𝑘 (∑ 𝑘=1

𝑖=0

𝑑ℓ𝑗 𝑑ℓ𝑖 | | 𝑦 − ℓ𝑗 (𝑥𝑘 ) 𝑟(𝑥𝑘 , 𝑦(𝑥𝑘 ))) = 0 𝑑𝑥 𝑥 𝑑𝑥 𝑥𝑘 𝑖

(32)

𝑘

The xk in Eq. (32) designate the quadrature base points, which differ from the nodal interpolation points for m > n. Consider quadrature with n interior points. Since the trial functions are polynomials of degree n+1, the diffusion term is of degree 2n. For a first order reaction, the source term is of degree 2n+2. An n point Gaussian quadrature is exact through degree 2n-1, so neither term would be integrated exactly. On the other hand, Lobatto quadrature with n interior points gives exact integration through degree 2n+1, so it gives exact integration for the diffusion term, but misses exact integration of a linear source term by one degree. This is the same level of discrepancy found when the moments method is approximated with Gaussian quadrature. Radau [23]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

quadrature is one degree less accurate than Lobatto quadrature, so Radau quadrature also integrates the diffusion term exactly. Using Lobatto quadrature with n interior points, Eq.(32) reduces to: 𝑛+1

𝛿𝑗,0 2𝐵𝑖0 𝑦0 + 𝛿𝑗,𝑛+1 2𝐵𝑖1 𝑦𝑛+1 + ∑ 𝐶𝑗𝑖 𝑦𝑖 − 𝑊𝑗  𝑟(𝑥𝑗 , 𝑦𝑗 ) = 0

(33)

𝑖=0

where: 𝑛+1

𝐶𝑗𝑖 = ∑ 𝑊𝑘 𝐴𝑘𝑗 𝐴𝑘𝑖 = 𝛿𝑗,𝑛+1 𝐴𝑛+1,𝑖 − 𝛿𝑗,0 𝐴0,𝑖 − 𝑊𝑗 𝐵𝑗𝑖

(34)

𝑘=0

Since Lobatto quadrature can perform the integration by parts exactly, it follows that the stiffness matrix, C, can be calculated by either expression above. The right expression is identical to Eq. (29). Given this relationship, it is clear that at the interior points, Eq. (33) is identical to Eq. (27) and equivalent to Eq. (18), i.e. collocation. To examine the different treatment of the boundary conditions, we compare Eqs. (19) and (33) at only one boundary, since the approximation is the same at the other boundary. Using the right expression of Eq. (34) for C, the boundary equation at j = 0 is: 𝑛+1

𝑛+1

2𝐵𝑖0 𝑦0 − ∑ 𝐴0𝑖 𝑦𝑖 − 𝑊0 (∑ 𝐵0𝑖 𝑦𝑖 + 𝑟(0, 𝑦0 )) = 0 𝑖=0

(35)

𝑖=0

Eq. (35) differs from (19) by the extra term on the right, which is the boundary quadrature weight multiplying the residual at the boundary, R(0,y) from Eq. (16). This procedure sets a weighted average of the two residuals to zero. It is still a form of collocation, since residuals are set to zero at the boundary. The boundary condition is not satisfied exactly, but the procedure drives both residuals to zero at an exponential rate. Eq. (35) is equally valid for Gauss points, since the quadrature weights are zero at the boundaries and it reduces to boundary collocation. The left expression of Eq. (34) is not valid for Gauss or Chebyshev points, because the quadrature is not accurate enough to perform the integration by parts exactly. For an infinite Bi number, Eq. (35) reduces to the Dirichlet condition, y0 = 0. For a first order reaction, the full Galerkin method requires a more accurate calculation of the source term in Eq. (32). The term gives a mass matrix, which is identical to Eq. (28) with the substitution of the trial functions, ℓ, for the reduced weighting functions in Eq. (32). A more accurate integration of the source term in Eq. (32), produces an equation identical in form to Eq. (30). For the Galerkin method, both the stiffness and mass matrices are symmetric. The mass matrix is full for the Galerkin method and is diagonal or lumped for collocation at Lobatto points. For more complex rate expressions, a full mass matrix adds calculations which do not normally improve the accuracy enough to warrant the extra complexity. [24]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

In summary, collocation at Lobatto quadrature base points is an accurate approximation of the Galerkin method, but only when a natural or weak boundary condition treatment is used. Collocation at Gauss points is an accurate approximation of the moments method. Collocation at Radau points is in between and reduces to boundary collocation at one boundary and a natural treatment at the other. The Clenshaw-Curtis quadrature is not accurate enough to give a good approximation to either the Galerkin or moments methods. Chebyshev points approximates MWR with the extra radical, 1/√(1 − 𝑥 )2 , in the weight function. The Chebyshev points and quadrature weights are intermediate between Gauss and Lobatto points and weights (see Fig. 2), so can also be justified on that basis. For Chebyshev points, we propose the use of a natural boundary condition treatment, i.e. Eq. (33), with C calculated by Eq. (29) or the right expression of Eq. (34). We will call this a natural or weak boundary condition treatment also even though the quadrature does not produce an accurate approximation of Eq. (32). Formulations using the stiffness matrix and a natural treatment of flux boundary conditions are referred to as weak formulations. Unlike with Gauss and Lobatto points, the stiffness matrix for Chebyshev points is not symmetric (except for n < 4). Consequently, almost twice the computations are required to solve matrix problems. Also note the natural boundary condition treatment causes boundary rate terms to appear in the approximation, Eq. (35). For a nonlinear rate expression, all n + 2 equations are nonlinear, whereas with Gauss points the two boundary equations are linear. These linear equations can be eliminated initially so only n nonlinear equations must be solved iteratively. The relationship between collocation at Lobatto points and the Galerkin method is what motivated Villadsen and Stewart (1967) to select Lobatto points. However, they considered only Dirichlet conditions. The natural boundary condition treatment for the Galerkin method is advocated by Finlayson and Scriven (1966). They point out this treatment sets a combination of boundary and interior residuals to zero, as in Eq. (35). The natural treatment is standard in finite element applications. The use of this procedure with collocation at Lobatto points was suggested by Young (1977) and is discussed more fully by Funaro (1992, pp. 143, 204). Canuto (1986) discussed an apparently similar procedure for Chebyshev and Lobatto points. Canuto et al. (2006) describe the procedure more clearly, but unfortunately rename it G-NI (Galerkin with Numerical Integration), which is likely to cause much confusion. Shen, et al. (2011) use the more appropriate designation of collocation in the weak form. However, all of these texts also describe the strong, boundary collocation treatment and do not claim a major benefit for a natural treatment. Consequently, most applications use boundary collocation. We see no reason to rename the method. The weak formulation is the standard, correct method, while the strong treatment is incorrect. We believe the examples below and in section 5 provide sufficient support for use of a natural or weak boundary condition treatment rather than boundary collocation and justify the words “correct” and “incorrect” used above.

[25]

4.3.1 Linear Source, Third Kind B.C., Galerkin/Moments 10 Fig. 10 shows the flux errors from Fig. 9 updated with results from full Galerkin and 10 moments methods and Lobatto and 10 Chebyshev collocation with a natural treatment of the boundary conditions. These 10 results are labeled “nat.” while the boundary 10 collocation results from Fig. 9 are labeled “bc”. 10 0

-2

Flux Error

-4

-6

-8

-10

10-12

Lobatto, nat. Gauss Cheb. nat. Lobatto, bc Cheb., bc Galerkin Moments

The improvement by using natural boundary conditions is impressive, especially with 10 2 4 6 8 10 12 14 n Lobatto points. With the natural boundary Fig. 10 Flux Error, Bi = 10  = 5, k = 1 condition treatment, the convergence rates for the three choices of points are similar to those in Fig. 6. Since the problem reduces to the Dirichlet problem for large Bi, the effect of this parameter was investigated by solving the problem with Bi = 2, 5, 10 and 50 for comparison. With the largest value, the condition approaches a Dirichlet condition. For the Lobatto cases in Fig. 10 with n = 8, the ratio of the error with boundary collocation relative to a natural treatment, is 1.4x103. This error ratio at n = 8 varies from 7.0x103 when Bi = 2 to 0.3x103 for Bi = 50. The convergence rates are similar, so the disparity grows with n. Clearly, the error with boundary collocation can be several orders of magnitude greater even for quite large values of Bi. The errors with Chebyshev points are reduced to a lesser extent by using a natural boundary condition treatment, but still almost an order of magnitude at n = 8. -14

As discussed above, the natural boundary condition treatment, Eq. (35), sets a weighted average of the boundary and interior residuals to zero. Neither residual will be identically zero, but they converge to zero at an exponential rate. Fig. 11 shows the convergence behavior of these two residuals as a function of n for Bi = 10. Graphs for the other values of Bi are similar. As shown in Eq. (35), the ratio of the two residuals is equal to the boundary quadrature weight, W0 and Wn+1 or approximately n-2 for Lobatto quadrature and half that for Clenshaw-Curtis quadrature (Chebyshev points).

10

2

10

0

R at boundary

10-2

Residual

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

10-4 10

-6

10-8 10-10 10

b.c.residual

Lobatto Chebyshev Gauss

-12

5 10 15 20 n Fig. 11 Boundary condition residual and interior residual at boundary, Bi = 10,  = 5, k = 1

We also note that in Fig. 10 the Galerkin and moments methods improve in a stair step fashion. When n is even, the results with Galerkin and Lobatto points are identical and the [26]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

results with moments and Gauss points are identical. This result is an artifact caused by the symmetry of the solution about x = 0.5. For example, with n = 4 the trial solution is a 5th order polynomial, but the highest order term has a zero coefficient because of the symmetry. Lobatto points normally miss exact integration of the Galerkin method by one degree, but since the problem symmetry knocks out the highest degree, the two methods agree. With n = 3 the polynomial is 4th order. Adding a 5th order term does not improve the Galerkin solution. With n = 3 collocation at Lobatto points does not agree with the Galerkin method because the integration is approximate. An analogous situation applies for the moments method and collocation at Gauss points. When the differential equation contains values of the dependent variable (y in this case), not exclusively its derivatives, this very special set of circumstances, i.e. linear, constant coefficient and symmetric solution, are essentially the only times when collocation at Lobatto points is identical to a Galerkin method and collocation at Gauss points is identical to the moments method. If the source r(x,y) in Eq. (10) is dependent only on x and is a polynomial of degree n or less, collocation at Gauss and Lobatto points are identical to moments and Galerkin methods, respectively. Other conditions can easily be evaluated by comparing the accuracy of the quadrature to the degree of polynomial which must be integrated. Usually, Gauss and Lobatto points give approximations to the moments and Galerkin methods, but often very good ones. Some texts incorrectly state that collocation at Gauss points replicates the Galerkin method for more general conditions [Boyd (2000), p. 89]. 4.4 Mass Conservation and Fluxes One is usually interested in the flux at the boundaries. For example, if a fluid were flowing on both sides of the slab we would want to know the rate of mass or heat transfer from the slab. For the symmetric problem, the transfer is quantified by a single normalized flux, η given by Eqs. (13) and (14). For a nonsymmetric problem one would generally want the breakdown of left and right side fluxes, while in multiple dimensions the flux profiles along the boundaries may be important. The divergence theorem, Eq. (14), or average energy equation, gives only the total flux. If the sum of the individual boundary fluxes (or integral in multiple dimensions) equals the total flux, this gives us greater confidence in the individual values. Conversely, if the fluxes do not achieve an overall balance, we have less confidence in the solution. The total of the fluxes will give the average rate provided the method conserves mass. The divergence theorem, Eq. (14), is nothing more than an overall balance (of mass, energy, etc.). It is derived by integrating Eq. (10) across the domain. In general, a method will be conservative if the integral of the residual is zero. If the MWR weights contain unity in some combination the method will be conservative. Moments and collocation at Gauss points are conservative. If formulated with monomials or Legendre polynomials the first weight function would be x0 or P0 = 1. In our formulation of section 4.2, ∑ ℓ∗𝑖 (𝑥 ) = 1. Since these methods are conservative, Eq. (14) is obeyed, so the [27]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

same value of the boundary flux is calculated by integration, Eq. (21) or by differentiation of the approximate solution, Eq. (22). This explains why the approximations in Fig. 6 are the same with either method of calculation for collocation at Gauss points. For the Galerkin method and collocation at Lobatto points, the sum of all the Lagrange interpolating polynomials is also unity. However, if we consider Dirichlet boundary conditions, the boundary values are directly substituted into Eq. (2) so the first and last interpolating polynomials are not used as weight functions. The method appears not to be conservative due to the leftover terms on the right side below: 1 1 𝑑𝑦 1 (36) | + ∫ 𝑟(𝑥, 𝑦) 𝑑𝑥 = ∫ [ℓ0 (𝑥 ) + ℓ𝑛+1 (𝑥 )]𝑅(𝑥, 𝒚)𝑑𝑥 𝑑𝑥 0 0 0 Using quadrature, Eq. (36) is approximated by: 1 𝑑𝑦 1 (37) | + ∫ 𝑟(𝑥, 𝑦) 𝑑𝑥 = 𝑊0 𝑅(0, 𝒚) + 𝑊𝑛+1 𝑅(1, 𝒚) 𝑑𝑥 0 0 The two terms on the right side are those needed to correct the fluxes in the natural boundary condition treatment, Eq. (35), so the correction makes the method conservative. To be consistent with the Galerkin method, the flux on the left side should be approximated by: 𝑛+1

𝑛+1

𝑛+1

𝑑𝑦 | = ∑ 𝐴0𝑖 𝑦𝑖 + 𝑊0 (∑ 𝐵0𝑖 𝑦𝑖 + 𝑟(0, 𝑦0 )) = ∑ 𝐴0𝑖 𝑦𝑖 + 𝑊0 𝑅(0, 𝒚) 𝑑𝑥 𝑥=0 𝑖=0

𝑖=0

(38)

𝑖=0

The same equation written with the stiffness matrix is: 𝑛+1

𝑑𝑦 | = − ∑ 𝐶0𝑖 𝑦𝑖 + 𝑊0 𝑟(0, 𝑦0 ) 𝑑𝑥 𝑥=0

(39)

𝑖=0

Analogous equations apply at the right boundary. When computed from these equations, the two fluxes will be consistent with Eq. (14), the average energy equation or divergence theorem. For the full Galerkin method, the fluxes must be calculated with the equivalent expression using the mass matrix, D. For the Robin boundary conditions, the expressions are consistent with Eq. (10b), so it would be far simpler to calculate fluxes by multiplying 2Bi by the boundary value, y0 or yn+1. However, that calculation is subject to roundoff errors when Bi is large. Note that Eqs. (38) and (39) are also valid with Gauss points, since the boundary quadrature weights are zero. This equivalence is useful when writing one computer code which will work with all types of points, provided one is not averse to a few unnecessary calculations. Using these equations with Chebyshev points will also make that method conservative, which provides some additional justification for using a natural treatment with Chebyshev points. The relationship shown above, between the overall balance and natural boundary condition treatment was noted by Finlayson and Scriven (1966). For flux boundary conditions, Ferguson (1971) suggested replacing boundary collocation with the integral of the rate, see Eq. (14), to

[28]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

enforce mass conservation and improve accuracy. This integral procedure has several limitations. The natural treatment is equivalent, more general and easier to implement.

5. Parabolic Problem - Falling Liquid Film In their original article, Villadsen and Stewart (1967) considered two parabolic problems: (1) the transient version of the catalyst pellet problem and (2) eigenvalues for the classic Graetz problem. As an example of a parabolic equation we consider a convective mass transfer problem. The problem is one of mass transfer from a liquid film flowing down a wall. Bird, et al. (1960, p. 538) describe the problem and present a penetration solution valid near the inlet. The problem was also treated by Villadsen and Michelsen (1978) using collocation at the base points of Radau quadrature. We wish to supplement their results by considering other choices of points, Gauss, Radau, Lobatto and Chebyshev points. The problem has a no flux (Neuman) condition at one end, so it provides an opportunity to further test different methods for treating boundary flux conditions. Consider laminar flow of a liquid film down a solid wall. Let the gas-liquid interface be at x = 0 and the solid wall at x = 1. The film is in laminar flow, so the velocity profile is parabolic, 0 at x = 1 and a maximum at x = 0. Like most convective heat and mass transfer problems, the no slip condition at the wall creates a singularity there. Assume the entering liquid is at one composition and the gas-liquid interface is at a different composition. Since the wall is solid, assume no flux there. The dimensionless governing equations for the problem are: 𝜕𝑦 𝜕 2 𝑦 (40) (1 − 𝑥 2 ) = 𝜕𝑧 𝜕𝑥 2 with: 𝜕𝑦 (40a) (1, 𝑧) = 0 and 𝑦(𝑥, 0) = 1 𝑦(0, 𝑧) = 0, 𝜕𝑥 The heat transfer analog of the problem is obvious. This model can be compared to a simpler lumped parameter model where the transfer is governed by a mass transfer coefficient: 2 𝑑𝑦̅ (41) = −𝑆ℎ 𝑦̅ 3 𝑑𝑧 where Sh is the Sherwood number. For heat transfer, the analogous parameter would be the Nusselt number, a dimensionless heat transfer coefficient. 𝑦̅ is the mixing cup average composition given by: 1

𝑦̅(𝑧) =

∫0 (1 − 𝑥 2 ) 𝑦 𝑑𝑥 1

∫0 (1 − 𝑥 2 ) 𝑑𝑥

 =

3 1 ∫ (1 − 𝑥 2 ) 𝑦(𝑥, 𝑧) 𝑑𝑥 2 0

(42)

By using the average energy equation or divergence theorem, the Sherwood number can be calculated from the solution by two different methods: −2 𝑑𝑦̅  1  𝜕𝑦 𝑆ℎ = = (0, 𝑧) (43) 3𝑦̅ 𝑑𝑧 𝑦̅ 𝜕𝑥

[29]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

The Sherwood number is another normalized flux quantity, like the effectiveness factor, η, in the catalyst pellet problem of Section 4. The conventional method of solution by the orthogonal collocation or pseudospectral method is to use the differentiation matrix to convert the problem to a coupled set of ordinary differential equations: 𝑛+1

(1 −

𝑥𝑗2 )

𝑑𝑦𝑗 (𝑧) − ∑ 𝐵𝑗𝑖 𝑦𝑖 (𝑧) = 0 𝑑𝑧

(44)

𝑖=0

for j = 1,…,n, with 𝑛+1

𝑦𝑗 (0) = 1, 𝑦0 (𝑧) = 0, and ∑ 𝐴𝑛+1,𝑖 𝑦𝑖 (𝑧) = 0

(44a)

𝑖=0

As before, A and B give approximations to the first and second derivatives as defined in Eq. (20). Eq. (44a) approximates the no flux condition at x = 1 using boundary collocation as suggested in most texts on the orthogonal collocation or pseudospectral method. Based on the results in section 4 for the catalyst pellet problem, we can expect the boundary collocation approach to be inferior to a natural treatment, except for Gauss or Radau-left points since they give a quadrature weight of zero at the no flux boundary, Wn+1 = 0. The base points and weights for Gauss, Radau-left, Lobatto and Clenshaw-Curtis quadrature (Chebyshev points) are compared in Fig. 2. The Radau-left points used by Villadsen and Michelsen are skewed away from x = 0 and are more densely spaced near x = 1. This choice is somewhat counterintuitive, since for this problem large gradients occur near x = 0 for small z and the point spacing is less dense in that area. The Radau-right points are the mirror image of the Radau-left points and have a greater density of points near x = 0. Although Villadsen and Michelsen did not state their reasoning, most likely they selected Radau-left points to achieve the greatest quadrature accuracy with Wn+1 = 0. As discussed above in section 4, an alternative weak formulation can be constructed using the stiffness matrix. For Eq. (40) the weak formulation is: 𝑛+1

(1 −

𝑥𝑗2 )𝑊𝑗

𝑑𝑦𝑗 + ∑ 𝐶𝑗𝑖 𝑦𝑖 = 0 𝑑𝑧

(45)

𝑖=1

for j = 1,…,n+1 where C is the stiffness matrix given by Eq. (29). In Eq. (45) we have used the boundary conditions, Eqs. (40a), which is not obvious since a zero was substituted. By examining the definition of C it is clear that at the interior points Eq. (45) is equal to Eq. (44) multiplied by the quadrature weights, Wj, so the equations are equivalent. They differ only at the boundary point, n+1, where Eq. (45) is equal to:

[30]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

𝑛+1

𝑛+1

∑ 𝐴𝑛+1,𝑖 𝑦𝑖 − 𝑊𝑛+1 ∑ 𝐵𝑛+1,𝑖 𝑦𝑖 = 0 𝑖=0

(46)

𝑖=0

Eq. (46) differs from Eq. (44a) by the additional term which is again the interior residual evaluated at the boundary multiplied by the boundary weight, Wn+1. This is the same relationship as before, see Eq. (35). With this natural treatment of the boundary condition, we can expect the residual of the boundary condition to converge along with the residual of the differential equation as in Fig. 11. So, the no flux boundary condition is embedded in Eq. (45). The boundary weight is zero for Gauss and Radau-left points, so the two formulations, Eqs. (44a) and (46), are equivalent for those cases. The stiffness matrix, C, is symmetric for Gauss, Radau (left and right) and Lobatto points and for Chebyshev points when n < 4. The weak form, Eq. (45), with Lobatto points approximates the Galerkin method. The stiffness matrix is exact, but the mass matrix is approximated by the diagonal one on the left side of Eq. (45). Due to the parabolic velocity profile, Lobatto quadrature misses exact integration of the Galerkin mass matrix by 3 degrees. A full Galerkin method can be determined by replacing the diagonal approximation with a more accurate full mass matrix. For the conventional formulation values, y0 and yn+1, are eliminated using the boundary conditions so Eq. (44) is reduced to: 𝑛

(1 −

𝑥𝑗2 )

𝑑𝑦𝑗 − ∑ 𝐵̂𝑗𝑖 𝑦𝑖 = 0 𝑑𝑧

(47)

𝑖=1

where: 𝐴𝑛+1,𝑖 𝐴𝑛+1,𝑛+1 Similarly, elimination of boundary values from Eq. (45) gives: 𝐵̂𝑗𝑖 = 𝐵𝑗𝑖 − 𝐵𝑗,𝑛+1

𝑛

𝑑𝑦𝑗 (1 − 𝑥𝑗2 )𝑊𝑗 + ∑ 𝐶̂𝑗𝑖 𝑦𝑖 = 0 𝑑𝑧

(48)

𝑖=1

where: 𝐶̂𝑗𝑖 = 𝐶𝑗𝑖 − 𝐶𝑗,𝑛+1

𝐶𝑛+1,𝑖 𝐶𝑛+1,𝑛+1

for j = 1,…,n 5.1 Continuous Solutions in z Eqs. (47) or (48) are a coupled set of ordinary differential equations. They would normally be solved numerically by one of many methods for initial value problems, e.g. Runge-Kutta or others [Hairer, et al. (1993), Hairer and Wanner (1996)]. Solutions of this type are often called the Method of Lines, referring to the lines traced out in z by the collocation points in x. However, to gain greater insight, we will solve the problem analytically. The analytical solution can be found by assuming a solution of the form 𝑦 = 𝑒 −𝜆𝑧 and then computing the eigenvalues and eigenvectors of the generalized eigenproblem: [31]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

̂ 𝒚 = 𝜆𝑾 ̂𝒚 (49) 𝑪 ̂ is replaced by −𝑩 ̂ for the first formulation, and 𝑾 ̂ is the diagonal matrix In Eq. (49), 𝑪 multiplying dy/dz in either formulation. If the problem is solved by a full Galerkin method the dense mass matrix prevents the elimination of the boundary point, yn+1, so the coefficients of dyn+1/dz are nonzero, C is used in ̂ and the full mass matrix is substituted for 𝑾 ̂. Consequently, since i and j = 1,…n + place of 𝑪 1 there is an additional eigenvalue with the Galerkin method. Some would call this a spurious eigenvalue. A dense mass matrix creates two complications for numerical solutions of the problem: (1) because of the full mass matrix a linear algebra problem must be solved even for an explicit method and (2) the extra eigenvalue is usually large making the problem stiff so an explicit method is not appropriate anyway. The continuous time solution to the problem is [Franklin (1968)]: 𝑛

𝑦𝑖 =

𝑛

𝑛

−1 ( ) −𝜆𝑘 𝑧 ∑ ∑ 𝑣𝑖𝑘 𝑣𝑘𝑗 𝑦𝑗 0  𝑒 𝑘=1 𝑗=1

= ∑ 𝑎𝑖𝑘 𝑒 −𝜆𝑘 𝑧

(50)

𝑘=1

where vik, are the eigenvectors. Coefficients 𝑎𝑛+1,𝑘 for the boundary value yn+1 are calculated by substitution of Eq. (50) into Eq. (44a) or (46). Using Eq. (42) the mixing cup composition is: 𝑛

𝑛

𝑛

3 𝑦̅ = ∑ ∑ 𝑊𝑖 (1 − 𝑥𝑖2 ) 𝑎𝑖𝑘 𝑒 −𝜆𝑘𝑧 = ∑ 𝑐𝑘 𝑒 −𝜆𝑘 𝑧 2 𝑘=1 𝑖=1

(51)

𝑘=1

Determination of the Sherwood number requires calculation of the boundary flux. As shown in Eq. (43), the boundary flux can be calculated from either 𝜕𝑦/𝜕𝑥|𝑥=0 or from 𝑑𝑦̅/𝑑𝑧. If calculated correctly as discussed in the section 4.4, it makes no difference which method is used. The coefficients of the boundary flux are conveniently given by the weak form of the approximation: 𝑛+1

𝑛

𝑛

𝑛

𝜕𝑦 | = − ∑ 𝐶0,𝑖 𝑦𝑖 = − ∑ ∑  𝐶0,𝑖 𝑎𝑖𝑘 𝑒 −𝜆𝑘 𝑧 = ∑ 𝑑𝑘 𝑒 −𝜆𝑘𝑧 𝜕𝑥 𝑥=0 𝑖=1

𝑘=1 𝑖=1

(52)

𝑘=1

Eq. (52) gives the same boundary flux as the differentiation of Eq. (51). The Sherwood number is then: 𝑆ℎ =

−2 ∑𝑛𝑘=1 𝜆𝑘 𝑐𝑘 𝑒 −𝜆𝑘 𝑧 ∑𝑛𝑘=1 𝑑𝑘 𝑒 −𝜆𝑘 𝑧 = 𝑛 ∑𝑘=1 𝑐𝑘 𝑒 −𝜆𝑘 𝑧 3 ∑𝑛𝑘=1 𝑐𝑘 𝑒 −𝜆𝑘 𝑧

(53)

The asymptotic Sherwood number applies for large z: −2𝜆1 𝑆ℎ∞ = lim 𝑆ℎ = 𝑧→∞ 3

(54)

The analytical solution to the problem is like Eq. (51) except that it is an infinite series. The collocation solution approximates a truncation of the infinite series. In the analytical solution, all

[32]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

the eigenvalues are not only real and negative but also the coefficients, c, are all positive. The solution is a series of different modes each decaying at a different rate. The statements above apply to the analytical solution, but what about the numerical solutions? Whether an approximation produces real eigenvalues and positive values for c, is one test for ̂ and 𝑾 ̂ in Eq. (49) are symmetric and the suitability of an approximate method. When 𝑪 positive definite the eigenvalues are all real and positive. The weak formulation, Eqs. (45) and (48), produces symmetric positive definite matrices for all but Chebyshev points with n > 3. As explained above, the conventional formulation, Eqs. (44) and (47), and weak formulation are equivalent when the boundary quadrature weight is zero, Wn+1 = 0, which is true for Gauss and Radau left points. Boundary collocation makes the matrix asymmetric for Lobatto and Radau right points and creates additional asymmetry for Chebyshev points. The consequences of this asymmetry will be discussed shortly. An eigenvalue problem like this can be solved more efficiently when matrices are symmetric. When the matrices are not symmetric there is no advantage to treating the eigenproblem in ̂, to give the generalized form, Eq. (49), so one should divide through by diagonal matrix, 𝑾 identity matrix on the right hand side. To test for possible complex eigenvalues and negative coefficients, c, all the nonsymmetric eigenproblems were solved for n ≤ 60. For the weak formulation, Chebyshev points produce only real eigenvalues and positive coefficients, c. For the conventional formulation, Eqs. (44) and (47), Radau-right points give complex eigenvalues for n > 4, Lobatto points give complex values for n = 4-6 and n > 11, Chebyshev points give complex values for n = 28, 45, 50 and 58. The coefficients, c, are negative for many cases when the eigenvalues are real. The propensity to produce complex eigenvalues and negative coefficients is directly related to the magnitude of the boundary quadrature weights, Wn+1, and hence the degree of difference between Eqs. (44a) and (46). For large n, the Chebyshev (Clenshaw-Curtis) boundary weights are half the Lobatto and Radau boundary weights. On the surface, these results appear to be different from previous work which has analyzed approximations to the constant coefficient heat equation. Gottlieb and Lustman (1983) found that Chebyshev points produce only real eigenvalues for the heat equation with general boundary conditions that include those in Eq. (40a). The analysis of the heat equation by Canuto, et al. (1988, p.407) proved real eigenvalues for Lobatto and Chebyshev points with Neumann and Dirichlet boundary conditions. However, they consider a weak formulation with a natural boundary condition treatment, like that used here. Funaro (1992, pp. 143, 204) also considered eigenvalues with a natural boundary condition treatment. Unfortunately, these previous works are frequently cited to support the suitability of the conventional formulation using boundary collocation. It is obviously incorrect to generalize the results to a different boundary condition treatment, and in this case to a problem with variable coefficients [33]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

To reconcile these differences, we did some limited testing of the constant coefficient heat equation with the more common boundary collocation treatment. We found that with boundary collocation, Chebyshev points produced only real negative eigenvalues, but the coefficients, c, are frequently negative. For Lobatto points, the eigenvalues were complex for n > 6, so with Lobatto points the conventional formulation with boundary collocation is unsuitable even for the constant coefficient problem. The accuracy of our results were double checked with two completely independent calculations: (1) in Fortran using the LAPack code DGEEV [Anderson, et al. (1999)] and (2) in Octave using the eig function and code built upon that of Trefethen (2000). The results agreed to within the limits of machine precision (about 13 digits). Tables 2 and 3 summarize this discussion concerning the suitability of the various formulations. Based on these results, we conclude that Lobatto and Radau-right points are unsuitable using the conventional formulation with boundary collocation but are suitable using the weak formulation. We conclude the conventional formulation with Chebyshev points (excluding n = 28, 45, 50 and 58) is marginally acceptable. The other formulations are suitable, especially those that produce symmetric matrix problems. Table 2 Conventional Formulation, Eq. (44) Equivalent Points Eigenvalues To Eq. (45) Gauss yes real Chebyshev no complex Lobatto no complex Radau Left yes real Radau Right no complex

Table 3 Weak Formulation, Eq. (45) Points Symmetric Eigenvalues Gauss yes real Chebyshev no real Lobatto yes real Radau Left yes real Radau Right yes real

Figs. 12 shows the first few eigenvalues for n = 6. The notation “bc” (boundary collocation) designates the conventional formulation of the Chebyshev method, i.e. Eq. (44) and (44a). Fig. 13 shows the convergence of the first and fourth eigenvalues. The convergence of the coefficients mirrors that of the eigenvalues, so are not presented. The relative convergence rates of the various points are like those for the catalyst pellet problem. However, in this case we have added the Radau points, which tend to have accuracy between that of Lobatto and Gauss points. Chebyshev points with boundary collocation give the worst convergence, which is disconcerting since it is probably the most popular formulation. Although no formal analysis was performed, these calculations appear to reach the limits of roundoff with about 13 digits of accuracy. Figs. 12 and 13, show that the first n/2 eigenvalues are usually accurate and the higher eigenvalues are greater than the actual ones.

[34]

10

104

10

Exact Chebyshev bc Chebyshev Radau Right Radau Left Lobatto Gauss Galerkin

3

10

0

Lobatto Gauss Radau L Radau R Chebyshev Chebyshev bc Galerkin

-2

10-4

error

10

eigenvalue

2

10-6

4

10-8

1 10-10

10

1

1

2

3

4 5 n Fig. 12 Eigenvalues calculated with n = 6

6

7

10-12

5

10

15

20

25

n Fig. 13 Error in first and fourth eigenvalues

Villadsen and Michelsen give asymptotic relationships for the eigenvalues and coefficients of the analytical solution. From information given there and the approximate solutions the following asymptotic relationships were determined: 𝜆𝑘 → (4𝑘 − 1.672)2 𝑐𝑘 → 3.820/𝜆𝑘

(55)

These relationships give values within 0.1% for k = 6 to 51 of the values calculated with 100 Lobatto points. 6 10 Fig. 12 shows that the higher eigenvalues Lobatto are greater than the actual ones. The Gauss 5 10 Chebyshev magnitude of the largest eigenvalues is of Radau Left interest when explicit numerical integration Radau Right 4 10 Galerkin methods are considered. For these 2 (4n-1.672) methods, the maximum step size is 3 10 inversely proportional to the magnitude of 2 the maximum eigenvalue. Previous work 10 has found that asymptotically, the 101 maximum eigenvalues for the constant coefficient heat equation increase with n4 0 10 5 10 15 [Canuto, et al. (2006)]. Fig. 14 shows the n Fig. 14 Maximum eigenvalues vs. number of points first 15 maximum eigenvalues for this problem. For those shown, the maximum values are up to two orders of magnitude larger than the actual ones. Larger eigenvalues have been calculated and those for n = 30 - 100 fit to the following equation: 𝑛 𝑏 (56) 𝜆𝑚𝑎𝑥 → 𝑎 ( ) 30 Values for Eq. (56) are given in Table 4. It appears that if still larger values were fit, the true asymptote would have an exponent of 6. The additional power of 2 relative to the constant ̂, given the point spacing is O(n-2) coefficient problem is to be expected from the norm of 𝑾

Max Eigenvalue

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

[35]

near the boundaries. As a result, the eigenvalues for large n are much larger than Lobatto Gauss Chebyshev Galerkin 6 6 6 6 a 4.05x10 18.3x10 7.98x10 62.5x10 those for the constant coefficient heat 5.822 5.932 5.869 5.690 b equation. The eigenvalues of Radau points were not fit, since for large n the values for Radau left points follow those for Gauss points and the values for Radau right points follow those for the Lobatto points. Table 4 Parameters of Eq. (56)

For this problem, explicit numerical integration methods in z are only practical for relatively small n. If for larger n, the step size is limited by stability considerations, then relative to Lobatto points, Chebyshev and Gauss points require 2 and 4.5 times as many steps, respectively.

If one is interested in accurate solutions for very small z then larger n is required. The behavior at small z is much like that of an infinite series analytical solution when the series is truncated as it must be for practical calculations. Fig. 16 shows the profiles at z = 0.01 for n = 3. For the solutions in Fig. 16, 1% is a typical error in average composition. Contrary to intuition, the Radau-right points are not appreciably more accurate even though the points are more densely spaced near x = 0.

Sh

yavg

1 20 Although Eq. (56) is interesting, it is of little 18 0.8 16 Lobatto practical consequence if only a few points are 14 Gauss 0.6 12 required to give the desired accuracy. The Chebyshev 10 Radau Left real power of these methods is that they Radau Right 0.4 8 frequently produce good accuracy with only a exact 6 few points. For example, Fig. 15 shows y̅ and Sh as a function of axial position for the 0.2 4 various point choices with n = 2. All the higher terms damp out by z = 0.10, so the solution is asymptotic from there on. The 2 0 0.1 0.2 0.3 0.4 slope and intercept of the asymptotic part is z determined by the first eigenvalue and Fig. 15 Average composition and Sh vs z, n = 2 coefficient, respectively. Since one would normally be interested in transferring all or most of a component, the figure shows the portion of the curve for approximately 90 percent transfer which is reached at z = 0.4. Except for Gauss points, all the methods give good results for z > 0.05.

1.2

1

0.8

y

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Lobatto n = 3 Gauss n = 3 Chebyshev n = 3 Radau Lft n = 3 Radau Rgt n = 3 exact

0.6

0.4

0.2

0

0

0.2

0.4

0.6

0.8

x Fig. 16 Calculated profiles at z = 0.01, n = 3

[36]

1

Sh

60 Larger n is required to determine accurate 50 flux or Sh values at small z. Fig. 17 shows 40 Lobatto n = 2 the calculation of Sh for Lobatto points with 30 Lobatto n = 4 Lobatto n = 6 different n. From a penetration solution valid 20 Lobatto n = 8 for small z, the Sherwood number is Lobatto n =10 exact proportional to 1/√𝑧. Any polynomial 10 solution will eventually breakdown as one moves toward the discontinuous initial profile, so all the different points display the general behavior seen in Fig. 17. However, there are differences. At small z, some tend 10 10 10 10 10 z to level out at larger Sh and have a larger Fig. 17 Sh vs z, Lobatto points maximum positive deviation after crossing over at larger z. For a range of n, the maximum positive deviations are 4 to 5 percent for Lobatto and Radau-left points, 25 to 35 percent for Gauss and Radau-right points and 10 to 15 percent for Chebyshev points. This behavior can be seen in Fig. 15 to some extent. -4

-3

-2

-1

0

0

10 Fig. 18 shows the overall error in the Lobatto calculated flux versus n. The errors in y̅ looks Gauss 10 Radau L similar, but at a given n are 5-10 times Radau R 10 smaller at z = 0.10 and 10-50 times smaller at Chebyshev Chebyshev bc z = 0.01. Fig. 18 looks like Fig. 13, but the 10 curves are more tightly clustered for the z = 0.01 various points. For z > 0.01 and n = 6, the 10 error is typically less than 1% in the flux and 10 0.05% in y̅ . The errors in y̅ are typically less z = 0.10 than 0.1% for n = 3 and z > 0.1. These 10 5 10 15 20 25 results indicate that the behavior of the n Fig. 18 Flux error at z = 0.01 and 0.10 eigenvalues (Eq. (56), Table 3) at large n is of no practical importance unless one is only interested in the behavior at very small z. For these conditions, a finite element method is better suited. -2

-4

error flux

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

-6

-8

-10

-12

30

5.2 Summary Remarks - Global Methods and Fluxes The conventional formulation using a strong, boundary collocation treatment for flux conditions, Eqs. (19) and (44a), produces poor approximations except for cases when the quadrature weight is zero at the flux boundary (see Appendix A). With Lobatto and Chebyshev points, the results are bad enough for the boundary value problem in section 4, the parabolic problem in section 5, and other problems [Ferguson (1971), Finlayson (1972), Elnashaie and Cresswell (1973), Young (2018)], that they can be called wrong. Their use for parabolic problems often produces nonphysical complex eigenvalues. Even when the eigenvalues are not complex, it [37]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

gives the slowest rate of convergence. Errors are typically several orders of magnitude greater, see Figs.10,13 and 18. Some problems show significant errors in both the internal profiles and fluxes, while for others only the fluxes show large errors. The weak formulation is the correct treatment for Lobatto and Chebyshev points, e.g. Eqs. (33), (35), (45) or (46). In addition to greater accuracy, the weak formulation with natural boundary conditions is more efficient due to symmetric matrices for all but Chebyshev points. For flux boundary conditions, boundary collocation is the correct approach with Gauss points (or Radau points with a weight of zero at the flux boundary). If one prefers the more intuitive boundary collocation treatment, then one should follow the long-standing recommendation and use Gauss points [Young and Finlayson (1976), Michelsen and Villadsen (1981)]. Collocation at Gauss points works quite well and is the most popular choice in OC applications. This conclusion is contrary to the oft stated claim that nonzero boundary weights somehow facilitate the approximation of boundary conditions. If the boundary conditions are correctly treated, collocation at Lobatto points is usually slightly more accurate than collocation at Gauss points, i.e. equivalent to about one increment of n. Collocation at Lobatto points is somewhat better suited for explicit time stepping methods due to smaller more accurate eigenvalues, but if large n is required, explicit methods are unsuitable in any case. Collocation at Chebyshev points usually produces somewhat less accurate flux calculations and can be less efficient because matrix problems are always nonsymmetric. Collocation at Chebyshev points appears not to be competitive with the other choices unless the problem requires large n and is amenable to FFT calculations.

6. Multidimensional Applications The global orthogonal collocation methods described in the forgoing examples are readily extended to two dimensions using tensor product trial functions, Eq. (3). An example for the Poisson problem can be found in Villadsen and Stewart (1967). This author’s first application of the method was in 1971 for a class assignment to solve the problem in section 4 for finite cylindrical catalyst pellets, i.e. r-z coordinates. This problem was later the subject of a technical paper [Sorensen, et al. (1973)]. Numerous other examples are available in the literature [e.g. Michelsen and Villadsen (1972,1981), Finlayson (1974, 2003), texts listed in section 3.4]. Global methods have also been applied to distorted domains, which can be mapped to rectangles. This approach was used by Young and Finlayson (1976) and Orzag (1980), and is fully discussed in Canuto, et al. (2006). Although these methods permit some geometric flexibility, finite element based methods have much greater flexibility.

7. Orthogonal Collocation - Finite Elements The examples above used global polynomial trial functions, i.e. a single polynomial represents the solution. Although the global orthogonal collocation method works quite well for many problems, it has difficulties when the solution cannot be accurately represented by a single low [38]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

order polynomial, e.g. the catalyst pellet problem of section 4 with large φ or the falling film problem in section 5 for small z. In addition, finite element methods have greater flexibility for representing irregularly shaped domains. There was an explosion of interest in finite elements in the 1970’s [Zienkiewicz (1971), Strang and Fix (1973), Hughes (1987)]. Their success in structural mechanics led to interest in other disciplines. It is not surprising that soon after the Villadsen and Stewart article, collocation methods were extended to finite element trial functions using the same idea of collocating at quadrature points. There were several related developments in this area that were made independently. A finite element method differs from a global method by partitioning the domain and using a different polynomial in each subdomain or element. If a single element is used, the two approaches are the same, so the finite element approach is more general. With a finite element approximation, one also has some flexibility in the degree of continuity at the boundary or interfaces between elements. The notation Cn is used to describe the continuity, where n denotes the highest derivative which is continuous. To develop a finite element approximation in one dimension, we first lay out a grid: 𝑥0 < 𝑥1 < ⋯ < 𝑥𝑘 < ⋯ < 𝑥𝑛𝑒 (57) and define element size as Δ𝑥𝑘 = 𝑥𝑘 − 𝑥𝑘−1 . Within each element a local coordinate is used, 𝜉 = (𝑥 − 𝑥𝑘−1 )/Δ𝑥𝑘 . The nodes within each element are at the polynomial roots or quadrature points, ξi, so a specific node is given by 𝑥𝑖𝑘 = 𝑥𝑘−1 + Δ𝑥𝑘 𝜉𝑖 . The double and single subscript notation are related by 𝑥𝑘 = 𝑥𝑛+1,𝑘 = 𝑥0,𝑘+1 . We are interested in methods which are directly related to the global collocation methods discussed above and will subdivide these methods into those which are C1, i.e. continuous first derivatives, and those that are C0, i.e. simple continuity. With finite elements, the interface between elements is an internal boundary. The approximate solution and the flux should be continuous at these internal boundaries. The discussion of flux boundary condition treatment above applies directly to the approximations used for these internal boundary conditions. For a Galerkin method, the interelement flux conditions are usually treated in weak form like natural boundary conditions. If the conditions at the element interfaces are strongly or exactly enforced, then: 𝑦𝑛+1,𝑘 = 𝑦0,𝑘+1 and 𝑛+1

𝑛+1

𝑖=0

𝑖=0

1 1 ∑ 𝐴𝑛+1,𝑖  𝑦𝑖𝑘 = ∑ 𝐴0𝑖  𝑦𝑖,𝑘+1 Δ𝑥𝑘 Δ𝑥𝑘+1

(58)

If there is a material change at the interface with an abrupt change in say thermal conductivity, then Eq. (58) would be modified to reflect continuity of the flux. This strong enforcement of interface conditions is equivalent to boundary collocation. An alternative to the flux condition in Eq. (58) is to use trial functions with built in continuity. For cubic trial functions, Hermite cubic polynomials could be used [Hildebrand (1987), p. 282]. For higher order polynomials, a grid [39]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

like that in Fig. 19 could be used, where the Gauss Pts. X Single Nodes function and its derivative are specified at Elements, double node the element boundaries (double nodes) and X X X X X X X X X X X X the function values are interpolated at the 1 Fig. 19 Four C Quartic Elements single nodes. With built in continuity, the number of equations and unknowns is reduced by ne – 1 relative to explicit enforcement with the side conditions, Eq. (58). There are other tradeoffs between the two approaches. These elements are easily extended to multiple dimensions. 7.1 Quadrature in Finite Element Methods From the forgoing discussion it is clear orthogonal collocation is successful because it is an accurate approximation to integrated MWR methods, while simpler to apply. Following Villadsen and Stewart’s (1967) idea, the key is to place the collocation points at quadrature base points to achieve an accurate approximation to the integrals in either the moments method using Gauss points or the Galerkin method using Lobatto points. Recall from sections 4.2 and 4.3, that with orthogonal collocation, the integration is one degree shy of that required for a full moments or Galerkin method. For finite elements, is the same level of accuracy sufficient? To answer this question, we review the integration accuracy requirements for finite element methods. For a general purpose code, analytical integrals are cumbersome to implement, especially for distorted isoparametric elements, so most codes use numerical integration. Early applications used Gaussian quadrature, since it gives the highest accuracy for a given number of function evaluations. These quadrature calculations are time consuming, especially for nonlinear timedependent problems which required many such calculations. Collocation seemed like an obvious way to circumvent these problems, since integration is not required. Many studies have compared the computational cost of conventional and collocation-like methods [Young (1977), Gray and Van Genuchten (1978), Houstis, et al. (1978), Finlayson (1979), Frind and Pinder (1979), Dyksen, et al. (1984), Melenk, et al. (2001), Schillinger, et al. (2015)]. The computational cost depends on the trial functions selected, the quadrature method and the number of quadrature points used. An obvious question is – how accurate does the integration have to be? To describe the integration issues, consider a general mass matrix, like Eq. (28): 1

𝐷𝑘𝑗 = ∫ 𝑤𝑘 (𝑥 )𝜓𝑗 (𝑥)𝑑𝑥

(59)

0

where the wk are the weight functions and the ψj are the trial functions. We say the method employs full integration if the quadrature is accurate enough to integrate D exactly. With this definition, the integration may still be approximate since the equation could have nonlinear terms or variable coefficients. We call it reduced integration if the integration is not accurate enough to integrate D exactly.

[40]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Several studies have investigated the question of integration requirements for Galerkin finite element methods [Ciarlet and Raviart (1972), Fix (1972), Strang and Fix (1973), Raviart (1973), Fried (1974), Ciarlet (1978)]. Although there are some side issues related to problem smoothness, these studies show that for most problems the maximum rate of convergence can be achieved with less than full integration. Exact integration of all but the highest degree terms in Eq. (59) is sufficient to maintain the maximum convergence rate. For example, cubic elements can produce a 4th order convergence rate, i.e. the same as for interpolation of a function. For a Galerkin method, full integration would require an exact integration of all terms through 6th degree, since both terms in Eq. (59) are cubic. However, the 4th order convergence rate is maintained when only terms through 5 th degree are integrated exactly. 5th degree integration accuracy is achieved with 3 Gauss points or 4 Lobatto points. Although the theoretical analyses cited above apply only for Galerkin methods, it is obvious that the same rule applies to the moments method. A moments method with cubic trial functions would use discontinuous linear weight functions, so full integration would require exact integration through 4th degree. To achieve a 4th order convergence rate only terms through 3rd degree must be integrated exactly, requiring 2 Gauss points or 3 Lobatto points. This illustration is for cubics, but for higher order trial functions the quadrature keeps pace with the required integration accuracy. If n is incremented by one to increase the degree of the trial function, the integrand in Eq. (59) and the quadrature accuracy both increase by two, so the integration is one degree less than full integration for any n. For time dependent finite element problems, a related issue is the form of the mass matrix. When the problem is amenable to explicit time stepping methods, a diagonal or lumped mass matrix is preferred. If the mass matrix is not diagonal, the algebraic problem for each time step is like that for an implicit method. Recall from sections 4.2 and 4.3 that global orthogonal collocation produces a lumped mass matrix due to the node placement and quadrature. Fried and Malkus (1975) found the same result could be achieved for C0 Lagrange trial functions by using Lobatto quadrature and placing the nodal interpolation points at the quadrature points. Furthermore, this procedure produced no loss of convergence rate. Here again, the integration accuracy is one degree less than that needed for full integration. Apparently, Fried and Malkus did not consider integrating the stiffness matrix with the same reduced quadrature scheme as was done in later work. 7.2 C1 Collocation at Gauss Points Section 4.2 shows that collocation at Gauss points approximates the moments method. The moments method extends naturally to finite elements with residual weighting by discontinuous C-1 piecewise polynomials. The interface conditions must be strongly enforced as in Eq. (58) or with a basis having built in C1 continuity, as in Fig. 19. If the integrals in the moments method are approximated by n point Gaussian quadrature, full integration is missed by one degree. From the discussion in section 7.1, this integration is accurate enough to maintain the maximum possible convergence rate, and it reduces to collocation at Gauss points. Collocation [41]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

at Gauss points was described by DeBoor and Swartz (1973) and Douglas and Dupont (1973). Carey and Finlayson (1975) applied the method using Lagrange trial functions with the conditions of Eq. (58) satisfied explicitly. Several texts describe the method [Davis (1984), Lapidus and Pinder (1999), Finlayson (2003)]. We also note the method of moments is related to the H1 Galerkin method [Douglas, et al. (1974)]. The H1 Galerkin method is not a Galerkin method in the classic sense, because the trial and weight functions are different. The H1 Galerkin method weights the residual by the second derivative of the trial functions, so it is a classic moments method as described by N.M. Krylov (1926), see section 3.2. When the trial functions are C1, the weight functions are C-1, discontinuous piecewise polynomials. However, the weight functions for the moments method span only one element, which contributes to improved efficiency. Since moments and the H1 Galerkin methods are equivalent, C1 Orthogonal Collocation at Gauss points approximates the H1 Galerkin method. 7.3 C0 Collocation at Lobatto Points Collocation at Lobatto points has also been extended to finite elements. First, there is the Hybrid-Collocation-Galerkin method [Diaz (1975), Dunn and Wheeler (1976), Wheeler (1977)]. It is a C0 procedure using the Lagrange trial functions shown in Fig. 20. Collocation is applied at the internal Lobatto points. At the interface nodes, a Galerkin-like integral condition is used with weighting by the linear hat functions, top of Fig. 20. This approach is reminiscent of Ferguson’s (1971) integral approach (see section 4.4) for treating flux boundary conditions in global collocation. Another approach, usually considered to be fundamentally different, is described by Gray (1977), Young (1977, 1981), Hennart (1982), and later Leyk (1986, 1997). Their method uses a plain vanilla Galerkin method with Lagrange trial functions like Fig. 20 but with integration performed using Lobatto quadrature. The interpolation points are coincident with the quadrature points as in global orthogonal collocation and the finite element procedure of Fried and Malkus (1975). However, Fried and Malkus used Lobatto quadrature only for the mass matrix, while here all terms are integrated with the same reduced quadrature procedure. The quadrature misses full integration by one degree, so based on the discussion in section 7.1, it is accurate enough to maintain the maximum possible convergence rate. This procedure was called the Lobatto-Galerkin method by Young. Gray and Hennart considered only the quadratic case with integration by Simpson’s rule, while Young considered trial functions Fig. 20 Finite element Lagrange trial of any order. This method was later renamed the hp functions

[42]

[from Young (1981)]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Spectral Element method [Maday and Patera (1989), Vosse and Minev (1996), Canuto, et al. (2007), Karniadakis and Sherwood (2013)]. Apparently, these authors were unaware of the earlier work since none is referenced in a preponderance of the spectral literature. Maday and Patera are normally credited with discovery of the method. This method is a straight forward extension to finite elements of the ideas in Section 4.3. There we showed that at the interior points, the weak form, Eq. (33), is equivalent to collocation, Eq. (18). The same relationship holds at the interior points with finite elements. We will also show that at the element interfaces, the hybrid-collocation-Galerkin method and Lobatto-Galerkin or hp spectral element methods are equivalent. However, the Lobatto-Galerkin or hp spectral element method gives a matrix structure which is better suited for calculations. Once again consider the catalyst pellet problem. Convert Eq. (33) to an element k: 𝑛+1

𝑑𝑦 𝑑𝑦 1 ∑ 𝐶𝑗𝑖 𝑦𝑖𝑘 − Δ𝑥𝑘 𝑊𝑗  𝑟(𝑥𝑗𝑘 , 𝑦𝑗𝑘 ) = 0 𝛿𝑗,0 | − 𝛿𝑗,𝑛+1 | + 𝑑𝑥 𝑥0,𝑘 𝑑𝑥 𝑥𝑛+1,𝑘 Δ𝑥𝑘

(60)

𝑖=0

At the interior nodes, the first two terms of Eq. (60) are zero and by substitution of Eq. (34) it is: 𝑛+1

𝑊𝑗 ∑ 𝐵𝑗𝑖 𝑦𝑖𝑘 + Δ𝑥𝑘 𝑊𝑗  𝑟(𝑥𝑗𝑘 , 𝑦𝑗𝑘 ) = Δ𝑥𝑘 𝑊𝑗 𝑅𝑗𝑘 = 0 Δ𝑥𝑘

(61)

𝑖=0

for j = 1,…,n. Rjk designates the residual at node j in element k, so Eq. (60) is clearly equivalent to collocation at the interior nodes, just like in global methods, see Eq. (18). At the boundary node, xk, the weight function spans both elements (see Fig. 20), so the equation for the last node of element k and the first node of k+1 must be combined: 𝑛+1

𝑑𝑦 1 ∑ 𝐶𝑛+1,𝑖 𝑦𝑖𝑘 − Δ𝑥𝑘 𝑊𝑛+1  𝑟(𝑥𝑛+1,𝑘 , 𝑦𝑛+1,𝑘 ) − | + 𝑑𝑥 𝑥𝑛+1,𝑘 Δ𝑥𝑘 𝑖=0 𝑛+1

(62)

𝑑𝑦 1 ∑ 𝐶0,𝑖 𝑦𝑖,𝑘+1 − Δ𝑥𝑘+1 𝑊0  𝑟(𝑥0,𝑘+1 , 𝑦0,,𝑘+1 ) = 0 + | + 𝑑𝑥 𝑥𝑜,𝑘+1 Δ𝑥𝑘+1 𝑖=0

Using the standard natural boundary condition treatment at the interface, the flux terms cancel and we are left with: 𝑛+1

𝑛+1

𝑖=0

𝑖=0

1 1 ∑ 𝐶𝑛+1,𝑖 𝑦𝑖𝑘 + ∑ 𝐶0,𝑖 𝑦𝑖,𝑘+1 − (Δ𝑥𝑘 𝑊𝑛+1 + Δ𝑥𝑘+1 𝑊0 ) 𝑟(𝑥𝑛+1,𝑘 , 𝑦𝑛+1,𝑘 ) = 0 Δ𝑥𝑘 Δ𝑥𝑘+1

(63)

The conditions at the interface, Eq. (58), are strongly imposed for the values, but not the flux. As with boundary flux conditions in the global methods, Eq. (63) is equivalent to setting a weighted average of residuals to zero. This relationship is made evident by the substitution of Eq. (34) into Eq. (63): 𝑛+1

𝑛+1

1 1 ( ∑ 𝐴𝑛+1,𝑖 𝑦𝑖𝑘 − ∑ 𝐴0,𝑖 𝑦𝑖,𝑘+1 ) − Δ𝑥𝑘 𝑊𝑛+1 𝑅𝑛+1,𝑘 − Δ𝑥𝑘+1 𝑊0 𝑅0,𝑘+1 = 0 Δ𝑥𝑘 Δ𝑥𝑘+1 𝑖=0

𝑖=0

[43]

(64)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

The term within the large parenthesis is the residual of the interface flux condition, Eq. (58), while the terms on the right are the interior residuals at the interface, i.e. R as defined in Eq. (61). The value of the residual on both sides of the interface appear, since the residual is discontinuous. This result is analogous to Eqs. (35) and (46). The method is still basically collocation, but in this case a weighted average of residuals is set to zero. The residual of the interface condition and both interior residuals will be driven to zero as the grid (either n or ne) is refined. Since Eq. (60) is equivalent to collocation at the interior nodes, the only difference between this procedure and the hybrid-collocation-Galerkin procedure is the treatment of the element interface. The hybrid method uses the linear hat functions to span the interface while the Galerkin method uses the endpoint Lagrange polynomials shown in Fig. 20, resulting in Eq. (63). We will show that since the linear basis is included within the Lagrange basis, the hybrid method is equivalent. For a linear source function, the residual is a polynomial of degree n + 1, so the linear weighting of the hybrid method gives integrals of degree n + 2. Lobatto quadrature, exact for 2n + 1, is more than sufficient, so it is used for all integrals. We make use of the following relationships which come from the treatment of the second order term in the hybrid method. Integration by parts gives: 1 𝑑 2 ℓ𝑖 𝑑ℓ𝑖 1 𝑑ℓ𝑖 𝑑ℓ𝑖 ∫ 𝜉 2 𝑑𝜉 = 𝜉 | −∫ | 𝑑𝜉 = − 𝛿𝑛+1,𝑖 + 𝛿0,𝑖 𝑑𝜉 𝑑𝜉 0 𝑑𝜉 𝑛+1 0 0 𝑑𝜉 1

𝑛+1

or

𝑛+1

(65)

∑ 𝜉𝑗 𝑊𝑗 𝐵𝑗𝑖 = 𝐴𝑛+1,𝑖 − ∑ 𝑊𝑗 𝐴𝑗𝑖 = 𝐴𝑛+1,𝑖 − 𝛿𝑛+1,𝑖 + 𝛿0,𝑖 𝑗=0

𝑗=0

Eq. (65) is for the left half of the hat function. An analogous equation applies to the right half. For the hybrid method, the interface equation in weak form, i.e. after integration by parts, is: 𝑛+1

𝑛+1

𝑖=0

𝑖=0

1 1 1 1 1 𝑑ℓ𝑖 1 𝑑ℓ𝑖 ∑ 𝑦𝑖𝑘 ∫ ∑ 𝑦𝑖,𝑘+1 ∫ 𝑑𝜉 − 𝑑𝜉 − Δ𝑥𝑘 ∫ 𝜉𝑟𝑘 𝑑𝜉 − Δ𝑥𝑘+1 ∫ (1 − 𝜉)𝑟𝑘+1 𝑑𝜉 = 0 Δ𝑥𝑘 Δ𝑥𝑘+1 0 𝑑𝜉 0 𝑑𝜉 0 0

(66)

where rk designates the source function in element k. The integrals can be simplified with Eq. (65); however, to show equivalence, we undo the integration by parts reflected in Eq. (66). After substitution from Eq. (65) for the left element and its analog for the right one, Eq. (66) becomes: 𝑛+1

𝑛+1

𝑛+1

𝑛+1

𝑖=0

𝑗=0

1 1 ∑ 𝑦𝑖𝑘 (𝐴𝑛+1,𝑖 − ∑ 𝜉𝑗 𝑊𝑗 𝐵𝑗𝑖 ) − ∑ 𝑦𝑖,𝑘+1 (𝐴0,𝑖 + ∑(1 − 𝜉𝑗 )𝑊𝑗 𝐵𝑗𝑖 ) Δ𝑥𝑘 Δ𝑥𝑘+1 𝑖=0

𝑗=0

𝑛+1

𝑛+1

−Δ𝑥𝑘 ∑ 𝜉𝑗 𝑊𝑗 𝑟𝑗𝑘 − Δ𝑥𝑘+1 ∑(1 − 𝜉𝑗 )𝑊𝑗 𝑟𝑗,𝑘+1 = 0 𝑗=0

𝑗=0

After collecting the residual terms in Eq. (67) we have: [44]

(67)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

(

𝑛+1

𝑛+1

𝑛+1

𝑛+1

𝑖=0

𝑖=0

𝑗=0

𝑗=0

1 1 ∑ 𝐴𝑛+1,𝑖 𝑦𝑖𝑘 − ∑ 𝐴0,𝑖 𝑦𝑖,𝑘+1 ) − Δ𝑥𝑘 ∑ 𝜉𝑗 𝑊𝑗 𝑅𝑗𝑘 − Δ𝑥𝑘+1 ∑(1 − 𝜉𝑗 )𝑊𝑗 𝑅𝑗,𝑘+1 = 0 Δ𝑥𝑘 Δ𝑥𝑘+1

(68)

Since the interior residuals are zero by Eq. (61), we are left with only interface values, so Eq. (68) is identical to Eq. (64) above. The hybrid-collocation-Galerkin method [Diaz (1975), Dunn and Wheeler (1976), Wheeler (1977)] is equivalent to the methods of Gray (1977), Young (1977), Hennart (1982) and Leyk (1986, 1997), Maday and Patera (1989), Vosse and Minev (1996), Karniadakis and Sherwood (2013). Previously, these methods were considered different. Given this equivalence the convergence and super-convergence properties (discussed below) for the hybrid method apply equally to all. Although this demonstration of equivalence is for a specific problem, it generalizes to other equations, because the linear hat functions are included within the Lagrange basis. Since this one basic procedure has gone by several different names, we ask what is a logical name for it? Is it orthogonal collocation or a Galerkin method with specific quadrature and nodal locations? The relationship between the Galerkin method and orthogonal collocation at Lobatto points has been known since Villadsen and Stewart (1967) and the method is obviously a natural extension. It was called a Lobatto-Galerkin method [Young (1981)], a name which caused subsequent work to reference it as just another run of the mill finite element paper. The other names hp Spectral Element method and G-NI or SEM-NI (Galerkin or Spectral Element Method with Numerical Integration) also fail to highlight the method’s most important feature, viz. that quadratures are not required. Whatever name is used it should include the word “collocation”, which is synonymous with “no integration required”. The method is equivalent to collocation. Even the element interface condition corresponds to setting residuals to zero, see Eq. (64). By all rights it should be called the C0 Orthogonal Collocation Finite Element Method, while the alternate collocation at Gauss points should be called the C1 Orthogonal Collocation Finite Element Method. These names include the important word “collocation” and honor the original one given by Villadsen and Stewart. These names will be used for the remainder of this article. 7.4 Convergence Rate and Efficiency Several theoretical convergence results have been proven for the Orthogonal Collocation Finite Element methods. Dupont (1976) has summarized the results for the C0 Galerkin method and the H1 Galerkin method which is equivalent to moments. The convergence and superconvergence rates for the orthogonal collocation methods are the same as for the finite element methods they approximate. The rates for C 1 Orthogonal Collocation are the same as for the H1 Galerkin method or moments, while those for the C0 Orthogonal Collocation method are the same as for a C0 Galerkin method with Lagrange trial functions. The overall convergence rate, e.g. L2 norm, is O(Δxn+2) which is optimal since n+1 is the polynomial degree. However, they also exhibit some superconvergence properties. C1 Orthogonal Collocation converges with O(Δx2n+2) for both the solution and its first derivative, i.e. fluxes, at element boundary nodes [DeBoor and Swartz (1973), Douglas and Dupoint (1973)]. The C0 method exhibits the same overall convergence and superconvergence rates with less stringent [45]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

conditions on smoothness [Diaz (1975), Dunn and Wheeler (1976), Wheeler (1977)]. It should be no surprise that to achieve this rate the flux must be calculated by the method described in section 4.4, i.e. using Eq. (60) evaluated at the element boundary. It has also been shown that for the C0 method the solution at the Lobatto points within the element is O(Δxn+3) [Nakao (1984)]. These estimates are based on certain assumptions which will not apply in all cases, but they indicate the potential accuracy of these methods. The convergence and superconvergence results have been demonstrated in some of the references above and by Carey, et al. (1981). Due to the equivalence demonstrated, these convergence rates apply to all the C0 collocation methods discussed in section 7.3, regardless of the name given. These convergence rates are in terms of element size, Δx. In hp finite element terminology, this refinement is call h refinement. The other alternative is to increase the degree of the polynomial trial function, i.e. p refinement. Maday and Patera (1989) analyzed this alternative and confirmed exponential convergence is achieved with the C0 Orthogonal Collocation or hp spectral element method. Actually, an exponential convergence rate is implied by the O(Δxn+2), since for increasing n, the rate is not proportional to any fixed power of Δx. Apparently, the superconvergence properties of this method have been overlooked in the spectral literature. The improved efficiency of collocation at Gauss points relative to a conventional Galerkin method is well documented [Houstis, et al. (1978), Finlayson (1979), Frind and Pinder (1979)]. A Hermite cubic or higher order grid like Fig. 19 would likely be more efficient than explicitly imposing continuity with Eq. (58), especially in multiple dimensions. These trial functions reduce the dimension of the problem substantially, but there is some loss of efficiency since the nodes and collocation points no longer coincide. Placing single nodes at Gauss points when possible, as in Fig. 19, should realize some economy since at least some of the nodes and collocation points are coincident. Collocation at Gauss points saves some computations relative to integrated moments, but the savings are not nearly as great as with global methods. These methods are easily extended to multiple dimensions. A two-dimensional extension of C1 Orthogonal Collocation was described by Prenter and Russell (1976). The method has also been extended to distorted grids and has found many applications in water resource modeling [Pinder, et al. (1978), Frind and Pinder (1979), Lapidus and Pinder (1999)]. However, the element distortion is more restrictive than with a conventional C0 isoparametric finite element method. When the same trial functions are used, both finite element moments and collocation at Gauss points provide substantial savings over the Galerkin [46]

Fig. 21 Domains of influence, - - - - Moments or Collocation, ______ Galerkin from Young (1977)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

method, especially in multiple dimensions. Fig. 21 shows the domains of influence for Hermite cubic trial functions in two dimensions. Since the weight functions are local to each element, the approximating matrix has only 16 nonzero entries per row for collocation or moments rather than 36 for the Galerkin method. One study found collocation to be just as accurate as the Galerkin method [Frind and Pinder (1979)]. Another comparison for a convective-diffusion problem found that a Galerkin method is more accurate than collocation but only slightly more accurate than moments [Young (1977)]. When the computational complexity was considered, moments was the most attractive of these methods. The moments method warrants more consideration than it has received in the past. The C0 orthogonal collocation method also extends easily to multiple dimensions using tensor product Lagrange trial functions. It has the following attractive computational features: 1. no quadratures to perform 2. nodes and quadrature points coincide 3. lumped mass or capacitance matrix 4. symmetric matrices for self-adjoint operators 5. can be used with conventional isoparametric elements 6. computational molecules with star structure for rectangular grids Many of these features are lost in the hybrid formulation. Since the nodes and quadrature points coincide, no interpolations are required to calculate coefficients. The lumped mass matrix permits the use of explicit time stepping methods. More efficient solution techniques can be used for symmetric matrix problems. Fig. 22 shows the computational molecules for grids of rectangular elements when no cross terms appear in the differential equation. The linear case is equivalent to the 5-point finite difference method. With quadratic, cubic and quartic elements the average number of points in the difference equations are 7, 9 and 11, respectively. The cubic method is a true fourth order method with the same number of neighboring nodes as a 9-point finite difference or linear finite element method. Several studies have found this procedure to be very efficient relative to conventional Galerkin finite element methods with full integration [Gray (1977), Young (1977), Melenk, et al. (2001), Schillinger, et al. (2015)].

Fig. 22 Difference stencils in 2D - - - - C0 Orthogonal Collocation ______ C0 Galerkin from Young (1977)

[47]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

In one study, the computational characteristics of the C0 orthogonal collocation method were found to be effective for solving convective-diffusion problems in two spatial dimensions with a combination of explicit and implicit time stepping [Young (1978,1981)]. Explicit time stepping was efficient due to lumping and to the small computational molecules. The implicit equations were solved with Line Successive Overrelaxation (LSOR) with an optimum relaxation parameter calculation [Varga (2009)]. Compared to the linear case, equivalent to 5-point finite differences, the iterations required per node are roughly 10, 20 and 25 percent greater for quadratic, cubic and quartic elements, respectively. Due to the difference stencils shown in Fig. 22, the computational effort per node with cubics was roughly equivalent to those for a 9point finite difference or linear finite element method. However, the accuracy per node is substantially greater and produces a method which is more efficient overall. The method was also compared to a cubic Hermite Galerkin method and it gave similar accuracy with far less calculation. This code was developed as a true hp finite element method and refinements were made both by increasing the number of elements and by increasing the degree of the trial functions. However, quartic and higher order elements were no more efficient than cubics. The C0 orthogonal collocation finite element method can be used with conventional isoparametric elements [Gray (1977), Gray and Van Genuchten (1978) ), Karniadakis and Sherwood (2013)] in two and three dimensions. It has also been extended to triangular and tetrahedral elements providing even greater flexibility [Dubiner (1991), Sherwin and Karniadakis (1995), Cohen, et al. (2001), Mulder (2001)]. All of these developments are based on Villadsen and Stewart’s idea of placing collocation points at quadrature base points.

Summary Orthogonal Collocation, Pseudospectral and Differential Quadrature methods are the same, so this review and its conclusions applies to all. Lanczos’ 80 year old idea was to collocate at the roots of Chebyshev orthogonal polynomials. His idea was refined by Villadsen and Stewart’s 51-year-old idea to collocate at Lobatto or Gaussian quadrature base points, also the roots of orthogonal polynomials. The idea of collocating at quadrature base points has proven useful for extending the method to finite element approximations. This retrospective: 1. Reviews development of the method in the context of numerical integration of the Methods of Weighted Residuals; 2. Shows equivalence of: moments and tau methods, accurately approximated by collocation at Gauss points; 3. Illustrates a simple modified formulation which gives more efficient symmetric matrix problems for self-adjoint operators with Gauss, Lobatto or Radau points; 4. Provides compelling support for the correct, weak or natural formulation of flux boundary conditions with Lobatto and Chebyshev points;

[48]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

5.

6.

Shows how the idea of collocating at quadrature points extends naturally to: a. C1 Orthogonal Collocation Finite Element method – collocating at Gauss points b. C0 Orthogonal Collocation Finite Element method – collocating at Lobatto points Demonstrates equivalence of the C0 Orthogonal Collocation Finite Element method with: a. Hybrid-Collocation-Galerkin method of Diaz (1975), Dunn and Wheeler (1976), and Wheeler (1977) b. the methods of Gray (1977), Young (1977, 1981), Hennart (1982), and Leyk (1986, 1997) c. hp Spectral Element method of Maday and Patera (1989), Vosse and Minev (1996), Canuto, et al. (2007), Karniadakis and Sherwood (2013)

Notation a – trial function coefficients, solution coefficients Eq. (1)

A – first derivative differentiation matrix defined by Eq. (20) B – 2nd derivative differentiation matrix defined by Eq. (20) Bi – Biot number, dimensionless external mass transfer coefficient c – coefficients of mixing cup average composition, Eq. (51) C – stiffness matrix defined by Eqs. (29) and (34) d – coefficients of flux, Eq. (52) D – mass matrix defined by Eq. (28) f – generic function g – generic function ℓ - Lagrange interpolating polynomial defined by Eq. (2) ℓ* - Lagrange polynomial through interior points only, Eq. (25) m – number of quadrature points n – number of interior points ne – number of elements P – Jacobi polynomial r – reaction rate, source function r̂ - normalized source function R – differential equation residual Sh – Sherwood number, dimensionless internal mass transfer coefficient, Eq. (43) v – eigenfunctions w – weight or test function in MWR W – quadrature weight x – coordinate y – dimensionless composition y* - exact solution y̅ - mixing cup average composition, Eq. (42) z – axial coordinate Greek Letters [49]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

α– parameter in Jacobi polynomial weight function, Eq. (4) β – parameter in Jacobi polynomial weight function, Eq. (4) δ – Dirac delta or Kronecker delta ϵp – definition of Lp error, Eq. (23) ϵη – error in η, Eq. (24) η – effectiveness factor, normalized flux, Eqs. (13) or (14) λ – eigenvalue ξ – local element coordinate τ – coefficient in residual, Eq. (9) φ – Thiele modulus ψ – generic trial or basis function

References Ames, W. F., Nonlinear partial differential equations in engineering, Academic Press (1965) Anderson, E., Bai, Z., Bischof, C., Blackford, S., Demmel, J., Dongarra, J., Du Croz, J., Greenbaum, A., Hammarling, S., McKenney, A., Sorensen, D., LAPACK Users' Guide (Third ed.) SIAM, Philadelphia, PA (1999). Bellman, R. and Casti, J.: “Differential Quadrature and Long-Term Integration,” J. Math. Anal. Appl., 134, 235-238 (1971). Bellman, R., Kashef B.G., and Casti, J.: “Differential Quadrature: A Technique for the Rapid Solution of Nonlinear Partial Differential Equations,” J. Comp. Physics 10, 40-52 (1972). Bellman, R., Methods of Nonlinear Analysis, 2 Academic Press, New York (1973). Bellomo, N.: “Nonlinear Models and Problems in Applied Sciences from Differential Quadrature to Generalized Collocation Methods,” Math. Comput. Modelling, 26, 4, 13-34 (1997). Bernardi, C. and Maday, Y.: “Spectral Methods”, in Handbook of Numerical Analysis Vol. 5, in Techniques of Scientific Computing (Part 2), Ciarlet & Lions (ed.) (1997). Bert, C.W. and Malik, M.: “Differential Quadrature Method in Computational Mechanics: A Review,” Appl. Mech. Review, 49, 1 (1996). Bird, R.B., Stewart, W.E. and Lightfoot, E.N., Transport Phenomena, John Wiley and Sons, (1960). Boyd, J.P., Chebyshev and Fourier Spectral Methods, 2nd ed., Dover Publications (2000). Bubnov, I.G., “Report on the works of Professor Timoshenko which were awarded the Zhuranskyi Prize,” Symposium of the Institute of Communication Engineers, No. 81, All Union Special Planning Office (1913). Canuto, C. and Quarteroni, A.: “Spectral and Pseudo-spectral Methods for Parabolic Problems with Nonperiodic Problems Boundary Conditions,” Calcolo, 18, 3 197-217 (1981). Canuto, C., Boundary Conditions in Chebyshev and Legendre Applications, NASA Contract Report 172469 (1986) [50]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Canuto, C., Hussaini, M., Quarteroni, A. and Zang, T., Spectral Methods in Fluid Dynamics, Springer, Berlin (1988) Canuto, C., Hussaini, M., Quarteroni, A. and Zang, T., Spectral Methods: Fundamentals in Single Domains, Springer, Berlin (2006) Canuto, C., Hussaini, M., Quarteroni, A. and Zang, T., Spectral Methods Evolution to Complex Geometries and Applications to Fluid Dynamics, Springer, Berlin (2007) Carey, G. and Finlayson, B.A., “Orthogonal Collocation on Finite Elements”, Chem. Engr. Sci. 30, 587-596 (1975). Carey, G., Humphrey, D., and Wheeler, M. F., “Galerkin and Collocation-Galerkin Methods with Superconvergence and Optimal Fluxes,” International Journal for Numerical Methods in Engineering, 17, 6, 939-950 (1981). Ciarlet, P.G. and Raviat, P.A.: “The Combined Effect of Curved Boundaries and Numerical Integration in Isoparametric Finite Element Methods,” in The Mathematical Foundations of the Finite Element Method with Applications to Partial Differential Equations, A.K. Aziz, ed., Academic Press, New York, pp. 409–474 (1972). Ciarlet, P.G. The Finite Element Method for Elliptic Problems, North-Holland, Amsterdam, (1978). Clenshaw, C.W. and Curtis, A.R., “A Method for Numerical Integration and an Automatic Computer,” Numerische Mathematik, 2, 197 (1960). Clenshaw, C. W. and Norton, H. J., "The Solution of Nonlinear Ordinary Differential Equations in Chebyshev Series," Comp. J., 6, 88-92 (1963). Cohen, G., Joly, P. , Roberts, J. E. and Tordjman, N.: “Higher Order Triangular Finite Elements With Mass Lumping for the Wave Equation,” SIAM J. Numer. Anal., 38, 6, 2047-2078 (2001). Collatz, L., The Numerical Treatment of Differential Equations, Springer, 1960 Courant, R., remarks on “weighted averages of the residual”, in discussion following - Biezeno, C.B. “Graphical and Numerical Methods for Solving Stress Problems,” Proc. 1 st Intl. Congress Appl. Mech., Delft (1924), p. 17. Crandall, S. H., Engineering Analysis, McGraw-Hill, (1956). Davis, M.E.: Numerical Methods and Modeling for Chemical Engineers, John Wiley, New York (1984) DeBoor, C. and Swartz, B., “Collocation at Gaussian Points”, SIAM J. Numer. Anal., 10 582– 606 (1973). Diaz, J., A Hybrid Collocation-Galerkin Method for Two-point Boundary Value Problems Using Continuous Piecewise Polynomial Spaces, Ph.D. Thesis, Rice Univ. (1975) Douglas, J., Jr., Dupont, T., “A Finite Element Collocation Method for Quasilinear Parabolic Equations”, Math. Comp., 27 17-28 (1973). [51]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Douglas, J., Jr., Dupont, T. and Wheeler, M. F., “H1-Galerkin Methods for the Laplace and Heat Equations”, Math. Aspects of Finite Elements in Partial Differential Equations, Carl De Boor, Ed., Academic Press, 383 (1974). Dubiner, M.: “Spectral Methods on Triangles and Other Domains,” J. of Sci. Comp., 6, 4, 345– 390 (1991). Dunn, R., and Wheeler, M. F., “Some Collocation-Galerkin Methods for Two-Point Boundary Value Problems,” SIAM J. Numer. Anal., 13, 5, 720-733 (1976). Dupont, T., “A Unified Theory of Superconvergence for Galerkin Methods for Two-Point Boundary Problems”, SIAM J. Numer. Anal., 13, 3, 362–368 (1976). Dyksen, W.R., Houstis, E. N., Lynch, R. E. and Rice, J. R.: “The Performance of the Collocation and Galerkin Methods with Hermite Bi-cubics,” SIAM J. Num. Anal., 21 4, 695–715 (1984). El Daou, M., Ortiz, E.L., Samara, H.: “A Unified Approach to the Tau Method and Chebyshev Series Expansion Techniques,” Computers Math. Applic., 25, 3, 73-82 (1993). El Daou, M., Ortiz, E.L.: "The Tau Method as an Analytic Tool in the Discussion of Equivalence Results Across Numerical Methods," Computing, 60, 365–376 (1998). Elnashaie, S.S.E. and Cresswell, D.L., “On Dynamic Modelling of Porous Catalyst Particles,” Chem. Engr. Sci. 28, 1387 (1973). Faddeeva, V.N., Computational Methods of Linear Algebra, Dover Publications, New York (1959). Ferguson, N.B., Orthogonal Collocation as a Method of Analysis in Chemical Reaction Engineering, Ph.D. thesis, University of Washington, Seattle, WA (1971). Ferguson, N.B. and Finlayson, B.A., “Error Bounds for Approximate Solutions to Nonlinear Ordinary Differential Equations,” Amer. Inst. Chem. Engr. J., 18 (5), 1053-1059 (1972). Finlayson, B.A. and Scriven, L.E., “The Method of Weighted Residuals – A Review,” Appl. Mech. Review. 19, 9, 735-748 (1966). Finlayson, B.A., The Method of Weighted Residuals and Variational Principles, Academic Press, New York, NY (1972). Finlayson, B.A., “Orthogonal Collocation in Chemical Reaction Engineering,” Cat. Rev. SciEng. 10, 69-138 (1974). Finlayson, B.A.: “Orthogonal Collocation on Finite Elements – Progress and Potential,” Advances in Computer Methods for Partial Differential Equations-III, R. Vichnevetsky (ed.), IMACS(AICA), Rutgers U., New Brunswick, N.J., 107-113 (1979) Finlayson, B.A., Nonlinear Analysis in Chemical Engineering, Ravenna Park Publ., Seattle (2003).

[52]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Fix, G.J. “Effect of Quadrature Errors in the FInite Element Approximation of Steady State, Eigenvalue and Parabolic Problems,” in The Mathematical Foundations of the Finite Element Method with Applications to Partial Differential Equations, A.K. Aziz (ed.), Academic Press, pp. 525–556 (1972). Franklin, J., Matrix Theory, Prentice-Hall, Englewood Cliffs, NJ (1968). Frazer, R. A, Jones, W. P., and Skan, S. W.: "Approximation to Functions and to the Solutions of Differential Equations," Great Britain Aero. Res. Council, London, Report and Memo No. 1799, (1937). Fried, I.: “Numerical Integration in the Finite Element Method,” Computers & Structures, 4, 5, 921–932 (1974). Fried, I. and Malkus, D.S.: “Finite Element Mass Matrix Lumping by Numerical Integration with No Convergence Rate Loss,” Intl. J. Solids and Structures, 11, 4, 461–466 (1975). Frind, E.O. and Pinder, G.F.: “A Collocation Finite Element Method for Potential Problems in Irregular Domains,” Int. J. Num. Methods Engr., 14 681-701, (1979) Fujita, H.: ‘Diffusion with a Sharp Moving Boundary,” J. Chem. Phys. 21, 700-705 (1953). Funaro, D., Polynomial Approximation of Differential Equations, Lecture notes in Physics, m8, Springer, Heidelberg (1992). Galerkin, B.G.: “Series Solution of Some Problems Related to Elastic Equilibrium of Rods and Plates,” Vestn. Inzhenerov Tech., 19, 897-908 (1915). Golub, G.H. and Welsch, J.H., “Calculation of Gaussian Quadrature Rules”, Mathematics of Computation, 23, 221-230 (1969). Gottlieb, D. and Lustman, L., “The Dufort-Frankel Method for Parabolic Initial Boundary Value Problems”, Computers and Fluids, 11, 2, 107-120 (1983). Gottlieb, D. and Orzag, S.A.: Numerical Analysis of Spectral Methods: Theory and Applications, SIAM, Philadelphia, PA (1977). Gray, W.G., “An Efficient Finite Element Scheme for Two-Dimensional Surface Water Computations”, Finite Elements in Water Resources, W.G. Gray, G.F. Pinder and C.A. Brebbia (ed.), Pentech Press, London (1977). Gray, W.G. and Van Genuchten, M.T.: “Economical Alternatives to Guassian Quadrature Over Isoparameteric Quadralaterals,” Intl. J. Num. Meth. Eng. 12 1478-1484 (1978). Hairer, E., Nörsett, S.P. and Wanner, G., Solving Ordinary Differential Equations I: Nonstiff Problems, 2nd ed., Springer Verlag, Berlin (1993). Hairer, E. and Wanner, G., Solving Ordinary Differential Equations II: Stiff and DifferentialAlgebraic Problems, 2nd ed., Springer Verlag, Berlin (1996). Hennart, J.-P.: “Topics in Finite Element Discretization of Parabolic Evolution Problems,” Lecture Notes in Math. 909, Springer-Verlag, Berlin, Heidelberg, (1982). [Cohen, et al. (2001) [53]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

list an earlier reference which could not be located: Hennart, J.-P., Sainz, E., and Villegas, M.: “On the Efficient Use of the Finite Element Method in Static Neutron Diffusion Calculations,” Computational Methods Nuclear Engrg., 1, 3.87 (1979)] Hildebrand, F.B., Introduction to Numerical Analysis, 2nd ed., Dover Publications, New York, NY (1987). Houstis, E.N., Lynch, R.E., Rice, J.R. and Papatheodorou, T.S.: “Evaluation of Numerical Methods for Elliptic Partial Differential Equations,” J. Comp. Phys., 27, 3, 323–350 (1978). Hughes, T.J.R., The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Prentice-Hall (1987). Johnson, D.: Chebyshev Polynomials in the Spectral tau Method and Applications to Eigenvalue Problems, NASA Contractor Report No. 198451 (1996). Kantorovich, L.V.: “On an Approximation Method for the Solution of a Partial Differential Equation”, Dokl. Akad. Nauk SSSR, 2, 532-536 (1934) Kantorovich, L.V. and Krylov, V.I, Approximate Methods in Higher Analysis, Interscience (1958). Karniadakis, G. and Sherwin, S.: Spectral/hp Element Methods for Computational Fluid Dynamics, Oxford Univ. Press (2013). Karpilovskaya, E.B.: “Convergence of the Collocation Method for Certain Boundary Value Problems of Mathematical Physics,” Sibirsk. Matem. Zh., 4, 3, 632-640 (1963). Kravchuk, M.P., “On Krylov’s Method in the Theory of Approximate Integration of Differential Equations,” Tr. Fiz. – Mat. Viddilu VUAN, 5, 2, 12-33 (1926) [reference from Lucka and Lucka (1992)] Kravchuk, M.P.: “Application of the Method of moments to the Solution of Linear Differential and Integral Equations,” (in Ukrainian), Kiev. Soobshch. Akad. Nauk UkSSR, 1, 168 (1932); [referenced in Kantorovich and Krylov (1958), Mikhlin (1964), Shuleshko (1959)]. Křížek, M. and Neittaanmäki, P., “Bibliography of Superconvergence”, in Finite Element Method, Superconvergence, Post-Processing and a Posteriori Estimates, 196, 315-348 Marcel Dekker, New York (1998). Krylov, N.M., “On a Method of Approximate Integration for Which the Ritz Method as Well as the Method of Least Squares are Special Cases,” C.R. Acad. Sci. Paris, 182, 676-678 (1926) (reference from Lucka and Lucka, 1992) Krylov, N.M., Collected Works in 3 volumes (in Ukranian) , Izd. Akad. Nauk. UkrSSR, Kiev, vol. 1-3, (1961) [reference from Lucka and Lucka, 1992]. Krylov, V.I., Approximate Calculation of Integrals, Macmillan, New York, NY (1962), Dover (2012). Lanczos, C., “Trigonometric Interpolation of Empirical and Analytical Functions”, J. Math. Phys., 17, 123-199 (1938). [54]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Lanczos, C., Applied Analysis, Prentice-Hall, Englewood Cliffs, NJ (1956). Lanczos, C., “Legendre Versus Chebyshev Polynomials,” in Topics in Numerical Analysis, Proceedings of the Royal Irish Academy Conference on Numerical Analysis 1972, J.J.H. Miller, ed., Academic Press, London, pp. 191–201 (1973). Lapidus, L. and Pinder, G.F.: Numerical Solution of Partial Differential Equations in Science and Engineering, John Wiley, New York (1999) Leyk, Z., “A Co-Collocation-Like Method for Boundary Value Problems,” Numerische Mathametik, 49, 39-53 (1986). Leyk, Z., “A Co-Collocation-Like Method for Elliptic Equations on Rectangular Regions,” J. Austral. Math. Soc. Ser. B , 38, 368-387 (1997). Lucka, T.F. and Lucka, A.Y.: “Development of Direct Methods in Mathematical Physics in the Works of M.P. Kravchuk,” translated from Ukrain. Math. J., 44, 7, 931-939 (1992). Maday, Y. and Patera, A. T.: “Spectral Element Methods for the Incompressible Navier-Stokes Equations” In State-of-the-Art Surveys on Computational Mechanics, A.K. Noor, editor, ASME, New York (1989). Melenk, J.M., Gerdes, K. and Schwab, C.: “Fully Discrete hp-Finite Elements: Fast Quadrature,” Comp Meth. in Appl. Mech. and Eng., 190, 32, 4339–4364 (2001). Michelsen, M.L. and Villadsen, J.V.: "A Convenient Computational Procedure for Collocation Constants", Chemical Engineering J, 4, 64-68 (1972). Michelsen, M.L. and Villadsen, J.V.: “Polynomial Solution of Differential Equations”, in Foundations of Computer-Aided Chemical Process Design, Proceedings of an International Conference (pp. 341-368), Mah, R.S.H. and Seider, W.D. ed. (1981). Mikhlin, S. G., Variational Methods in Mathematical Physics, Pergamon/Macmillan (1964). Mulder, W.A.: “Higher-Order Mass-Lumped Finite Elements for the Wave Equation,” J. Comput. Acoustics, 9, 2, 671-680 (2001). Nakao, M., “Superconvergence Estimates at Jacobi Points of the Collocation-Galerkin Method for Two Point Boundary Value Problems”, J. Information Processing, 7 31-34 (1984). Nielson, K.L.: Methods in Numerical Analysis, MacMillan, NY, pp. 150-4 (1956). Ortiz, E.L.: “Step by Step tau Method – Part I. Piecewise Polynomial Approximations,” Comput. and Math. Applic., 1, 381-392 (1975). Orzag, S.A., “Comparison of Pseudo Spectral and Spectral Approximation,” Studies in Applied Mathematics 51, 3, 253-259 (1972). Orszag, S.A.: “Spectral Methods for Problems in Complex Geometries,” J. of Comp. Physics, 37 1, 70–92 (1980). Peyret, R.: Spectral Methods for Incompressible Viscous Flows, Springer (2002).

[55]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Picone, M.: “Sul Metodo Delle Minime Potenze Ponderate e Sul Metodo di Ritz per il Calcolo Approssimato nei Problemi Della Fisica-Matematica,” Rend. Circ. Mat. Palermo, 52, 225-253 (1928). Pinder, G.F., Frind, E.O. and Celia, M A.: “Groundwater Flow Simulation Using Collocation Finite Elements,” Proc. Second Int. Conf. Finite Elements in Water Resources, (eds Brebbia et al.), 1.171-1.185, Pentech Press, London, (1978). Prenter, P.M. and Russell, R.D.: ”Orthogonal Collocation for Elliptic Partial Differential Equations”, SIAM J. Numer. Anal., 13, 923 (1976). Quan, J.R. and Chang, C.T.: “New Insights in Solving Distributed System Equations by the Quadrature Method – I. Analysis,” Comput. Chem. Engng., 13, 7, 779-788 (1989). Raviart, P.A.: “The Use of Numerical Integration in Finite Element Methods for Parabolic Equations,” in Topics in Numerical Analysis, Proceedings of the Royal Irish Academy Conference on Numerical Analysis 1972, J.J.H. Miller, ed., Academic Press, London, pp. 233– 264 (1973). Ritz, W.: “Über eine neue Methode zur Lösung gewisser Variationsprobleme der mathematischen”, Physik, J. reine angew. Mathematik 135 1-61 (1908), also collected works, Paris (1911). Runge, M.: “Über Empirische Functionen und die Interpolation Zwischen Aequidistanten Ordinaten,” Zeitschrift für Math. und Physik tome XLVI (1901). Schillinger, D., Evans, J.A., Frischmann, F., Hiemstra, R.R., Hsu, M.-C., and Hughes, T.J.R.: “A Collocated C0 Finite Element Method: Reduced Quadrature Perspective, Cost Comparison with Standard Finite Elements, and Explicit Structural Dynamics,” Intl. J. Num. Meth. Engrg., 102, no. 3-4, 576-631 (2015). Shen,J., Tang,T. and Wang,L.-L., Spectral Methods: Algorithms, Analysis and Applications, Springer-Verlag, Series in Computation Mathematics, vol. 41, Berlin (2011). Sherwin, S.J. and Karniadakis, G.E.: “A New Triangular and Tetrahedral Basis for High-Order (hp) Finite Element Methods,” Intl. J. Num. Meth. in Engr., 38, 22, 3775–3802 (1995). Shu, C., Differential Quadrature and Its Application in Engineering, Springer, London, 2000 Shuleshko, P., “A New Method of Solving Boundary-Value Problems of Mathematical Physics, Generalizations of Previous Methods of Solving Boundary-Value Problems of Mathematical Physics,” Austral. J. Appl. Sci. 10, 1-7, 8-16 (1959). Slater, J. C.: "Electronic Energy Bands in Metal," Phys. Rev., 45, 794-801 (1934). Sorensen, J.P., Guertin, E.W., and Stewart, W.E.: “Computational Models for Cylindrical Catalyst Particles,” Amer. Inst. Chem. Eng. J., 19, 5, 969-975 (1973). Strang, G. and Fix, G.J.: An Analysis of the Finite Element Method, Prentice-Hall, 1973 Trefethen, L.N., Spectral Methods in Matlab, SIAM, Philadelphia, PA (2000). [56]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Trefethen, L.N. "Is Gauss Quadrature Better than Clenshaw-Curtis?". SIAM Review 50, 67–87 (2008). Van de Vosse, F.N. and Minev, P.D.: Spectral elements methods: Theory and applications, EUT Report 96-W-001, Eidenhoven Univ. Tech., The Netherlands (1996) 96-W-001, Eindhoven University of Technology, 1996. Varga, R., Matrix Iterative Analysis, Springer (2009). Villadsen, J.V. and Stewart, W.E., “Solution of Boundary-Value Problems by Orthogonal Collocation,” Chem. Engr. Sci. 22, 1483-1501 (1967). Villadsen J., Selected Approximation Methods for Chemical Engineering Problems, Inst. for Kemiteknik Numer. Inst. Danmarks Tekniske Hojskole (1970). Villadsen, J.V. and Michelsen, M.L., Solution of Differential Equation Models by Polynomial Approximation, Prentice-Hall, Englewood Cliffs, NJ (1978). Wheeler, M.F.: “A C0-Collocation-Finite Element Method for Two-Point Boundary Value and One Space Dimension Parabolic Problems,” SIAM J. Numer. Anal., 14, 1, 71-90 (1977). Wright, K., "Chebyshev Collocation Methods for Equations," Comp. J., 6, 358-365 (1964). Yamada, H., Rept. Res. Inst. Fluid Engng., Kyushu Univ. 3, 29 (1947); referenced by H. Fujita (1953). Also: “A method of approximate integration of the laminar boundary-layer equation,” Rept. Res. Inst. Fluid Engng., Kyushu Univ. 6, 87-98 (1950). Young, L.C. and Finlayson, B.A., “Mathematical Models of the Monolith Catalytic Converter: Part I. Development of Model and Application of Orthogonal Collocation,” Amer. Inst. Chem. Engr. J., 22, 2, 331-343 (1976). Young, L.C., “A Preliminary Comparison of Finite Element Methods for Reservoir Simulation,” Advances in Computer Methods for Partial Differential Equations-II, R. Vichnevetsky (ed.), IMACS(AICA), Rutgers U., New Brunswick, N.J., 307-320 (1977) Young, L.C., “A Finite-Element Method for Reservoir Simulation,” Soc. Petr. Engrs. J. 21(1) 115-128, (Feb. 1981), paper SPE 7413 presented Oct. (1978). Young, L.C., Orthogonal Collocation Revisited, http://tildentechnologies.com/Numerics (2018). Zienkiewicz, O.C., The Finite Element Method in Engineering Science, McGraw-Hill (1971).

Appendix A – Galerkin Method with Flux Boundary Conditions Consider a problem approximated by Lagrange interpolating polynomials as in Eq. (2). The trial functions are substituted into the differential equation to form the residual as in Eq. (7). Assume the boundary conditions are: 𝑦 = 0 at 𝑥 = 0 𝑑𝑦 + 𝛼𝑦 = 0 at 𝑥 = 1 𝑑𝑥 The usual procedure in orthogonal collocation, pseudospectral or differential quadrature methods is to satisfy both of these conditions exactly, i.e. use boundary collocation [Finlayson [57]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

(1972), p. 101; Villadsen and Michelsen (1978), p. 137; Bert and Malik (1996); Bellomo (1997); Trefethen (2000), p. 137; Boyd (2000), p. 111, Peyret (2002), p. 59]: 𝑦0 = 0 𝑛+1

∑ 𝐴𝑛+1,𝑖 𝑦𝑖 + 𝛼𝑦𝑛+1 = 0 𝑖=0

where A is the first derivative matrix defined in Eq. (20). After elimination of the boundary values, the approximation is: 𝑛

𝑦(𝑥 ) = ∑[ℓ𝑖 (𝑥 ) + 𝑏𝑖 ℓ𝑛+1 (𝑥 )] 𝑦𝑖

(A.1)

𝑖=1

where 𝑏𝑖 = −𝐴𝑛+1,𝑖 /(𝐴𝑛+1,𝑛+1 + 𝛼). When the boundary conditions are satisfied exactly, the Galerkin method weights the residual, R, by the trial functions which obey the boundary conditions, i.e. the functions in the brackets of Eq. (A.1). If the Galerkin method integrals are approximated by quadrature using the n + 2 interior and boundary points, the result is: 𝑛+1

∑ 𝑊𝑘 𝑅(𝑥𝑘 , 𝒚)[ℓ𝑖 (𝑥𝑘 ) + 𝑏𝑖 ℓ𝑛+1 (𝑥𝑘 )] = 𝑊𝑖 𝑅(𝑥𝑖 , 𝒚) + 𝑊𝑛+1 𝑏𝑖 𝑅 (1, 𝒚)𝑏𝑖 = 0

(A.2)

𝑘=0

Eq. (A.2) is not equivalent to a collocation method because of the second term involving the residual at the boundary. Collocation sets only the first interior residual term to zero, so the nonzero second term seriously violates the Galerkin method criteria. The boundary weight, Wn+1, is approximately n-2 for Lobatto quadrature, so this inconsistency creates an error over and above that caused by the approximate quadrature. Numerical experiments suggest the error introduced is significant. Since the residual at the boundary decreases exponentially with n, the method with boundary collocation converges exponentially, but usually at a much slower rate, especially for fluxes. Boundary collocation also causes errors in the overall mass or energy balance (see Section 4.4). Boundary collocation works well with Gauss points or Radau points with Wn+1 = 0 (or W0 = 0 if the flux condition is at x = 0). However, Lobatto or Chebyshev points produce greater accuracy with a natural treatment of the boundary conditions as described in Section 4.3. Many articles and texts on pseudospectral and differential quadrature methods state that having a nonzero boundary weight somehow facilitates the approximation of boundary conditions. If one uses boundary collocation, the exact opposite is true.

Appendix B – Links to Available Software Software useful for implementing MWR and collocation can be found at the following internet sites: Funaro (1992) FORTRAN 77 (80 programs, 3500 lines): morespace.unimore.it/danielefunaro/routines/ Don, Somonoff and Costa Fortran 90 PseudoPack high performance software: www.cfm.brown.edu/people/wsdon/PseudoPack.htm Weideman and Reddy (2000) Matlab Differentiation Suite: [58]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

appliedmaths.sun.ac.za/~weideman/research/differ.html Trefethen (2000) software from book: people.maths.ox.ac.uk/trefethen/spectral.html Shen, Tang and Wang (2011) approximately 50 Matlab functions from book: www.ntu.edu.sg/home/lilian/book.htm Gautschi’s (2005) OPQ Suite, orthogonal polynomials and quadrature: www.cs.purdue.edu/archives/2002/wxg/codes/OPQ.html Burkardt (2015), quadrature, numerous other codes in Fortran 90, C++ and Matlab: www.people.sc.fsu.edu/~jburkardt/ Karniadakis and Sherwood (2013), code in C/C++ oriented toward spectral elements www.nektar.info/2nd_edition Finlayson (2003) companion to books, Matlab, Excel, Fortran, Phython, and more www.chemecomp.com Young (2018), Matlab/Octave, Fortran 90 and C++ : www.tildentechnologies.com/Numerics

[59]