Bivariate Lagrange interpolation at the node points of non-degenerate Lissajous curves
Wolfgang Erb, Christian Kaethner, Mandy Ahlborg, Thorsten M. Buzug
BBivariate Lagrange interpolation at the node pointsof non-degenerate Lissajous curves
Wolfgang Erb · Christian Kaethner · Mandy Ahlborg · Thorsten M. Buzug
September 26, 2014
Abstract
Motivated by an application in Magnetic Particle Imaging, westudy bivariate Lagrange interpolation at the node points of Lissajous curves.The resulting theory is a generalization of the polynomial interpolation the-ory developed for a node set known as Padua points. With appropriately de-fined polynomial spaces, we will show that the node points of non-degenerateLissajous curves allow unique interpolation and can be used for quadraturerules in the bivariate setting. An explicit formula for the Lagrange polynomi-als allows to compute the interpolating polynomial with a simple algorithmicscheme. Compared to the already established schemes of the Padua and Xupoints, the numerical results for the proposed scheme show similar approxi-mation errors and a similar growth of the Lebesgue constant.
Keywords bivariate Lagrange interpolation · quadrature formulas · Chebyshev polynomials · Lissajous curves
A challenging task for multivariate polynomial interpolation is the construc-tion of a suitable set of node points. Depending on the application, these nodepoints should provide a series of favorable properties including a unique inter-polation in given polynomial spaces, a slow growth of the Lebesgue constantand simple algorithmic schemes that compute the interpolating polynomial.The construction of suitable point sets for multivariate interpolation has a
W. Erb, C. Kaethner, M. Ahlborg, T.M. BuzugUniversität zu LübeckRatzeburger Allee 16023562 LübeckE-mail: [email protected]: [email protected]: [email protected]: [email protected] a r X i v : . [ m a t h . NA ] N ov W. Erb et al. long-standing history. For an overview, we refer to the survey articles [13,14]and the references therein. Examples of remarkable constructions in the bi-variate setting are the point sets introduced by Morrow and Patterson [22],Xu [24], as well as some generalizations of them [17]. A modification of theMorrow-Patterson points, introduced as Padua points [7], is particularly in-teresting for the purposes of this article.In some applications, the given data points are lying on subtrajectories ofthe euclidean space. In this case, aside from the above mentioned favorableproperties, it is mandatory that the node points are part of these trajectories.Lissajous curves are particularly interesting examples for us, as they are usedas a sampling path in a young medical imaging technology called MagneticParticle Imaging (MPI) [15].In MPI, the distribution of superparamagnetic iron oxide nanoparticles isreconstructed by measuring the magnetic response of the particles. The mea-surement process is based on the combination of various magnetic fields thatgenerate and move a magnetic field free point through a region of interest.Although different trajectories are possible, this movement is typically per-formed in form of a Lissajous curve [19]. The reconstruction of the particledensity from the data on the Lissajous trajectory is currently done in a verystraight forward way, either by solving a system of linear equations based ona pre-measured system matrix or directly from the measurement data [16].By using multivariate polynomial interpolation on the nodes of the samplingpath, i.e. the Lissajous curve, we hope to obtain a further improvement in thereconstruction process.Of the node points mentioned above, the Padua points, as described in[3], are the ones with the strongest relation to Lissajous curves. They can becharacterized as the node points of a particular degenerate Lissajous figure.Moreover, they satisfy a series of remarkable properties: they can be describedas an affine variety of a polynomial ideal [5], they form a particular Chebyshevlattice [11] and they allow a unique interpolation in the space Π n of bivariatepolynomials of degree n [3]. Furthermore, a simple formula for the Lagrangepolynomials is available and the Lebesgue constants are growing slowly as O (cid:0) log n (cid:1) [3].The aim of this article is to develop, similar to the Padua points, an inter-polation theory for node points on Lissajous curves. To this end, we extendthe generating curve approach as presented in [3] to particular families ofLissajous curves in [ − , . In this article, we will focus on the node pointsof non-degenerate Lissajous curves, which are important for the application inMPI [18]. Not all of the above mentioned properties of the Padua points willbe carried over to the node points of Lissajous figures. However, the resultingtheory will have some interesting resemblences, not only to the theory of thePadua points, but also to the Xu points.We start our investigation by characterizing the node points Lisa n,p of non-degenerate Lissajous curves. Based on the node points
Lisa n,p , we will derivesuitable quadrature formulas for integration with product Chebyshev weightfunctions. Next, we will provide the main theoretical results on bivariate in- agrange interpolation at Lissajous nodes 3 terpolation based on the points
Lisa n,p . We will show that the points
Lisa n,p allow unique interpolation in a properly defined space Π Ln,p of bivariate poly-nomials. Further, we will derive a formula for the fundamental polynomials ofLagrange interpolation. This explicit formula allows to compute the interpo-lating polynomial with a simple algorithmic scheme similar to the one of thePadua points [8]. We conclude this article with some numerical tests for thenew bivariate interpolating schemes. Compared to the established interpolat-ing schemes of the Padua and Xu points, the novel interpolation schemes showsimilar approximation errors and a similar growth of the Lebesgue constant.
In this article, we consider π -periodic Lissajous curves of the form γ n,p : R → R , γ n,p ( t ) = (cid:16) sin( nt ) , sin(( n + p ) t ) (cid:17) , (1)where n and p are positive integers such that n and n + p are relatively prime.Based on the calculations in [1] (see also [20]), the Lissajous curve γ n,p isnon-degenerate if and only if p is odd. In this case, γ n,p : [0 , π ) → R is animmersed plane curve with precisely n ( n + p ) − n − p self-intersection points.In the following, we will always assume that p is odd and sample the Lissajouscurve γ n,p along the n ( n + p ) equidistant points t k := 2 πk n ( n + p ) , k = 1 , . . . , n ( n + p ) . In this way, we get the following set of Lissajous node points:
Lisa n,p := (cid:110) γ n,p ( t k ) : k = 1 , . . . , n ( n + p ) (cid:111) . (2)To characterize the set Lisa n,p , we divide γ n,p ( t k ) for the even and oddintegers k . For this decomposition, we use the fact that n and n + p arerelatively prime. Then, if n is odd, every odd integer k can be written as k = (2 i + 1) n + 2 j ( n + p ) with i, j ∈ Z . If k is even, we can write k = 2 in +(2 j +1)( n + p ) with i, j ∈ Z . If n is even, the same holds with the roles of n and n + p switched. In this way, we get the decomposition Lisa n,p = Lisa b n,p ∪ Lisa w n,p with the sets Lisa b n,p := (cid:26) γ n,p (cid:18) (2 i + 1) n + 2 j ( n + p )4 n ( n + p ) 2 π (cid:19) : i, j ∈ Z (cid:27) , (3) Lisa w n,p := (cid:26) γ n,p (cid:18) in + (2 j + 1)( n + p )4 n ( n + p ) 2 π (cid:19) : i, j ∈ Z (cid:27) . (4)Two examples of Lissajous curves γ n,p with the corresponding node points Lisa b n,p and Lisa w n,p are illustrated in Figure 1. To get a compact representation W. Erb et al. − − . . − − . . (a) Lissajous figure γ , , | Lisa , | = 17 . − − . . − − . . (b) Lissajous figure γ , , | Lisa , | = 27 . Fig. 1
Illustration of non-degenerate Lissajous curves γ n,p . The node points Lisa n,p of γ n,p are arranged on two different grids (black, white) corresponing to the sets Lisa b n,p and Lisa w n,p . of Lisa b n,p and Lisa w n,p , we use the following notation for the Gauß-Lobattopoints: z nk := cos (cid:18) kπn (cid:19) , n ∈ N , k ∈ Z . (5)Then, evaluating the points (3) and (4) explicitly for the Lissajous curve (1),we get the following characterization: Lisa b n,p = (cid:26) ( − i + j (cid:16) z n + p )(2 i +1) p , z n jp (cid:17) : i = 0 , . . . , n + p − j = 0 , . . . , n (cid:27) , (6) Lisa w n,p = (cid:26) ( − i + j (cid:16) z n + p )2 ip , z n (2 j +1) p (cid:17) : i = 0 , . . . , n + pj = 0 , . . . , n − (cid:27) . (7)Since p is assumed to be odd and relatively prime to n , p is relativelyprime to n as well as to n + p ) . Therefore, by rearranging the points, wecan drop the number p in the lower indices of the Gauß-Lobatto points in(6) and (7). Due to the point symmetry of the Lissajous curve γ n,p , the term ( − i + j which preceeds the points in (6) and (7) can also be dropped by furtherrearrangement. This leads to the following simple characterization of the pointsets Lisa b n,p and Lisa w n,p : Lisa b n,p = (cid:26)(cid:16) z n + p )2 i (cid:48) +1 , z n j (cid:48) (cid:17) : i (cid:48) = 0 , . . . , n + p − j (cid:48) = 0 , . . . , n (cid:27) , (8) Lisa w n,p = (cid:26)(cid:16) z n + p )2 i (cid:48) , z n j (cid:48) +1 (cid:17) : i (cid:48) = 0 , . . . , n + pj (cid:48) = 0 , . . . , n − (cid:27) . (9)With this characterization, we can also divide the points Lisa n,p into thesets
Lisa int n,p and
Lisa out n,p denoting the points lying in the interior and on the agrange interpolation at Lissajous nodes 5 boundary of the square [ − , respectively. We have Lisa int n,p := (cid:26)(cid:16) z n + p )2 i (cid:48) +1 , z n j (cid:48) (cid:17) : i (cid:48) = 0 , . . . , n + p − j (cid:48) = 1 , . . . , n − (cid:27) ∪ (cid:26)(cid:16) z n + p )2 i (cid:48) , z n j (cid:48) +1 (cid:17) : i (cid:48) = 1 , . . . , n + p − j (cid:48) = 0 , . . . , n − (cid:27) , Lisa out n,p := (cid:110)(cid:16) z n + p )2 i (cid:48) +1 , ± (cid:17) : i (cid:48) = 0 , . . . , n + p − (cid:111) ∪ (cid:110)(cid:16) ± , z n j (cid:48) +1 (cid:17) : j (cid:48) = 0 , . . . , n − (cid:111) . From the representation of the Lisa points in (8) and (9), it is possible tocount the number of points in the different sets. They are listed in Table 1.
Table 1
Cardinality of the different Lisa sets.Set Number of elements
Lisa n,p n ( n + p ) + 2 n + p Lisa b n,p ( n + 1)( n + p )Lisa w n,p n ( n + p + 1)Lisa int n,p n ( n + p ) − n − p Lisa out n,p n + 2 p From the representation in (3) and (4) and its identification in (6) and (7),we can deduce that γ n,p (cid:18) (2 i + 1) n + 2 j ( n + p )4 n ( n + p ) 2 π (cid:19) = γ n,p (cid:18) (2 i + 1) n − j ( n + p )4 n ( n + p ) 2 π (cid:19) ,γ n,p (cid:18) in + (2 j + 1)( n + p )4 n ( n + p ) 2 π (cid:19) = γ n,p (cid:18) − in + (2 j + 1)( n + p )4 n ( n + p ) 2 π (cid:19) holds for all i, j ∈ Z . Moreover, in (3) and (4) the boundary points are repre-sented by j ∈ n Z and i ∈ ( n + p ) Z , respectively. Thus, for interior points in Lisa b n,p ∩ Lisa int n,p , i.e. all points in (3) satisfying j (cid:54) = n Z , there exist at leasttwo different ≤ k, k (cid:48) ≤ n ( n + p ) in (2) that represent the same point. Thesame holds for all interior points in the second set Lisa w n,p .Therefore, all points in Lisa int n,p are self-intersection points of the Lissajouscurve γ n,p . Since | Lisa int n,p | = 2 n ( n + p ) − n − p corresponds to the total numberof self-intersection points of a non-degenerate Lissajous curve (see [1]), we canconclude that Lisa int n,p is precisely the set of all self-intersection points of theLissajous curve γ n,p . Finally, since | Lisa int n,p | + | Lisa out n,p | = 4 n ( n + p ) , we canalso conclude that there are exactly two different ≤ k, k (cid:48) ≤ n ( n + p ) thatrepresent the same point in Lisa int n,p and that every point in
Lisa out n,p is describedby exactly one ≤ k ≤ n ( n + p ) in (2). W. Erb et al.
In order to identify the different integers k in (2) that describe the samepoint A ∈
Lisa n,p , we introduce for k, k (cid:48) ∈ Z the equivalence relation k Lisa n,p ∼ k (cid:48) ⇔ γ n,p ( t k ) = γ n,p ( t k (cid:48) ) . We say that k ∈ Z belongs to the equivalence class [ A ] , A ∈
Lisa n,p , if γ n,p ( t k ) = A . Therefore, by the above argumentation, there is exactly one ≤ k ≤ n ( n + p ) in the equivalence class [ A ] if A ∈
Lisa out n,p and exactly twoif
A ∈
Lisa int n,p . Remark 1
There are some remarkable relations between the Lisa, Padua andXu points. In formal terms, if p = 0 in the characterization (8) and (9) of theLisa points, the points Lisa n, correspond with the even Xu points XU n asdefined in [24]. Moreover, if p = in (8) and (9), we obtain the even Paduapoints PD n of the second family (see [8] and (20), (21) in Section 6) with aslight adjustment in the range of the indices. A further comparison of thesethree point sets in terms of numerical simulations is given in the last sectionof this article. Finally we would like to add that the Lisa points, similarly tothe Padua points, can be considered as two-dimensional Chebyshev lattices ofrank 1 (see [11]). In this section, we study quadrature rules for bivariate integration defined bypoint evaluations at the points
Lisa n,p . As underlying polynomial spaces in R , we consider Π n = span { T i ( x ) T j ( y ) : i + j ≤ n } , where T i ( x ) denotes the Chebyshev polynomial T i ( x ) = cos( i arccos x ) of thefirst kind. It is well-known (cf. [24]) that { T i ( x ) T j ( y ) : i + j ≤ n } is anorthogonal basis of Π n with respect to the inner product (cid:104) f, g (cid:105) := 1 π (cid:90) − (cid:90) − f ( x, y ) g ( x, y ) 1 √ − x (cid:112) − y d x d y. (10)The corresponding orthonormal basis is given by { ˆ T i ( x ) ˆ T j ( y ) : i + j ≤ n } ,where ˆ T i ( x ) = (cid:26) , if i = 0 , √ T i ( x ) , if i (cid:54) = 0 . Using the trajectory γ n,p , it is possible to reduce a double integral of the formused in (10) into a single integral for a large class of bivariate polynomials. Lemma 1
For all polynomials P ∈ Π n +4 p − with (cid:104) P, T n + p ) ( x ) T n ( y ) (cid:105) = 0 ,the following formula holds: π (cid:90) − (cid:90) − P ( x, y ) 1 √ − x (cid:112) − y d x d y = 12 π (cid:90) π P ( γ n,p ( t ))d t. (11) agrange interpolation at Lissajous nodes 7 Proof
We check (11) for all basis polynomials T i ( x ) T j ( y ) in the space Π n +4 p − .For the left hand side of (11) we get the value if ( i, j ) = (0 , and oth-erwise. For the right hand side of (11) we get also if ( i, j ) = (0 , . For ( i, j ) (cid:54) = (0 , we get for P ( x, y ) = T i ( x ) T j ( y ) the expression π (cid:90) π P ( γ n,p ( t ))d t = 12 π (cid:90) π T i (sin( nt )) T j (sin(( n + p ) t ))d t = 12 π (cid:90) π cos (cid:16) int − i π (cid:17) cos (cid:16) j ( n + p ) t − j π (cid:17) d t. We now determine for which indices ( i, j ) this integral is different from . Thisis only the case if in = j ( n + p ) and i − j is even. Since the numbers n and n + p are relatively prime, this can only be the case if i = k ( n + p ) , j = nk and k ∈ N is an even number. We see that the smallest possible value for k is k = 2 andthe second smallest is k = 4 . Furthermore, the sum of the respective indices isgiven by i + j = (2 n + p ) k . Therefore, we can conclude that for all indices ( i, j ) satisfying i + j = 1 , . . . , n + 2 p − and i + j = 4 n + 2 p + 1 , . . . , n + 4 p − the right hand side of (11) vanishes. If i + j = 4 n + 2 p , the above integral isnonzero only if i = 2( n + p ) and j = 2 n . (cid:117)(cid:116) To get a quadrature formula supported on the points
Lisa n,p , we define asuitable polynomial subspace Π Qn,p = span { T i ( x ) T j ( y ) : ( i, j ) ∈ Γ Qn,p } with the index set Γ Qn,p ⊂ N given by Γ Qn,p := (cid:110) ( i, j ) : i + j ≤ n − (cid:111) ∪ p − (cid:91) m =0 (cid:26) ( i, j ) : i + j = 4 n + m, j < n (4 p − m ) p (cid:27) . Note that the particular index ( i, j ) = (2( n + p ) , n ) is not included in Γ Qn,p andthat Lemma 1 is applicable for all polynomials P ∈ Π Qn,p . An example of theindex set Γ Qn,p is shown in Figure 2. Clearly, the polynomial space Π Qn,p satisfies Π n − ⊂ Π Qn,p ⊂ Π n +4 p − and the dimension of Π Qn,p can be computed as dim Π Qn,p = | Γ Qn,p | = 8 n ( n + p ) + 4 n − p −
1) = 4( | Lisa n,p | − n − p ) − p − . For points
A ∈
Lisa n,p , we define the quadrature weights w A := n ( n + p ) , if A ∈
Lisa out n,p , n ( n + p ) , if A ∈
Lisa int n,p . Then, we get the following quadrature rule based on the node set
Lisa n,p : W. Erb et al. ij Fig. 2
Illustration of the index set Γ Q , with black and white bullets. We have | Γ Q , | = 56 black and white bullets. The black bullets correspond to indices describing the polynomialspace Π n − . The black cross is not contained in Γ Q , . It corresponds to the special index ( i, j ) = (6 , appearing in Lemma 1. Theorem 1
For all P ∈ Π Qn,p the quadrature formula π (cid:90) − (cid:90) − P ( x, y ) 1 √ − x (cid:112) − y d x d y = (cid:88) A∈ Lisa n,p w A P ( A ) (12) is exact.Proof For all trigonometric π -periodic polynomials q of degree less than n ( n + p ) , the following composite trapezoidal quadrature rule is exact: π (cid:90) π q ( t ) dt = 14 n ( n + p ) n ( n + p ) (cid:88) k =1 q ( t k ) . Since Π Qn,p ⊂ Π n +4 p − and Π Qn,p ⊥ T n + p ) ( x ) T n ( y ) , we have by Lemma 1the identity π (cid:90) − (cid:90) − P ( x, y ) 1 √ − x (cid:112) − y d x d y = 12 π (cid:90) π P ( γ n,p ( t ))d t. Thus, if we show that for P ∈ Π Qn,p the trigonometric polynomial P ( γ n,p ( t )) is of degree less than n ( n + p ) , we immediately get the quadrature formula π (cid:90) − (cid:90) − P ( x, y ) 1 √ − x (cid:112) − y d x d y = 14 n ( n + p ) n ( n + p ) (cid:88) k =1 P ( γ n,p ( t k ))= (cid:88) A∈ Lisa n,p w A P ( A ) . agrange interpolation at Lissajous nodes 9 To finish the proof we consider the representation of the polynomial P ∈ Π Qn,p in the orthogonal basis { T i ( x ) T j ( y ) : ( i, j ) ∈ Γ Qn,p } and get P ( γ n,p ( t )) = (cid:88) ( i,j ) ∈ Γ Qn,p a ij T i (sin nt ) T j (sin( n + p ) t )= (cid:88) ( i,j ) ∈ Γ Qn,p a ij cos (cid:0) int − i π (cid:1) cos (cid:0) j ( n + p ) t − j π (cid:1) for some coefficients a ij ∈ R . In order for the trigonometric polynomials inthis formula to have a degree less than n ( n + p ) , the indices ( i, j ) have tosatisfy the condition ( i + j ) n + jp < n ( n + p ) . In the case that i + j < n , we have ( i + j ) n + jp ≤ ( i + j ) n + 4 np < n ( n + p ) and the condition above is satisfied.In the case that i + j = 4 n + m with ≤ m ≤ p − , we have ( i + j ) n + jp =4 n + mn + jp < n ( n + p ) and the condition above is satisfied for all j with j < n (4 p − m ) p . By definition, this condition is exactly satisfied for all indices ( i, j ) ∈ Γ Qn,p and therefore for all polynomials P ∈ Π Qn,p . (cid:117)(cid:116) Remark 2
Lemma 1 and Theorem 1 are generalizations of corresponding re-sults proven in [3] for the Padua points. An analogous formula also existsfor the Xu points (see [22,24]). Furthermore, the cardinality of the Xu points XU n is known to be minimal for exact integration of bivariate polynomialsin Π n − with respect to a product Chebyshev weight function (see [21,24]).Since | Lisa n,p | > | XU n | = 2 n ( n + 1) , this is not the case for the Lisa points.On the other hand, as illustrated in Figure 2, the space Π Qn,p , for which (12) isexact, shows a remarkable asymmetry. As for multivariate interpolation, theconstruction of suitable nodes for cubature rules has a long history. For anoverview, we refer to the survey article [10].
Given the quadrature formulas of the last section, we now investigate bivariateinterpolation at the points
Lisa n,p . The corresponding interpolation problemcan be formulated as follows: for given function values f A ∈ R , A ∈
Lisa n,p ,we want to find a unique bivariate interpolating polynomial L n,p f such that L n,p f ( A ) = f A for all A ∈
Lisa n,p . (13)To set this problem correctly, we have to fix an underlying interpolation space.This space is linked to Π Qn,p and defined as Π Ln,p := span { T i ( x ) T j ( y ) : ( i, j ) ∈ Γ Ln,p } on the index set Γ Ln,p := (cid:110) ( i, j ) : i + j ≤ n (cid:111) ∪ p − (cid:91) m =1 (cid:26) ( i, j ) : i + j = 2 n + m, j < n (2 p − m ) p (cid:27) . ij (a) Index set Γ L , with | Γ L , | = 17 . ij (b) Index set Γ L , with | Γ L , | = 27 . Fig. 3
Illustration of the index sets Γ Ln,p for n = 2 , p = 1 and n = 2 , p = 3 . The blackbullets correspond to indices describing the polynomial space Π n . Examples of sets Γ Ln,p with different values of p are given in Figure 3. Thereproducing kernel K Ln,p : R × R → R of the polynomial space Π Ln,p is givenas K Ln,p ( x A , y A ; x B , y B ) = (cid:88) ( i,j ) ∈ Γ Ln,p ˆ T i ( x A ) ˆ T i ( x B ) ˆ T j ( y A ) ˆ T j ( y B ) . It is straightforward to check that the kernel K Ln,p has the reproducing property (cid:104)
P, K
Ln,p ( x, y ; · ) (cid:105) = P ( x, y ) for all polynomials P ∈ Π Ln,p . We have Π n ⊂ Π Ln,p ⊂ Π n + p ) − . The dimen-sion of the polynomial space Π Ln,p is given as dim( Π Ln,p ) = | Γ Ln,p | = (2 n + 1)(2 n + 2)2 + p − (cid:88) m =1 (cid:24) n (2 p − m ) p (cid:25) = 2 n + n (2 p + 2) + p = | Lisa n,p | . Therefore, the dimension dim( Π Ln,p ) of the polynomial space Π Ln,p correspondsprecisely to the number of distinct points in
Lisa n,p .Soon, we will deduce a formula for the fundamental polynomials of La-grange interpolation with respect to the points in
Lisa n,p and show that the agrange interpolation at Lissajous nodes 11 interpolation problem (13) has a unique solution. To this end, we investigatean isomorphism between the polynomial space Π Ln,p and the subspace Π trig ,L n ( n + p ) := (cid:26) q ∈ Π trig2 n ( n + p ) : q ( t k ) = q ( t k (cid:48) ) for all k, k (cid:48) with k Lisa n,p ∼ k (cid:48) (cid:27) (14)of π -periodic trigonometric polynomials Π trig2 n ( n + p ) := q ( t ) = n ( n + p ) (cid:88) m =0 a m cos( mt ) + n ( n + p ) − (cid:88) m =1 b m sin( mt ) : a m , b m ∈ R . Theorem 2
The operator E γ : Π Ln,p → Π trig ,L n ( n + p ) , E γ P ( t ) = P ( γ n,p ( t )) , t ∈ [0 , π ] , defines an isometric isomorphism from the space (cid:0) Π Ln,p , (cid:104)· , ·(cid:105) (cid:1) onto the space Π trig ,L n ( n + p ) equipped with the inner product (cid:104) q , q (cid:105) = 12 π (cid:90) π q ( t ) q ( t )d t .Proof The system (cid:110) ˆ T i ( x ) ˆ T j ( y ) : ( i, j ) ∈ Γ Ln,p (cid:111) forms an orthonormal basis ofthe space Π Ln,p . The image e ij ( t ) := E γ (cid:16) ˆ T i ( x ) ˆ T j ( y ) (cid:17) ( t ) , ( i, j ) ∈ Γ Ln,p , of this basis under the linear operator E γ is given by e ij ( t ) = , if ( i, j ) = (0 , , √ (cid:0) int − i π (cid:1) , if j = 0 , i < n + p ) , √ (cid:0) j ( n + p ) t − j π (cid:1) , if i = 0 , j ≤ n, (cid:0) int − i π (cid:1) cos (cid:0) j ( n + p ) t − j π (cid:1) , otherwise . (15)For ( i, j ) ∈ Γ Ln,p , j (cid:54) = 2 n , the functions e ij ( t ) are trigonometric polynomialsof degree less than n ( n + p ) . The only trigonometric polynomial of exactdegree n ( n + p ) is precisely e , n . By the definition of the operator E γ , thevalues E γ P ( t k ) and E γ P ( t k (cid:48) ) , t k (cid:54) = t k (cid:48) coincide if γ n,p ( t k ) = γ n,p ( t k (cid:48) ) is a self-intersection point of γ n,p . This is precisely encoded in the constraints given in(14). We can conclude that E γ maps Π Ln,p into the space Π trig ,L n ( n + p ) .For polyonomials P , P ∈ Π Ln,p , the product polynomial P P is an elementof the space Π n +4 p − and satisfies (cid:104) P P , T n + p ) ( x ) T n ( y ) (cid:105) = 0 . Therefore,by Lemma 1, the set (cid:8) e ij : ( i, j ) ∈ Γ Ln,p (cid:9) is an orthonormal system in Π trig ,L n ( n + p ) ,and thus, E γ is an isometric embedding from Π Ln,p into Π trig ,L n ( n + p ) . Now, if we can show that the dimensions of Π Ln,p and Π trig ,L n ( n + p ) coincide,the proof is finished. To this end, we consider in Π trig2 n ( n + p ) the Dirichlet kernel D n ( n + p ) ( t ) := 1 + cos(2 n ( n + p ) t ) + 2 n ( n + p ) − (cid:88) k =1 cos( kt )4 n ( n + p ) = sin(2 n ( n + p ) t ) cos t n ( n + p ) sin t . It is well known that the trigonometric polynomials D k n ( n + p ) ( t ) := D n ( n + p ) ( t − t k ) , k = 1 , . . . , n ( n + p ) , are precisely the Lagrange polynomials in the space Π trig2 n ( n + p ) with respect tothe points t k , k = 1 , . . . , n ( n + p ) , i.e. D k n ( n + p ) ( t k (cid:48) ) = δ k,k (cid:48) , ≤ k, k (cid:48) ≤ n ( n + p ) . In general, the polynomials D k n ( n + p ) do not lie in the subspace Π trig ,L n ( n + p ) .However, we can define a basis for Π trig ,L n ( n + p ) by using the linear combinations l A ( t ) := (cid:88) k =1 ,..., n ( n + p ): k ∈ [ A ] D k n ( n + p ) ( t ) , A ∈
Lisa n,p . (16)Clearly, the polynomials l A are elements of Π trig ,L n ( n + p ) , and l A ( t k ) is equal toone if k ∈ [ A ] and zero if k / ∈ [ A ] . Also, by the orthogonality of the functions D k n ( n + p ) , we have (cid:104) l A , l B (cid:105) = 0 if A (cid:54) = B . Therefore, the system { l A : A ∈
Lisa n,p } forms an orthogonal basis of Π trig ,L n ( n + p ) and dim( Π trig ,L n ( n + p ) ) = | Lisa n,p | .This corresponds exactly with the dimension of the space Π Ln,p . (cid:117)(cid:116) Theorem 3
For A = ( x A , y A ) ∈ Lisa n,p , the polynomials L A := E − γ l A havethe representation L A ( x, y ) = w A (cid:18) K Ln,p ( x, y ; x A , y A ) −
12 ˆ T n ( y ) ˆ T n ( y A ) (cid:19) (17) and are the fundamental polynomials of Lagrange interpolation in the space Π Ln,p on the point set
Lisa n,p , i.e. L A ( B ) = δ A , B , A , B ∈
Lisa n,p . The interpolation problem (13) has a unique solution in Π Ln,p and the interpo-lating polynomial L n,p f is given by L n,p f ( x, y ) = (cid:88) A∈ Lisa n,p f A L A ( x, y ) . agrange interpolation at Lissajous nodes 13 Proof
From the definition (16) of the trigonometric polynomials l A and themapping E γ it follows immediately that the polynomials L A = E − γ l A satisfy L A ( B ) = δ A , B for B ∈
Lisa n,p . Moreover, since the trigonometric polynomials { l A : A ∈
Lisa n,p } form an orthogonal basis of the space Π trig ,L n ( n + p ) , Theorem2 implies that the polynomials { L A : A ∈
Lisa n,p } form an orthogonal basisof Lagrange polynomials for the space Π Ln,p as well.It remains to prove (17). To this end, we compute the decomposition of thepolynomials l A in the basis e ij given in (15) and use the inverse of the operator E γ to obtain (17). The proof will be given only for A ∈
Lisa b n,p having therepresentation A = ( x A , y A ) = γ n,p (cid:18) (2 r + 1) n + 2 s ( n + p )4 n ( n + p ) 2 π (cid:19) = ( − r + s (cid:16) z n + p )(2 r +1) p , z n sp (cid:17) . We first suppose that
A ∈
Lisa b n,p is an interor point such that the two points k, k (cid:48) ∈ [ A ] ∩ [1 , n ( n + p )] , k (cid:54) = k (cid:48) that represent the same A are given as k = (2 r + 1) n + 2 s ( n + p ) mod 4 n ( n + p ) ,k (cid:48) = (2 r + 1) n − s ( n + p ) mod 4 n ( n + p ) . Using simple trigonometric transformations, the basis function l A can be writ-ten as l A ( t ) = D k n ( n + p ) ( t ) + D k (cid:48) n ( n + p ) ( t )= 24 n ( n + p ) r + 1) nπ ) cos(2 n ( n + p ) t )+ 2 n ( n + p ) − (cid:88) m =1 cos (cid:0) smπ n (cid:1)(cid:16) cos (cid:16) (2 r +1) mπ n + p ) (cid:17) cos( mt ) + sin (cid:16) (2 r +1) mπ n + p ) (cid:17) sin( mt ) (cid:17) . Now, using the explicit expression (15) of the basis polynomials e ij and com-paring the coefficients in the decomposition of l A , we get the following formulafor the inner product (cid:104) l A , e ij (cid:105) = π (cid:82) π l A ( t ) e ij ( t )d t , ( i, j ) ∈ Γ Ln,p : (cid:104) l A , e ij (cid:105) = n ( n + p ) , if ( i, j ) = (0 , , √ n ( n + p ) , if i = 0 , j = 2 n, √ − ( r + s ) j n ( n + p ) cos (cid:0) j spπ n (cid:1) , if i = 0 , j < n, √ − ( r + s ) i n ( n + p ) cos (cid:16) i (2 r +1) pπ n + p ) (cid:17) , if i (cid:54) = 0 , j = 0 , − ( r + s )( i + j ) n ( n + p ) cos (cid:16) i (2 r +1) pπ n + p ) (cid:17) cos (cid:0) j spπ n (cid:1) , otherwise, = 24 n ( n + p ) (cid:40) ˆ T n ( y A ) , if i = 0 , j = 2 n, ˆ T i ( x A ) ˆ T j ( y A ) , if ( i, j ) ∈ Γ Ln,p \ (0 , n ) . Therefore, l A ( t ) can be decomposed as l A ( t ) = ˆ T n ( y A )4 n ( n + p ) e , n ( t ) + (cid:88) ( i,j ) ∈ ΓLn,pj (cid:54) =2 n T i ( x A ) ˆ T j ( y A )4 n ( n + p ) e ij ( t ) . Now, using the inverse mapping E − γ together with the definition of the repro-ducing kernel K Ln,p , we can conclude: L A ( x, y ) = E − γ l A ( x, y ) = w A (cid:18) K Ln,p ( x, y ; x A , y A ) −
12 ˆ T n ( y ) ˆ T n ( y A ) (cid:19) . If A ∈
Lisa b n,p is a point on the boundary of the square [ − , , the number k can be represented as k = (2 r + 1) n and the basis function l A is given as l A ( t ) = D k n ( n + p ) ( t ) = 14 n ( n + p ) r + 1) nπ ) cos(2 n ( n + p ) t )+ 2 n ( n + p ) − (cid:88) m =1 cos (cid:18) r + 12( n + p ) mπ (cid:19) cos( mt ) + sin (cid:18) r + 12( n + p ) mπ (cid:19) sin( mt ) . Now, similar calculations to the above yield (17) with the half sized weightfunction w A = n ( n + p ) . Finally, for all points A ∈
Lisa w n,p , (17) can be obtainedby analogous calculations using the representation (4) instead of (3). (cid:117)(cid:116) Remark 3 (17) has a remarkable resemblence to the Lagrange polynomials ofthe Padua points. For the Padua points, the analog statement of Theorem 3can be proved very elegantly by using ideal theory (cf. [5]). This approach was,however, not successful for the more general Lissajous nodes. Here, we had touse the isomorphism E γ and Theorem 2 instead. In view of Theorem 3, the solution to the interpolation problem (13) in Π Ln,p is given as L n,p f ( x, y ) = (cid:88) A∈ Lisa n,p w A f A (cid:18) K Ln,p ( x, y ; x A , y A ) −
12 ˆ T n ( y ) ˆ T n ( y A ) (cid:19) . The representation of the polynomial L n,p f ( x, y ) in the orthonormal Cheby-shev basis { ˆ T i ( x ) ˆ T j ( y ) : ( i, j ) ∈ Γ Ln,p } can now be written as L n,p f ( x, y ) = (cid:88) ( i,j ) ∈ Γ Ln,p c i,j ˆ T i ( x ) ˆ T j ( y ) agrange interpolation at Lissajous nodes 15 with the coefficients c i,j = (cid:104)L n,p f, ˆ T i ( x ) ˆ T j ( y ) (cid:105) given by c i,j = (cid:88) A∈ Lisa n,p w A f A ˆ T i ( x A ) ˆ T j ( y A ) , if ( i, j ) ∈ Γ Ln,p \ (0 , n ) , (cid:88) A∈ Lisa n,p w A f A ˆ T n ( y A ) , if ( i, j ) = (0 , n ) . Using a matrix formulation, this identity can be written more compactly. Weintroduce the coefficient matrix C n,p = ( c ij ) ∈ R n + p ) × n +1 by c ij = (cid:26) (cid:104)L n,p f, ˆ T i ( x ) ˆ T j ( y ) , if ( i, j ) ∈ Γ Ln,p , , otherwise . Next, we define the diagonal matrix D f (Lisa n,p ) = diag ( w A f A , A ∈
Lisa n,p ) ∈ R | Lisa n,p |×|
Lisa n,p | . Further, for a general finite set
X ⊂ R of points, we introduce the matrices T x ( X ) = ˆ T ( x B ) · · · ... · · · ˆ T n + p ) − ( x B ) (cid:124) (cid:123)(cid:122) (cid:125) B∈X ∈ R n + p ) ×|X | , T y ( X ) = ˆ T ( y B ) · · · ... · · · ˆ T n ( y B ) (cid:124) (cid:123)(cid:122) (cid:125) B∈X ∈ R n +1 ×|X | . Finally, we define the mask M n,p = ( m ij ) ∈ R n + p ) × n +1 by m i,j = , if ( i, j ) ∈ Γ Ln,p \ (0 , n ) , / , if ( i, j ) = (0 , n ) , , if ( i, j ) / ∈ Γ Ln,p . Now, the coefficient matrix C n,p of the interpolating polynomial can be com-puted as C n,p = (cid:0) T x (Lisa n,p ) D f (Lisa n,p ) T y (Lisa n,p ) T (cid:1) (cid:12) M n,p , (18)where (cid:12) denotes pointwise multiplication of the matrix entries. For an arbi-trary point B ⊂ R , the evaluation L n,p f ( B ) of the interpolation polynomial L n,p f at B is then given by L n,p f ( B ) = T x ( B ) T C n,p T y ( B ) . (19) Remark 4
The matrix formulation in (18) and (19) is almost identical to theformulation of the interpolating scheme of the Padua points given in [8]. Thisis due to the similarity in the representation (17) of the Lagrange polynomials.The main difference between the schemes lies in the form of the mask M n,p .The mask M n,p for the Lisa points has an asymmetric structure determinedby the index set Γ Ln,p , whereas the matrix is an upper left triangular matrixfor the Padua points. Two examples of such a structure are given in Figure 3.
Based on the results derived in the last sections, we perform numerical simula-tions on the behaviour of the
Lisa points ( LS ) in comparison to some alreadyestablished point sets. Unless explicitly mentioned, we assume p = 1 for allnumerical simulations of the Lisa points. For the comparison point sets, ourfocus is on the Xu points (XU) [24] and the Padua points (PD) [7]. Based onthe Chebychev-Lobatto points given by (5), and in correspondance to (8) and(9), the odd Xu points XU n +1 are defined as the union of the setsXU b2 n +1 = (cid:8) ( z n +12 i , z n +12 j ) : 0 ≤ i ≤ n, ≤ j ≤ n (cid:9) , XU w2 n +1 = (cid:8) ( z n +12 i +1 , z n +12 j +1 ) : 0 ≤ i ≤ n, ≤ j ≤ n (cid:9) , with the cardinality | XU n +1 | = 2( n + 1) . In turn, the even Padua points PD n (2nd family) are defined as the union of the setsPD b2 n = (cid:8) ( z n +12 i +1 , z n j ) : 0 ≤ i ≤ n, ≤ j ≤ n (cid:9) , (20)PD w2 n = (cid:8) ( z n +12 i , z n j +1 ) : 0 ≤ i ≤ n, ≤ j ≤ n − (cid:9) . (21)The cardinality can be calculated as | PD n | = ( n +1)(2 n +1) . The distributionsof the Lisa , Xu and Padua points are shown for small degrees of n in Figure 4.The point sets are introduced in such a way that an equally chosen n resultsin a similar cardinality.The stability of the mapping f → L n,p f is evaluated by means of thegrowth of the Lebesgue constant. Here, we calculate the values of the Lebesgueconstant Λ LS n,p = max B∈ [ − , (cid:88) A∈ Lisa n,p | L A ( B ) | of the Lisa points up to a degree of n = 60 . We compare them with the least-squares fitting of the Lebesgue constant for the Padua and the Xu points. Asshown in [7], it holds for the Padua points that Λ PD n = ( π log(2 n + 1) + 1 . and as presented for the Xu points in [2] that Λ XU n +1 = ( π log(2 n + 2)) . Fig-ure 5(a) indicates that the asymptotic growth of Λ LS n, corresponds to the order O (cid:0) log n (cid:1) of the Lebesgue constant Λ PD2 n . In Figure 5(b) it is shown, how avariation of the parameter p of the Lisa points changes the growth of theLebesgue constant. Here, we consider p = { , , , } and excluded each entryfor n and p not being relatively prime. In total, these numerical evaluations agrange interpolation at Lissajous nodes 17 − − . . − − . .
51 LSXUPD (a) Point sets for n = 2 . − − . . − − . .
51 LSXUPD (b) Point sets for n = 5 . Fig. 4
Visualizations of the
Lisa ( LS ), Xu (XU) and Padua (PD) point sets. suggest the conjecture that the Lebesgue constant of the Lisa points is of thesame order O (cid:0) log n (cid:1) as the Lebesgue constant of the Padua and Xu points(see [4]). n L e b e s gu ec on s t a n t LSXUPD (a)
Lisa n, , XU n +1 and PD n point sets. n L e b e s gu ec on s t a n t p = p = p = p = (b) Lisa point sets for p ∈ { , , , } . Fig. 5
Lebesgue constants up to a degree of n = 60 for the Lisa points in comparison tothe least-squares fitting of the Lebesgue constant of the Xu and Padua points.
For a further evaluation of the
Lisa points, we perform numerical interpo-lations with the Xu, Padua and
Lisa points on the Franke-Renka-Brown testset [12,23]. In order to simulate the Xu points as well as the Padua points, thenumerical algorithms presented in [6,9] are used. The maximum interpolationerrors are computed on a uniform grid of × points defined in a region Ω = [0 , × [0 , . As mentioned above, the degree n is defined to result in asimilar total number of points, i.e. a similar cardinality. For our simulationswe take n ∈ { , , , } . The results are shown in Table 2–4. It can be seenthat the maximum interpolation error of all three point sets shows a similar behaviour in terms of degree n , with respect to the chosen test function. Interms of the Lisa point sets, we evaluated the behaviour of p in addition tothe aforementioned comparisons. We can state that the influence of varying p ,with respect to the maximum interpolation error and the nodes used for theevaluation, is almost negligible. Table 2
Interpolation errors for the points
Lisa n, . n F F F F F F F F F F Table 3
Interpolation errors for the points XU n +1 . n F F F F F F F F F F Table 4
Interpolation errors for the Padua points PD n . n F F F F F F F F F F Acknowledgements
The authors gratefully acknowledge the financial support of the Ger-man Federal Ministry of Education and Research (BMBF, grant number 13N11090), theGerman Research Foundation (DFG, grant number BU 1436/9-1 and ER 777/1-1), theEuropean Union and the State Schleswig-Holstein (EFRE, grant number 122-10-004).agrange interpolation at Lissajous nodes 19
References
1. Bogle, M.G.V., Hearst, J.E., Jones, V.F.R., Stoilov, L.: Lissajous knots. J. Knot TheoryRamifications (2), 121–140 (1994)2. Bos, L., Caliari, M., De Marchi, S., Vianello, M.: A Numerical Study of the Xu Poly-nomial Interpolation Formula in Two Variables. Computing (3), 311–324 (2006)3. Bos, L., Caliari, M., De Marchi, S., Vianello, M., Xu, Y.: Bivariate Lagrange interpo-lation at the Padua points: the generating curve approach. J. Approx. Theory (1),15–25 (2006)4. Bos, L., De Marchi, S., Vianello, M.: On the Lebesgue constant for the Xu interpolationformula. J. Approx. Theory (2), 134–141 (2006)5. Bos, L., De Marchi, S., Vianello, M., Xu, Y.: Bivariate Lagrange interpolation at thePadua points: The ideal theory approach. Numer. Math. (1), 43–57 (2007)6. Caliari, M., De Marchi, S., Sommariva, A., Vianello, M.: Padua2DM: fast interpolationand cubature at the Padua points in Matlab/Octave. Numer. Algorithms (1), 45–60(2011)7. Caliari, M., De Marchi, S., Vianello, M.: Bivariate polynomial interpolation on thesquare at new nodal sets. Appl. Math. Comput. (2), 261–274 (2005)8. Caliari, M., De Marchi, S., Vianello, M.: Bivariate Lagrange interpolation at the Paduapoints: Computational aspects. J. Comput. Appl. Math. (2), 284–292 (2008)9. Caliari, M., Vianello, M., De Marchi, S., Montagna, R.: Hyper2d: a numerical code forhyperinterpolation on rectangles. Appl. Math. Comput. (2), 1138–1147 (2006)10. Cools, R.: Constructing cubature formulae: the science behind the art. Acta Numerica , 1–54 (1997)11. Cools, R., Poppe, K.: Chebyshev lattices, a unifying framework for cubature with Cheby-shev weight function. BIT (2), 275–288 (2011)12. Franke, R.: Scattered data interpolation: Tests of some methods. Math. Comp. (157),181–200 (1982)13. Gasca, M., Sauer, T.: On the history of multivariate polynomial interpolation. J. Com-put. Appl. Math. (1-2), 23–35 (2000)14. Gasca, M., Sauer, T.: Polynomial interpolation in several variables. Adv. Comput.Math. (4), 377–410 (2000)15. Gleich, B., Weizenecker, J.: Tomographic imaging using the nonlinear response of mag-netic particles. Nature (7046), 1214–1217 (2005)16. Grüttner, M., Knopp, T., Franke, J., Heidenreich, M., Rahmer, J., Halkola, A., Kaeth-ner, C., Borgert, J., Buzug, T.M.: On the formulation of the image reconstruction prob-lem in magnetic particle imaging. Biomed. Tech. / Biomed. Eng. (6), 583–591 (2013)17. Harris, L.A.: Bivariate Lagrange interpolation at the Geronimus nodes. Contemp. Math. , 135–147 (2013)18. Kaethner, C., Ahlborg, M., Bringout, G., Weber, M., Buzug, T.M.: Axially elongatedfield-free point data acquisition in magnetic particle imaging. IEEE Trans. Med. Imag.(2014). Accepted for publication19. Knopp, T., Biederer, S., Sattel, T.F., Weizenecker, J., Gleich, B., Borgert, J., Buzug,T.M.: Trajectory analysis for magnetic particle imaging. Phys. Med. Biol. (2), 385–3971 (2009)20. Lamm, C.: There are infinitely many Lissajous knots. Manuscr. Math. (1), 29–37(1997)21. Möller, H.M.: Kubaturformeln mit minimaler Knotenzahl. Numer. Math. , 185–200(1976)22. Morrow, C.R., Patterson, T.N.L.: Construction of algebraic cubature rules using poly-nomial ideal theory. SIAM J. Numer. Anal. , 953–976 (1978)23. Renka, R.J., Brown, R.: Algorithm 792: Accuracy tests of acm algorithms for interpo-lation of scattered data in the plane. ACM Trans. Math. Softw. (1), 78–94 (1999)24. Xu, Y.: Lagrange interpolation on Chebyshev points of two variables. J. Approx. Theory87