PUBLICATIONS DE L’INSTITUT MATHÉMATIQUE Nouvelle série, tome 100(114) (2016), 141–162 DOI: 10.2298/PIM1614141O

BIVARIATE GENERALIZED BERNSTEIN OPERATORS AND THEIR APPLICATION TO FREDHOLM INTEGRAL EQUATIONS

Donatella Occorsio and Maria Grazia Russo

Abstract. We introduce and study the sequence of bivariate Generalized Bernstein operators {B_{m,s} }_{m,s} , m, s ∈ ℕ,

B_{m,s} = I − (I − B_{m} )^{s} , B_{m}^{i} = B_{m} (B_{m}^{i−1} ),

where B_{m} is the bivariate Bernstein operator. These operators generalize the ones introduced and studied independently in the univariate case by Mastroianni and Occorsio [Rend. Accad. Sci. Fis. Mat. Napoli 44 (4) (1977), 151– 169] and by Micchelli [J. Approx. Theory 8 (1973), 1–18] (see also Felbecker [Manuscripta Math. 29 (1979), 229–246]). As well as in the one-dimesional case, for m fixed the sequence {B_{m,s} (f)}_{s} can be successfully employed in order to approximate “very smooth” functions f by reusing the same data points f( \frac{i}{m} , \frac{j}{m} ), i = 0, 1, …, m, j = 0, 1, …, m, since the rate of convergence improves as s increases. A stable and convergent cubature rule on the square [0, 1]^{2}, based on the polynomials B_{m,s} (f) is constructed. Moreover, a Nyström method based on the above mentioned cubature rule is proposed for the numerical solution of two-dimensional Fredholm integral equations on [0, 1]^{2}. The method is numerically stable, convergent and the involved linear systems are well conditioned. Some algorithm details are given in order to compute the entries of the linear systems with a reduced time complexity. Moreover the procedure can be significantly simplified in the case of equations having centrosymmetric kernels. Finally, some numerical examples are provided in order to illustrate the accuracy of the cubature formula and the computational efficiency of the Nyström method.

1. Introduction

The Generalized Bernstein polynomials B_{m,s} (F) of a continuous function F in [0, 1] were introduced in [15] (see also [16, 7]) and defined by

B_{m,s} (F; x) = ∑_{i=1}^{s} \binom{s}{i} (−1)^{i−1} B_{m}^{i}(F; x),

2010 Mathematics Subject Classification: 41A10, 41A63, 65D32, 65R20. Key words and phrases: iterated Bernstein polynomials, multivariate polynomial approxima-

tion, cubature formula, Nyström method.

Communicated by Gradimir Milovanović.

where B_{m}^{i} = B_{m}(B_{m}^{i−1} ), i = 1, …, s, s ∈ ℕ, and B_{m,1} (F) ≡ B_{m}(F) is the ordinary m-th Bernstein polynomial. Therefore B_{m}(F) approximates the function F by using the values of F on equidistant points of the interval [0, 1].

A remarkable property shared by the sequence {B_{m,s} (F)}_{m} is the improvement of the rate of convergence to the function F, as well as the smoothness of F increases (see [16]). This means that the sequence does not suffer of the saturation phenomena that occurs in the case of the classical Bernstein polynomials.

In the present paper we introduce the bivariate Generalized Bernstein (shortly GB) operator B_{m,s} , defined as the tensor product of two univariate operators B_{m,s} . The sequence {B_{m,s} (f)}_{m} uniformly converges to f, for any continuous function in the square S = [0, 1] × [0, 1]. Moreover for m fixed {B_{m,s} (f)}_{s} → L_{m,m} (f), where L_{m,m} (f) denotes the bivariate Lagrange polynomial interpolating f at the grid of (m + 1)^{2} equally spaced points in the square S.

We first prove that, similarly to the univariate case, the rate of convergence of the sequence {B_{m,s} (f)}_{m} to the function f improves, as well as the smoothness of the function increases. Moreover by incrementing the parameter s, with m fixed, the most relevant consequence is an acceleration of the approximation process by reusing the same values of the function f, computed for a fixed m.

An application of B_{m,s} is obtained by replacing f with B_{m,s} (f) in integrals of the type ∫_{S} f. In this way we obtain a cubature rule on the grid of equally spaced points in the square S. We prove that this formula, which can be easily implemented, is stable and convergent. As a consequence of the contribution of s in speeding up the convergence, we will discuss the connection between m and s in order to obtain the maximum rate of convergence of the cubature formula with the minimal number of data. Hence, for any m, a maximum value will be determined up to which incrementing s.

Finally we give a possible application of the aforesaid cubature rule in a Nyström method for approximating the solution of integral equations of the type

(1.1)  f(x, y) − μ ∫_{S} f(z, t)k(x, y, z, t)dz dt = g(x, y), (x, y) ∈ S,

where μ ∈ ℝ, k, g are given functions and f is the unknown function. This kind of equations is of interest in engineering areas, Computer Graphics, Mathematical Physics etc., where many problems can be modeled by one or two dimensional Fredholm integral equations of the second kind.

About the problem in two dimensions there are few results in the literature. The proposed numerical strategies make use of collocation, Galerkin or Nyström methods based on piecewise approximating polynomials [1, 9, 11] or Monte Carlo methods [10]. Recently Nyström methods based on the global polynomial approximation using the zeros of orthogonal polynomials have been proposed in [19] and [14] for equations defined on squares and triangles, respectively. Both of the involved procedures give very good results, especially when the kernel and the known function are smooth inside the domain and with possible singularities on the boundaries. However several times in the practice, the kernel and the right–hand side in (1.1) are only pointwise given on equally spaced grids. For this type of problems, if on one hand the numerical methods based on piecewise polynomials are cheap, on the other they provide a low order of convergence, even if the involved functions are very smooth.

We will prove that, under suitable conditions, the Nyström method based on the GB polynomials leads to a linear system which is uniquely solvable and wellconditioned too. In addition we will give error estimates in some Sobolev-type spaces and we will show that the rate of convergence of the method follows the degree of smoothness of the known functions.

Moreover we will discuss how the computational effort can be reduced when the kernel k is a centrosymmetric function, i.e.,

k(x, y, z, t) = k(1 − x, 1 − y, 1 − z, 1 − t).

Here, since the matrices of the linear systems are centrosymmetric too, the solution can be carried out with a reduction of time complexity exceeding the 90%.

The outline of the paper is as follows. Section 2 contains some notation and preliminary results. In Section 3 we define the bivariate Generalized Bernstein operator B_{m,s} proving some convergence results and studying the behavior of the sequence {B_{m,s} } w.r.t. both the parameters m and s. Section 4 is devoted to the cubature rule based on B_{m,s} , giving the computational details and some numerical tests. In Section 5 some results about the Nyström method approximating the solution of the integral equation (1.1) are given. A particular attention is paid to the cases of kernel functions leading to matrices with suitable symmetry properties and the consequent reduction of the computational cost. Some numerical tests are given showing the efficiency of the proposed method. Finally Section 6 contains the proofs.

2. Notations and preliminary results

All along the paper the notation C ≠ C(a, b, c, …) will mean that the positive constant C is independent of the parameters a, b, c, …. On the contrary C = C(a, b, c, …) will highlight that C depends on a, b, c, ….

Moreover by N_{0}^{m} we will indicate the set N_{0}^{m} = {0, 1, …, m} and therefore N_{0}^{m} × N_{0}^{m} = {(i, j)}_{i=0,1,...,m; j=0,1,...,m} .

From now on ℙ_{m,m} will denote the space of all bivariate algebraic polynomials of degree at most m in each variable separately.

Setting S := [0, 1] × [0, 1], C(S) will indicate the space of continuous functions in two variables, equipped with the uniform norm on the square S:

∥f∥_{∞} = max_{(x,y)∈S} |f(x, y)|.

Now set φ_{1}(x) = \sqrt{x(1 − x)} , φ_{2}(y) = \sqrt{y(1 − y)} and denote by f_{x} and f_{y} the function f(x, y) as a function of the only variable y or x respectively.

