Averaged cubature schemes on the real positive semiaxis

Stratified cubature rules are proposed to approximate double integrals defined on the real positive semiaxis. In particular, anti-Gauss cubature formulae are introduced and averaged cubature schemes are developed. Some of their appropriate modifications are also studied. Several numerical experiments are given to testify the performance of all the formulae.


Introduction
This paper deals with the numerical approximation of double integrals in which the integrand function may increase as x, y → ∞ and/or x,y → 0. Then, we consider where f is a known function defined on the domain [0, ∞) × [0, ∞) and w(x, y) = x y e −(x+y) , , > −1 is a weight function which includes eventual algebraic singularities. Let us note that w turns out to be the product of two Laguerre weights w (x) = x e −x and w (y) = x e −y .
The numerical treatment of (1) can be approached in two ways [1,2]. The first procedure consists of approximating each integral by using well-known quadrature rules. This is an "indirect" technique that takes the advantage of the fact that the univariate rules have been more treated and explored than the multivariate ones. The second approach consists of constructing cubature schemes from scratch. In the bivariate case, an example is given in [3] where the nodes are zeros of suitable bivariate orthogonal polynomials (see also [4]).
In this paper, we focus on the "indirect" procedure. The optimal cubature formula is the following Gauss-Laguerre cubature rule [5] based on m × n nodes where {x k } m k=1 and {y j } n j=1 are the zeros of the Laguerre polynomials which are orthogonal with respect to the weights w α and w β , respectively, { k } m k=1 and { j } n j=1 are the corresponding coefficients of the Gauss-Laguerre formula, and the term R m,n (f) denotes the related cubature error of G m,n . The rule is obtained as the tensor product of two univariate Gauss-Laguerre quadrature formulae. Hence, it can be easily constructed [5] and R m,n (f) = 0 if f ∈ ℙ 2m−1,2n−1 , where ℙ 2m−1,2n−1 denotes the set of all algebraic bivariate polynomials of degree at most 2m − 1 in the variable x and 2n − 1 in the variable y.
One can associate to the cubature rule (2), the truncated version introduced in [6]: Here, the upper limit of summations κ and ι are the indeces determined by with 1 , 2 ∈ (0,1) fixed, and R m,n (f ) identifies the error. As shown in [6] (see also [7,8]), the cubature error of (3) has the same magnitude as the error of (2) and is convenient especially when the function f in (1) is bounded or does not increase too fast, so that the last terms of the sums in (2) can be neglected. Of course, in order to construct formula (3) all the m + n weights and nodes of (2) have to be computed but the number of function evaluations is reduced. Moreover, an application of the truncated formula (3) to integral equations implies a considerable computational saving (see [6][7][8][9]).
The aim of this paper is to give estimates of the magnitude of the errors R m,n (f ) and R m,n (f ) . We remark that in [  x = min{x j |x j ≥ 4m 1 }, y = min{y j |y j ≥ 4n 2 }, nature of the integrand. Specifically, they state the order of convergence depending on the smoothness properties of the function f, providing a lower bound which include an unknown constant independent of m, n, and f. In this paper, we provide numerical estimates of the error for each fixed m and n so that one can determine the number of points m and n needed to approximate the integral with a prescribed accuracy. The estimate does not depend on unknown constants and it is not asymptotic. So, it does not provide a decay law depending on the smoothness of f. Neverthless, it allows us to develop new cubature rules, which turns out to have several advantages in terms of accuracy and computational cost. In the univariate case, an interesting approach to estimate the error of the m-point Gauss formula G m consists of using another quadrature rule Q l with l > m nodes and degree of exactness higher than 2m − 1. A popular choice of Q l is the Gauss-Kronrod rule [10] which has 2m + 1 nodes, including the m nodes of G m , and has degree of exactness at least 3m + 1. However, in some cases, among which the Laguerre one, the formula fails to exist: the nodes are not real and distinct and the weights are not positive [10][11][12][13]. Another possible choice is given by "stratified" quadrature formulae (see, e.g., [14][15][16][17][18]). They are linear combinations of two rules, usually the m-point Gauss rule G m and a new rule U ℓ which involves a number of points ℓ greater than m and has degree of exactness greater than 2m − 1. A first example of a "stratified" rule was given by Laurie [17] in 1996 who proposed the so-called averaged Gaussian quadrature formula G L 2m+1 : Basically, it is the average between the m-point Gauss rule G m and the m + 1-point anti-Gauss formula G A m+1 which is such that where I denotes the corresponding univariate integral of the form (1). Formula G L 2m+1 has degree of exactness at least 2m + 1, involves 2m + 1 real and distinct zeros and has positive weights. All these properties hold true also in the Laguerre case for each α > − 1, α being the parameter of the Laguerre weight w α appearing in the integral I [17]. In 2002, Ehrich [16] proposed the optimal stratified extension for the Gauss-Laguerre and the Gauss-Hermite formulae. Specifically, he proved that, if we handle with the Laguerre weight w α , the formula G L * 2m+1 of the highest possible degree of exactness 2m + 2 is unique and is obtained for However, this formula has positive nodes only for α ≥ 1, whereas it has one negative node for − 1 < α < 1. In 2007, Spalević [18] developed the generalized averaged formula G S 2m+1 which is the unique stratified formula with the highest possible degree of exactness and, in the Laguerre and the Hermite cases, coincides with G L * 2m+1 . In 2016, Djukić, Reichel, and Spalević, in an attempt to ensure that all nodes be internal, proposed in [15] the truncated generalized averaged Gauss quadrature rule obtained by removing the last r < m rows and columns from the Jacobi matrix of G S 2m+1 . In particular, taking r = m − 1 one gets the truncated rule G T m+2 that in the Laguerre case has positive zeros for each α > − 1. Let us mention that here the term "truncated" does not have the same meaning as in formula (3). In fact, in (3) the truncation refers to the sequence of nodes whereas in G T m+2 the truncation refers to the Jacobi matrix. Consequently, in this last case, all the computed nodes are only the necessary ones. To avoid confusion, from now on, we will name G T m+2 "reduced" generalized averaged Gauss rule.
In the multidimensional case, in [19] the authors consider integrals defined on bounded intervals [a,b] and estimate the error of the corresponding Gauss cubature rule by using the generalized averaged formula G S 2m+1 .
In order to estimate the error R m,n , in this paper, we develop a "stratified" cubature formula which is obtained as an average of the cubature formula (2) and the anti-Gauss cubature formula. According to our knowledge, this last formula has never been proposed to approximate double integrals and, even if we get it as tensor product of the univariate anti-Gauss rule, some new useful properties are proved. Moreover, its stability and convergence is studied in suitable weighted spaces equipped with the uniform norm. In addition, we write the "reduced" generalized averaged Gauss cubature rule which is the bivariate version of the formula proposed in [15].
Then, we give two truncated rules aiming at estimating the error R m,n (f ) . Here, what we develop are the truncated anti-Gauss cubature rule and the "reduced" truncated generalized averaged Gauss cubature rule. Basically, we consider the two developed cubature schemes, truncate the sequences of nodes and prove that these rules provide an error which is of the same magnitude of the original cubature formula by only using a fraction of nodes.
The paper is structured in five main sections. Section 2 focuses on the estimation of the error R m,n , Section 3 pays attention to R m,n , and Section 4 gathers several numerical examples to show and compare the performance of the four considered cubature formulae. Section 5 collects the proofs of the theoretical results and, finally, Section 6 draws some conclusions and briefly sketches possible research developments.

On the error of the Gauss-Laguerre cubature formula
The aim of this section is to estimate the error R m,n by introducing "stratified" cubature formulae or averaged schemes of the form where η is a given real parameter, G m,n (f) is as in (2), and U m+ 1,n+ 1 is a new cubature formula.

The averaged Gauss-Laguerre cubature rule
A possible approach for constructing the new rule U m+ 1,n+ 1 is to use the well-known 1D anti-Gauss formula [17] to approximate each integral of (1). In this way, we get a cubature scheme which we will refer to as anti-Gauss cubature formula G A where a j = 2j + + 1 , j ≥ 0, b 0 = Γ(1 + ) , b j = j(j + ) , j ≥ 1. Similarly for p j . Then, the nodes {x k } m+1 k=1 and {ỹ j } n+1 j=1 are the zeros of the polynomials p m+1 and p n+1 , respectively, defined as and, consequently, are the eigenvalues of the following two matrices: The coefficients ̃k and ̃j are also related to the matrices J m+1 and J n+1 . Indeed, they are determined by the first component of the associated normalized real eigenvectors.
Next Proposition contains the main properties of G A m+1,n+1 (f ).

Proposition 1
The following properties hold true. Similarly, for each fixed j,

Let us mention that the first property implies
This suggests that a "stratified" cubature rule of the form (6) can be obtained with = 1∕2 and U m+1,n+1 = G A m+1,n+1 ∕2 , that is From now on, we will refer to the above formula as the averaged Gauss-Laguerre cubature rule.
Let us also note that G L 2m+1,2n+1 requires the computation of 2(m + n + 1) nodes and weights and if we use the algorithm devised by Golub and Welsch in [20] their computation has a cost equal to 2c(m 2 + n 2 ) + 2O(m) + 2O(n) , that is one half of the cost of the cubature Gauss rule G 2m,2n which is 4c(m 2 + n 2 ) + 2O(m) + 2O(n).
Moreover, the first property of Proposition 1 suggests that the cubature error for . This allows us to reach a certain accuracy with a number of evaluation of the integrand function which is smaller than that required by G 2m,2n . Furthermore, we can use G L 2m+1,2n+1 (f ) to provide an accurate approximation of the cubature error (2), that is Let us mention that R [1] m,n (f ) provides an accurate numerical estimates for R m,n (f); for each fixed m and n, see Table 8 for a numerical evidence.
Finally, we underline that the inequalities proved in Property 4 are Possé-Chebyshev-Markov-Stieltjes type inequalities [21, p. 33] and they allow us to state the following corollary.
Corollary 1 For a fixed quadrature node x k and ỹ k and for each differentiable function g such that g (j) > 0 for j = 0,...,2m + 1 one has The Corollary is essential to prove the stability of the formula G A m+1,n+1 in suitable weighted spaces. In spaces of continuous functions C ≡ C([0, ∞]) , the stability is trivial.
Let us now investigate on the stability and convergence properties of the anti-Gauss cubature formula in suitable weighted spaces C u . This setting is useful when we want to consider functions f(x,y) that may tend to infinity with algebraic growth as x,y → 0 + and with an exponential growth as x, y → ∞ . Fix the weight functions and, setting we define the weighted space C u as the set of all continuous functions on (0, ∞) × (0, ∞) such that We endow the space C u with the norm For smoother functions, we introduce the Sobolev-type space of index 1 ≤ r ∈ ℕ where 1 (x) = √ x , 2 (y) = √ y and f y means that the function f is considered as a function of the only variable x. Similarly for f x .
In C u , we define the error of best polynomial approximation by means of bivariate polynomials in ℙ m,n as (see, e.g., [ where C is a positive constant independent of m, n and f.
Theorem 1 For any f ∈ C u if α > γ 1 − 1 and β > γ 2 − 1, the anti-Gauss cubature formula is stable and where C is a positive constant independent of m, n, and f.
By the previous theorem, we can deduce the following estimate for the error of the averaged rule (13) Corollary 2 For any f ∈ C u , if α > γ 1 − 1 and β > γ 2 − 1, then

The reduced generalized averaged Gauss cubature rule
In this subsection, we want to extend to the bivariate case the one-dimensional quadrature rule G S 2m+1 having the highest possible degree of exactness 2m + 2. Such rule has several formulations [16,18,22] and has been applied to matrix functions approximation with applications to complex networks analysis [23,24].
For our aims, we will consider the formulation presented in [18] and we will refer to it as the generalized averaged Gauss cubature rule G S 2m+1,2n+1 . It is obtained as tensor product of G S 2m+1 and G S 2n+1 and reads as Formula (18) retains nice properties that can be automatically deducted by the corresponding univariate rule [16,18,22]. We summarize them in the following proposition.

Proposition 2
The generalized averaged Gauss cubature rule (18) has the following properties.
The nodes are all real. Specifically, when they are ordered in increasing order, they are such that where x k are the Gauss nodes and x * k are the zeros of the following polynomial i.e., the eigenvalues of the following matrix The same holds for the nodes ȳ j with β and n in place of α and m, respectively. • The weights are all positive. In details, assuming that they are ordered in increasing order, one has where * k can be computed from the first components of the associated normalized real eigenvectors of the matrix J * m+1 .
According to the last property, the introduced cubature rule (18) has nodes outside the domain (0, ∞) . Therefore, in order to overcome this problem, we follow what has been proposed in [15] for the univariate case, that is, we consider the Jacobi matrix J 2m+1 given in (18) p n+1 (y) = (y − a n−1 )p n+1 (y) − b n+1 p n (y), The above Proposition suggests that we can use G T m+2,n+2 (f ) to estimate the cubature error R m,n (f) by computing the following difference: 3 On the error of the truncated Gauss-Laguerre cubature rule In this section, we develop truncated versions of both the rules previously introduced, preserving the same accuracy of the original cubature error. We aim at reducing the number of function evaluations, producing at the same time an estimate for the error of the truncated Gauss-Laguerre rule (3).
Generally, truncated rules are very useful in the numerical solution of integral equations [6,8,9]. In fact, collocation methods lead to square linear systems whose order depends on the number of cubature points. If we approximate the integral operator by using a cubature scheme with m × n nodes, then the system will be of order mn. Moreover, if the kernel and the right-hand side of the equation are not smooth, to obtain a good accuracy, we need to take m,n ≫ 1, so that the system becomes very large. Truncated rules use only a fraction of the nodes in the cubature scheme, reaching the same precision by solving a smaller linear system, with a significant saving in computing time.

The truncated averaged Gauss-Laguerre cubature rule
Let us consider the rule (7). Among all the computed nodes {x k } m+1 k=1 and {ỹ j } n+1 j=1 let us select the node x and the node ỹ such that where 1 , 2 ∈ (0, 1) are two fixed parameters.
Then, we define the truncated anti-Gauss Laguerre cubature formula as Denoting next proposition shows that, for a certain class of polynomials, formula (22) gives an error opposite in sign of the error of formula (3).

, it follows
According to Proposition 4, one has This identity suggest us to define the truncated averaged Gauss-Laguerre cubature rule as Such formulae can be used to estimate the cubature error of (3) Let us note that R [1] m,n (f ) provides an accurate numerical estimate for R m,n (f ) ; for each fixed m and n, see Table 9 for a numerical evidence, which also confirms the lower bound given in [6].

The truncated reduced generalized averaged Gauss cubature rule
Similarly to what has been done for the anti-Gauss rule, let us consider now the formula (19) and introduce the nodes where 1 , 2 ∈ (0, 1) are two fixed parameters. Then, we define truncated reduced generalized averaged Gauss cubature rule as Next proposition shows that the error given by the above formula is of the same order as the one given by the original formula (19). Then, in order to use fewer points, we can use this formula to estimate the cubature error R m,n (f ) , i.e., [1] m,n (f ).
Proposition 5 Let f ∈ C u where u is given by (15) with γ 1 < α + 1 and γ 2 < β + 1. Then, where C and c are positive constants independent of m and n.

Numerical examples
The aim of this section is to show and compare the performance of the four cubature rules applied to several integrals of the form (1).
In each example, for increasing value of m and n, we compute the relative errors of the Gauss and anti-Gauss cubature rule m,n , A m+1,n+1 , together with the corresponding relative error L 2m+1,2n+1 furnished by the averaged formula (13). Then, we show the relative errors of the associated truncated rules ̊m ,n , ̊A m+1,n+1 and ̊L m+1,n+1 to highlight the advantage of such formulae in terms of accuracy and evaluation of functions. In addition, we test the performance of the reduced generalized averaged Gauss rule and its truncated version, by showing the relative errors T m+2,m+2 and ̊T m+2,n+2 .
In each example, we also compute the Estimated Order of Convergence (EOC), i.e., where here err m,n is the absolute error of the cubature rule we are dealing with.
All the computations were performed on an Intel Xeon E-2244G system with 16Gb RAM, running MATLAB 9.10. The software developed is only prototypal, but it is available from the authors upon request.
Example 1 Let us test the convergence of the considered cubature rules to the following integral with f (x, y) = sin(x + y)x 3 y and w(x,y) = e −x−y . The integrand function is very smooth and then we expect a fast convergence. The numerical tests are performed with high-precision arithmetic. Specifically, to make the round-off errors introduced during the computations negligible, the computations are achieved with 110 significant decimal digits. By choosing m = n, Table 1 contains the relative errors of the [2] m,n (f ).
EOC m,n = log err m,n err 2m,2n log 2 , classical Gauss-Laguerre formula G m,m (f), the anti-Gauss rule G A m,m (f ) and the corresponding averaged rule G L 2m+1,2m+1 (f ) , together with the corresponding EOC m,m . By comparing the second and fourth column, we can see that the anti-Gauss formula furnishes an error opposite in sign and with the same magnitude of the Gauss rule. Consequently, if we approximate the integral by using the averaged rule, then we get a better approximation as highlighted by the sixth column.
In Table 2, we collect the results we get by approximating the integral with the truncated version of the Gauss, anti-Gauss rule, and the corresponding averaged formula, by fixing 1 = 2 = 0.4. The truncated anti-Gauss rule still retains the same properties of its complete formulation, that is gives an error opposite in sign and with the same magnitude of the truncated Gauss-rule. This implies that if we use the truncated averaged rule we get the same accuracy of the complete rule but with fewer points. For instance, we get an error of the order 10 − 9 with 27 × 2 = 54 nodes (instead of 33 × 2 = 66 nodes) and 365 evaluations of functions (instead of 545 evaluations). Table 3 contains the relative errors we get with the reduced generalized averaged formula G T m+2,m+2 and its truncated version G T m+2,m+2 . They are smaller than the errors given by the single rules G m,m and G A m+1,m+1 but the averaged rule G L 2m+1,2m+1 is more precise. By the last column, we can see that, also in this case, the truncated form gives the same accuracy of the complete rule with fewer points.   In Table 4, we want to compare the performance of the reduced generalized averaged formula G T m+2,m+2 with the complete generalized averaged formula from which it derives G S 2m+1,2m+1 to show that the first formula reaches the same accuracy of the second one even if the number of points is reduced. To this end, to make internal the complete formula, we set the parameters of the weight function as α = 2 and β = 1 and f (x, y) = sin (x + y).

Example 2
The aim of this test is to show that in the two-dimensional case (and even more in the multidimensional one) the averaged rules are more competitive with respect to the simple Gauss formula. Let us consider the following integral which is of the form (1) with The function f is smooth and, for each fixed y, is bounded and does not increase too fast, so we need to use a large value of m and n to get a good accuracy [7]. The exact value of the integral is not known. Since the use of a high precision considerably increases the computation times when m is large, we have computed the exact value numerically with Wolfram Mathematica. We perform the computation with 32 significant decimal digits.
By inspecting Table 5 we can deduce, for instance, that to get an error of the order 10 − 8 we have two options: using the averaged rule G L 2m+1,2n+1 with m = n = 64 which w(x, y) = e −(x+y) f (x, y) = 1 1 + y + 2x e x∕4 (x − 2) 2 + 1 .  requires the computation of 2(2m + 1) = 258 nodes and 64 2 + 65 2 = 8321 evaluation of functions, or using the simple Gauss-Laguerre formula G m,m with m = n = 128. However, in the latter case, we have to compute 256 nodes and 16384 evaluation of functions.
If we move to the truncated averaged rules, the number of function evaluations is further reduced. As showed in Table 6, to reach an error of order 10 − 8 , by using the truncated averaged rule G L 2m+1,2m+1 , we only need to compute the integrand function in 36 2 + 37 2 = 2665 points. On the contrary, the computational cost of the truncated Gauss G m,m rule is 71 2 = 5041 evaluations. Table 7 contains the results concerning the reduced generalized averaged rule from which we can see that G T m+2,m+2 is less accurate than the rule G L 2m+1,2m+1 .

Example 3 Let us consider an integral in which the integrand function is not smooth as in the previous examples, that is
where we can fix w α (x) = x − 1/10 e −x and w β (y) = y − 1/5 e −y . Here, we compute the tests with double arithmetic precision and since the exact value is not known, we consider as exact the approximated value obtained by the Gauss cubature formula with m = 1024. The integrand function f (x, y) = |y−1| 5∕2 25+x 3 +y 3 is smooth with respect to the variable x but belongs to the Sobolev space W 2 with respect to the variable y. Then,  Tables 8, 9, and 10. Once again, the computational saving of the truncated averaged rule is confirmed.
In this example, we have also add a column in each table to show that the proposed formulae can be used to estimate the relative error given by the classical Gauss cubature formulae G m,m and G m,m as emphasized in (14), (21), (23), and (24).  Table 9 Numerical results for Example 3 with 1 = 2 = 0.2  are the Gauss and anti-Gauss rule applied on the polynomial g. Using the fact that these two rules give the opposite errors on all polynomials of degree at most 2m + 1, we conclude that Property 2 It is inherited by the one-variable anti-Gauss rule [17,26]. Property 3 The positivity follows by the fact that the j th weight is the squared of the first component of the eigenvector associated to the j th eigenvalue. Let us now give the proof of the expression of the weights k , by following essentially the one given in [26] for bounded intervals. The identity related to the weights j can be proved in the same way. For each fixed k, let us consider the polynomial of degree 2m Since G m (f k ) = 0, the univariate averaged rule (4) gives the identity On the other hand, for each k, there exist a polynomial q m− 1,k of degree m − 1 such that where ∥⋅∥ denotes the usual norm in the Hilbert space L 2 (ℝ + ) . Then, by combining (25) with (26)  Consequently, by the assumption on α and β, we deduce the boundedness of the integrals and then the stability of the formula. About the convergence, denoted by P ∈ ℙ 2m−1,2n−1 we have Then, taking infimum on ℙ 2m−1,2n−1 , by assumptions we get the thesis.
Proof of Corollary 2 By definition, taking Theorem 1 into account and being |R m,n (f )| ≤ CE 2m−1,2n−1 (f ) u one get Proof of Proposition 3 Both assertions follow by the properties of the corresponding univariate quadrature rule G T m+2 [14].

Proof of Proposition 4
By the first assertion of Proposition 1 we can write Therefore, we have Let us now prove that S 1 and S 2 are negligible so that they cannot influence the sign of the first term at the right-hand side.
Let u(x,y) = u 1 (x)u 2 (y) with u i (x) = (1 + x) i x i e −x∕2 , where η i ≥ 0 for each i = 1,2 and γ 1 < α + 1, γ 2 < β + 1. Then, setting = min{ 1 , 2 } and M = min{m, n} , we have S 1 = sign(S 1 )|S 1 |, where [9, formula 4.2] with C and c two positive constants independent of M and p, and with c depending on . Hence, S 1 approaches zero and then whatever its sign, it does not change the sign of R m,n (f ) . A similar result can be proven for the sum S 2 so that we have the assertion.

Proof of Proposition 6
Let P M the polynomial of best approximation for f, we can write p(x k ,ỹ j ).
By virtue of the exactness of formula G T m+2,n+2 for constants, we have Moreover, by proceeding as done in Proposition 4 for the sums S 1 and S 2 , we get

Conclusions and research perspectives
In this paper, we have proposed four cubature rules to approximate double integrals defined on the interval [0, +∞) . The first one is the averaged Gauss cubature formula which is an average of the Gauss cubature formula and the anti-Gauss cubature rule, introduced here for the first time. Then, we have extended to the bivariate case the formula proposed in [15] by writing the reduced generalized averaged Gauss cubature rule. Finally, we have investigated on the corresponding truncated versions obtained by removing "large" nodes. Our numerical results confirm the theoretical properties of these schemes. In the case of the averaged rule, the formula does not only provide a good estimate for the remainder terms R m,m (f ) and R m,m (f ) but gives a good approximation of the integral with a reduced number of evaluation of functions with respect to the one given by the Gauss cubature formula G 2m,2n . To this advantage, the reduced computational cost for computing nodes and weights is added, compared to the related cost required for the construction of G 2m,2n . If we move on the truncated versions, the number of evaluation of function is streactly reduced.
All the above formulae are constructed as tensor product of the corresponding univariate formulae and then are based on the zeros of polynomials orthogonal with respect to a weight function of the type u(x) = x α e −x .
Regarding research perspectives, we want to construct these formulae by considering the more general weight u(x) = x e −x [27][28][29]. In this case, the coefficients of the recurrence relation are not known in a closed form and then must be evaluated numerically. The Mathematica package "Orthogonal Polynomial" [30,31] computes such coefficients (and then evaluates nodes and weights) for the Gauss and anti-Gauss rule but not for the reduced generalized averaged formula. Then, we want to develop a software package, written in Matlab, which includes the computation of the recurrence coefficients of all the cubature formulae investigated here.
Our package will also include solvers for the numerical solution of integral equations defined on the set (0, ∞) × (0, ∞) based on the cubature schemes presented in this paper [32]. Their properties will imply a considerable computational saving in the numerical solution of the linear system that will have to be solved.