Applied Mathematics and Computation 183 (2006) 885–889 www.elsevier.com/locate/amc Superconvergence and the computation of generalized turning points for B.V.P’s Basem S. Attili * UAE University, Department of Mathematics and Computer Science, P.O. Box 17551, Al-Ain, United Arab Emirates Abstract We consider projection methods in general and collocation at gauss points in particular for the numerical computation of generalized turning points for two point B.V.P’s. We discretize the original problem using collocation to obtain a finite dimensional problem. We will employ a direct method for the characterization and computation of simple turning points. Proper extensions of this direct method will be used for the computation of bifurcation points and cubic turning points. We comment also on the superconvergence which will be clear from the numerical examples.  2006 Elsevier Inc. All rights reserved. Keywords: Generalized turning points; Projection methods; Collocation; Superconvergence 1. Introduction We will consider the numerical implementation of the accurate location of generalized turning points for two point boundary value problems. Such problem has received much attention recently. It is a well known fact that the application of conventional numerical methods leads to difficulties near such singularities, see for example, the proceedings edited by Kupper, Mittelman and Weber [10], Keller and Antman [7], Keller [8] and the review by Stakgold [13]. The purpose of this paper is to present a general discussion on the projection methods applied to such singular points and in particular to consider the application of collocation at gauss points, to give an idea about the systems to be solved and to show that collocation proved to be highly accurate and superconverges to the singular point. Related work on the discretization can be found in, for example, Attili [1,4] where finite differences was used, Weber [14] had used multiple shooting for bifurcation problems only, Brezzi, Rappaz and Raviart [5] considered Rayleigh-Ritz method for simple turning points and Kikuchi [9] used finite elements. The problem under consideration is of the form F ðy; kÞ ¼ 0; * Corresponding author. E-mail address: b.attili@uaeu.ac.ae 0096-3003/$ - see front matter  2006 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2006.06.036 ð1Þ 886 B.S. Attili / Applied Mathematics and Computation 183 (2006) 885–889 where F:Rn · R ! Rn ; and F is a C3–function. Here F(y, k) = 0 is the finite dimensional discretization of the two point B.V.P through, say, the use of collocation. We assume the codimension of Range(Fy) is no more than one on Rn and Fk 62 Range(Fy) when the codim. = 1. It follows then (Attili [1,3]) that the solution manifold F1(0) is a p-dimensional (p = 1) and the singular manifold is an n + 1  p dimensional. As a result, N(Fy) = Span{/0} and RangeðF y Þ ¼ fy 2 Rn ; wT0 y ¼ 0g. To characterize the singular manifold, choose T, r 2 Rn so that the matrix   Fy r A¼ ; ð2Þ TT 0 is nonsingular near (y0, k0) the singular point. This can be the case at the simple turning point if we choose r 62 Range(Fy) and TT/0 = 1. Next we solve the systems    
v 0 A ¼ ; ð3Þ uT gT  A ¼ ð0; 1Þ g 1 for u,v and g. Lemma 1. Both systems in (1.3) are uniquely solvable and (a) g = uTFyv, (b) g0 ¼ uT F 0y v  ðuT r0 gT þ gT T 0T vÞ, where the prime denotes differentiation in y and k. The proof is straight forward and can be seen from the discussion presented before the statement of the lemma. Note that in applying this procedure, we normally choose r and TT to be constants which reduces (b) of the lemma to gy = uTFyyv and gk = uTFykv. Now to numerically compute the generalized turning point, we solve the extended system F ðy; kÞ ¼ 0; gðy; kÞ ¼ 0: ð4Þ Note also that solving this system means finding the intersection between the solution manifold represented by F(y, k) = 0 and the singular manifold represented by g(y, k) = 0. Newtons method can be used to solve the system in (4). To find the gradient of g, we use the formula gy  uT ðF y ðy þ v; kÞ  F y ðy; kÞÞ=; ð5Þ which is based on the exact formula for derivative of g. The system in (4) is regular if and only if ry gðy 0 ; k0 Þv; ð6Þ is of full rank, see Attili [1,3]. To consider the bifurcation case, we will assume Fk 2 Range(Fy). As a result N(Fy, Fk) will be spanned by two vectors and hence we will have the following lemma Lemma 2. Let F ~y ¼ ðF y ; F k Þ; ~y ¼ ðy; kÞ and if Fy has a one dimensional null space spanned by /0 and wT0 F k ¼ 0. Then N(Fy) = Span{U0, U1} where U0 = (/0, 0) and U1 ¼ ð^y ; cÞ where ^y 2 RangeðF y Þ and c 2 R such that F y ^y þ cF k ¼ 0. This lemma means that v and g will each have two components. In general, in the case of bifurcation the system to be solved is not (1) but rather a perturbed system to unfold the singularity; that is, F(y, k) + gr = 0; g 2 R and , r 62 Range(Fy, Fk) and g will be identically zero at the bifurcation point. Also, to guarantee the unique solvability of the extended system, a similar condition like the one in (6) must hold and it will be a symmetric 2 · 2 matrix, see Attili [1,2]. In the cubic turning point case, the problem will be a two parameter one; that is, F ðy; k; lÞ ¼ 0: ð7Þ To compute such singular points will require the extension twice to produce F ðy; k; lÞ ¼ 0; gðy; k; lÞ ¼ 0; g~ðy; k; lÞ ¼ 0; ð8Þ since the cubic turning point (y0, k0, l0) for F(y, k, l) = 0 with respect to one parameter corresponds to a simple turning point for (4) with respect to the other one. More details on these ideas can be found in Attili [1,4] and Spence and Werner [12]. B.S. Attili / Applied Mathematics and Computation 183 (2006) 885–889 887 The outline of the remaining sections of this paper will be as follows. In Section 2, we will present the discretization of the B.V.P’s through projection methods, together with some convergence results. Finally, in Section 3 we will consider the details of the collocation and numerical examples. 2. Discretization using projection methods Let us consider the application of the results obtained in the previous section to a simple problem which is sufficient to illustrate the results, y 00 þ f ðy; kÞ ¼ 0; a 6 x 6 b; yðaÞ ¼ yðbÞ ¼ 0: ð9Þ To find v and g the problem to be solved is of the form Lv þ gr ¼ v00 þ f ðy; kÞv þ gr ¼ 0; vðaÞ ¼ vðbÞ ¼ 1 T T v ¼ vð1=2Þ ¼ 1; ð10Þ where L:C2[a, b] \ B ! C[a, b] and B = {y:y(a) = y(b) = 0}. While to find u and g, we solve L u þ gT  ¼ 0; u r ¼ 1; ð11Þ where L* is the adjoint of L. Now, to solve (9)–(11) numerically, we must deal with the finite dimensional forms of such systems; that is, the discretized form. To do so through the projection methods, we proceed n n as follows. Let fhi gi¼1 be a set of C2–functions satisfying the boundary conditions and fli gi¼1 be a set of linear 00 functionals on C[a, b] such that the matrix {łi(ki)} is nonsingular where ki is such that k i ¼ hi . Define Pn to be a n n 2 projection operator: Pn that is, a linear operator which maps R ! R and idempotent ðP n ¼ P n Þ. Also, Pnt = w if and only if w ¼ i¼1 ci k i and li(w) = li(t); i = 1, 2, . . . , n. Here {hi} can be chosen as piecewise polynomial splines, for more details, one can refer to Reddien [11] and De Boor and Swartz [6]. With the choice of the 00 linear functionals to be point interpolator, then v00 + fv = gr will P be approximated by v + Pnfvn = gPnr (projection approximation) where vn = Pnv. Now taking vn ¼ ai hi , then to approximate vn,un and gn, with   Ln P n r An ¼ ; ð12Þ T 0 we solve the systems         fai g 0 fbi g 0 An ¼ and ATn ¼ : gn 1 gn 1 ð13Þ To determine P n un we need the following lemma Lemma 3. If {bi} is as given in (2.5), then P n un ¼ P bi li . The use of this lemma is in the computation of gn since gn ¼ hun ; Ln vn i ¼ hun ; P n Lvn i ¼ hP n un ; Lvn i: ð14Þ This equation implies that using collocation (to be explained later), gn is the weighted residual for PnLvn at the collocation points. Now to show that the systems in (13) are uniquely solvable; that is, A1 n exists, the following lemma guarantees that and is a straight forward result of the definition of An in (12) Lemma 4. There exists a neighborhood D of (y0, k0) so that An(x) ! A(x) uniformly on D as n ! 1. 1 1 A direct result of this lemma is that A1 uniformly on D as n ! 1. n exists for n large enough and An ! A The following error estimates follow easily from this lemma. 888 B.S. Attili / Applied Mathematics and Computation 183 (2006) 885–889 Lemma 5. For all (y, k) in some neighborhood D of (y0, k0) there exists a constant M such that kv  vn k þ kg  gn k 6 MkðI  P n Þvk and ku  un k þ kg0  g0n k ! 0 as n ! 1, where the prime denotes differentiation. T T T T 1 1 The proof is straightforward since ðvn ; gn ÞT ¼ A1 n ð0; 1Þ ; ðv; gÞ ¼ A ð0; 1Þ , ðun ; g n Þ ¼ An ð0; 1Þ and 1 T * (u, g) = A (0, 1). Similar results were given in Reddien [11] and for more on the proofs the reader can refer to it. Now let us consider the convergence of the critical parameter k. To do so, expand F(yn, kn) about (y0, k0) the turning point, we obtain, F ðy n ; kn Þ ¼ F ðy 0 ; k0 Þ þ F y ðy 0 ; k0 Þðy n  y 0 Þ þ F k ðy 0 ; k0 Þðkn  k0 Þ þ Oðkðy n ; kn Þ  ðy 0 ; k0 ÞkÞ: Left multiplying by u0 ð15Þ and simplifying, we obtain u F ðy n ; kn Þ þ Oðkðy n ; kn Þ  ðy 0 ; k0 ÞkÞ kn  k0 ¼ 0 u0 F k ðy 0 ; k0 Þ ð16Þ which will imply 2 jkn  k0 j 6 MðkðI  P n Þu0 k  kF ðy n ; kn Þk þ kðy n ; kn Þ  ðy 0 ; k0 Þk Þ ð17Þ since u0 F n ðy; kÞ ¼ u0 ðI  P n ÞF ðy n ; kn Þ. Similar results were given for Raleigh-Ritz method and simple turning points, for more details refer to Brezzi, Rappaz and Raviart [5]. A very important note here is that if kðI  P n Þu0 k ! 0 as n ! 1 this will imply superconvergence. The following theorem gives error estimates for the simple turning point case. Theorem 6. Let (y0, k0) be a simple turning point for (2.1) and g(y0, k0) = 0. Then there exists an integer N and a neighborhood D of (y0, k0) in C2 · R such that the solutions (yn, kn) to PnF = 0 and gn = 0 exist and are unique in D for all n P N. Moreover, there exists a constant M such that kðy n ; kn Þ  ðy 0 ; k0 Þk 6 MðkP n y 000  y 000 k þ kP n v000  v000 kÞ; where v0 is the solution to (2.2). Note that the norm on the left is taken in C2[a, b] · R and on the right in C[a, b]. The case of bifurcation is similar to that of turning points. It will not be considered here and some related work can be found in Weber [15]. 3. Collocation at gauss points and numerical examples In this section, we will describe the collocation at Gaussian points. For that reason let us work with a more general two point boundary value problem of the form m1 X ek ðs; kÞuðkÞ ðsÞ ¼ f ðs; kÞ; a 6 s 6 b; k 2 R; ð18Þ LðuÞ ¼ uðmÞ ðsÞ þ k¼0 subject to the m– linearly independent homogenous boundary conditions m1 X ½aik uðkÞ ðaÞ þ bik uðkÞ ðbÞ ¼ 0; 1 6 i 6 m; ð19Þ k¼0 where aik and bik are constants. To describe our collocation scheme we start by defining Pn:a = s0 < s1 <  < sn = b to be a partition of [a, b], P(Pn, k, p): the family of Cp[a, b] functions which are polynomials of degree k on each subinterval of Pn and P 0 (Pn, k, p): the subfamily of P(Pn, k, p) consisting of those functions which satisfy Eq. (19). dn Now after selecting a partition Pn we use it to choose the collocation points ðsi Þi¼1 ; that is, we choose d points per subinterval by picking 1 6 r1 < r2 <    < rd 6 1, as the zeros of the d–th Legendre polynomial (the points used in standard Gauss-quadrature rule). Then set B.S. Attili / Applied Mathematics and Computation 183 (2006) 885–889 sidþm ¼ siþ2 þ siþ1 þ rm ðsiþ2  siþ1 Þ ; 2 m ¼ 1; 2; . . . ; d; 889 i ¼ 0; 1; . . . ; n  1: This choice of collocation points gives an approximation of optimal order. Such phenomenon is called superconvergence, see DeBoor and Swartz [6]. Note that we have nd + 1 collocation points. We seek a function un(s) in P 0 (Pn, m + d, m) which satisfies the two point boundary value problem exactly at dn collocation points ðsi Þi¼1 . One of the natural representations of un(s) is of the form un ðsÞ ¼ c1i ðs  siþ1 Þ mþd þ c2i ðs  siþ1 Þ mþd1 þ    þ cimþd ðs  siþ1 Þ þ cimþdþ1 ; for s 2 [si, si+1]; 0 6 i 6 n  1. This choice means that un(s) is uniquely determined by the n(m + d + 1) coefficients cji ; i ¼ 0; 1; 2; . . . ; n  1 and j = 1, 2, . . . , m + d + 1. These coefficients are calculated from the m boundary conditions, nd + 1 collocation conditions and (m + 1)(n  1) continuity conditions which add up to n(m + d + 1) conditions. This will produce a system of n(m + d + 1) equations in n(m + d + 1) unknowns which is used in the systems in (13) to obtain u, v and g. Then the extended system is formed which when solved gives the singular point. To obtain the partials of g the difference formula given in (5) is used. For numerical experimentation, we used Mu þ k expðuÞ ¼ 0; uð0Þ ¼ uð1Þ ¼ 0; which has a simple turning point at the critical parameter k = 3.513830714. Using collocation at gauss points with quintic splines with meshes of size h = 1/2,1/4 and 1/8 the values of k obtained were 3.5137066, 3.5138306 and 3.5138307. It is clear that such method is highly accurate and the error is O(h6). To illustrate the bifurcation case. Consider the problem y 00 þ 3:9 exp½1=ð1 þ kyÞ ¼ 0; yð0Þ ¼ yð1Þ ¼ 0: Again superconvergence was obtained by using collocation at gauss points with quartic splines. Computing the quantity jkh  k0j for step sizes h = 1/4, 1/5 and 1/6, the results obtained were respectively .1592–4, .3772–5 and .1313–5 with estimated order of convergence O(h6). As can be seen from these examples that accuracy obtained is very high with relatively large step size h. This implies that the extended systems which are solved at each iteration are very small in size. Which in turn reduces the amount of work needed to compute such singular points accurately. References [1] B.S. Attili, Numerical computation of cusp singularities for operator equations, J. Dynam. Syst. Geo. Theory 2 (2004) 57–64. [2] B. Attili, Continuation method for computing certain singular points for index-1 differential algebraic equations, Appl. Math. Comput. 175 (2006) 1026–1038. [3] B. Attili, A direct method for the characterization and computation of bifurcation points with corank 2, Computing 47 (1992) 149– 159. [4] B. Attili, Efficient finite differences for parameter-dependent singularities in two-point boundary value problems, Comput. Meth. Appl. Mech. Eng. 190 (2001) 5543–5549. [5] F. Brezzi, T. Rappaz, P. Raviart, Finite dimensional approximation of nonlinear problems. Part I. branches of nonsingular solutions, Numer. Math. 36 (1981) 1–28. [6] C. De Boor, B. Swartz, Collocation at Gaussian points, SIAM J. Numer. Anal. 10 (1973) 582–606. [7] J. Keller, B. Antman, Bifurcation Theory and Nonlinear Eigenvalue Problems, Benjamin, New York, 1969. [8] J. Keller, Numerical solution of bifurcation and nonlinear eigenvalue problems, in: P.H. Rabinwitz (Ed.), Application of Bifurcation Theory, Academic Press, New York, 1977, pp. 359–384. [9] F. Kikuchi, An iterative finite element scheme for bifurcation analysis of semi-linear elliptic equations, Institute of Space and Aeronautical Science, University of Tokyo, Report No. 542, vol. 41, no.6, 1976. [10] T. Kupper, H.D. Mittelmann, H. Weber, Numerical Methods for Bifurcation Problems, Isnm 70, Birkhauser, Basel, 1984. [11] G.W. Reddien, Projection methods for two point boundary value problems, SIAM J. Rev. 22 (1980) 156–171. [12] A. Spence, B. Werner, Nonsimple turning points and cusps, SIAM J. Numer. Anal. 2 (1982) 413–427. [13] I. Stakgold, Branching solutions of nonlinear equations, SIAM Rev. 13 (1971) 283–332. [14] H. Weber, Numerical treatment of bifurcation problems for ordinary differential equations, Numer. Math. 32 (1979) 17–29. [15] H. Weber, The numerical approximation of secondary bifurcation problems, in: Numerical Solution of Nonlinear Equations, Lecture Notes in Math. 878, Springer Verlag, Berlin, 1981, pp. 407–427.