For smoother functions, i.e., for functions having some partial derivatives which can be discontinuous on the boundaries of S, we introduce the following Sobolev– type space

W_{r} = {f ∈ C(S) : M_{r}(f) := max{∥f_{y}^{(r)} φ_{1}^{r} ∥_{∞} , ∥f_{x}^{(r)} φ_{2}^{r} ∥_{∞} } < ∞}, r ⩾ 1,

where the superscript (r) denotes the rth derivative of the one-dimensional function f_{y} or f_{x}. W_{r} will be equipped with the norm ∥f∥_{W_{r}} = ∥f∥_{∞} + M_{r}(f).

Finally, following [14] we introduce a modulus of smoothness on C(S). Recalling the definition of the φ modulus of smoothness of Ditzian and Totik [6] for a one-variable continuous function F:

ω_{φ}^{k}(F, t) = sup_{h⩽t} max_{x∈[4h^{2} k^{2} ,1−4h^{2} k^{2} ]} |Δ_{hφ}^{k} F (x)| + inf_{P ∈ℙ_{k−1}} max_{x∈[0,4t^{2} k^{2} ]} |F (x) − P (x)|

+ inf_{P ∈ℙ_{k−1}} max_{x∈[1−4t^{2} k^{2} ,1]} |F (x) − P (x)|,

where φ(x) = \sqrt{x(1 − x)} and

Δ_{hφ}^{k} F (x) = ∑_{i=0}^{k} (−1)^{i−1} \binom{k}{i} F ( x + hφ(x)( \frac{k}{2} − i) ) ,

we set

(2.1)  ω_{S}^{k}(f, t) = max{ sup_{x∈[0,1]} ω_{φ}^{k}(f_{x}, t), sup_{y∈[0,1]} ω_{φ}^{k}(f_{y}, t)} .
It is well known [6] that if in addition max_{x∈[0,1]} |F^{(r)} (x)φ^{r}(x)| < ∞, then
(2.2)  ω_{φ}^{k}(F, t) ⩽ Ct^{r} max_{x∈[0,1]} |F^{(r)} (x)φ^{r}(x)|, k ⩾ r.

Therefore, if f ∈ W_{r}, by using (2.2) it results

(2.3)  ω_{S}^{k}(f, t) ⩽ Ct^{r}M_{r}(f), k ≧ r.

Now we recall some basic facts about the univariate Generalized Bernstein operator (see [8, 18] and the reference therein).

For any continuous function F in [0, 1] let B_{m}(F) be the m-th Bernstein polynomial

B_{m}(F; x) = ∑_{k=0}^{m} p_{m,k} (x)F (t_{k}), p_{m,k} (x) = \binom{m}{k} x^{k}(1 − x)^{m−k} ,

where t_{k} = \frac{k}{m} , k ∈ N_{0}^{m}. The Generalized Bernstein operators \{B_{m,s} \}_{s} are defined as

B_{m,s} = I − (I − B_{m})^{s}, B_{m}^{1} = B_{m}, B_{m}^{i} = B_{m}(B_{m}^{i−1} ), i = 2, …, s, s ∈ ℕ,

where I denotes the identity operator on the space of continuous functions. From the definition it follows that

(2.4)  B_{m,s} (F; x) = ∑_{j=0}^{m} p_{m,j}^{(s)} (x)F (t_{j}),

where

(2.5)  p_{m,j}^{(s)} (x) = ∑_{i=1}^{s} \binom{s}{i} (−1)^{i−1} B_{m}^{i−1} (p_{m,j} ; x).

Setting

p_{m}^{(s)} (x) := [p_{m,0}^{(s)} (x), p_{m,1}^{(s)} (x), …, p_{m,m}^{(s)} (x)]^{T} and p_{m}(x) := [p_{m,0} (x), …, p_{m,m} (x)]^{T},

the following vectorial expression holds true [17]

(2.6)  p_{m}^{(s)} (x)^{T} = p_{m}(x)^{T}C_{m,s} ,

where C_{m,s} ∈ ℝ^{(m+1)×(m+1)} is defined as

(2.7)  C_{m,s} = I_{m} + (I_{m} − A_{m} ) + … + (I_{m} − A_{m} )^{s−1}
= A_{m} ^{−1} [I_{m} − (I_{m} − A_{m} )^{s}] = [I_{m} − (I_{m} − A_{m} )^{s}]A_{m} ^{−1}

being A_{m} ∈ ℝ^{(m+1)×(m+1)} the matrix

(A_{m} )_{i,j} = p_{m,j} (t_{i}), (i, j) ∈ N_{0}^{m} × N_{0}^{m},

and I_{m} the identity matrix of order m + 1.

By induction on s, the following recurrence relation holds true

(2.8)  C_{m,2^{p}} = C_{m,2^{p−1}} + (I_{m} − A_{m} )^{2^{p−1}} C_{m,2^{p−1}} ,

which allows a fast construction of the subsequence {B_{m,2^{p}} }_{p=1,2,...} , by means of the identity

(2.9)  B_{m,2^{p}} (f; x) = 2B_{m,2^{p−1}} (f; x) − B_{m,2^{p−1}}^{2} (f; x).

Finally we recall that the eigenvalues of the matrix A_{m} are given by [4]

(2.10)  λ_{0}^{(m)} = λ_{1}^{(m)} = 1, λ_{i}^{(m)} = ∏_{j=1}^{i} ( 1 − \frac{j − 1}{m} ) , i = 2, 3, …, m.

Therefore, by (2.7), the eigenvalues μ_{k}^{(m,s)} of the matrix C_{m,s} are

(2.11)  μ_{k}^{(m,s)} = \frac{1 − (1 − λ_{k}^{(m)} )^{s}}{λ_{k}^{(m)}} , k ∈ N_{0}^{m}.

3. The bivariate GB operator

From now on denote by P_{ij} = (t_{i}, t_{j}), (i, j) ∈ N_{0}^{m} ×N_{0}^{m}, with t_{k} = \frac{k}{m} , k ∈ N_{0}^{m}. By B_{m,s} given in (2.4), we define the bivariate Generalized Bernstein operator

B_{m,s} as

B_{m,s} := B_{m,s} ⊗ B_{m,s} : C(S) → ℙ_{m,m} , m, s ∈ ℕ.

Let f ∈ C(S). For s = 1, B_{m,s} (f) reduces to the bivariate Bernstein polynomial (see for instance [4]),

B_{m}(f; x, y) = ∑_{i=0}^{m} ∑_{j=0}^{m} p_{m,i} (x)p_{m,j} (y)f(P_{ij} ).

Revisiting some properties holding true in the univariate case, it is not hard to deduce that

B_{m,s} (f; 0, 0) = f(P_{00}), B_{m,s} (f; 0, 1) = f(P_{01}),

B_{m,s} (f; 1, 0) = f(P_{10}), B_{m,s} (f; 1, 1) = f(P_{11}),

∑_{i=0}^{m} ∑_{j=0}^{m} p_{m,i}^{(s)} (x)p_{m,j}^{(s)} (y) = 1,

B_{m,s} (Q_{1}; x, y) = Q_{1}(x, y), Q_{1} ∈ ℙ_{1,1},

that means that B_{m,s} (f) interpolates f at the corners of the square S and preserves the bivariate polynomials of degree 1 in each variable separately.

By definition, and taking into account (2.4), the polynomial B_{m,s} (f; x, y) can be expressed as

(3.1)  B_{m,s} (f; x, y) = ∑_{i=0}^{m} ∑_{j=0}^{m} p_{m,i}^{(s)} (x)p_{m,j}^{(s)} (y)f(P_{ij} ),

with {p_{m,k}^{(s)} }_{k∈N_{0}^{m}} defined in (2.5). By (2.6), B_{m,s} (f) can be also represented in the Bernstein basis, i.e.,

(3.2)  B_{m,s} (f; x, y) = p_{m}(x)^{T}C_{m,s} F_{m}C_{m,s}^{T} p_{m}(y),

where the entries of the matrix F_{m} ∈ ℝ^{(m+1)×(m+1)} are

(3.3)  F_{m}(i, j) = f(P_{ij} ), (i, j) ∈ N_{0}^{m} × N_{0}^{m}.

