Constructing Cubic Splines on the Sphere Except where reference is made to the work of others, the work described in this thesis is my own or was done in collaboration with my advisory committee. This thesis does not include proprietary or classi ed information. Mosavverul Hassan Certi cate of Approval: Narendra Kumar Govil Professor Department of Mathematics University of Montana Amnon J. Meir, Chair Professor Mathematics and Statistics Bertram Zinner Associate Professor Mathematics and Statistics George T. Flowers Acting Dean Graduate School Constructing Cubic Splines on the Sphere Mosavverul Hassan A Thesis Submitted to the Graduate Faculty of Auburn University in Partial Ful llment of the Requirements for the Degree of Master of Science Auburn, Alabama August 10, 2009 Constructing Cubic Splines on the Sphere Mosavverul Hassan Permission is granted to Auburn University to make copies of this thesis at its discretion, upon the request of individuals or institutions and at their expense. The author reserves all publication rights. Signature of Author Date of Graduation iii Thesis Abstract Constructing Cubic Splines on the Sphere Mosavverul Hassan Master of Science, August 10, 2009 (M.Sc., I.I.T. Guwahati{India, 2006) (B.Sc., Ranchi University, 2002) 47 Typed Pages Directed by Amnon J. Meir A method to approximate functions de ned on a sphere using tensor product cu- bic B-splines is presented here. The method is based on decomposing the sphere into six identical patches obtained by radially projecting the six faces of a circumscribed cube onto the spherical surface. The theory of univariate splines has been general- ized in di erent forms to functions of several variables. Among these extensions the tensor product splines are the easiest to handle. Although the tensor product splines are restricted to rectangular domains rendering their applicability limited they are extremely e cient compared to other surface approximation techniques which are far more complicated and hence computationally less attractive. iv Acknowledgments I would like to acknowledge and thank my professor, Dr. A.J. Meir for his continuous help and support he provided me throughout the thesis work. His vast knowledge, experience and patience helped me explore and bring my work to its conclusion. His constant encouragement motivated me to enrich myself with the scienti c acumen necessary for the present work. I would also like to express my gratitude to my committee members Dr. Narendra Kumar Govil and Dr. Bertram Zinner for their advice and support. v Style manual or journal used Journal of Approximation Theory (together with the style known as \aums"). Bibliograpy follows van Leunen?s A Handbook for Scholars. Computer software used The document preparation package TEX (speci cally LATEX) together with the departmental style- le aums.sty. vi Table of Contents List of Figures viii 1 Introduction 1 1.1 Objective . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1.1 Spline approximation and its signi cance . . . . . . . . . . . . 2 1.2 Spline Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2.1 B-spline Representation . . . . . . . . . . . . . . . . . . . . . . 5 1.2.2 Tensor Product Splines . . . . . . . . . . . . . . . . . . . . . . 9 1.2.3 Error Estimates . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2 Radial Projection 16 2.1 Characteristics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 3 Analysis 19 3.1 Univariate Cubic Spline Interpolation . . . . . . . . . . . . . . . . . . 19 3.1.1 Radial Projection: The One Dimensional Case . . . . . . . . . 21 3.1.2 Periodic Splines on a Square . . . . . . . . . . . . . . . . . . . 25 3.2 B-spline representation on a Square . . . . . . . . . . . . . . . . . . . . 27 3.3 Tensor Product Splines . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 4 Conclusion 35 Bibliography 36 Appendices 37 A Notations 38 A.0.1 One dimensional case . . . . . . . . . . . . . . . . . . . . . . . 38 A.0.2 Bivariate case . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 vii List of Figures 3.1 Approximation of the function f( ) = sin ;N = 28 . . . . . . . . . . . 24 3.2 Approximation of the function f( ) = sin ;N = 60 . . . . . . . . . . . 24 3.3 Approximation of the function f( ) = sin cos ;N = 28 . . . . . . . . 24 3.4 Approximation of the function f( ) = sin cos ;N = 60 . . . . . . . . 24 3.5 Approximation of the function f( ) = sin3 ;h = 2:5 10 1 . . . . . . 26 3.6 Approximation of the function f( ) = sin3 ;h = 6:25 10 2 . . . . . 26 3.7 Approximation of a function f =2C1[a;b];h = 1 . . . . . . . . . . . . . 26 3.8 Approximation of a function f =2C1[a;b];h = 1:5625 10 2 . . . . . . 26 3.9 Approximation of a function f( ) = sin 4 ;h = 6:25 10 2 . . . . . . 29 3.10 Approximation of a function f( ) = sin 4 ;h = 1:56 10 2 . . . . . . 29 3.11 Function f(x;y) = x6y6 . . . . . . . . . . . . . . . . . . . . . . . . . . 32 3.12 Approximation of a function f(x;y) = x6y6;h = 1 . . . . . . . . . . . 32 3.13 Approximation of a function f(x;y) = x6y6;h = 1:25 10 1 . . . . . . 33 3.14 Mesh on the sphere Sr . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 3.15 Mesh on the cube Bd . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 viii Chapter 1 Introduction 1.1 Objective Cubic splines on a spherical domain may be used for the approximation of func- tions de ned on such domains which serve as a tool for modeling and analyzing of physical processes. However evaluation of these functions explicitly is sometimes dif- cult due to the limitations of the underlying processes. In such cases it becomes pertinent to form an approximation of the de ned function. We also run into cases where the evaluation of the exact function is computationally expensive forcing us to use approximation techniques. The problem of functional approximation may be broadly classi ed into two categories. The rst involving problems where the exact function is unknown and approximation technique is based on function value at cer- tain set of discrete points. The second class of problems is related to physical process modeling. These usually involve operator equations. Our aim is to look into a suitable set of approximations A and develop means to select an appropriate approximation. Since we want the function to be approximated by some member of the approxima- tion set A we need to device a method to select this member. Usually this is done by choosing an approximation member such that the error is within a certain factor of the least error that can be achieved. Functions arising from physical processes are generally smooth implying an obvious need for the approximants to be su ciently smooth. Functional approximations on the sphere are of research interest since many geophysical applications including oceanography, climate modeling and modeling of earth?s gravitational potential involve large amount of data on the surface of the 1 earth (basic model is a sphere) or on the satellite orbit (approximately a spherical manifold). 1.1.1 Spline approximation and its signi cance It appears that an obvious choice for function approximation is the polynomial pm2Pm because of its relative smoothness and easy manipulation on a digital com- puter. However it turns out that interpolating polynomials do not always converge to the function being interpolated [3]. The following theorem justi es this. Theorem 1.1. Let [a;b] be xed and suppose that for each k > 1;tk1;tk2;:::;tkk is a collection of points in [a;b]: Then there exists a function f2C[a;b] such that jj(f Lmf)jj1 !1 as m!1 (1.1) where Lmf is the unique polynomial of order k interpolating f at tk1;tk2;:::;tkk: This leads to the approximation using smooth piecewise polynomials i.e splines. 1.2 Spline Theory We address here the one dimensional case of approximating a given function using univariate splines. Mathematically we may represent it as f(xi) = sm(xi); where xi denote the nodal points. That is to say that the constructed spline agrees with the function values at a certain set of points termed here as the nodal points. We would want this spline to possess a certain degree of smoothness. We then represent the spline on any partition say as a linear combination of the basis elements of the linear space to which it belongs. These basis elements are generally called the B-splines. 2 De nition 1.1 (Piecewise Polynomials [3]). Let a = x0 < x1 < x2 < < xn 1 < xn = b; and write = fxign0: The set partitions the interval [a;b] into n subintervals, Ii = [xi;xi+1); i = 0;1;:::;n 2; and In 1 = [xn 1;xn]: Given a positive integer m, let PPm( ) = ff : there exists polynomials p0;p1;:::;pn 1 in Pm with f(x) = pi(x) for x2Ii;i = 0;1;:::;n 1g; where Pm =fp(x) : p(x) = mX i=1 cixi 1; c1;c2;:::;cm; x realg: We call PPm( ) the space of piecewise polynomials of order m with knots x1;x2;:::; xn 1: Switching from the approximation of a given function by a polynomial to approx- imation using piecewise polynomials provides us with a degree of exibility. However piecewise polynomials are not necessarily smooth. To maintain exibility and at the same time allow a certain degree of global smoothness we now de ne a class of functions known as polynomial splines. De nition 1.2. Let be a partition of the interval [a;b] as in De nition 1.1, and let m be a positive integer. Let Sm( ) = PPm( )\Cm 2[a;b]; where PPm( ) is the space of piecewise polynomials de ned in (1.1). We call Sm( ) the space of polynomial splines of order m (degree m 1) with respect to . Polynomial splines spaces are nite dimensional linear spaces. The dimension of this space of splines is dim(Sm( )) = n+m 1: The above de nition clearly implies that any polynomial on of degree 6m 1 is a spline function of degree m 1 on 3 . In general a spline of degree m 1 is represented by di erent polynomials in each interval Ii; i = 0;1;:::;n 1: This may give rise to dicontinuities in its (m 1) th derivatives at the internal nodes x1;x2;:::;xn 1. The nodes for which this actually happens are called active nodes. Let sm2Sm( ) be a spline of degree m 1 de ned on the partition . Let us denote the restriction of this spline function as smj[xi;xi+1] where smj[xi;xi+1] = m 1X j=0 sji(x xi)j; if x2[xi;xi+1] so we have mn coe cients to determine. Again we have the continuity conditions at the internal nodes. Each internal node has m 1 continuity conditions which amounts to (n 1)(m 1) conditions. We therefore have mn (n 1)(m 1) = m + n 1 coe cients to determine. Since we are talking about an interpolatory spline we have smj[xi;xi+1](xi) = f(xi) fi for i = 0;:::;n where the n+1 function values are known. We now have (m+n 1) (n+1) = m 2 coe cients still unaccounted for. To be more precise we still need m 2 conditions to determine the spline completely. This leads to imposing further constraints i.e conditions for periodicity or conditions for the spline to be natural. Mathematically they are represented as 1. Periodic splines, if slm(a) = slm(b); l = 0;:::;m 2: (1.2) 2. Natural splines, if for m = 2p 1; with l> 2 sp+jm (a) = sp+jm (b) = 0; j = 0;1:::;p 2 (1.3) 4 We will be dealing with cubic periodic splines throughout our work and unless otherwise mentioned m = 4: 1.2.1 B-spline Representation Before we de ne the B-Spline representation for a spline sm2Sm( ) and what the B-splines themselves are we de ne the concept of divided di erence since B-splines can be de ned in terms of the divided di erence. De nition 1.3 (Divided Di erence [2]). The n-th divided di erence of a function f at the points x0;x1;:::;xn (which are assumed to be distinct) is the leading coe cient (i.e the coe cient of xn) of the unique polynomial pn+1(x) of degree n which satis es pn+1(xi) = f(xi); i = 0;1;:::;n: Mathematically the n-th divided di erence is denoted as f[x0;x1;:::;xn] = nX i=0 f(xi) !0n+1(xi); (1.4) where !n+1(x) = nY i=0 (x xi): We now de ne the B-splines in terms of divided di erence De nition 1.4 (Normalized B-splines). The normalized B-splines of degree m 1 relative to the distinct nodes xi;xi+1;:::;xi+m is de ned as Bi;m(x) = (xi+m xi)g[xi;:::;xi+m]: (1.5) where g(t) = (t x)m 1+ = 8> >< >>: (t x)m 1 if x6t; 0 otherwise. 5 To nd an explicit expression for the normalized B-splines we state here the the Uniqueness Theorem for an interpolating polynomial which forms a basis for for this explicit expression. Theorem 1.2. Given n + 1 distinct points x0;x1;:::;xn and n + 1 corresponding values f(x0);f(x1);:::;f(xn) there exists a unique polynomial pn+12Pn+1 such that pn+1(xi) = f(xi) for i = 0;1;:::;n: The uniqueness of the interpolating polynomial provides for the comparison be- tween the Lagrange?s form of the interpolating polynomial and the Newton?s divided di erence formula. This comparison along with the notion that the divided di erence is the coe cient of xn in the interpolating polynomial yields the explicit represen- tation for the n-th divided di erence as de ned in equation (1.4). Using equation (1.5) in the expression for normalized B-splines we arrive at the following explicit representation for B-splines Bi;m(x) = (xi+m xi) mX j=0 (xj+i x)m 1+ mY l=0l6=j (xj+i xl+i) : (1.6) We note here that the m-th order normalized B-spline have active nodes xi;xi+1;:::;xm and vanish outside the interval [xi;xi+m]: The term normalized for B-splines has been introduced since B-splines can have varied sizes depending on the location of the nodal points, for example Qi;1(x) = 8> >< >>: 1 xi+1 xi; xi 6x >< >>: 1 if xi 6x 0 for x2(xi;xi+m): In view of Theorem 1.4 we can now de ne any spline sm(x) 2Sm( ) uniquely as a linear combination of these basis elements i.e sm(x) = n 1X m+1 ciBi;m(x): (1.13) The real numbers ci are called the B-spline coe cients of sm: The nodes in equation (1.11) are generally chosen as periodic or coincident. For periodicity we must have x i = xn i b+a; xn+i = xi +b a: (1.14) Using equations (1.14) and (1.6) we have the following condition for periodicity B i;m(x) = Bn i;m(x+b a); i = 1;:::;m 1 (1.15) Since the B-splines vanish outside their local support i.e [xi;xi+m] the condition de- ned in equation (1.15) will be satis ed if and only if c i = cn i; i = 1;:::;m 1: (1.16) 8 1.2.2 Tensor Product Splines In this section we introduce the Tensor Product Splines as an extension of the univariate B-spline representation. We require here only a few concepts concern- ing tensor products of vector spaces to de ne the tensor product polynomial spline. Mathematically we denote the tensor product of two vector spaces U and V as U V: Let us x a eld say C and let U and V be vector spaces de ned over this eld. Then we can have De nition 1.6 (Bilinear Mapping). A function f from U V to the vector space P is said to be bilinear if it is linear in each of the two variables when the other is kept xed. The space of all such bilinear functions forms a vector space over the eld C under addition and scalar multiplication. If U is a linear space of functions de ned on some set X and V; a linear space of functions de ned on some set Y into R then for each u2U and v2V the rule w(x;y) = u(x)v(y); 8(x;y)2X Y (1.17) de nes a function on X Y called the tensor product of u with v [1] and is denoted as u v: Further the set of all nite linear combinations of the functions on X Y of the form u v is called the tensor product of U with V. Hence U V = ( nX i=1 i(ui vi) : i2R;ui2U;vi2V; (1.18) i = 1;:::;n; for some n ) (1.19) It can be veri ed that the tensor product de ned above is bilinear, i.e., the map U V !U V : (u;v)7 !u v (1.20) 9 is linear in each argument. ( 1u1 + 2u2) v = 1(u1 v) + 2(u2 v) u ( 1v1 + 2v2) = 1(u v1) + 2(u v2): (1.21) Remark 1.1. The tensor product discussed above forms a linear space of functions de ned on X Y and its dimension is given by the following proposition. Proposition 1.1 (Tensor Product Splines). If U and V are some vector spaces de ned over a eld C, then dim(U V) = (dimU)(dimV) As mentioned earlier the tensor product we use here is simply an extension from the univariate case to the bivariate case and hence we will be considering two sets of partitions one along the horizontal and one along the vertical. De nition 1.7 (Tensor product Splines). Consider the strictly increasing se- quences a = x0 0 in x and l 1 > 0 in y with respect to the partition along the horizontal and the vertical as de ned in equations (1.22)-(1.23) if the following conditions are satis ed 10 1. On each subrectangle Ri;j = [xi;xi+1] [yi;yi+1];s(x;y) is given by a polynomial of degree m 1 in x and l 1 in y: sjRi;j 2Pm Pl; i = 0;1;:::;n 1; j = 0;1;:::;p 1: 2. The function s(x;y) and all its partial derivatives @i+js(x;y) @xi@yj 2C(R); i = 0;1;:::;m 2; j = 0;1;:::;l 2: For our particular casem = l = 4 and the dimension of the vector space S( 1; 2) of all functions satisfying the above two conditions is dim(S( 1; 2)) = (m + n 1) (l + p 1) where i; i = 1;2 are the partitions along the horizontal and the vertical and are de ned as 1 = fxign0 and 2 = fyjgp0 respectively. We can de ne an extended partition as in De nition 1.5 and because of the tensor product nature of the vector space it is possible to work with the one dimensional B-spline basis. This gives us the unique tensor product representation of a spline s(x;y) in terms of its basis functions s(x;y) = n 1X i= m+1 p 1X j= l+1 ci;jBi;m(x)Bj;l(y) (1.24) where Bi;m(x) and Bj;l(y) are the normalized B-splines de ned on the partitions 1 and 2 respectively. For the approximating function to be periodic in both x(horizontal) and y(vertical) direction we have @is(a;y) @xi = @is(b;y) @xi ; i = 0;:::;m 2; c6y 6d @js(x;c) @yj = @is(x;d) @yj ; j = 0;:::;l 2; a6x6b: (1.25) 11 1.2.3 Error Estimates We present here a priori error bounds for the interpolation procedures introduced in Sections 1.2.1 and 1.2.2. The derivation of the interpolation error, f Sf and its derivatives in L2 norm and the L1 norm are given. Here Sf denotes the approximation to the function f. It is observed that if the function f is su ciently smooth, Sf is a fourth order approximation to f in both L2 and L1. Theorem 1.5 (Variational Problem). Let and f ff0;f1;:::;fn;f0a;f0bg be given and V fw 2 PC22(I) j w(xi) = f(i); 0 6 i 6 n and Dw(xi) = f0i; i = 0 and ng. The variational problem of nding the functions p2V which minimize jjD2wjj22 over all w2V has the unique solution Sf: The function p2V is a solution of the variational problem if and only if D2p;D2 2 = 0 for all 2V0 fw2PC22(I) jw(xi) = 0;0 6 i 6 n;and Dw(xi) = 0; i = 0 and ng: By de nition we have jjD2p+ D2 jj22 = (D2p;D2p)2 + 2(D2p;D2 )2 + (D2 ;D2 )2: (1.26) By the orthogonality condition we have jjD2p+ D2 jj22 = (D2p;D2p)2 + (D2 ;D2 )2; (1.27) and this would gives us the following corollary Corollary 1.1 (First integral relation). If f2PC22(I); then jjD2fjj22 =jjD2Sfjj22 +jjD2Sf D2fjj22: (1.28) 12 Theorem 1.6 (Preliminary Result). If f2PC22(I); then jjD2(f Sf )jj2 6 jjD2fjj2 (1.29) jjD(f Sf )jj2 6 2 1hjjD2fjj2 (1.30) jj(f Sf )jj2 6 2 2h2jjD2fjj2 (1.31) Inequality (1.29) follows directly from Corollary 1.1. We can see from the above theorem that Sf is a second-order approximation to f. Intuitively we may assume that a smoother function will result in a higher order of convergence and hence we have the following theorem. Theorem 1.7. If f2PC14 ; then jjf Sfjj16 5584h4jjD4fjj1: (1.32) Moreover if f2C5(I) and is a uniform partition, then jjf Sfjj16h4( 1384jjD4fjj1 + 1240hjjD5fjj1) (1.33) Proof. Since f Sf = f Sh + Sh Sf where Sh is the cubic Hermite interpolate of f, we have jjf Sfjj16jjf Shjj1 +jjSh Sfjj1: (1.34) Now jjf Shjj16 1384h4jjD4fjj1 (1.35) 13 and Sh Sf = nX i=1 e0ih0i(x); hence jjSh Sfjj1 = jj nX i=1 e0ih0i(x)jj1 6 jje0jj1 n nX i=1 h0i(x) o ; (1.36) where jje0jj1 max16i6nje0ij6 124h3jjD4fjj1: (1.37) Now, since jh0ij+jh0i+1j6 h4 for all x2[xi;xi+1] and 0 6i6n; (1.38) we have from equation (1.37) jj nX i=1 e0ih0i(x)jj16 196h4jjD4fjj1; (1.39) which gives us the required result jjf Sfjj16 1384h4jjD4fjj1 + 196h4jjD4fjj1 = 5584h4jjD4fjj1: (1.40) We proceed now to the error bounds of the bivariate interpolation procedure and nd that as in case of the the univariate splines the approximation Sf for a su ciently smooth function is fourth order accurate in both the L2 and the L1-norm. Theorem 1.8. If f2PC24(U); then jjf Sfjj2 6 4 4 h4jjD4xfjj2 +h2k2jjD2xD2yfjj2 +k4jjD4yfjj2 : (1.41) 14 The above theorem gives us the error bound in L2-norm. For the error bound in the L1-norm we have the following theorem Theorem 1.9. If f2PC14 (U); then jjf Sfjj16 5384h4jjD4xfjj1 + 49h2k2jjD2xD2yfjj1 + 5384k4jjD4yfjj1; (1.42) which clearly shows that the approximating spline is fourth order accurate. 15 Chapter 2 Radial Projection The radial projection method which we describe here (see [4], p. 24) is a method to radially project the points on the surface of the cube onto the sphere. Since the tensor product B-splines are restricted to rectangular domains all calculations will essentially be done on the cube. The terminology used here will run as follows: The surface of the cube will be termed as the box Bd centered at the origin and of side length 2d: We will denote the sphere centered at the origin with radius r as Sr: Mathematically we may represent the box and the sphere as follows Bd = fxjx(x1;x2;:::;xn)2Rn;jjxjj1 = dg Sr = faja(a1;a2;:::;an)2Rn;jjajj2 = rg The radial projection from the box to the sphere is de ned as a mapping P : Bd !Sr given by P(x) = r xjjxjj = a; where as the inverse mapping from the sphere to the box is given by P 1(a) = d ajjajj 1 = x: 16 2.1 Characteristics The radial projection P is a one-one mapping from the box Bd to the sphere Sr: We mention below some related properties. The following lemma shows that the mapping P and its inverse P 1 are both Lipschitz continuous. Lemma 2.1. The radial projection P and its inverse P 1 satisfy the inequalities jjP(x) P(y)jj 6 2rjjxjj jjx yjj ; (2.1) jjP 1(a) P 1(b)jj1 6 2djjajj 1 jja bjj1 : (2.2) Proof. jjP(x) P(y)jj = jj rjjxjjx rjjyjjyjj = rjjxjjjjyjj n jj jjyjjx jjxjjy jj o = rjjxjjjjyjj n jj jjyjj(x y) y(jjyjj jjxjj) jj o 6 rjjxjjjjyjj n jjyjjjjx yjj+jjyjj jjyjj jjxjj o 6 2rjjxjj n jjx yjj o In a similiar way we can prove inequality (2.2) Corollary 2.1. The radial projection P and its inverse P 1; are globally Lipscitz continuous mappings, that is, jjP(x) P(y)jj6 2rdjjx yjj; (2.3) and jjP 1(a) P 1(b)jj16 2dnr jja bjj1: (2.4) 17 This follows directly from the lemma 2.1 by observing that for anyx2Rn;jjxjj16 x 6pnjjxjj1 and that for any a2 Sr; rpn 6jjajj1 6 r: For any x2 Bd; we have x2Bd; d6x6dpn: 18 Chapter 3 Analysis In this chapter we discuss the construction of Tensor Product Splines as a natural extension of the B-spline representation of a spline. We also estimate the di erence (f sm(x)) where f is the function de ned at the nodal points and sm(x) is the approximating spline. Univariate spline representation is analyzed in Section 3.1. In Section 3.2 we analyze the B-spline representation of a spline and in Section 3.3 extend the construction to Bivariate splines and analyze it. 3.1 Univariate Cubic Spline Interpolation Let us consider a partition of an interval [a;b] as de ned in De nition 1.1 with a x0 = xn b and the corresponding function evaluations at the nodal points fi; i = 0;1;:::;n 1: Our aim here will be to develop an e cient method to construct a periodic cubic spline interpolating the function values at the distinct nodal points [5]. Since the degree of the spline is m 1 = 3 the spline must be twice continuously di erentiable i.e the second order derivative must be continuous. We introduce here the following notations fi = s4(xi); m = s04(xi); and M = s004(xi); i = 0;1;:::;n 1: Due to the periodic consideration we have fn+j = fj and Mn+j = Mj for j = 0;1: Since s4;i 12P4; s004;i 1 is linear and s004;i 1(x) = Mi 1xi xh i +Mix xi 1h i for x2[xi 1;xi] (3.1) 19 where hi = xi xi 1; i = 1;:::;n: Integrating (3.1) twice we get s4;i 1(x) = Mi 1(xi x) 3 6hi +Mi (x xi 1)3 6hi +Ci 1(x xi 1) + ~Ci 1 (3.2) and the constants Ci 1 and ~Ci 1 are determined by imposing the end point values s4(xi 1) = fi 1 and s4(xi) = fi:This gives us, for i = 1;:::;n ~Ci 1 = fi 1 Mi 1h2i 6 ; Ci 1 = fi fi 1 hi hi 6 (Mi Mi 1): (3.3) Imposing the continuity of the rst derivatives at xi; we get s04(x i ) = hi6 Mi 1 + hi3 Mi + fi fi 1h i = hi+13 Mi hi+16 Mi+1 + fi+1 fih i+1 = s04(x+i ); where s04(x i ) = limt!0s04(x t): This gives us the following linear system also know as the M-continuity system iMi 1 + 2Mi + iMi+1 = di; i = 1;:::;n (3.4) where i = hih i +hi+1 ; i = hi+1h i +hi+1 di = 6h i +hi+1 (fi+1 fih i+1 fi fi 1h i ) i = 1;:::;n: It is clear from (3.2) that the only unknows are M0;M1;:::;Mn 1: So our task of nding a periodic cubic spline representation interpolating the given function values 20 now reduces to solving the linear system (3.4) of n equations and n unknowns. This construction of the spline produces a system tridiagonal in nature. In matrix notation it is represented as 0 BB BB BB BB BB @ 2 n 1 0 ::: n 1 n 2 2 n 2 ... 0 ... ... ... 0 ... 0 1 2 1 n 0 0 n 2 1 CC CC CC CC CC A 0 BB BB BB BB BB @ Mn 1 Mn 2 ... M1 M0 1 CC CC CC CC CC A = 0 BB BB BB BB BB @ dn 1 dn 2 ... d1 d0 1 CC CC CC CC CC A (3.5) and can be easily solved on a computer using existing techniques. 3.1.1 Radial Projection: The One Dimensional Case In the above discussion we have considered the spline to be periodic in nature as we want to apply the construction on a circle C with radius r = 1 such that the rst and last nodal point coincide i.e x0 = xn. The radial projection of the nodal points on the square onto a unit circle involves basic geometry. Let the four corners, of the square under consideration be P1(1;1);P2( 1;1);P3( 1; 1);P4(1; 1) in this order. Let P(x;y) be any point on the side of the square, say the side joining P1 and P2: Let P0(x0;y0) be the radial projection of the point P(x;y) on the circle. Then we have the relation tan = y 0 x0 = y x (3.6) where is the angle between the X-axis and the vector joining the point P to the origin. Since the point P0 lies on the circle we have y02 +x02 = 1: (3.7) 21 Solving for x0 and y0 we have x0 = xpx2 +y2 y0 = ypx2 +y2 Taking points systematically anticlockwise on the square starting at P1 and then projecting them onto the circle we have points on the circle which we now label as xi; i = 0;1;:::;n 1: This provides us with the setting required to apply the univariate spline constructed in Section 3.1. However we are still required to nd the interval lengths hi; i = 0;1;:::;n 1 between the nodal points. To do this we make use of the inner product of two vectors. Let ~xi and ~xi 1 be two vectors then cos i = <~xi;~xi 1 >jj~x ijjjj~xi 1jj ; (3.8) where i as the angle between two vectors labelled here as ~xi and ~xi 1: Hence the arc length or the interval length hi = r i: Table 3.1: The following table shows the error and the observed rate of convergence for various step sizes for the function f( ) = sin( ). Step Size Error Order of Convergence 2:243994752564133e 01 1:521985704066554e 04 1:047197551196593e 01 7:326907390286181e 06 3:980413480161499e+ 00 5:067084925144792e 02 4:029381312631907e 07 3:995561255766796e+ 00 2:493327502848618e 02 2:703779530789965e 08 3:809570269779505e+ 00 1:236847501412775e 02 1:903539716356657e 09 3:785053304137925e+ 00 The experimentally observed order of convergence is de ned as p = ln (Erj+1Erj ) ln (hj+1hj ) ; 22 where Erj denotes the error at the jth re nement and is de ned as Erj =jjhx(f s4(x))jj2 ; and hx represents the step size for the intermediate points taken to test the spline function developed. 23 0 1 2 3 4 5 6 7 !1 0 1 Plot of the function 0 1 2 3 4 5 6 7 !1 0 1 Plot of the Spline Function 0 1 2 3 4 5 6 7 !1 0 1 x 10 !3 Plot of the difference Figure 3.1: Approximation of the function f( ) = sin ;N = 28 0 1 2 3 4 5 6 7 !1 0 1 Plot of the function 0 1 2 3 4 5 6 7 !1 0 1 Plot of the Spline Function 0 1 2 3 4 5 6 7 !1 0 1 2 x 10 !4 Plot of the difference Figure 3.2: Approximation of the function f( ) = sin ;N = 60 0 1 2 3 4 5 6 7 !0.5 0 0.5 Plot of the function 0 1 2 3 4 5 6 7 !0.5 0 0.5 1 Plot of the Spline Function 0 1 2 3 4 5 6 7 !1 0 1 x 10 !3 Plot of the difference Figure 3.3: Approximation of the function f( ) = sin cos ;N = 28 0 1 2 3 4 5 6 7 !0.5 0 0.5 Plot of the function 0 1 2 3 4 5 6 7 !0.5 0 0.5 1 Plot of the Spline Function 0 1 2 3 4 5 6 7 !1 0 1 2 x 10 !4 Plot of the difference Figure 3.4: Approximation of the function f( ) = sin cos ;N = 60 24 3.1.2 Periodic Splines on a Square De ning periodic splines on a square is based again on the construction given in Section 3.1. However in this case the nodal points on the square are not projected on a circle. Hence with equally spaced nodes on the edges of the square i.e hi = hi+1 we have from equation (3.4) i = 12; i = 12 di = 3h2 i (fi+1 2fi +fi 1); i = 1;:::;n: Table 3.2: The following table shows the error and the observed rate of convergence for various step sizes for the function f( ) = sin3( ). Step Size Error Order of Convergence 1:000000000000000e+ 00 2:993908237928666e 02 5:000000000000000e 01 3:738564612692465e 03 3:001473631942863e+ 00 2:500000000000000e 01 1:044464521796866e 04 5:161649073884993e+ 00 1:250000000000000e 01 4:952085825195140e 06 4:398583359488427e+ 00 6:250000000000000e 02 2:868804616210702e 07 4:109514698407482e+ 00 3:125000000000000e 02 1:758372110067624e 08 4:028137400984083e+ 00 1:562500000000000e 02 1:093598905987019e 09 4:007084798744635e+ 00 Table 3.3: The following table shows the error and the observed rate of convergence for various step sizes for a function f =2C1[a;b] Step Size Error Order of Convergence 1:000000000000000e+ 00 8:728715607973290e 03 5:000000000000000e 01 8:728715607973281e 03 1:441541926716714e 15 2:500000000000000e 01 2:822924332662890e 03 1:628578924802774e+ 00 1:250000000000000e 01 9:724452119338235e 04 1:537501583136462e+ 00 6:250000000000000e 02 3:437050954186293e 04 1:500445730420506e+ 00 3:125000000000000e 02 1:215162841762690e 04 1:500021580057984e+ 00 1:562500000000000e 02 4:295235806309275e 05 1:500340417834547e+ 00 25 0 1 2 3 4 5 6 7 8 !1 0 1 Number of Points Function f 0 1 2 3 4 5 6 7 8 !1 0 1 Number of Points Spline 0 1 2 3 4 5 6 7 8 !2 0 2 x 10 !3 Number of Points Function f ! spline Figure 3.5: Approximation of the function f( ) = sin3 ;h = 2:5 10 1 0 1 2 3 4 5 6 7 8 !1 0 1 Number of Points Function f 0 1 2 3 4 5 6 7 8 !1 0 1 Number of Points Spline 0 1 2 3 4 5 6 7 8 !1 0 1 x 10 !5 Number of Points Function f ! spline Figure 3.6: Approximation of the function f( ) = sin3 ;h = 6:25 10 2 It is evident from the Figure 3.8 below that the function f =2C1[a;b] and hence the observed order of convergence. In view of the approximation power of splines [3] we may expect that the order of approximation attainable will increase with the smoothness of the class of functions F being approximated. However this is true only up to a limit. In fact if F\Pm =;; then the maximal order of convergence possible for the class F is m; no matter how smooth F is assumed. In our case m = 4: This is known as the saturation result. 0 1 2 3 4 5 6 7 8 0 0.5 1 Number of Points Function f 0 1 2 3 4 5 6 7 8 0 0.5 1 Number of Points Spline 0 1 2 3 4 5 6 7 8 !0.1 0 0.1 Number of Points Function f ! spline Figure 3.7: Approximation of a function f =2C1[a;b];h = 1 0 1 2 3 4 5 6 7 8 0 0.5 1 Number of Points Function f 0 1 2 3 4 5 6 7 8 0 0.5 1 Number of Points Spline 0 1 2 3 4 5 6 7 8 !4 !2 0 2 4 x 10 !3 Number of Points Function f ! spline Figure 3.8: Approximation of a function f =2C1[a;b];h = 1:5625 10 2 26 3.2 B-spline representation on a Square Our aim here is to construct a cubic spline which is represented as a linear combination of the B splines Bi;m(x) as in Section 1.4. We consider the extended partition =fxign+m 1 m+1 of the interval [a;b] as de ned in the De nition 1.5. Using the periodicity condition available to us through the equations (1.14),(1.15) and (1.16) we have a system of linear equations available to us which we can represent in a matrix form. We denote the matrix of all basis functions as B: Using the properties of the B-spline and noting that the spline function must agree with the function values at the nodal points xi; i = 0;1;:::;n 1 we arrive at the following matrix representation of the linear system, BC = F; which implies C = B 1F: B = 0 BB BB BB BB BB BB BB BB BB BB BB BB BB B@ 0 ::: ::: ::: ::: ::: 0 B0 3 B0 2 B0 1 B10 0 ... ... ... ... ... 0 B1 2 B1 1 B20 B21 0 ... ... ... ... ... 0 B2 1 B30 B31 B32 0 ... ... ... ... ... 0 0 B41 B42 B43 0 ... ... ... ... 0 0 0 ... ... ... ... ... ... ... 0 ... 0 0 ... ... ... ... ... ... ... 0 ::: ::: ::: Bn 3n 6 Bn 3n 5 Bn 3n 4 0 0 0 0 ::: ::: ::: ::: Bn 2n 4 Bn 2n 3 0 0 0 ::: ::: ::: ::: ::: Bn 1n 3 Bn 1n 2 0 1 CC CC CC CC CC CC CC CC CC CC CC CC CC CA (3.9) The matrix B has as its entries B-splines evaluated at the nodal points where Bji represents Bi;m(xj); i = m+ 1; m+ 2;:::;n 1 and j = 0;1;:::;n 1: We have 27 the zero entries in the above matrix because the B-spline Bi;m(x) vanishes outside its local support [xi;xi+m]: The vector of all coe cients is represented as C and the vector of all function values at the nodal points is represented as F, then C = 0 BB BB BB BB BB BB BB BB BB BB BB BB BB B@ C0 C1 C2 C3 C4 C5 ... Cn 3 Cn 2 Cn 1 1 CC CC CC CC CC CC CC CC CC CC CC CC CC CA and F = 0 BB BB BB BB BB BB BB BB BB BB BB BB BB B@ f0 f1 f2 f3 f4 f5 ... fn 3 fn 2 fn 1 1 CC CC CC CC CC CC CC CC CC CC CC CC CC CA (3.10) The problem of representing the cubic spline as a linear combination of the basis functions is numerically equivalent to the problem of nding the B-spline coe cients, i.e., the vector C: Since the matrix B and the vector F is completely known we can easily solve for C on a digital computer. Table 3.4: The following table shows the error and the observed rate of convergence for various step sizes for the function f( ) = sin( 4 ). Step Size Error Order of Convergence 1:000000000000000e+ 00 3:412826951770349e 03 5:000000000000000e 01 9:534612985020590e 05 5:161649073886882e+ 00 2:500000000000000e 01 4:520615188600750e 06 4:398583359486799e+ 00 1:250000000000000e 01 2:618848335825963e 07 4:109514698295172e+ 00 6:250000000000000e 02 1:605166782410018e 08 4:028137400685653e+ 00 3:125000000000000e 02 9:983146507216004e 10 4:007084797322440e+ 00 28 As we are approximating a function de ned on a square of edge length 2, the approximating spline will have a periodicity of b a which in our case is 8: Hence the factor of 4 in the function f( ) = sin 4 : 0 1 2 3 4 5 6 7 8 !1 0 1 Number of Points Function f 0 1 2 3 4 5 6 7 8 !1 0 1 Number of Points Spline 0 1 2 3 4 5 6 7 8 !5 0 5 x 10 !7 Number of Points Function f ! spline Figure 3.9: Approximation of a function f( ) = sin 4 ;h = 6:25 10 2 0 1 2 3 4 5 6 7 8 !1 0 1 Number of Points Function f 0 1 2 3 4 5 6 7 8 !1 0 1 Number of Points Spline 0 1 2 3 4 5 6 7 8 !2 0 2 x 10 !9 Number of Points Function f ! spline Figure 3.10: Approximation of a function f( ) = sin 4 ;h = 1:56 10 2 Table 3.5: The following table shows the error and the observed rate of convergence for various step sizes for the function f =2C1[a;b] Step Size Error Order of Convergence 1:000000000000000e+ 00 7:968190728250351e 03 5:000000000000000e 01 7:968190728250389e 03 7:047538308392802e 15 2:500000000000000e 01 2:576965563447752e 03 1:628578922185994e+ 00 1:250000000000000e 01 8:877169929217066e 04 1:537501539184818e+ 00 6:250000000000000e 02 3:137585531721901e 04 1:500445025990139e+ 00 3:125000000000000e 02 1:109295998711029e 04 1:500010409963090e+ 00 1:562500000000000e 02 3:921497851505102e 05 1:500167662913683e+ 00 3.3 Tensor Product Splines In this section we will construct the tensor product spline rst on a square patch and then extend this notion to the cube treating each face of the cube as a square 29 patch. This will give us a method to approximate any function posed on the sphere with the help of the mapping P and P 1 from the surface of the cube to the surface of the sphere and back respectively. From equation (1.24) we have the unique tensor product representation of a spline de ned over a rectangle [a;b] [c;d] with respect to the extended partition 1 and 2 given by De nition 1.5 in the horizontal and vertical directions respectively. Keeping in mind that we do not want periodicity on a single patch but rather across the four faces of the cube we do not implement the periodic conditions for the single patch which we create. The cube in question is the box Bd as de ned in Chapter 2. It is centered at the origin with d = 1: It is clear from the above setup that that the total number of coe cients that need to be determined are (n+3) (p+3) since m = l = 4: However we only have the function values at the nodal points (xi;yj); i = 0;1;:::;n; j = 0;1;:::;p: Hence we still need 2n+ 2p+ 8 conditions to have a unique spline representation on the patch. These extra conditions are given as restrictions on the derivatives of the spline at the boundary and the corner points of the grid formed by the partitions 1 and 2: We give here a brief list of the boundary conditions generally associated with the tensor product cubic splines. Boundary conditions of the rst type @s @x(xi;yj) = f x ij; i = 0;n; j = 0;1;:::;p @s @y(xi;yj) = f y ij; i = 0;1;:::;n; j = 0;p @2s @x@y(xi;yj) = f xy ij ; i = 0;n; j = 0;p (3.11) The total number of rst boundary conditions here is 2n+ 2p+ 8: Boundary conditions of the second type 30 @2s @x2(xi;yj) = f x ij; i = 0;n; j = 0;1;:::;p @2s @y2(xi;yj) = f y ij; i = 0;1;:::;n; j = 0;p @4s @x2@y2(xi;yj) = f xy ij ; i = 0;n; j = 0;p (3.12) The total number of second boundary conditions is again 2n+ 2p+ 8: Boundary conditions of the third type Boundary conditions of the third type are called periodic boundary conditions. Peri- odicity with respect to the horizontal variable must be Px = b a and with respect to the vertical variable must be Py = d c: For our particular case we must have Px = Py = 8 s(x0;yj) = s(xn;yj); j = 0;1;:::;p s(xi;y0) = s(xi;yp); i = 0;1;:::;n @ks @xk(x0;yj) = @ks @xk(xn;yj); j = 0;1;:::;p; k = 1;2 @ls @yl(xi;y0) = @ls @yl(xi;yp); i = 0;1;:::;n; l = 1;2 @2ks @xk@yk(x0;yj) = @ks @xk@yk(xn;yj); j = 0;1;:::;p; k = 1;2 @2ks @xk@yk(xi;y0) = @2ks @xk@yk(xi;yp); i = 0;1;:::;n; k = 1;2 (3.13) Applying the boundary conditions of the rst type in order to determine the coe cients we end up with a system of linear equations which in matrix form may be represented as MCN = F C = M 1FN 1 (3.14) 31 where M(p+3) (p+3) is the matrix of B-splines in the vertical direction, N(n+3) (n+3) is the matrix of B-splines in the horizontal direction, C(p+3) (n+3) is the coe cient matrix and F(p+3) (n+3) is the matrix of all function values at the nodal points in the grid along with its derivatives at the boundary and corner points. Table 3.6: The following table shows the error and the observed rate of convergence for various step sizes for the function f(x;y) = x6y6 Step Size Error Order of Convergence 1:000000000000000e+ 00 7:060256060377759e+ 00 5:000000000000000e 01 9:550664958381095e 01 2:886047419513821e+ 00 2:500000000000000e 01 5:683782321068801e 02 4:070677975228285e+ 00 1:250000000000000e 01 3:408501271897639e 03 4:059641876435491e+ 00 6:250000000000000e 02 2:095067852178469e 04 4:024068647573170e+ 00 !1 !0.5 0 0.5 1 !1 !0.5 0 0.5 1 0 0.2 0.4 0.6 0.8 1 X ! Axis Graph of the Function values at Intermediate Points Y !Axis Function Values 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Figure 3.11: Function f(x;y) = x6y6 !1 !0.5 0 0.5 1 !1 !0.5 0 0.5 1 !0.4 !0.2 0 0.2 0.4 0.6 0.8 1 X ! Axis Graph of the Spline Function values at intermediate Points Y !Axis Function Values !0.2 0 0.2 0.4 0.6 0.8 Figure 3.12: Approximation of a function f(x;y) = x6y6;h = 1 To obtain the the tensor product spline representation on the cube we note that spline function must be periodic with periodicity along x-direction Px, y-direction Py and z-direction Pz and Px = Py = Pz = 8. Even though the tensor products are de ned on individual patches, the spline function must be periodic across the four adjacent faces of the cube. Let us consider here the tensor product de ned on a 32 !1 !0.5 0 0.5 1 !1 !0.5 0 0.5 1 0 0.2 0.4 0.6 0.8 1 X ! Axis Graph of the Spline Function values at intermediate Points Y !Axis Function Values 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Figure 3.13: Approximation of a function f(x;y) = x6y6;h = 1:25 10 1 single patch. Then to obtain periodicity along the horizontal we need to wrap the ctitious nodes u i;i = m+1;:::; 1 and un+i;i = 1;:::;m 1 along the horizontal so that they coincide with the nodal points on the adjacent sides. For periodicity along the vertical we do the same, i.e., wrap the nodes v i;i = m + 1;:::; 1 and vp+i;i = 1;:::;m 1 along the vertical so that they coincide with nodal points on the adjacent sides. To do this we must have the following u i = ui hw; i = 1;2;3 un+i = un +hw; i = 1;2;3 (3.15) and v i = vi hw; i = 1;2;3 vp+i = vp +hw; i = 1;2;3 (3.16) where ui denotes the nodal points along the horizontal, vj denotes the nodal points along the vertical and hw is the step size of the adjacent side. 33 Table 3.7: The following table shows the error and the observed rate of convergence for various step sizes for the function f(x;y;z) = sin(xyz) Step Size Error Order of Convergence 1:000000000000000e+ 00 1:670920218139921e 01 5:000000000000000e 01 7:961222685614932e 03 4:391509023007923e+ 00 2:500000000000000e 01 4:648371361131741e 04 4:098192780859205e+ 00 1:250000000000000e 01 2:836252157206563e 05 4:034667624804363e+ 00 6:250000000000000e 02 1:759114473732576e 06 4:011064527357291e+ 00 3:125000000000000e 02 1:097374360545641e 07 4:002721689943840e+ 00 !1 !0.5 0 0.5 1 !0.5 0 0.5 !0.8 !0.6 !0.4 !0.2 0 0.2 0.4 0.6 0.8 x y z Figure 3.14: Mesh on the sphere Sr !1 !0.5 0 0.5 1 !1 !0.5 0 0.5 1 !1 !0.5 0 0.5 1 x y z Figure 3.15: Mesh on the cube Bd 34 Chapter 4 Conclusion We developed and analyzed a method to approximate functions posed on the sphere. The method described here is based on tensor product splines and because of the tensor nature of the resulting space many algebraic properties of univariate polynomial splines can easily be carried over. Although the tensor product splines have a restricted use they are computationally advantageous due to their easy imple- mentation on a digital computer. 35 Bibliography [1] Carl de Boor, \A Practical Guide to Splines," Applied Mathematical Sciences, Vol. 27. Springer, New York, 1978. [2] Paul Dierckx, \Curve and Surface Fitting with Splines," Oxford Science Publi- cations, 1993. [3] Larry L. Schumaker, \Spline Functions:Basic theory," Wiley, New York, 1981. [4] Necibe Tuncer, \A Novel Finite Element Discretization of Domains with Spheroidal Geometry," Ph.D Thesis, Auburn University, AL, May 2007. [5] Al o Quarteroni, Ricardo Sacco, Fausto Saleri, \Numerical Mathematics," Springer, 2000. [6] Martin H. Schultz, \Spline Analysis," Prentice-Hall, Inc., 1973. [7] M. G. Cox, \The Numerical evaluation of B-Splines," Journal of the Institute of Mathematical Applications, 10 (1972), 134-149. [8] A. Ralston, \A rst course in Numerical Analysis," Mc-Graw-Hill, New York, 1965. [9] C. Ronchi, R. Iaconco, P. S. Paolucci \The \cubed-sphere": A New Method for the Solution of Partial Di erential Equations in Spherical Geometry," J. Comput. Phy., 124 (1996), 93-114. [10] C. R. Trass, \Smooth Approximation of Data on the Sphere," Computing, 38 (1987), 177-184. [11] Amnon J. Meir, Necibe Tuncer, \Radially Projected Finite Elements," SIAM J. Sci. Comput., 31 (2007), 2368-2385. 36 Appendices 37 Appendix A Notations A.0.1 One dimensional case I [a;b] fxja6x6bg: For each nonnegative integer t and for each p; 1 6 p 61; PCpt(I) will denote the set of all real valued functions f(x) such that 1. f(x) is t 1 times continuously di erentiable, 2. there exists xi; 0 6i6n 1; with a = x0