Example 2.1: Series expansion at a singularity 

>

 

(2.1.1)

 

>

 

>

 

Form the integral on and apply numerical integration. 

>

 

(2.1.2)

 

>

 

(2.1.3)

 

Evaluating numerically to 25 digits yields the following result. 

>

 

(2.1.4)

 

Note that the integrand has a logarithmic singularity at x=0 as the following series expansion shows. 

>

 

(2.1.5)

 

The hybrid symbolic-numeric technique (automated in Maple) includes the concept of term-by-term integration of such a series expansion, employing symbolic integration for the individual terms. 

>

 

Example 2.2: Subtracting off a singularity 

>

 

>

 

(2.2.1)

 

>

 

(2.2.2)

 

>

 

Note that the graph of this integrand goes to at . However, this is an integrable singularity. It behaves like near . 

>

 

>

 

Note: The technique of subtracting off a singularity illustrated below, takes place automatically within Maple's numerical integration routines. 

Suppose it is desired to compute the result to 25 digits of accuracy. 

>

 

(2.2.3)

 

The generalized series expansion of    at    takes the following form. 

>

 

(2.2.4)

 

The non-regular part is 

>

 

(2.2.5)

 

The new expression 

>

 

(2.2.6)

 

is analytic on the interval [0, 1]. Thus it can be integrated easily by the default numerical integration method. 

>

 

(2.2.7)

 

>

 

(2.2.8)

 

>

 

Integrating    is easy because it has the indefinite integral 

>

 

(2.2.9)

 

and its definite integral can therefore be computed symbolically. 

>

 

(2.2.10)

 

Finally, summing the two values, we obtain the value for the original definite integration problem. 

>

 

(2.2.11)

 

>

 

Example 2.3: Algebraic transformation of variables 

>

 

>

 

(2.3.1)

 

>

 

(2.3.2)

 

Note: The change of variables illustrated below takes place automatically within Maple's numerical integration routines. 

The generalized series expansion of at is of the form 

>

 

(2.3.3)

 

Applying the change of variables  (i.e., ) yields 

>

 

(2.3.4)

 

>

 

The new integrand is analytic on the interval of integration. 

Thus it can be integrated easily by the default numerical integration method, even at high accuracy. 

Suppose that the result is desired to 50 digits of accuracy. 

>

 

(2.3.5)

 

>

 

Tensor Products 

• A tensor product is a finite sum of terms where each term is a product of univariate functions:

• There are many examples in mathematics:

• The rank of a tensor product is the minimum number of terms among all the equivalent representations.
  

• It is a natural tensor product when the factors of each term can be derived from the original function by a finite linear combination of linear functionals.

The Splitting Operator 

• For a bivariate function  and a point where  , the splitting operator is defined by
     .
  

• The point is called a splitting point.
  

• Note: splits    into a rank-one tensor product.
  

•  interpolates   on the lines and , which means that    and  .
  

• These properties allow us to generate a natural tensor product series simply by iterating the splitting operator while varying the choice of splitting point.

Tensor Product Series Expansions (a la Newton Interpolation Series) 

• Not every bivariate function    can be expressed as a tensor product of finite rank.

• For an approximation in   where    is continuous, we can uniformly approximate    by a rank-tensor product:

• Define the remainder :   

• The following simple algorithm generates a natural tensor product series expansion    of a given function    on .
  

• Refining terminology introduced by Chapman in his doctoral thesis, we call    the Geddes-Newton series expansion of    to    terms with respect to the splitting points   .

• Comment on computational cost: Note that Algorithm NTPS in its most general form as presented above has one expensive step. Specifically, step 2(a) calls for a two-dimensional search of the region for the point where the numerical maximum is attained for the absolute remainder .

• Achieving efficient computation: It turns out to be sufficient to replace step 2(a) by two one-dimensional searches (or a single one-dimensional search in the case of symmetric functions). Moreover, for the one-dimensional searches it is sufficient to sample a discrete set of points defined by the midpoints of previous splitting points. (This is based on computational experience; some proofs are needed here.)

Convergence of the Series (Informal) 

• Note similarities with the case of univariate polynomial interpolation via the Newton interpolation series.

• In both cases, addition of a new term in the series corresponds to the addition of a new interpolation point.

• However, in the case of Geddes-Newton series, the interpolation property holds not only at the new point but at all points on two lines passing through that point.