Setting

(3.4)  G_{m,s} = C_{m,s} F_{m}C_{m,s}^{T} ,

by (3.2) it results

(3.5)  B_{m,s} (f; x, y) = p_{m}(x)^{T}G_{m,s} p_{m}(y),

i.e., the polynomial B_{m,s} (f) can be seen as the bivariate Bernstein polynomial of a continuous function g such that g(P_{ij} ) = (G_{m,s} )(i, j), (i, j) ∈ N_{0}^{m} × N_{0}^{m}.

By (2.9) the following recurrence relation can be deduced

B_{m,2^{p}} (f; x, y) = 2B_{m,2^{p−1}} (f; x, y) − B_{m,2^{p−1}}^{2} (f; x, y).

Therefore, for m fixed, we have a fast algorithm in order to construct the subsequence B_{m,2^{p}} (f), p = 1, 2, ….

Now we want to state some results about the convergence of the sequence {B_{m,s} (f)}_{m,s} , discussing either the case m → ∞ and s fixed, and also the other case s → ∞ and m fixed. Both situations show interesting aspects we want to highlight. Consider s fixed at first. By using a result in [8] we get

Theorem 3.1. Let f ∈ C(S). For any fixed s it results

(3.6)  (∀m)  ∥B_{m,s} (f)∥_{∞} ⩽ 2^{2s} ∥f∥_{∞} .

Moreover for m sufficiently large (say m > m_{0})

(3.7)  ∥f − B_{m,s} (f)∥_{∞} ⩽ C\{ω_{S}^{2s} ( f, \frac{1}{\sqrt{m}}) + \frac{∥f∥_{∞}}{m^{s}} \},

where C is a positive constant depending on s and independent of f and m.

Corollary 3.1. If f ∈ W_{r} and s ⩾ r/2, there holds

(3.8)  ∥f − B_{m,s} (f)∥_{∞} ⩽ C \frac{M_{r}(f)}{m^{r/2}} , C ≠ C(m, f), C = C(s).

Remark 3.1. Inequality (3.6) states the uniform boundedness of the operator sequence {B_{m,s} } w.r.t. m. As in the univariate case, {B_{m,s} (f)} comes out to be a sequence of so-called “quasi-interpolant” polynomials, whose rate of convergence, according to (3.7), improves as well as the smoothness of the function increases. Therefore, when f ∈ W_{2r}, for m → ∞, the sequence {B_{m,r} (f)}_{m} approximates f with order O(1/m^{r}), being r the smallest value of the parameter s for obtaining this order of convergence.

Now we explore the behavior of the polynomial sequence {B_{m,s} (f)}_{s} with respect to s, when m is fixed.

In what follows L_{m,m} (f) will denote the bivariate Lagrange polynomial interpolating f at the equally spaced grid points {P_{ij} }_{(i,j)∈N_{0}^{m} ×N_{0}^{m}} , i.e., L_{m,m} (f) ∈ ℙ_{m,m} and

L_{m,m} (f; P_{ij} ) = f(P_{ij} ), (i, j) ∈ N_{0}^{m} × N_{0}^{m}.

Theorem 3.2. Let m ∈ ℕ be fixed and assume s ⩾ 1. With λ_{m}^{(m)} and μ_{m}^{(m,s)} defined in (2.10) and (2.11), for any function f ∈ C(S) and for any grid point P_{ij} , (i, j) ∈ N_{0}^{m} × N_{0}^{m}, we have

(3.9)  |f(P_{ij} ) − B_{m,s} (f; P_{ij} )| ⩽ C∥f∥_{∞} (1 − λ_{m}^{(m)} )^{2s} ,

and

(3.10)  ∥B_{m,s} (f)∥_{∞} ⩽ C∥f∥_{∞} (μ_{m}^{(m,s)} )^{2},

where in both cases C ≠ C(s, f), C = C(m). Moreover, for any function f ∈ C(S),

(3.11)  lim_{s→+∞} B_{m,s} (f; x, y) = L_{m,m} (f; x, y)

holds uniformly w.r.t. (x, y) ∈ S.

Remark 3.2. Since (1 − λ_{m}^{(m)} )^{2s} = O(m^{−2s} ), (3.9) assures a very fast convergence at the grid points P_{ij} , (i, j) ∈ N_{0}^{m} × N_{0}^{m}. This kind of “point wise super convergence” can be useful in those contexts in which the quality of the process depends only on the behavior at the grid points. The Computer Aided Geometric Design, for instance, is a field, where this property is especially advisable.

Remark 3.3. Taking into account that

lim_{s}(μ_{m}^{(m,s)} )^{2} = ( \frac{1}{λ_{m}^{(m)}} )^{2},

(3.10) represents a bound for the operator norm when m is fixed and s → ∞, i.e.,

(3.12)  sup_{s} sup_{f≠0} \frac{∥B_{m,s} (f)∥_{∞}}{∥f∥_{∞}} ⩽ C( \frac{1}{λ_{m}^{(m)}} )^{2} ∼ \frac{e^{2m}}{2πm} < ∞.

This bound seems to be very pessimistic for “small" values of s. On the other hand (3.12) is not surprising, since {B_{m,s} (f)}_{s} uniformly converges to the Lagrange polynomial L_{m,m} (f) and the Lebesgue constants of the polynomial interpolation at equidistant points diverge exponentially (see for instance [13]).

4. The GB cubature rule

The above introduced GB operator can be usefully employed in the numerical cubature. Indeed for integrals of the type ∫_{S} f(x, y)dxdy it is possible to deduce the following cubature rule,

(4.1) ∫_{S} f(x, y)dxdy = ∫_{S} B_{m,s} (f; x, y)dxdy + R_{m,s} (f) =: Σ_{m,s} (f) + R_{m,s} (f).

By (3.5) and taking into account that

∫_{0}^{1} p_{m,k} (t)dt = \frac{1}{m + 1} , k ∈ N_{0}^{m},

it is not hard to prove that

(4.2)  Σ_{m,s} (f) = \frac{1}{(m + 1)^{2}} ∑_{i=0}^{m} ∑_{j=0}^{m} G_{m,s} (i, j),

where G_{m,s} is the matrix defined in (3.4).

An equivalent and more useful expression from the computational point of view, is

(4.3)  Σ_{m,s} (f) = ∑_{i=0}^{m} ∑_{j=0}^{m} D_{i,j}^{(s)} f(t_{i}, t_{j}),

where, for any (i, j) ∈ N_{0}^{m} × N_{0}^{m},

D_{i,j}^{(s)} = \frac{1}{m + 1} ( ∑_{r=0}^{m} C_{m,s} (r, i)) \frac{1}{m + 1} ( ∑_{k=0}^{m} C_{m,s} (k, j)) =: D_{i}^{(s)} D_{j}^{(s)} .

Now we prove that the rule is numerically stable and convergent and that the rate of convergence improves as well as the smoothness of the integrating function increases.

Theorem 4.1. With the notation used in (4.1)–(4.3) and for any f ∈ C(S), the cubature formula is convergent

(4.4)  |R_{m,s} (f)| ⩽ C\{ ω_{S}^{2s} ( f, \frac{1}{\sqrt{m}}) + \frac{∥f∥_{∞}}{m^{s}} \} , C ≠ C(m, f), C = C(s)

and numerically stable, i.e.,

(4.5)  sup_{m} ∑_{i=0}^{m} ∑_{j=0}^{m} |D_{i,j}^{(s)} | < ∞.

Moreover, for any f ∈ W_{r} and with s ⩾ \frac{r}{2}, the following estimate holds true

(4.6)  |R_{m,s} (f)| ⩽ \frac{C}{m^{r/2}} M_{r}(f), C ≠ C(m, f), C = C(s).

Remark 4.1. We want to discuss on the possible choices of m and s in order to obtain the maximum rate of convergence with the minimal number of data. For the sake of simplicity assume f ∈ W_{2s}, for s large enough, and define the sequences

(c_{s})_{s→∞} := ( \frac{M_{2s} (f)}{m^{s}} )_{s→∞} , (d_{m})_{m→∞} := ( \frac{M_{2s} (f)}{m^{s}} )_{m→∞} ,

