Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 8s (2024) 724 https://internationalpubls.com Function Quadrature Using Quartic Spline Interpolation Kallur V Vijayakumar 1*, Hariprasad A S 2, Shilpa R3 1 Kallur V Vijayakumar, BMS Institute of Technology and Management, Bengaluru-560064. kallurvijayakumar@bmsit.in 2Hariprasad A S, Rajarajeswari College of Engineering, Bengaluru-560074. ashariprasad@yahoo.co.in 3Shilpa R, Research Scholar, Rajarajeswari College of Engineering, Bengaluru-560074. Shilpa8906@gmail.com* corresponding author Article History: Received: 20-04-2024 Revised: 10-06-2024 Accepted: 24-06-2024 Abstract: This research paper explores the application of quartic splines for numerical integration, particularly in situations with curved boundaries in two-dimensional space. The paper dives into the concept of quartic splines, which are piecewise polynomial functions where each piece is a polynomial of the fourth degree. These splines are designed to ensure continuity of the function itself, its first, second, and third derivatives across the boundaries between the pieces. The core technique involves using quartic splines to approximate the curved boundary of a two-dimensional domain. This allows for the calculation of integrals defined over that domain by transforming it into a simpler shape (often a rectangle) where traditional numerical integration methods like Gauss-Legendre quadrature can be applied effectively. The paper validates the method by applying it to various curved domains and comparing the obtained integral values with known analytical solutions (if available). The absolute errors are calculated and tabulated, demonstrating the accuracy of the approach. The paper showcases the efficiency of the method by solving numerical integration problems from existing literature. These examples serve as evidence that the method is effective and can handle real-world scenarios. A significant advantage of this method is its ability to handle situations where only coordinate data for the boundary is available. Quartic splines can accurately reconstruct the curved boundary from these coordinates, enabling the integration process. Keywords: Numerical Integration, Quadrature, Spline, curved boundary, spline. Introduction Numerical integration, the art of approximating definite integrals through computational methods, finds widespread application in various scientific and engineering disciplines. Among different techniques, quadrature using quartic splines emerges as a robust and accurate approach, offering numerous benefits over earlier traditional approaches. In this comprehensive discussion, we delve into the theoretical foundations, practical implementation, and performance analysis of quartic spline quadrature, equipping you with the knowledge to confidently apply it to your own integration problems. Definite integrals represent the area under a curve within a specified interval. However, analytically solving integrals becomes cumbersome or impossible for many functions. Numerical integration bridges this gap by approximating the area through weighted evaluations of the function at specific points, known as nodes. The objective is to find a formula that expresses the integral as a Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 8s (2024) 725 https://internationalpubls.com linear combination of these function values, minimizing the error between the approximation and the true integral. Spline functions offer flexibility in approximating complex curves. Quartic splines, polynomials of degree four, can closely fit diverse function behaviors through piecewise construction. By partitioning the interval of integration into subintervals, we can represent the integral as the sum of integrals over each subinterval, which we further approximate using quartic splines. Splines are piecewise polynomial functions, meaning they consist of multiple polynomial segments joined smoothly at specific points called knots. This construction allows them to capture intricate curve shapes that single polynomials might struggle to represent accurately. Quartic splines, in particular, use polynomials of degree four, providing a balance between flexibility and computational efficiency. They can model a wide variety of function behaviors with relatively few knots. A quartic spline is a type of piecewise polynomial function that is used for interpolation or approximation of data points over an interval. Specifically, it is defined on an interval [a, b] and is characterized by the following properties: Continuity: A quartic spline is continuous at each knot, which are the distinct points in the interval [a, b]. This means that the function and its derivatives up to the fourth order are continuous across the entire interval. Piecewise Quartic Polynomials: Within each subinterval determined by the knots, the quartic spline is represented by a quartic (degree-4) polynomial. The use of piecewise quartic polynomials allows the spline to flexibly capture and approximate the behavior of the underlying data within each interval. Interpolatory Conditions: The quartic spline satisfies interpolatory conditions, meaning that it accurately passes through specified data points. In other words, the function coincides with given data values at the knots, ensuring that the spline interpolates the data over the entire interval. Uniqueness: In certain conditions, there exists a unique quartic spline that fulfills the continuity, piecewise quartic polynomial, and interpolatory conditions. This uniqueness ensures that the spline provides a distinct representation of the data within the specified interval. In summary, a quartic spline is a smooth and continuous piecewise polynomial function of degree 4, designed to interpolate or approximate data points over an interval [a, b]. Quartic splines use polynomials of degree four, providing a balance between flexibility and computational efficiency. They can model a wide variety of function behaviors with relatively few knots. Knot Placement: We strategically place knots, points demarcating the subintervals, to influence the spline's shape and accuracy. Different strategies, such as equispaced, Chebyshev, or adaptive knot placement, exist, each with its own strengths and weaknesses. Within each subinterval, we construct a quartic spline that interpolates or passes through the function values at specific points within the interval. These points, typically chosen as function extrema or points where the integral is accurately known, further enhance the approximation. We integrate the constructed quartic splines over their respective subintervals and express the overall integral as a weighted sum of these sub- integrals. The weights depend on the chosen knots and spline construction method. Partitioning the Interval: We divide the interval of integration [a, b] into n subintervals using knots xi: 𝑎 = 𝑥0 < 𝑥1 <. . … . . < 𝑥𝑛−1 < 𝑥𝑛 = 𝑏. Constructing Quartic Splines: Within each subinterval [𝑥𝑖 , 𝑥𝑖+1], we construct a quartic spline 𝑆𝑖(𝑥) that interpolates or approximates the function 𝑓(𝑥) at specific points within that subinterval. We approximate the integral of 𝑓(𝑥) over each subinterval Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 8s (2024) 726 https://internationalpubls.com using the integral of the corresponding quartic spline: ∫ 𝑓(𝑥)𝑑𝑥 𝑥𝑖+1 𝑥𝑖 ≈ ∫ 𝑆𝑖(𝑥)𝑑𝑥 𝑥𝑖+1 𝑥𝑖 . The overall integral is then approximated as the sum of these sub-integrals: ∫ 𝑓(𝑥)𝑑𝑥 𝑏 𝑎 ≈ ∑ ∫ 𝑆𝑖(𝑥)𝑑𝑥 𝑥𝑖+1 𝑥𝑖 𝑛−1 𝑖=0 . Quartic Splines can closely fit diverse function behaviors, leading to more accurate approximations than lower-degree polynomials. They ensure continuity and smoothness of the approximation, even at knot points, making them suitable for applications requiring continuous derivatives. The number of knots and their placement can be adjusted to refine the approximation in critical regions, offering flexibility for different problem complexities. The integro cubic spline methods over a uniform partition were studied in [5,23], but their error orders are lower. Later, an integro quintic spline approach over a uniform partition was discussed in [6]. In The present study focus is on creating spline functions that interpolate a given function's first derivatives at specific points (knots) and match the second derivatives between those knots. The derivation of these quartic splines in a standard form. An interesting approach is taken, Instead of solving for the spline function values at the knots, the second derivatives are chosen to coincide with the initial knots' values. References [12] and [22] likely provide more information on this technique. A method for numerically integrating over curved two-dimensional domains. The theoretical foundation relies on several concepts: Green's theorem for boundary integration [2] (likely referring to a specific reference on vector calculus), Gauss-Legendre quadrature rule, a popular technique for numerical integration [1], Parametric relations describing the boundary curve (either exact form or obtained through spline interpolation polynomials) - references [3], [4], [11], [14], [15], [17], [19], might be relevant for these methods. The application of the method on a curved domain shaped like a lunar model. The integration of challenging functions over this curved domain is performed using both explicit parametric relations and parametric quartic spline interpolating polynomials. This highlights the advantage of splines when only coordinate data is available along the boundary. Another application of the numerical integration scheme. Here, it's used to compute the geometric moments of a general ellipse. Overall, this chapter presents a valuable approach for numerical integration over complex, curved boundaries. By leveraging quartic splines, the method can handle scenarios where only coordinate data is available. Quartic spline interpolant Definition: Let f be a function defined on the interval [a, b] and let the knots   1 0 N i i x + = 0 1 2 1 .... N N a x x x x x b + =      = be the ( 2)N + distinct points. Note that ' i x s divide [a, b] into (N+1) subintervals. A function S is said to be a quartic spline on the interval [a, b], is , , andS S S S   are continuous in [a, b] and S is a polynomial of degree  4 in each knot interval 1 [ , ] i i x x − . We shall now present the quasi – Hermite lacunary quartic spline which interpolates the first derivatives of a given function at the knots and the second derivatives between them as proposed in [12, 22]. Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 8s (2024) 727 https://internationalpubls.com Theorem: Given the real numbers ( ) ( 0,1,2,.., 1), i i f x f i N = = + ( ) ( 0,1,2,.., ), i i x h f i N + + = = 0 0 ( ) ,f x f= and 1 1 ( ) N N f f x + + = , then there exists a unique quartic spline S(x) such that (1) On each subinterval 1 [ , ], 0,1,2,.., , i i x x i N + = S(x) coincides with the quartic polynomial 2 3 4( ) ( ) i i i i i i S x S u a bu c u d u e u= = + + + + , 1 ( ) / , ( ) 0 1 i i i u x x h h x x u + = − = −   (1) (2) , andS S S  satisfy the interpolatory conditions. 0 0 1 1 , ( 0,1,2,.., 1) , ( 0,1,2,.., ) , i i i i N N S f i N S f i N S f S f  + + + +  = = +  = = = = (2) Now, we shall prove the above theorem. The proof is constructive: it gives on a algorithm to compute the spline. Proof: We consider the quartic spline ( ) i S x over the interval 1 [ , ] i i x x + defined in eqn (1). On satisfying the interpolatory conditions of eqn (2), we obtained 1 1 1 1 2 2 ( ) (0) ( ) (1) ( ) (0) / ( ) (1) / 2 / 3 / 4 / ( ) (2 6 12 ) / i i i i i i i i i i i i i i i i i i i i i i i i i i i i S x S a f S x S a b c d e f S x S b h f S x S b h c h d h e h f S x h c d e h f     + + + + + = = = = = + + + + =  = = =  = = + + + =  + = + + = (3) Solving the above linear equations in eqn (3), we obtain 2 2 1 2 2 2 1 2 2 2 2 2 1 1 1 , , ( 36 24 ) (36 24 ) ( 24 24 6 ) ( 12 6 ) (24 8) ( 24 8) (12 6) (12 2) 2 ( 12 6) (12 6) ( 6 4) i i i i i i i i i i i i i i i i i i i i i a f b hf A c f f hf hf h f A d f f hf hf h f A e f f hf hf                   + + + + + + + = = = − + + −   + − + − + − + +   = − + − + + − + − −  = − + + − + − + + 2 1 2 ( 6 2) , 2(6 6 1) i h f A    + + − + + = − + (4) On substituting for ai, bi, ci, di and ei in eqn (1), we obtain the following alternative expression Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 8s (2024) 728 https://internationalpubls.com 2 2 2 2 2 2 1 2 2 2 2 2 2 2 2 2 1 3(2 1) 3(2 1) ( ) (1 ) (1 2 ) (1 ) + (3 2 ) (1 ) (3 2) (3 1) (1 ) + (1 ) (1 ) (1 ) (1 ) 2 i i i i i i S u f u u u u f u u u u u u hf u u u u hf u u u u h f      + + + − −    = − + − − − + −         − − −      − − − − − + − +             (5) where, 1 2 ( ) / , ( ), (6 6 1) i i i i u x x h h x x   + = − = −   = − + (6) It is clear that Si(u) is not defined for 0= , i.e., when (3 3) / 6 =  . Now using the continuity of the second derivatives over the intervals [xi-1, xi] and [xi, xi+1], we obtain the following relation as in [22].   1 1 2 1 1 2 1 ( 1)(3 1) (1 2 ) (3 2) ( 1)(2 1) (8 8 1) (2 1) 2 , 1,2,.., i i i i i i i i f f f h f f f h f f i N             − + − + − + + − − − + − + −    = − − + − + + −   + − = (7) The above system of N linear equations will determine f1, f2,..,fN, since f0 and fN+1,   1 0 N i i f + =  and   0 N i i f + =  are already known. We see that (3 3) / 6   and for some simple values of 0,1/ 3.2 / 3,1 = , the coefficient matrix in the linear system of eqn (7) is nonsingular. We shall now consider the simple case =0 only. Case =0 The following explicit expression over the interval 1 [ , ] i i x x + represent the quartic spline interpolant, 2 3 4( ) ,.... i i i i i i S u a bu c u d u e u= + + + + (8) where 1 2 2 1 1 2 1 1 ( ) / , ( ), , , , 4 4 3 , 3 3 2 / 2 i i i i i i i i i i i i i i i i i i i i i u x x h h x x a f b hf c h f d f f hf hf h f e f f hf hf h f + + + + + = − = −   = = =  = − + − − −   = − + + + (9) and,   1 N i i f = can be evaluated by the below said recurrence relation for the quartic spline from eqn(5.2.7), 2 1 1 1 1 ( / 2)( ) ( /12)( ) ( 1,2,..., ) i i i i i i f f h f f h f f i N − − + −    = + + + + = (10) by forward recurrence computing. Numerical integration over curved domains in 2-space Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 8s (2024) 729 https://internationalpubls.com Consider the integral [44, 3] ( , ) xy def II f x y dxdy  =  (11) where, xy  is the curved domain. Let’s denote the closed boundary of the domain xy  , by xy  . We discretise this boundary xy  into N curved arcs by points (( , ), 1,2,.., 1). k k x y k N= + Let the i-th arc be joined by the vertices ( , ) i i i V x y= and 1 1 1 ( , ) i i i V x y + + + = and denoted as , 1i i  +  . Thus, we can now represent the closed boundary xy  as a sum of N curved arcs: , 1 1 N xy i i i   + =  =  (12) Here, each , 1i i  +  signifies a specific curved arc bounded by the vertices ( , ) i i i V x y= and 1 1 1 ( , ) i i i V x y + + + = . Hence, by the use of the boundary integration formula of Green’s theorem for the closed boundary xy  and the relation of eqn (12), we obtain: , 1 1 I I = ( , ) ( , ) xy i i N i f x y dxdy x y dy   + =  =   (13) where ( , ) ( , ) x x y f u y du   =  (14) in which  is fixed, i x  and arbitrary. Let us now choose the replacement 2 2 x x u s  − +    = +        (15) using for u from eqn (15) into eqn (14), we have 1 1 ( , ) , 2 2 2 x x x x y f s y ds    −  − − +       = +               (16) Further by using ( , )x y from eqn (16) into eqn (13), we obtain , 1 , 1 1 1 1 1 1 1 I I ( , ) ˆ ˆ ˆ( ) ( ) ( ) ˆ , ( ) 2 2 2 ˆ ˆ ˆ( ) ( ) ( ) ˆ , ( ) 2 2 2 i i i i i i N i N i i i i i yN i i i i i y x y dy x t x t x t f s y t ds dy x t x t x t f s y t ds         + + + =  =  − = − =    − − +      = +                 − − +      = +                     1 1 ( ) i y t dt        (17) Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 8s (2024) 730 https://internationalpubls.com where, ˆ ( ) i x t and ˆ ( ) i y t refer to the appropriate parametric form over the curved boundary arc joining the points ( , ) i i x y and 1 1 ( , ) i i x y + + Now letting 1 1 2 2 i i i i y y y y t T+ + − +    = +        (18) and substituting in equation (17), we finally obtain : 1 1 1 1 1 1 ( ) ( ) ( ) , ( ) ( ) 2 2 2 2 N i i i i i i i i X T X T X T y y I f s Y T ds Y T dT    + = − −   − − + −        = +                      (19) where 1 1 1 1 ˆ( ) , 2 2 ˆ( ) 2 2 i i i i i i i i i i i i y y y y X T x T y y y y Y T y T + + + +  − +    = +           − +    = +          (20) We shall now apply the Gauss Legendre Quadrature rules of order n and m along the directions s and T respectively. This gives 1 1 1 1 ( ) ( ) ( ) ( ) , ( ) 2 2 2 2 m m mN n m n m m n mi k i i i k i k j k i k j i k i j k X T y y X T X T I W W Y T f s Y T    + = = =       − − − +  = +                 (21) From the above formula, it is quite obvious that m>n , for a bivariate polynomial function of degree less than or equal to (2n-1). The parametric form of ( )( ), ( ) i i X T Y T must be of degree greater than or equal to two for a curved boundary, Hence in equation (21), we must have 𝑚 ≅ (𝑛 + 𝑝) for the curved boundary, as p–point Gauss Legendre quadrature rule is necessary to integrate the following product in equation (21). ( )( )( ) / 2 ( )m m i k i k X T Y T − a polynomial of degree, say, (2p-1). On the other hand, we have chosen n–th order Gauss Legendre Quadrature rule to integrate the bivariate function ( )( ) ( )' ' 1 ( ) / 2 ( ) / 2, ( )m m i k i k i k f X T s X T Y T − + + as polynomial degree (2n-1) in s in equation (21). Application example: a lunar model We shall now consider the curved domain in the shape of a lunar model whose boundary is composed of two circular arcs. The outer circular arc satisfies the equation: ( ) ( ) 2 2 1/ 2 1/ 2 1/ 4x y− + − = and the inner circular arc satisfies the equation 2 2 1/ 4x y+ = . This is shown in figure: Fig 1. Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 8s (2024) 731 https://internationalpubls.com Figure 1: Lunar Model The numerical scheme developed in the previous section relies on having suitable parametric equations that describe the curved boundary of the domain xy  , as specified by equations (17) to (21). In some scenarios, the boundary of the curved domain can be precisely defined by mathematical relations. In such cases, it's feasible to transform these mathematical relations into explicit parametric equations. However, there are situations where the boundary curve is described by a set of discrete coordinate points, making it necessary to generate parametric equations through alternative methods. One such method involves utilizing spline interpolating polynomials of varying orders (e.g., 2nd, 3rd, 4th order, etc.) to generate parametric equations for the curved boundary. This approach proves to be advantageous even when dealing with a boundary curve represented by discrete coordinate points. Now, let's delve into the application of these two methods in the context of a lunar model, as discussed previously. Exact Mathematical Relations Method: If the boundary of the lunar model can be precisely defined by mathematical relations, we can derive explicit parametric equations directly from these relations. This method is straightforward and offers accurate representations of the boundary. However, it relies on the availability of exact mathematical descriptions of the boundary. Spline Interpolating Polynomials Method: In cases where the boundary of the lunar model is not described by exact mathematical relations but rather by a set of discrete coordinate points, spline interpolating polynomials can be employed to generate parametric equations. Spline interpolation allows us to construct a smooth curve that passes through these discrete points, providing a continuous representation of the boundary. The choice of spline order (e.g., 2nd, 3rd, 4th order) depends on the desired level of smoothness and accuracy required for the boundary representation. By applying these two methods, we can effectively obtain parametric equations for the curved boundary of the lunar model, enabling the implementation of the numerical scheme discussed earlier. This ensures accurate simulation and analysis of the lunar environment for the intended application. Explicit form of Parametric Equations The following parametric relations can be immediately obtained for the lunar model describe above. (1) Outer circular arc Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 8s (2024) 732 https://internationalpubls.com The outer circular arc is the boundary curve over the fourth first and second quadrant and its boundary is described by the equation 2 2( 1/ 2) ( 1/ 2) 1/ 4x y− + − = . Hence the parametric equations.      ( ) 1 1 1 1 cos , sin 2 2 2 2 0, / 2 , / 2, , 3 / 2,2 x y            = + = +         (22) (2) Inner circular arc The inner circular arc is the boundary curve over the third quadrant and its boundary is described by the equation 2 2 1/ 4x y+ = and hence the parametric equations are 1 1 cos , sin ,( / 2 0) 2 2 x y   = =   (23) Quartic spline interpolants as parametric equations The parametric equations described above can be represented by two integral functions. F() and G() and the associated boundary conditions on F() and G() over each curved boundaries in the four quadrants (see fig 1). We have used F() and G() to generate the parametric equations for (x(), y()). We also impose the conditions ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ), , 1,2 p p p p x F y G p   = = = . We shall now explain the method of finding F(), G() and their relations with x(), y() the parametric equations. Once x() and y() are found in terms of the quartic spline integral representations as in [22] for F() and G(), we can implement the numerical scheme described in eqns (11) – (21). We now present the integral function representations for the curved arcs in each of the four quadrants. (i) Over the first quadrant [0, / 2]  , we define : ( ) ( ) ( ) 0 0 cos sin , 2 2 (0) 0, ( / 2) 1/ 2, 0, ( / 2) t t F d G d F F G G            = =        = = − =   so that we have ( ) ( ) 1x F = + and ( ) ( ) 1/ 2y G = + The quartic splines for (x(), y()) over the first quadrant,  0, / 2  can be then computed by the MATLAB function : quartic_spline (M, a, b, N) in which M refers to function type, a, b refer to end points a = 0, b = /2 for the present case and (N+1) refers to the number of subdivisions made in [a, b]. Similar reasoning is applicable to the curved arcs of the remaining quadrants viz second third and fourth. (ii) Over the second quadrant, [ / 2, ]   ,we define : Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 8s (2024) 733 https://internationalpubls.com / 2 / 2 cos sin ( ) , ( ) 2 2 ( / 2) 0, ( ) 1/ 2, ( / 2) 0, ( ) 1/ 2 t t F d G d F F G G               = =        = = − = = −   so that ( ) ( ), ( ) ( ) 1/ 2x F y G   = = + . (iii) Over the third quadrant, [ / 2,0]  and we define : /2 /2 cos sin ( ) , ( ) 2 2 ( / 2) 0, (0) 1/ 2, ( / 2) 0, (0) 1/ 2 t t F d G d F F G G             = =        = = = = −   so that we have ( ) ( ), ( ) ( ) 1/ 2x F y G   = = + Here, we may note that it is required to generate a circular arc in first quadrant of the circle 2 2 1/ 4x y+ = . (iv) Over the fourth quadrant, [3 / 2,2 ]   and we define : 3 / 2 3 / 2 cos 1 sin ( ) , ( ) 2 2 (3 / 2) 0, (2 ) 0 and (3 / 2) 0, (2 ) 1/ 2 t t F d G d F F G G           +    = =        = = = =   so that ( ) ( ) 1/ 2, ( ) ( )x F y G   = + = It is important to note here that as shown in Fig1, the outer circular arcs covers the first, second and fourth quadrants of the circle. 2 2( 1/ 2) ( 1/ 2) 1/ 4x y− + − = Examples We Consider the Following Integrals: 𝐼𝐼𝑖 = ∬ 𝑓𝑖(𝑥, 𝑦) ⬚ 𝜋𝑥𝑦 𝑑𝑥𝑑𝑦, (𝑖 = 1, 2, 3, 4, 5, 6, 7) (24) The values if the integrals 𝐼𝐼𝑖 (𝑖 = 1, 2, 3, 4, 5, 6, 7) with a relative error around 10−14, are 𝐼𝐼1 = 𝐼Ω0 (𝑓1) = 0.20307626985342 𝐼𝐼2 = 𝐼Ω0 (𝑓2) = 0.20646770293563 𝐼𝐼3 = 𝐼Ω0 (𝑓3) = 638.55743274702 𝐼𝐼4 = 𝐼Ω0 (𝑓4) = 0.57263720432530 𝐼𝐼5 = 𝐼Ω0 (𝑓5) = 0.03137185199242 𝐼𝐼6 = 𝐼Ω0 (𝑓6) = 0.0062895812195655 𝐼𝐼7 = 𝐼Ω0 (𝑓7) = 0.6426990816987241 Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 8s (2024) 734 https://internationalpubls.com where xy  is the curved domain described in fig 1. The integrands ( , ) i f x y are the following [18].   2 2 2 1 2 2 2 2 3 1 3 1 1 ( , ) exp (9 2) (9 2) exp (9 1) (9 1) 4 4 4 49 10 1 1 1 exp (9 7) (9 3) exp ((9 4) (9 7) 2 4 5 f x y x y x y x y y y     = − − + − + − + − +          + − − + − − − − + −    (25) 2 2 2 ( , ) ( 0.5) ( 0.5)f x y x y= − + − (26) 19 3 ( , ) ( )f x y x y= + (27)  2 2 4 ( , ) exp (( 0.5) ( 0.5) )f x y x y= − − + − (28)  2 2 5 ( , ) exp 100(( 0.5) ( 0.5) )f x y x y= − − + − (29) 6 ( , ) cos(20( ))f x y x y= + (30) 7 ( , ) 1f x y = (31) We have implemented the numerical scheme outlined in the above by leveraging the explicit parametric equations to compute the integrals specified by Equation (24). To facilitate this computation, we've developed two MATLAB programs. These programs utilize the precise representation of the boundary curve of the lunar model provided by the parametric equations elucidated in Equations (24). Notably, even with as few as four points strategically placed along the curved boundary of the lunar model, employing Gauss-Legendre quadrature rules of order (32, 36) or (32, 40) yields near-exact results. A glimpse of the output generated is showcased in Table-1. Moving forward, we address scenarios where the boundary curve of the lunar model is characterized by a set of discrete coordinate points. In such instances, we employ this discrete information to derive parametric equations in the form of spline interpolating polynomials. Specifically, we explore the application of quartic splines, a topic thoroughly discussed in this chapter. By transitioning to spline interpolating polynomials, we ensure a continuous and smooth representation of the boundary curve, even when described by discrete points. This strategy enables us to maintain the accuracy and reliability of our numerical computations, crucial for robust analysis and simulation of the lunar environment. Numbe r of Points( np) I1 I2 I3 I4 I5 I6 I7 1 637.86297 06 0.2064222 98 0.5725849 29 0.0313720 08 0.2105005 52 0.6426250 25 0.0063188 13 2 638.39271 0.2064559 0.5726234 0.0313718 0.2105027 0.6426797 0.0062962 Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 8s (2024) 735 https://internationalpubls.com 55 95 33 82 43 04 11 3 638.50166 88 0.2064635 56 0.5726322 65 0.0313718 62 0.2105033 86 0.6426921 6 0.0062915 69 4 638.53381 84 0.2064658 95 0.5726350 33 0.0313718 56 0.2105036 14 0.6426960 48 0.0062903 88 5 638.54582 03 0.2064667 96 0.5726361 09 0.0313718 54 0.2105037 09 0.6426975 54 0.0062899 65 6 638.55108 42 0.2064672 0.5726365 94 0.0313718 53 0.2105037 54 0.6426982 32 0.0062897 86 7 638.55367 78 0.2064674 02 0.5726368 38 0.0313718 53 0.2105037 78 0.6426985 72 0.0062897 01 8 638.55507 29 0.2064675 12 0.5726369 72 0.0313718 52 0.2105037 91 0.6426987 58 0.0062896 55 9 638.55587 65 0.2064675 76 0.5726370 49 0.0313718 52 0.2105037 98 0.6426988 67 0.0062896 3 10 638.55636 55 0.2064676 16 0.5726370 97 0.0313718 52 0.2105038 03 0.6426989 33 0.0062896 14 100 638.55743 26 0.2064677 03 0.5726372 04 0.0313718 52 0.2105038 15 0.6426990 82 0.0062895 81 200 638.55743 27 0.2064677 03 0.5726372 04 0.0313718 52 0.2105038 15 0.6426990 82 0.0062895 81 300 638.55743 27 0.2064677 03 0.5726372 04 0.0313718 52 0.2105038 15 0.6426990 82 0.0062895 81 400 638.55743 27 0.2064677 03 0.5726372 04 0.0313718 52 0.2105038 15 0.6426990 82 0.0062895 81 500 638.55743 27 0.2064677 03 0.5726372 04 0.0313718 52 0.2105038 15 0.6426990 82 0.0062895 81 600 638.55743 27 0.2064677 03 0.5726372 04 0.0313718 52 0.2105038 15 0.6426990 82 0.0062895 81 700 638.55743 27 0.2064677 03 0.5726372 04 0.0313718 52 0.2105038 15 0.6426990 82 0.0062895 81 800 638.55743 27 0.2064677 03 0.5726372 04 0.0313718 52 0.2105038 15 0.6426990 82 0.0062895 81 900 638.55743 27 0.2064677 03 0.5726372 04 0.0313718 52 0.2105038 15 0.6426990 82 0.0062895 81 1000 638.55743 27 0.2064677 03 0.5726372 04 0.0313718 52 0.2105038 15 0.6426990 82 0.0062895 81 Table - 1 Spline Integration for 𝐼𝐼(𝑖) = ∬ 𝑓𝑖(𝑥, 𝑦) ⬚ 𝐷 𝑑𝑥𝑑𝑦, 𝑖 = 1(1)7, D: Domain of Lunar model I1 integral seems to converge to approximately 638.5574327 as the number of points Np increases. This suggests that I1 has reached a stable value, indicating convergence of the numerical method used to compute it. I2: Similarly, converges to a stable value of approximately 0.206467703 as Np Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 8s (2024) 736 https://internationalpubls.com increases. The convergence of I2 indicates that the numerical method accurately computes this integral. I3, I4, I5, I6, I7: These integrals also converge to stable values as Np increases, with each approaching approximately 0.572637204, 0.031371852, 0.210503815, 0.642699082, 0.572637204, 0.031371852, 0.210503815, 0.642699082, and 0.006289581 respectively. The convergence of all these integrals suggests the effectiveness and accuracy of the numerical method across different aspects of the problem. Overall, the trend observed in the table indicates that increasing the number of points Np in the numerical method leads to improved accuracy and convergence of the computed integrals. This convergence behavior is consistent with expectations from numerical methods, where increasing the resolution typically leads to more accurate results. Application example 2: general ellipse: Consider an ellipse in (x, y) plane with semi major axis of length a and semi minor axis of length b. The ellipse is then rotated anticlockwise by an angle of 𝛉 degrees and translated to coincide the centre of the ellipse to (x0, y0). The equation of the ellipse in the (x, y) plane is 2 2 2 2 1 x y a b + = . Let the new coordinate system with origin at (x0, y0) be denoted by ( , )x y  plane. The ellipse in ( , )x y  plane is then rotated by an angle 𝛉 in anticlockwise direction. Let us denote the points on the rotated ellipse with respect to ( , )x y  coordinate system as * *( , )x y . (i). The points on the ellipse 2 2 2 2 1 x y a b   + = are cos , sin ,0 2x a t y b t t  = =   (ii). The points * *( , )x y corresponding to the point ( , )x y  on the rotated ellipse are * *cos sin , sin cosx x y y x y      = − = + (iii). Points (x, y) corresponding to * *( , )x y must satisfy * * * * 0 0 * * 0 0 , . , , , x x y y i e x x x x y y y y x x x y y y  = =  = − = = − =  = + = + We shall now determine the points on the rotated ellipse. Let the points on the ellipse with respect to ( , )x y  should be {(a, 0), (0, b), (- a, 0), (0, - b)}. These points with respect to (x, y) system are { 0 0 0 0 ( cos , sin ),( sin , cos ),x a y b x b y b   + + − + 0 0 ( cos , sin ),x a y a − − 0 0 ( sin , cos )x b y b + − }. The remaining points on the ellipse can be determined in a similar manner. The geometric moments over the general ellipse described above is given by, p q pq A m x y dxdy=  (32) Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 8s (2024) 737 https://internationalpubls.com where A is the area of the ellipse. The compound moment values and the exact moment values are tabulated. We may note that the  - method, the rectangular grid method, trapezoidal integration method and contour integration method proposed also tabulated for the purpose of integration. Number of points on general ellipse np Moments Computed Values of Moments Theoretical Values 1.0e+006 * 498 m00 = 1.884955592153876e+003 m10 = 1.884955592153876e+004 m01 = 5.654866776461627e+004 m20 = 4.948008429403924e+005 m02 = 2.002765316663493e+006 m11 = 6.832964021557800e+005 fn=1 1.884955592022257e+003 fn=x 1.884955592161867e+004 fn=y 5.654866776206377e+004 fn=x^2 4.948008429123767e+005 fn=y^2 2.002765316586024e+006 fn=x*y 6.832964021556849e+005 0.00188495559202 0.01884955592162 0.05654866776206 0.49480084291238 2.00276531658602 0.68329640215568 998 m00 = 1.884955592153876e+003 m10 = 1.884955592153876e+004 m01 = 5.654866776461627e+004 m20 = 4.948008429403924e+005 m02 = 2.002765316663493e+006 m11 = 6.832964021557800e+005 fn=1 1.884955592145685e+003 fn=x 1.884955592154394e+004 fn=y 5.654866776445747e+004 fn=x^2 4.948008429386528e+005 fn=y^2 2.002765316658672e+006 fn=x*y 6.832964021557791e+005 1.884955592145685e+003 1.884955592154394e+004 5.654866776445747e+004 4.948008429386528e+005 2.002765316658672e+006 6.832964021557791e+005 1498 m00 = 1.884955592153876e+003 m10 = fn=1 1.884955592152255e+003 fn=x 1.884955592152255e+003 1.884955592153955e+004 5.654866776458480e+004 Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 8s (2024) 738 https://internationalpubls.com 1.884955592153876e+004 m01 = 5.654866776461627e+004 m20 = 4.948008429403924e+005 m02 = 2.002765316663493e+006 m11 = 6.832964021557800e+005 1.884955592153955e+004 fn=y 5.654866776458480e+004 fn=x^2 4.948008429400452e+005 fn=y^2 2.002765316662538e+006 fn=x*y 6.832964021557742e+005 4.948008429400452e+005 2.002765316662538e+006 6.832964021557742e+005 1998 m00 = 1.884955592153876e+003 m10 = 1.884955592153876e+004 m01 = 5.654866776461627e+004 m20 = 4.948008429403924e+005 m02 = 2.002765316663493e+006 m11 = 6.832964021557800e+005 fn=1 1.884955592153362e+003 fn=x 1.884955592153892e+004 fn=y 5.654866776460630e+004 fn=x^2 4.948008429402803e+005 fn=y^2 2.002765316663193e+006 fn=x*y 6.832964021557753e+005 1.884955592153362e+003 1.884955592153892e+004 5.654866776460630e+004 4.948008429402803e+005 2.002765316663193e+006 6.832964021557753e+005 Table - 2 Computed Values of Moments Conclusions we have presented quartic spline, which interpolates the first derivatives of a given function at the knots and the second derivatives between them. A simple approximation to function values is determined by choosing the second derivate to coincide with the initial knots. A numerical integration scheme Gauss Legendre quadrature rules and the explicit parametric relations along the curved boundary of the two dimensional Cartesian space is considered. The above numerical integration scheme is illustrated by computing several difficult integrals over the curved domain. The integral values with its absolute error of integrals are tabulated. Integration for 𝐼𝐼(𝑖) = ∬ 𝑓𝑖(𝑥, 𝑦) 𝐷 𝑑𝑥𝑑𝑦, 𝑖 = 1(1)4, D: Domain of as shown in figure 1 Acknowledgment: Both the authors acknowledge Late Dr. H T Rathod for his valuable suggestions, the support and guidance. Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 8s (2024) 739 https://internationalpubls.com References: [1] Anwar, M.N and El-Tarazi, M.N: Direct cubic spline with application to quadrature, communications in applied numerical methods, 5, pp 237 – 246 (1989). [2] Apostol T. M., Calculus, Vol. II, Second Edition Blaisdell (1969). [3] Elden, L., LeWittmeyer-Koch, H. B. Nielson Introduction to Numerical computation, Analysis and MATLABR illustrations, Overseas Press India Pvt, Ltd (2006). [4] Gustafasson, F. and Bergman, N., MATLAB for engineers explained, Springer International Edition (2003). [5] H. Behforooz, Approximation by integro cubic splines, Appl. Math. Comput. 175 (2006) 8–15. [6] H. Behforooz, Interpolation by integro quintic splines, Appl. Math. Comput. 216 (2010) 364–367 [7] H. T. Rathod and H. S. Govinda Rao, Integration of polynomials over n-dimensional linear polyhedra, Computers and Structures, Vol.65, No.6, pp.829-847 (1997). [8] H. T. Rathod and S. V. Hiremath, Boundary Integration of polynomials over an arbitrary linear hexahedron in Euclidean three dimensional space, Computer Methods in Applied Mechanics and Engineering, Vol.161, pp.155- 193 (1998). [9] H. T. Rathod, Md.Shafiqul.Islam, Bharath Rathod , K. Sugantha Devi , Finite element solution of Poisson Equation over Polygonal Domains using a novel auto mesh generation technique and an explicit integration scheme for linear convex quadrilaterals of cubic order Serendipity and Lagrange families, International Journal Of Engineering And Computer Science ISSN:2319-7242,Volume 7 Issue 1 January 2018, Page No. 23329-23482 [10] H.T. Rathod, A. S. Hariprasad, K.V.Vijayakumar, Numerical Integration over Curved Domains Using Convex Quadrangulations And Gauss Legendre Quadrature Rules International Journal Of Engineering And Computer Science, 2(11), 2013, 3290-3332 [11] Hanselman.D and B. Littlefield, Mastering MATLABR, Pearson Education Inc (2007). [12] John Mathews, Kurtis Fink, Numerical Methods Using Matlab, International Edition, 4th Edition, Pearson, Jan 2004. [13] Lo. S. H, Generating quadrilateral elements on plane and over curved surfaces, Comput. Stuct. 31(3) 421- 426(1989) [14] MATLAB Programming Version 7, The Maths Works INC, 2004. [15] MIcula C., Sanda Micula, Hand Book of Splines, Kluwer Academic Publishers (1999). [16] Moler, C. B., Numerical Computing with MATLAB, Society for Industrial and Applied Mathematics, Phildelphia (2004). [17] Otto S. R. and J. P. Denier., An Introduction to programming and numerical methods in MATLAB for beginners and experienced users, Cambridge University Press (2000). [18] Sallam S., W. Ameen Numerical solution of general nth-order differential equations via splines Applied Numerical Mathematics, 6(3), 1990, 225-238. [19] Shikin E. V., Alexander I. Plis, Handbook on splines for the user, CRC Press, 1995. [20] Sommariva A. and M. Vianello, Gauss-Green Cubature over spline curvilinear polygons, Applied Mathematics and Computation 183(2), 1098 - 1107 (2006) [21] Sommariva A. and M. Vianello, Gauss-Green Cubature over spline curvilinear polygons, Applied Mathematics and Computation Vol. 183 (2) pp 1098 - 1107 (2006) [22] Stroud A. H., Approximate calculation of multiple integrals, prentice Hall series in automated computation, Prentice Hall, Inc Englewood Cliffs N. J., (1971). [23] T. Zhanlav, R. Mijiddorj, The local integro cubic splines and their approximation properties, Appl. Math. Comput. 216 (2010) 2215–2219. [24] Zienkiewicz. O. C, Taylor. R. L and Zhu. J. Z, The finite element method, its basis and fundamentals, 6th Edn, Elsevier (2007).