• Intuitively, for any function    which is continuous in the unit square , if the splitting points are dense along the diagonal then it is reasonable to expect that    might converge uniformly to    on as .

• In other words, by choosing well distributed splitting points on the diagonal, we can cover the whole unit square with small grid cells.  Recall: The remainder vanishes on the boundary of each cell .

• Some formal convergence results are discussed later in this presentation.

A Brief Glance at a Proof Technique 

• The easiest case for proving convergence is when    is analytic in a sufficiently large region.

• The convergence theorem and its proof are similar to the classical case where   is a function of a single complex variable and denotes the interpolating polynomial.

• Recall the univariate case (Ref:  P.J. Davis, Interpolation and Approximation [Davis63]), as summarized below.

• Basic tool: A contour integral representation of   .

Let    be analytic on    and suppose that the interpolation points satisfy   .   Let    be the polynomial interpolating   at points .   ,  .   Then for all   , where   .   Conclusion:  If    then .

• For our case, we consider a bivariate function    of two complex variables and    denotes the Geddes-Newton series expansion of    to    terms with respect to the splitting points   .

• Notation:    and   .

• We derive a contour integration formula for the remainder   :  

• Here,    and   .

• Then we can conclude that  where the constant    depends on the geometry, and depends on separation constants (like before).

• If   and if the   boundary constant satisfies as then .

Example 4.1 

The tensor product series of rank generated in Example 3.1 for the symmetric function    is represented in the following matrix-vector form. 

(4.1.1)

 

(4.1.2)

 

(4.1.3)

 

(4.1.4)

 

(4.1.5)

 

(4.1.6)

 

(4.1.7)

 

Since is a symmetric function, we see that the vector is the same as with replaced by and the matrix is simply the transpose of . 

With the above definitions for the matrices and vectors, the tensor product series of Example 3.1 is as follows where, for illustration, we show and as well as . 

Rank 1 approximation: 

(4.1.8)

 

Rank 2 approximation: 

(4.1.9)

 

Rank 3 approximation: 

(4.1.10)

 

Note: The size of the expressions representing tensor product approximations grows exponentially as the number of terms increases, if we form the explicit expressions. However, by maintaining the representation in matrix-vector form we avoid this expression growth. In our applications, the expressions in and (which appear only in the vectors and ) will be evaluated to numerical values before performing the matrix-vector multiplications. 

Transformation to a Canonical 2-D Integral 

• An integral with non-constant limits can be converted to a new integral with constant limits via a change of variables.

• Specifically,

using the change of variables: 

;     

• Therefore limit attention to integrals on the unit square.
  

• Making the integrand symmetric: For an integral on the unit square, we can always make the integrand symmetric, meaning

. 

• Namely,    (symmetric part + anti-symmetric part) where

;     

• Integral of   over the unit square is , which gives us:

. 

• Therefore we always have a symmetric integrand on the unit square.

Computational Results in 2-D 

The results in the following table are for the 2-D integral: 

Tol	Time	Result	Pts	RelErr
	0.3	0.069551393139	15
	3.7	0.06955139313890799	21
	5.3	0.0695513931389079901727	26
	7.2	0.069551393138907990172712980	33
	14.2	0.06955139313890799017271297964487	38

• Tol = requested accuracy tolerance

• Time = CPU time in seconds

• Result = computed value of the integral (to the number of digits requested)

• Pts = number of splitting points

• RelErr = relative error, compared with an accurate reference value

Integration in High Dimensions: A Brief Summary 

• As presented in our ISSAC'05 paper, we have extended our hybrid symbolic-numeric integration method to families of analytic integrands with special algebraic structure over hyperrectangular regions in higher dimensions.
  

• The method is based on recursive dimension-halving via separation of variables.
  

• Multidimensional integrals are reduced to one-dimensional integrals using Geddes-Newton series expansions to generate accurate tensor-product approximations of the integrand.
  

• Test results show that our new adaptive integration algorithm is effective both with high accuracy and in high dimensions.

Example 6.1: Potential SpeedUp for a Bessel function 

Consider , the Bessel function of the second kind. 

(6.1.1)

 

 

trying DCUHRE (TOMS algorithm 698)

 

cuhre: evalhf callback failed, trying evalf callback

 

cuhre: result=-8.64489683576170798
	(6.1.2)

 

(6.1.3)

 