which have different speeds of convergence, depending on the ratios

\frac{c_{s+1}}{c_{s}} ∼ \frac{1}{m}, \frac{d_{m+1}}{d_{m}} = ( 1 − \frac{1}{m + 1} )^{s},

respectively.

Since m and s can be chosen independently, we fix m and we note that the theoretical error decreases faster till s ⩽ s̃ := log(1/m)/ log(1 − \frac{1}{m+1} ), i.e., when c_{s+1} /c_{s} ⩽ d_{m+1} /d_{m}. Therefore for s increasing, but s ⩽ s̃, the best performance is obtained. For s > s̃ the situation is quite reversed. Indeed the ratio d_{m+1} /d_{m} is less than the other and then the efficacy of incrementing s is reduced. Therefore, in order to accelerate the convergence, m has to be increased. So, for any m the value s̃ represents the maximum one for choosing s, being useless to overcome it. Table 1 contains the values of the thresholds s̃ = s̃(m), computed for some values of m.

Table 1. Tresholds for s

\begin{matrix}m & 16 & 32 & 64 & 128 & 256 & 512 & 1024 \\ s̃ & 46 & 113 & 269 & 624 & 1423 & 3198 & 7102 \end{matrix}

Remark 4.2. Finally, we want to discuss about the stability of the cubature formula w.r.t. s, i.e., when m is fixed and s → ∞. Setting

σ_{m,s} := ∑_{i=0}^{m} ∑_{j=0}^{m} |D_{i,j}^{(s)} |,

as a consequence of (3.11), lim_{s→∞} σ_{m,s} converges to the stability constant of the Newton-Cotes cubature formula related to the grid points {P_{ij} }, (i, j) ∈ N_{0}^{m} ×N_{0}^{m}.

Nevertheless for the values of s < s̃ that in practice are used, with s̃ as in Table 1, {σ_{m,s} }_{s} seems to have a good behaviour.

Fig. 1 shows the growth of σ_{m,s} for increasing values of s. For any fixed m, σ_{m,s} was computed for s ∈ [2, 2^{20}].

4.1. Algorithm details and performance of the cubature rule. We start with some details about the computation of B_{m,s} (f), analyzing at first the computational effort in the general case and showing how it can be drastically reduced for suitable choices of s.

As we have remarked before, B_{m,s} (f) can be regarded as the Bernstein polynomial of a function g(x, y) s.t. g(P_{ij} ) = G_{m,s} (i, j), (i, j) ∈ N_{0}^{m} × N_{0}^{m}. Therefore, if the matrix G_{m,s} is known, the computation of B_{m,s} (f) can be performed by means of a double de Casteljau scheme w.r.t. g(x, y). As it is well known, this algorithm

Figure 1.

computes the univariate polynomial B_{m,s} (F; x) with m^{2} long operations for any x and it is numerically stable. Therefore, for any fixed (x, y) ∈ S the computational effort will be O(m^{4}).

Now we focus our attention on the matrix G_{m,s} = C_{m,s} F_{m}C_{m,s}^{T} , the construction of which is essentially reduced to that of C_{m,s} . On the other hand, for any couple (m, s), the computation of C_{m,s} can be performed only once, since it does not depend on the function f.

Starting by (2.7) we observe that A_{m} is a centrosymmetric matrix, i.e., A_{m} = J_{m}A_{m}J_{m}, where J_{m} is the counter-identity matrix of order m + 1 (J_{m}(i, j) = δ_{i,m−j} , (i, j) ∈ N_{0}^{m} × N_{0}^{m}, being δ_{h,k} the Kronecker delta). Therefore, as remarked in [20], the construction of A_{m} can be performed in \frac{m^{3}}{2} long operations. Furthermore, since the set of centrosymmetric matrices is closed under product, if U, V are centrosymmetric matrices the product UV can be performed in O(m^{3}/4) long operations (see for instance [22, 3]). Thus, we can conclude that the centrosymmetric matrix C_{m,s} can be constructed in (s − 2)m^{3}/4 flops, instead of (s − 2)m^{3} flops, i.e., with a 75% reduction in time complexity.

A more significant reduction is obtained by choosing s = 2^{p}, i.e., considering the subsequence {C_{m,2^{p}} }_{p=0,1,...} . Indeed, by using (2.8), the matrix C_{m,s} can be determined by 2(log_{2} s − 1) products of centrosymmetric matrices and therefore in \frac{m^{3}}{2} (log_{2} s − 1) flops. For instance, for s = 256, the cost is 3.5m^{3} instead of 63.5m^{3} needed in the incremental scheme one-by-one.

Now we give some details about the computation of the coefficients of the cubature rule. Setting

D_{m}^{(s)} = \frac{1}{(m + 1)^{2}} [D_{0}^{(s)} , D_{1}^{(s)} , …, D_{m}^{(s)} ]^{T},

we rewrite the cubature rule in the following vectorial form

(4.7)  Σ_{m,s} (f) = D_{m}^{(s)} ^{T}F_{m}D_{m}^{(s)} .

In view of the centrosymmetry of C_{m,s} it follows D_{i}^{(s)} = D_{m−i}^{(s)} , i = 0, 1, …, m, (or equivalently the vector D_{m}^{(s)} is centrosymmetric, i.e., D_{m}^{(s)} = J_{m}D_{m}^{(s)} J_{m}), and therefore only D_{i}^{(s)} , i = 0, 1, …, [m/2], have to be computed. Moreover, if the integrand f(x, y) satisfies f(x, y) = f(1 − x, 1 − y) for each (x, y) ∈ S, then F_{m}, defined in (3.3), is a centrosymmetric matrix, so it is completely defined by means of (m + 1)^{2} /2 computations of the function f(x, y).

A partial reuse of the values of the function can be performed by constructing a sequence of the cubature rules with m = 2^{q}, q ⩾ 1, i.e., {Σ_{2^{q},s} }_{q=1,2,...} . Indeed, denoting by F_{m} the (m + 1)-order matrix in (4.7), the entries of the (2m + 1)-order matrix F_{2m} can be obtained with a reduction of 25 % of the evaluations of the function.

In addition we remark that if the integrand f satisfies f(x, y) = f(1 − y, 1 − x) for each (x, y) ∈ S, the computational reduction is drastic since F_{m} comes out to be a persymmetric matrix (i.e., F_{m}^{T} = J_{m}F_{m}J_{m}) and, therefore, it is completely defined by means of (m + 1)^{2} /4 computations of f.

We close the section showing the performance of cubature rule (4.2) by means of some numerical tests.

In the tables we report the approximating values of the integrals. For any choice of m we consider different values of s. Moreover an empty cell means that there is no improvement in the computation w.r.t. the value written in the corresponding left cell.

Here and in what follows, all the computations were performed in double arithmetic machine precision 2.2204 × 10^{−16} .

Example 4.1. Consider the integral

∫_{S} \frac{sin(x + y)}{(1 + xy)^{4}} dxdy.

Its value with 17 exact digits is 0.35054764241461881 .

Here f ∈ W_{r}, for any r. According to Remark 4.1, we expect that for fixed m and increasing values of s, chosen accordingly to Table 1, the performance of the cubature rule improves.

As Table 2 shows the machine precision is attained for m = 256 and s = 8, or for m = 128 and s = 16 or for m = s = 64.

Example 4.2. Now, we consider the integral

∫_{S} \frac{e^{x^{2} +y^{2}}}{(1 + x + y)^{6}} dx dy,

whose value with 17 exact digits is 0.05731445500095343 .

Also in this case f ∈ W_{r} for any r and the machine precision is already attained for m = 128, s = 32. We remark that by using (2.8), the global computational cost is 2m^{3} instead of 7.5 m^{3}. The numerical results are contained in Table 3.

Table 2. Numerical results in Example 1

m s = 8  s = 16  s = 32  s = 64
8 0.3505  0.3505  0.35054
16 0.350547  0.3505476  0.3505476  0.35054764
32 0.35054764  0.350547642  0.3505476424  0.35054764241
64 0.3505476424  0.350547642414  0.3505476424146 0.350547642414619
128 0.3505476424146  0.350547642414619
256 0.350547642414619

Table 3. Numerical results in Example 2

m s = 8  s = 16  s = 32
8 0.057  0.0573  0.0573
16 0.05731  0.05731  0.057314
32 0.057314  0.0573144  0.057314455
64 0.057314455  0.0573144550  0.05731445500
128 0.05731445500  0.05731445500095  0.057314455000953
256 0.0573144550009  0.057314455000953
512 0.057314455000953

Example 4.3. The value of the integral

∫_{S} \frac{(1 − xy)^{8.1}}{1 + x^{7}y^{8}} dx dy

with 17 exact digits is 0.31202047436387431 .

Since f ∈ W_{16}, according to the estimate in (4.6), the error behaves like O(m^{−8} ) for any s s.t. s ⩾ 16. We point out that the seminorm M_{16}(f) is at least of order 8×10^{9}. This means that it is impossible, in practice, to obtain the machine precision in the computation. Neverthless the rate of convergence tells that m ∼ 720 nodes occur at least in order to get 14 significant correct digits. As Table 4 shows, this precision is obtained with significantly lower number of points.

Moreover we stress that until s is less than the threshold values, the speed of convergence really improves for m fixed and large s (see in particular the results for m = 16, 32).

Example 4.4. Finally, we consider the integral

∫_{S} \frac{(1 − xy)^{2.1}}{1 + x^{7}y^{8}} dx dy,

whose value with 16 exact digits is 0.5998045286943496 .

Since f ∈ W_{4}, according to the estimate in (4.6), the error behaves like O(m^{−2} ) for any s s.t. s ⩾ 4. In this case it is possible to get the machine precision with a value of m that is far below the theoretical one, as Table 5 show.

Table 4. Numerical results in Example 3

m s = 8  s = 16  s = 32
8 0.312  0.3120  0.31202
16 0.31202  0.312020  0.3120204
32 0.3120204  0.312020474  0.31202047436
64 0.312020474  0.312020474363  0.31202047436387
128 0.31202047436  0.31202047436387
256 0.31202047436387
m s = 64  s = 128

8

16 0.312020474
32 0.31202047436 0.312020474363

Table 5. Numerical results in Example 4

m s = 8  s = 16  s = 32  s = 64
16 0.599804  0.5998046  0.5998045  0.5998045
32 0.5998045  0.59980452  0.599804528  0.599804528
64 0.5998045286  0.5998045286  0.5998045286  0.5998045286
128 0.59980452869  0.59980452869  0.59980452869  0.59980452869
256 0.599804528694 0.599804528694  0.599804528694  0.599804528694
512 0.5998045286943 0.59980452869434  0.59980452869434 0.59980452869434
1024 0.5998045286943 0.599804528694349
m s = 128  s = 256  s = 512  s = 1024  s = 2048

16

32 5.998045286
64 0.5998045286 0.5998045286 0.5998045286 0.5998045286 0.59980452869
128 0.59980452869 0.599804528694
256 0.599804528694 0.5998045286943

5. The Fredholm integral equation on the square

Denoting by

Kf(x, y) = μ ∫_{S} k(x, y, z, t)f(z, t)dz dt

(1.1) can be rewritten in operatorial form as

(5.1)  (I − K)f = g,

where I is the identity operator on C(S). Here and in the sequel we will denote k_{(z,t)} for meaning that the function of four variables k is considered as a function of the only pair (x, y).

Using standard arguments, it is possible to prove that if k(x, y, z, t) is continuous then K : C(S) → C(S) is compact and consequently the Fredholm Alternative holds true for (5.1) in C(S) (see for instance [1]). Moreover, if for some r ∈ ℕ

(5.2)  sup_{(z,t)∈S} ∥k_{(z,t)} ∥_{W_{r}} < +∞,

then Kf ∈ W_{r} for any f ∈ C(S).

Starting with the cubature rule (4.2), we can define the following discrete operator

K_{m}f(x, y) = μ ∑_{i=0}^{m} ∑_{j=0}^{m} D_{i,j}^{(s)} k(x, y, t_{i}, t_{j})f(t_{i}, t_{j}),

and consider the operator equation

(5.3)  (I − K_{m})f_{m} = g,

where f_{m} is unknown. Collocating on the pairs (t_{h}, t_{ℓ}), (h, ℓ) ∈ N_{0}^{m} × N_{0}^{m}, the quantities β_{ij} = f(P_{i,j} ), (i, j) ∈ N_{0}^{m} × N_{0}^{m}, come out to be the unknowns of the linear system

(5.4) β_{hℓ} − μ ∑_{i=0}^{m} ∑_{j=0}^{m} D_{i,j}^{(s)} k(t_{h}, t_{ℓ}, t_{i}, t_{j})β_{ij} = g(t_{h}, t_{ℓ}), (h, ℓ) ∈ N_{0}^{m} × N_{0}^{m}.

The matrix solution (β_{ij}^{∗} )_{i,j=0,1...,m} of this system, if it exists, allows us to construct the Nyström interpolant in two variables

(5.5)  f_{m}(x, y) = μ ∑_{i=0}^{m} ∑_{j=0}^{m} D_{i,j}^{(s)} k(x, y, t_{i}, t_{j})β_{ij}^{∗} + g(x, y),

which will approximate the unknown f. Now denote by Γ_{m,s} the coefficient matrix of system (5.4), which is a (m + 1) block matrix, the entries of which are matrices of order m + 1.

Denoting by cond(Γ_{m,s} ) the condition number in infinity norm of Γ_{m,s} , the following Theorem holds:

Theorem 5.1. Assume that k is continuous in S×S and that Ker{I−K} = {0} in C(S). Denote by f the unique solution of (5.3) in C(S) for a given g ∈ C(S). If in addition, for some r ∈ ℕ, k satisfies (5.2), g ∈ W_{r}, and

(5.6)  sup_{(x,y)∈S} ∥k_{(x,y)} ∥_{W_{r}} < +∞,

then, for m sufficiently large, system (5.4) is uniquely solvable and well conditioned too, since

cond(Γ_{m,s} ) ⩽ C, C ≠ C(m), C = C(s).

Moreover, for any s ⩾ r/2 there results

(5.7)  ∥f − f_{m}∥_{∞} ⩽ C \frac{∥f∥_{W_{r}}}{m^{r/2}} ,

where C ≠ C(m, f) and C = C(s).

5.1. Algorithm details and performance of the Nyström method. In what follows we propose some numerical tests. In the tables for each m we give the maximum relative error attained in the computation of f_{m}, defined as in (5.5), at the grid of equally spaced points [0 : 0.1 : 1] × [0 : 0.1 : 1].

All the linear systems were solved by the Gaussian elimination with partial pivoting, and therefore the main computational effort was of the order of m^{6}/3.

The computation were performed in double arithmetic machine precision. Moreover an empty cell means that there is no improvement in the computation w.r.t. the value written in the corresponding left or up cell.

Example 5.1. We consider the following equation

f(x, y) − \frac{1}{5} ∫_{S} f(z, t)e^{−(1+x)(1+z)−(1+y)(1+t)} dz dt

= 1 − \frac{e^{−2(2+x+y)} (5e^{1+x} − 1)(e^{1+y} − 1)}{(1 + x)(1 + y)} ,

where μ = 0.2, k(x, y, z, t) = e^{−(1+x)(1+z)−(1+y)(1+t)} , and

g(x, y) = 1 − \frac{e^{−2(2+x+y)} (5e^{1+x} − 1)(e^{1+y} − 1)}{(1 + x)(1 + y)} ,

as well as the exact solution f(x, y) = 1. Here, the known functions are very smooth and according with (5.7) we expect a fast convergence. The numerical results confirm our expectation.

Table 6. Numerical results in Example 5