The user information printed above identifies the cause of the slow computation: 

"evalhf callback failed" 

meaning that the function is unknown to the numerical library in C. Every function evaluation must be performed by calling back to Maple in software floating point mode rather than proceeding completely in hardware floating point. 

(6.1.4)

 

 

trying DCUHRE (TOMS algorithm 698)

 

cuhre: result=-.570570460032231774
	(6.1.5)

 

(6.1.6)

 

(6.1.7)

 

By developing appropriate approximations for as a function of two arguments, the speed of the numerical integration computation with integrand should be nearly as fast as the latter integral computation. 

Example 6.2: A function defined by a nonlinear BVP 

Consider the function defined by the following nonlinear boundary value problem on . 

(6.2.1)

 

(6.2.2)

 

Define the function by a procedure which numerically solves the BVP to obtain the function value. 

(6.2.3)

 

The computation is slow because each new evaluation point requires the numerical solution of a boundary value problem. 

Example 6.3: BesselY near z=0 

Consider again , the Bessel function of the second kind as in Example 6.1, and consider the series expansion first for the case when the order is a non-integer. 

(6.3.1)

 

In general for any order , the expansion about starts with a term in . 

(6.3.2)

 

Fact: For any order , as . 

In the overall approximation scheme for this function, evaluation at points where is close to will be handled by summing the series expansion. 

For our general bivariate approximation scheme based on tensor product series, we will restrict the region away from . In other words, the restriction is where, for example, we may choose . 

Still, polynomial approximation requires a high degree due to the nearby singularity. 

(6.3.3)

 

(6.3.4)

 

(6.3.5)

 

(6.3.6)

 

The command has computed the "infinite" Chebyshev series, truncated where the magnitude of the terms becomes smaller than the specified precision . Thus the degree obtained is a measure of the "approximation degree" required for the given function. 

Multiplication by yields a function which is finite at . 

(6.3.7)

 

(6.3.8)

 

(6.3.9)

 

This achieves some reduction in the degree required for polynomial approximation. 

However, the series expansion reveals a square root singularity at . For the case of integer order , the series expansion reveals a logarithmic singularity at . In both cases, rational approximation will be advantageous. 

Trying various rational degrees decreasing from degree , we find that degree is sufficient to achieve the desired accuracy . 

The resulting approximation for the function on the interval is of the form 

(6.3.10)

 

for numerator coefficients and denominator coefficients . 

The approximation for the original function on the interval is therefore 

. 

It can be verified that the maximum error is less than . 

Example 7.1: Tensor Product Series Approximation in a Region 

Consider the function , the Bessel function of the second kind, and suppose that we approximate this function in the region , in other words . This is a small region and no subdivisions of are required. 

The algorithm computes a tensor product series approximation based on the following distribution of splitting points. 

The matrix-vector representation of the tensor product series is as follows (carrying only 6 digits for this illustration). 

(7.1.1)

 

(7.1.2)

 

(7.1.3)

 

(7.1.4)

 

(7.1.5)

 

If one wanted the explict expression for the tensor product series it is given by 

However the expression is quite large and, more significantly, with the explicit expression it is impossible to evaluate to full accuracy in hardware floating point precision. The problem can be seen by the large diagonal elements appearing above, which produce a sum of terms with large values adding up to a relatively small final answer (a classical case of unstable numerical computation!). 

Fortunately, by maintaining the matrix-vector formulation we obtain a stable evaluation algorithm as will be seen. 

It can be verified that (when the digits are carried to hardware floating point precision): 

Example 7.2: Approximating the Univariate Functions 

Continuing with the approximation problem of Example 7.1, we must now compute approximations for each of the univariate functions appearing in the vectors and . 

(7.2.1)

 

(7.2.2)

 

The approach taken for each univariate function is to compute a Chebyshev series on the specified interval, truncated to accuracy . This determines the approximation degree required for the particular function and if this degree is too large (greater than degree 25) for any of the univariate functions then the corresponding interval is subdivided and the process is repeated. 

Since the degrees are not too large in this case, no subdivisions are required. Noting that rational approximations are more efficient than polynomial approximations for some functions, the algorithm will compute Chebyshev-Pade approximations if they lead to approximations of smaller total degree. 

For example, consider the entry and compute a Chebyshev series to the specified accuracy. Recall from Example 6.3 that for this function it is helpful to factor out a power of due to the singularity at and our algorithm will do so. 