m s = 16  s = 32  s = 64  s = 128  s = 256
5 0.14 × 10^{−6}  0.48 × 10^{−7}  0.46 × 10^{−7}
10 0.94 × 10^{−9}  0.11 × 10^{−9}  0.29 × 10^{−10}  0.95 × 10^{−11}  0.14 × 10^{−11}
15 0.21 × 10^{−10}  0.14 × 10^{−11}  0.10 × 10^{−12}  0.17 × 10^{−13}  0.19 × 10^{−14}
20 0.11 × 10^{−11}  0.31 × 10^{−13}  0.15 × 10^{−14}  0.22 × 10^{−15}
30 0.13 × 10^{−14}  0.22 × 10^{−15}

Example 5.2. We consider the following equation

f(x, y) − 0.3 ∫_{S} f(z, t) cos(xz) cos(yt)dz dt = e^{xy},

where μ = 0.3, k(x, y, z, t) = cos(xz) cos(yt), g(x, y) = e^{xy}. The kernel and the known function g are very smooth. Also in this case we expect a rapid convergence. The numerical results, given in Table 7, confirm our expectation.

Example 5.3. This example can be found in [21]:

f(x′, y′) − ∫_{−1}^{1} ∫_{−1}^{1} f(z′, t′)[x′siny′+ t′e^{z′}] dz′dt′ = x′e^{−y′} + 4x′siny′− \frac{7}{3}.

By the transformations x′ = 2x−1, y′ = 2y −1, t′ = 2t−1, z′ = 2z −1, the equation becomes

Table 7. Numerical results in Example 6

m s = 16  s = 32  s = 64  s = 128  s = 256
10 0.22 × 10^{−8}  0.12 × 10^{−9}  0.12 × 10^{−9}  0.36 × 10^{−10}  0.17 × 10^{−11}
15 0.92 × 10^{−9}  0.31 × 10^{−10}  0.51 × 10^{−11}  0.60 × 10^{−13}  0.50 × 10^{−13}
20 0.12 × 10^{−9}  0.20 × 10^{−11}  0.47 × 10^{−13}  0.41 × 10^{−14}  0.23 × 10^{−14}
30 0.59 × 10^{−11}  0.32 × 10^{−13}  0.45 × 10^{−15}
40 0.59 × 10^{−12}  0.58 × 10^{−15}

f(x, y) − 4 ∫_{S} f(z, t)[(2x − 1) sin(2y − 1) + (2t − 1)e^{2z−1} ]dz dt

= (2x − 1)e^{1−2y} + 4(2x − 1) sin(2y − 1) − \frac{7}{3}.

i.e., the equation is of the type (1.1) with

μ = 4, k(x, y, z, t) = (2x − 1) sin(2y − 1) + (2t − 1)e^{2z−1} ,

g(x, y) = (2x − 1)e^{1−2y} + 4(2x − 1) sin(2y − 1) − \frac{7}{3}.

In this case, the solution is known and it is f(x, y) = (2x − 1)e^{1−2y} − 1. The numerical results are shown in Table 8.

Table 8. Numerical results in Example 7

m s = 16  s = 32  s = 64
5 0.36 × 10^{−9}  0.97 × 10^{−14}
10 0.58 × 10^{−13} 0.15 × 10^{−13} 0.19 × 10^{−13}
20 0.19 × 10^{−13}
30 0.19 × 10^{−13}

Now we give some computational details that are useful in constructing the linear system. We start from the structure of the matrix Γ_{m,s} of linear system (5.4).

It is a block-matrix of order (m + 1) having the following expression

Γ_{m,s} = [ \begin{matrix}Γ^{(0,0)} & Γ^{(0,1)} & … & Γ^{(0,m)} \\ Γ^{(1,0)} & Γ^{(1,1)} & … & Γ^{(1,m)} \\ Γ^{(2,0)} & Γ^{(2,1)} & … & Γ^{(2,m)} \\ … & … & … & … \\ Γ^{(m,0)} & Γ^{(m,1)} & … & Γ^{(m,m)} \end{matrix} ]

The blocks Γ^{(h,l)} of order m + 1 are defined as

Γ^{(h,ℓ)} = δ_{h,ℓ}I_{m} − μK_{m}^{(h,ℓ)} U_{ℓ}, (h, ℓ) ∈ N_{0}^{m} × N_{0}^{m},

where U_{ℓ} = diag(D_{ℓ,0}^{(s)}, D_{ℓ,1}^{(s)}, …, D_{ℓ,m}^{(s)} ), I_{m} denotes the identity matrix of order

m + 1 and the entries of the matrix K_{m}^{(h,ℓ)} , are

K_{m}^{(h,ℓ)} (i, j) = k(x_{h}, y_{i}, x_{ℓ}, y_{j}), (i, j) ∈ N_{0}^{m} × N_{0}^{m}.

Setting

K_{m}^{∗} = [ \begin{matrix}K_{m}^{(0,0)} & K_{m}^{(0,1)} & … & K_{m}^{(0,m)} \\ K_{m}^{(1,0)} & K_{m}^{(1,1)} & … & K_{m}^{(1,m)} \\ K_{m}^{(2,0)} & K_{m}^{(2,1)} & … & K_{m}^{(2,m)} \\ … & … & … & … \\ K_{m}^{(m,0)} & K_{m}^{(m,1)} & … & K_{m}^{(m,m)} \end{matrix} ]

D_{m} = diag(D_{0}^{(s)} , D_{1}^{(s)} , …, D_{m}^{(s)} ),
U_{m}^{∗} = D_{m} ⊗ D_{m},  I_{m}^{∗} = I_{m} ⊗ I_{m},

the matrix Γ_{m,s} can be rewritten as follows

Γ_{m,s} = I_{m}^{∗} − μK_{m}^{∗}U_{m}^{∗}.

Now we discuss on how the global time complexity can be reduced for some choices of the kernel k given in (5.1). For the sake of brevity we assume m + 1 let be even (the case m + 1 odd follows by a little bit elaborate scheme).

Assuming the kernel k satisfies

(5.8)  k(x, y, z, t) = k(1 − x, 1 − y, 1 − z, 1 − t),

then the matrix Γ_{m,s} is centrosymmetric and therefore, in view of [12, Th. 1], the following orthogonal block similarity holds

P^{T}Γ_{m,s} P = [ \begin{matrix}F_{1} \\ & F_{2} \\ & & F_{3} \\ & & & F_{4} \end{matrix} ] ,

where P is a known orthogonal matrix and the matrices F_{i} ∈ ℝ^{\frac{(m+1)^{2}}{4} × \frac{(m+1)^{2}}{4}} are expressed in terms of suitable blocks of Γ_{m,s} . The complete and detailed expressions of P and F_{i}, i = 1, …, 4 can be found in [12].

Therefore, the solution of the linear system (5.4) is obtained by solving 4 independent linear systems of orders (m + 1)^{2}/4. Hence the final solution is computed with a a 93% time complexity reduction (as shown in the graphic).

If in addition the kernel k is “bisymmetric", i.e.,

(5.9)  k(x, y, z, t) = k(y, x, t, z),

then the matrix Γ_{m,s} is both symmetric and centrosymmetric and therefore it is a persymmetric matrix. In this case, Γ_{m,s} is completely defined by means of (m + 1)^{2}/2 symmetric blocks of order m + 1.

Example 5.4. We consider the following equation

f(x, y) − 0.4 ∫_{S} f(z, t)|x − z|^{4.5}|y − t|^{7.3}dz dt = e^{x+y} ,

where μ = 0.4, k(x, y, z, t) = |x − z|^{4.5}|y − t|^{7.3}, g(x, y) = e^{x+y} . In this case k satisfies both properties (5.8) and (5.9).

Figure 2. Time complexity comparison

Table 9. Numerical results in Example 8

m s = 8  s = 16  s = 32  s = 64
10 0.86 × 10^{−4}  0.61 × 10^{−4}  0.18 × 10^{−4}  0.36 × 10^{−5}
20 0.13 × 10^{−5}  0.13 × 10^{−5}  0.85 × 10^{−7}  0.57 × 10^{−8}
30 0.86 × 10^{−7}  0.84 × 10^{−7}  0.16 × 10^{−8}  0.46 × 10^{−9}
40 0.10 × 10^{−7}  0.10 × 10^{−7}  0.18 × 10^{−9}
50 0.19 × 10^{−8}  0.19 × 10^{−8}  0.65 × 10^{−10}
60 0.46 × 10^{−9}  0.46 × 10^{−9}  0.31 × 10^{−10}
70 0.13 × 10^{−9}  0.13 × 10^{−9}  0.66 × 10^{−11}
80 0.47 × 10^{−10}  0.46 × 10^{−10}  0.46 × 10^{−11}

Since the kernel k ∈ W_{4} the expected rate of convergence is O(1/m^{2}). However, as Table 9 evidences, 11 significant digits are obtained with m = 80, while the theoretical estimates assure only 8 s.d.

About the time complexity, we observe that the solution has been obtained for m = 80 by solving four linear systems of order of 1600.

6. The proofs

First of all we state some results concerning the univariate operator B_{m,s} defined in (2.4).

Lemma 6.1. For any continuous function F in [0, 1] it results

(6.1)  max_{x∈[0,1]} |B_{m,s} (F, x)| ⩽ (2^{s} − 1) max_{x∈[0,1]} |F (x)|,

i.e., B_{m,s} is a bounded operator on the space of continuous functions.

Proof. By (2.4) and (2.6)

|B_{m,s} (F; x)| ⩽ max_{x∈[0,1]} |F (x)| ∑_{j=0}^{m} |p_{m,j}^{(s)} (x)|

= max_{x∈[0,1]} |F (x)| ∑_{j=0}^{m} | ∑_{i=0}^{m} p_{m,i} (x)C_{m,s} (i, j)|

⩽ max_{x∈[0,1]} |F (x)|| ∑_{i=0}^{m} p_{m,i} (x) ∑_{j=0}^{m} C_{m,s} (i, j)Q_{j} |

⩽ max_{x∈[0,1]} |F (x)|∥C_{m,s} ∥,

where Q_{j} (x) := sgn{ ∑ _{i=0}^{m} p_{m,i} (x)C_{m,s} (i, j)} and ∥D∥ = max_{1⩽k⩽n} ∑ _{h=1}^{n} |D(k, h)|

denotes the infinity norm of the matrix D ∈ ℝ^{n×n} .

Taking into account that (see [18]).

(6.2)  ∥C_{m,s} ∥ ⩽ 2^{s} − 1,

this lemma follows. □

Now we recall a result in [8] about the convergence of the univariate operator B_{m,s} (see also [18]). It was proved that for any continuous function F in [0, 1] there holds

(6.3)  max_{x∈[0,1]} |F (x) − B_{m,s} (F, x)| ⩽ C\{ ω_{φ}^{2s} ( F, \frac{1}{\sqrt{m}}) + \frac{max_{x∈[0,1]} |F (x)|}{m^{s}}\} ,

where C is a positive constant depending on s and independent of F and m.

Proof of Theorem 3.1. By (3.1) and (2.6) we get

|B_{m,s} (f; x, y)| ⩽ ∥f∥_{∞} ∑_{i=0}^{m} |p_{m,i}^{(s)} (x)| ∑_{j=0}^{m} |p_{m,j}^{(s)} (y)|

= ∥f∥_{∞} | ∑_{j=0}^{m} ∑_{k=0}^{m} p_{m,k} (x)C_{m,s} (k, j)Q_{j} || ∑_{i=0}^{m} ∑_{r=0}^{m} p_{m,r} (y)C_{m,s} (r, i)T_{i} |,

where for any fixed x, y,

Q_{j} := sgn( ∑_{k=0}^{m} p_{m,k} (x)C_{m,s} (k, j)) , T_{i} := sgn( ∑_{r=0}^{m} p_{m,r} (y)C_{m,s} (r, i)) .

Setting d_{k} = ∑ _{j=0}^{m} C_{m,s} (k, j)Q_{j} , h_{r} = ∑ _{i=0}^{m} C_{m,s} (r, i)T_{i}, we can conclude that

(6.4) |B_{m,s} (f; x, y)| ⩽ ∥f∥_{∞} | ∑_{k=0}^{m} p_{m,k} (x)d_{k} || ∑_{r=0}^{m} p_{m,r} (y)h_{r} | ⩽ ∥f∥_{∞} ∥C_{m,s} ∥^{2}

and therefore (3.6) follows by (6.2).

In order to prove (3.7) we remark that we can consider the function f(x, y) as a function of y (freezing x) or viceversa. Therefore, for instance, adding and subtracting the univariate Generalized Bernstein polynomial B_{m,s} (f_{x}) defined in (2.4), we can write

|f(x, y) − B_{m,s} (f; x, y)| ⩽ |f_{x}(y) − B_{m,s} (f_{x}, y)| + |B_{m,s} (f_{x}, y) − B_{m,s} (f; x, y)|

= |f_{x}(y) − B_{m,s} (f_{x}; y)| + |B_{m,s} (f_{y} − B_{m,s} (f_{y}); x)| (6.5)

=: I_{1}(x, y) + I_{2}(x, y)

By using (6.3) it easily follows

(6.6)  I_{1}(x, y) ⩽ C\{ ω_{φ}^{2s} ( f_{x}, \frac{1}{\sqrt{m}}) + \frac{max_{y∈[0,1]} |f_{x}(y)|}{m^{s}} \} .

Moreover, using (6.1) and (6.3) once again, we get

(6.7)  I_{2}(x, y) ⩽ C\{ ω_{φ}^{2s} ( f_{y}, \frac{1}{\sqrt{m}}) + \frac{max_{x∈[0,1]} |f_{y}(x)|}{m^{s}} \} .

Therefore, (3.7) follows using (6.6)–(6.7) in (6.5) and taking the maximum on (x, y) ∈ S. □

Proof of Corollary 3.1. If f ∈ W_{r} (3.8) follows immediately by using (2.3) in (3.7). □

Proof of Theorem 3.2. First of all set

l_{m,k} (z) = ∏_{i=0,i≠k}^{m} \frac{(z − t_{i})}{(t_{k} − t_{i})} , and l_{m}(z) = [l_{m,0} (z), l_{m,1} (z), …, l_{m,m} (z)]^{T}.

Therefore the bivariate Lagrange interpolating polynomial at the grid points P_{ij} = (t_{i}, t_{j}), i, j = 0, 1, …, m, can be written as follows

L_{m,m} (f; x, y) = ∑_{i=0}^{m} ∑_{j=0}^{m} l_{m,i} (x)l_{m,j} (y)f(P_{ij} ) = l_{m}(y)^{T}F_{m}l_{m}(x) = l_{m}(x)^{T}F_{m}l_{m}(y),

where F_{m} is defined in (3.3). Recalling that [17]

(6.8)  p_{m}(z)^{T}A_{m}^{−1} = l_{m}(z)^{T},

in order to prove (3.9), start from the following identities:

f(P_{ij} ) − B_{m,s} (f; P_{ij} ) = L_{m,m} (f; P_{ij} ) − B_{m,s} (f; P_{ij} )
= l_{m}^{T}(t_{i})F_{m}l_{m}(t_{j}) − p_{m}(t_{i})^{T}C_{m,s} F_{m}C_{m,s}^{T} p_{m}(t_{j})
= p_{m}(t_{i})^{T}[A_{m}^{−1} F_{m}A_{m}^{−T} − C_{m,s} F_{m}C_{m,s}^{T} ]p_{m}(t_{j})
=: p_{m}(t_{i})^{T}Q_{m,s} p_{m}(t_{j}).

By (2.7)

Q_{m,s} = A_{m}^{−1} (I − A_{m})^{s}F_{m}[(I − A_{m})^{s}]^{T}A_{m}^{−T}

and consequently

f(P_{ij} ) − B_{m,s} (f; P_{ij} ) = p_{m}(t_{i})^{T}Q_{m,s} p_{m}(t_{j})
= l_{m}(t_{i})^{T}(I − A_{m})^{s}F_{m}((I − A_{m})^{s})^{T}l_{m}(t_{j})

= ((I − A_{m})^{s}F_{m}((I − A_{m})^{s})^{T})_{ij} .

Therefore, denoting by ρ(D) the spectral radius of a matrix D,