(7.2.3)

 

(7.2.4)

 

(7.2.5)

 

(7.2.6)

 

(7.2.7)

 

(7.2.8)

 

(7.2.9)

 

This polynomial approximation is of degree . By trying successive Chebyshev-Pade approximants of degrees decreasing until the maximum error is too large, the algorithm finds that the approximant (for total degree 10) achieves the desired accuracy. 

(7.2.10)

 

Rational approximation is advantageous over polynomial approximation for this function because, as discussed in Example 6.3, even though we have factored out a power of to improve the leading term of the series expansion about , it remains the case that has a logarithmic singularity at . The left endpoint of the interval of approximation is which means that the singularity is "nearby". 

(7.2.11)

 

Noting that the approximation has been computed for the function , the final approximation for the original function is 

In general, Chebyshev-Pade approximants will be used in this manner for the functions of appearing in the vector . In contrast, for the functions of in the vector the algorithm finds that rational approximation offers no advantage and the truncated Chebyshev series is used. 

For example, consider  

(7.2.12)

 

(7.2.13)

 

(7.2.14)

 

One finds that rational approximations offer no advantage for the univariate functions of meaning that these functions are "polynomial-like". 

Remark on Data Structures 

For each univariate function appearing in and in , we must store information about the corresponding polynomial or rational approximation. The most general form of approximation to be stored is the form appearing for the function considered above, which is of the form: 

where in that particular example and where 

(7.3.1)

 

Therefore for each entry in and in , we store a vector containing the following purely numerical information: 

Evaluation of the approximation at a point requires evaluation of the numerator and denominator polynomials in Chebyshev form. Since there are several such polynomials in Chebyshev form (the numerator and denominator in the approximations for each entry in , for example), it can be advantageous to use a vectorized form of the evaluation algorithm. 

Similarly, there would be approximations of the same form for the univariate functions of ν appearing in the vector . As noted above, in some cases the approximation will be a polynomial rather than a rational function, in which case the denominator is simply 1. 

Example 7.4: Stable Evaluation of the Matrix-Vector Form 

Suppose that we have stored polynomial or rational approximations for each univariate function appearing in the vectors and of the preceding Example. Given our overall approximation of the function in the region , suppose that we wish to evaluate at the point . 

Since we have stored accurate approximations for the univariate functions, we will evaluate the vectors and to the following values (for simplicity here, we just enter the correct numerical values). 

We will use here for illustration, but we would normally work with hardware floating point precision. Also, for this precision the number of terms required would be . 

(7.4.1)

 

(7.4.2)

 

Then the tensor product series approximation in the region has the following matrix-vector representation which can be evaluated to yield the desired numerical value. 

(7.4.3)

 

(7.4.4)

 

Note that has entries which grow large and this can cause unstable numerical evaluation unless we choose the correct order of evaluation. The bracketing above specifies the correct order of evaluation. 

(7.4.5)

 

The following is an example of unstable evaluation. 

Suppose that we first compute the product of the three matrices. 

(7.4.6)

 

(7.4.7)

 

(7.4.8)

 

This result has no correct digits! 

Fortunately, using the right-to-left evaluation of the matrix-vector form can be shown to be numerically stable. 

Example 7.5: Approximation in a Large Region 

Consider once again the Bessel function of the second kind, , but this time in a larger region . Our algorithm was run for this example and the region was subdivided based on the requirements: 

and  

where is the number of terms in a tensor product series expansion and refers to the degree of the polynomials appearing in the numerator or denominator of the univariate approximations. 

The algorithm produces a hash function such that, given an evaluation point the result of is an integer index into the array of numerical values defining the approximation on the subregion. Specifically, the hash function for this example is the following Maple expression indicating that the region was subdivided into 12 subregions. 

(7.5.1)

 

For example, for evaluation at the point we get 

(7.5.2)

 

indicating that we must evaluate the approximation stored for the subregion indexed by 2. 

The 12 subregions for the above example are pictured in the plot below. 

As can be seen, the individual subregions have index numbers ranging from 1 to 12. 

Definition of the function BVP_H (Solution of a nonlinear BVP) 

The function was defined in a preceding section by the following nonlinear boundary value problem on . 

(8.1.1)

 

(8.1.2)

 

Define the function by a procedure which numerically solves the BVP to obtain the function value.