|f(P_{ij} ) − B_{m,s} (f; P_{ij} )| ⩽ ∥(I − A_{m})^{s}∥_{∞}^{2} ∥f∥_{∞}
⩽ C∥f∥_{∞} ρ^{2}((I − A_{m})^{s}) = ∥f∥_{∞} (1 − λ_{m}^{(m)} )^{2s}

and (3.9) follows.

Now we want to prove (3.10). By (6.4) we get

|B_{m,s} (f; x, y)| ⩽ ∥f∥_{∞} ∥C_{m,s} ∥^{2} ⩽ C∥f∥_{∞} ρ^{2}(C_{m,s} )

with C = C(m). Therefore by (2.11), (3.10) follows.

Finally, taking into account (3.2) and lim_{s→∞} C_{m,s} = A_{m}^{−1} , by (6.8) we get

(∀f ∈ C(S)) lim_{s→∞} B_{m,s} (f; x, y) = l_{m}(y)^{T}F_{m}l_{m}(x) = L_{m,m} (f; x, y)

uniformly w.r.t. (x, y) ∈ S and (3.11) follows. □

Proof of Theorem 4.1. First we remark that (4.4) and (4.6) immediately follow by (3.7) and (3.8), respectively.

Moreover stability condition (4.5) of the cubature rule is obtained by applying the Uniform Boundedness Principle since (4.4) assures the convergence of the rule for any function in C(S) (see for instance [2]). □

Proof of Theorem 5.1. Holding (4.4), the Nyström method (5.3) is based on a cubature formula which is convergent for any continuous function. Hence, by standard arguments (see [1]) the method is stable and convergent in C(S) and well conditioned too. Moreover, it is also known that ∥f − f_{m}∥_{∞} ∼ ∥Kf − K_{m}f∥_{∞} .

On the other hand by (4.4) and taking into account (2.1), we get, since we are assuming s ⩾ r/2,

(6.9)  ∥Kf − K_{m}f∥_{∞} ⩽ C\{ \frac{∥fk_{(x,y)} ∥_{∞}}{m^{s}} + \frac{M_{r}(fk_{(x,y)} )}{\sqrt{m^{r}}} \} , C ≠ C(m, f).

Now assume for instance that M_{r}(fk_{(x,y)} ) = ∥\frac{∂^{(r)}}{∂z^{r}} (fk_{(x,y)} )φ^{r}∥_{∞} . In the same way we can proceed in the alternative case i.e., if M_{r}(fk_{(x,y)} ) = ∥\frac{∂^{(r)}}{∂t^{r}} (fk_{(x,y)} )φ^{r}∥_{∞} .

By (6.9) we get

∥Kf − K_{m}f∥_{∞} ⩽ \frac{C}{\sqrt{m^{r}}} ∑_{j=0}^{r} \binom{r}{j} ∥ \frac{∂^{(j)} f}{∂z^{j}} \frac{∂^{(r−j)} k_{(x,y)}}{∂z^{r−j}} φ^{r} ∥_{∞} ⩽ \frac{C}{\sqrt{m^{r}}} ∑_{j=0}^{r} a_{j} ∥f_{t}^{(j)} φ^{j} ∥_{∞}

where C ≠ C(m, f) and a_{j} = \binom{r}{j} sup_{(x,y)∈S} ∥\frac{∂^{(r−j)} k_{(x,y)}}{∂z^{r−j}} φ^{r−j} ∥_{∞} . Hence, since by (5.6) it follows that a_{j}, j = 0, …, r, are bounded, using [5, p. 310, Lemma 2.1], it results

∑_{j=0}^{r} a_{j} ∥f_{t}^{(j)} φ^{j} ∥_{∞} ⩽ C(∥f∥_{∞} + ∥f_{t}^{(r)} φ^{r}∥_{∞} ),

with C ≠ C(f), and then (5.7) follows. □

References

1. K. E. Atkinson, The Numerical Solution of Integral Equations of the second kind, Camb. Monogr. Appl. Comput. Math., 4. Cambridge University Press, Cambridge, 1997.

2. P. J. Davis, Interpolation & Approximation, Dover Publications, 1975.

3. I. T. Abu-Jeib, Algorithms for Centrosymmetric and Skew-Centrosymmetric Matrices, Missouri J. Math. Sci. 18(1) (2006), 1–8.

4. S. Cooper, S. Waldron, The eigenstructure of the Bernstein operator, J. Approx. Theory 105(1) (2000), 133–165.

5. Z. Ditzian, On interpolation of L_{p}[a, b] and weighted Sobolev spaces, Pacific J. Math. 90(2) (1980), 307–323.

6. Z. Ditzian, V. Totik, Moduli of smoothness, Springer Series in Computational Mathematics 9, Springer-Verlag, New York, 1987.

7. G. Felbecker, Linear kombinationen von iterierten Bernsteinoperatoren, Manuscripta Math. 29 (1979), 229–246.

8. H. H. Gonska, X. L. Zhou, Approximation theorems for the iterated Boolean sums of Bernstein operators, J. Comput. Appl. Math. 53 (1994), 21–31.

9. G. Han, R. Wang, Richardson extrapolation of iterated discrete Galerkin solution for twodimensional Fredholm integral equations, J. Comput. Appl. Math. 139(1) (2002), 49–63.

10. A. Keller, Instant Radiosity, Proceedings of the 24th annual conference on Computer graphics and interactive techniques (SIGGRAPH ’97). ACM Press/Addison-Wesley, New York, (1997), 49–56. doi:10.1145/258734.258769

11. F. Liang, F. R. Lin, A fast numerical solution method for two dimensional Fredholm integral equations of the second kind based on piecewise polynomial interpolation, Appl. Math. Comput. 216(10) (2010), 3073–3088.

12. Z. Liu, Y. Zhang, C. Ferreira, R. Ralha, On inverse eigenvalue problems for block Toeplitz matrices with Toeplitz blocks, Appl. Math. Comput. 216(6) (2010), 1819–1830.

13. G. Mastroianni, G.V. Milovanović, Interpolation Processes. Basic Theory and Applications, Springer Monogr. Math., Springer-Verlag, Berlin, (2009) xiv, 446 p. ISBN: 978-3-540-68346-9.

14. G. Mastroianni , G. Milovanović, D. Occorsio, A Nyström method for two variables Fredholm integral equations on triangles, Appl. Math. Comput. 219(14) (2013), 7653–7662.

15. G. Mastroianni, M. R. Occorsio, A generalization of the Bernstein operator, Rend. Accad. Sci. Fis. Mat. Napoli 44 (4) (1977), 151–169. (in Italian)

16. C. A. Micchelli, The saturation class and iterates of the Bernstein polynomials, Collection of articles dedicated to Isaac Jacob Schoenberg on his 70th birthday, I., J. Approx. Theory 8 (1973), 1–18.

17. D. Occorsio, A. C. Simoncelli, How to go from Bézier to Lagrange curves by means of generalized Bézier curves, Facta Univ. Ser. Math. Inform. (Niš) 11 (1996), 101–111.

18. D. Occorsio, Some new properties of Generalized Bernstein polynomials, Stud. Univ. Babeş-Bolyai Math. 56(3) (2011), 147–160.

19. D. Occorsio, M. G. Russo, Numerical methods for Fredholm integral equations on the square, Appl. Math. Comput. 218(5) (2011), 2318–2333.

20. D. Occorsio, M. G. Russo, A Nyström method for Fredholm integral equations based on equally spaced knots, Filomat, 28(1) (2014), 49–63.

21. A. Tari, Modified homotopy perturbation method for solving two-dimensional Fredholm integral equations, Int. J. Comput. Appl. Math. 5(5) (2010), 58–593.

22. Z. Peng, X. Hu, L. Zhang, The inverse problem of centrosymmetric matrices with a submatrix constraint, J. Comput. Math. 22(4) (2004), 535–544.

Department of Mathematics, Computer Science and Economics (Received 04 03 2015) University of Basilicata, Potenca (Revised 22 06 2015) Members of GNCS-INDAM, Italy

 donatella.occorsio@unibas.it
 mariagrazia.russo@unibas.it