Int. J. Anal. Appl. (2022), 20:60 Better Uniform Approximation by New Bivariate Bernstein Operators Asha Ram Gairola1, Suruchi Maindola1, Laxmi Rathour2,∗, Lakshmi Narayan Mishra3, Vishnu Narayan Mishra4 1Department of Mathematics, Doon University, Dehradun 248 001 Uttarakhand, India 2Ward No-16, Bhagatbandh, Anuppur 484 224, Madhya Pradesh, India 3Department of Mathematics, School of Advanced Sciences, Vellore Institute of Technology, Vellore 632 014, Tamil Nadu, India 4Department of Mathematics, Indira Gandhi National Tribal University, Amarkantak, Anuppur 484 887, Madhya Pradesh, India ∗Corresponding author: laxmirathour817@gmail.com, rathourlaxmi562@gmail.com Abstract. In this paper we introduce new bivariate Bernstein type operators BM,in (f ;x,y), i = 1,2,3. The rates of approximation by these operators are calculated and it is shown that the errors are signif- icantly smaller than those of ordinary bivariate Bernstein operators for sufficiently smooth functions. 1. Introduction In approximation theory the Bernstein operators are the most studied linear positive operators. Let, B[0,1] denote the class of bounded functions and f ∈ B[0,1], then the Bernstein operator Bn : B[0,1]→ C∞[0,1] of degree n with respect to f is defined by Bnf (x)= n∑ k=0 Pn,k(x)f ( k n ) where Pn,k(x)=   ( n k ) xk(1−x)n−k if k =0,1, ...,n, 0 if k < 0ork > n. Here C∞[0,1] is the class of infinitely differentiable functions on [0,1]. Received: Jun. 23, 2022. 2010 Mathematics Subject Classification. 41A35, 41A10, 41A25. Key words and phrases. linear operators; approximation by polynomials; rate of convergence. https://doi.org/10.28924/2291-8639-20-2022-60 ISSN: 2291-8639 © 2022 the author(s). https://doi.org/10.28924/2291-8639-20-2022-60 2 Int. J. Anal. Appl. (2022), 20:60 The operators Bnf (x) exhibit interesting approximation properties. We mention a few of them in the following. (1) Bnf (0)= f (0) and Bnf (1)= f (1) i.e. Bnf (x) has end point interpolation property, (2) The sequence Bnf (x) converges uniformly to f (x) whenever f is continuous on [0,1], (3) Bnf (x)≥ 0 for a positive function f (x), (4) The polynomials Bnf (x) are convex for a convex function f (x), (5) The Bernstein operators demonstrate simultaneous approximation property i.e. the derivative dr dtr ∣∣∣ t=x Bnf (t) converges to the corresponding derivative f (r)(x) whenever the later derivative exists. (6) The polynomials Bnf (x) provide a quadrature rule(see [1]) i.e. for f ∈ C[0,1],n ∈N ∫ 1 0 f (x)dx = lim n→∞ 1 n+1 n∑ k=0 f ( k n ) . For other interesting properties we refer to the monograph [2]. Inspite of extensive research work of more than one century, new and interesting results have been continuously appearing on the properties and applications of these polynomials, see [2] and the references therein. Although the operators Bnf (x) are easy to work with and show remarkable properties, they lack rapid convergence. In fact, the identity Bn(t2;x)= x2+ x(1−x) n implies that the order of approximation for the function, f (t)= t2 is not less than O(n−1) while the function is sufficiently smooth. In fact, it is known that the order of approximation is O(n−1/2) for the operators Bn(f ;x) and can not be improved as such. This fact is indicated by the following result due to Voronovskaya(see, [3]) Theorem 1.1. Let f (x) be bounded in [0,1] and suppose that the second derivative f ′′(x) exists at a certain point x of [0,1], then Bnf (x)− f (x)= 1 n ( x(1−x) 2 f ′′(x) ) +o ( 1 n ) . For the details see [4], pp. 102. In order to improve the degree of approximation by the Bernstein operators Arab et al. [5] have introduced kernel decomposition method recently. This method, then have been applied by Maria et al. in [6] and Gairola et al. in [7] for other sequence of operators. The aim of this paper is to extend this technique to modify the bivariate Bernstein operators and apply the modified operators for surface plotting. We show that the new bivariate operator approximate a surface better than the classical bivariate Bernstein operator (1.1). Let D denote the open square (0,1)× (0,1), and the closed square D := {(x,y) : 0≤ x ≤ 1,0≤ y ≤ 1}. Int. J. Anal. Appl. (2022), 20:60 3 The class B(D) consists of the functions f (x,y) defined and bounded on D. Then the bivariate Bernstein operators of degrees n,m with respect to f ∈ B(D) are defined by (cf. [8]) Bn,m(f ;x,y)= n∑ k=0 m∑ j=0 Pn,k(x)Pm,j(y)f ( k n , j m ) , (1.1) where Pn,k(x) = ( n k ) xk(1− x)n−k,k ∈ 0,n and Pm,j(y) = ( m j ) y j(1− y)m−j, j ∈ 0,m. We will need the class Cs(D),s ∈ N0 in order to establish estimates for the smooth functions. Cs(D) := { f ∈ C(D) : ∂rf ∂i∂r−i ∈ C(D), i ∈ 0, r, r ∈ 0,s } . Here, N0 :=N∪{0}. Approximation properties of multivariate Bernstein operators were studied by a number of re- searchers(see [8]- [16]) and they have established noteworthy results on degree of approximation, convergence properties etc. It is worth to mention the relevant work in [17]- [22]. In 1953 P.L. Butzer [23] discussed the simultaneous approximation properties of the operator 1.1. Moreover, Pop [24] obtained the rate of convergence in terms of Voronovskaja type asymptotic theorem and modulus of continuity for bivariate Bernstein operator Bm,n(f ;x,y) on the square D. In 1960 D.D. Stancu [25] defined another bivariate Bernstein operators on the triangle and derived the rate of con- vergence in terms of complete modulus of continuity. Subsequently, in 1992 Zhou [26] defined the two dimensional Bernstein-Durrmeyer operators. In [28], N. Deo introduced a certain modification of Bernstein operator and obtained direct results. The bivariate Bernstein operators have been also studied on a simplex by several authors(see [26]- [30]). By straightforward calculations it follows that Bn,m(x2+y2;x,y)= x2+y2+ x(1−x) n + y(1−y) m , so that as in the case of the univariate Berntsein operator Bn(f ;x), the degree of approximation is at most O ( 1 n , 1 m ) . Thus, this estimate can not be improved as such even for the sufficiently smooth function x2 + y2 on D. Although, these operators are easy to compute and demonstrate nice approximation properties, the slow rate of convergence has been a matter of concern. In order to achieve higher order of approximation we propose three bivariate Bernstein operators BM,rn,m(f ;x,y), r =1,2,3 as follows. 2. The first order operator BM,1n,m(f ;x,y) Let B(D) and C(D) denote the space of bounded functions and space of continuous functions f (x,y) on D with the defining norm ‖f‖ := sup (x,y)∈D |f (x,y)|. For f ∈ B(D), the operator BM,1nm (f ;x,y) is defined by BM,1n,m(f ;x,y)= n∑ k=0 m∑ j=0 PM,1 n,m;k,j (x,y)f ( k n , j m ) , (2.1) where PM,1 n,m;k,j (x,y)= PM,1 n,k (x)PM,1 m,j (y), PM,1 n,k (x)= a(x,n)Pn−1,k(x)+a(1−x,n)Pn−1,k−1(x), 4 Int. J. Anal. Appl. (2022), 20:60 PM,1 m,j (y)= b(y,m)Pm−1,j(y)+b(1−y,m)Pm−1,j−1(y) and a(x,n) = a1(n)x + a0(n), b(x,n) = b1(m)y + b0(m). Here ai(n),bi(m) n,m = 0,1,2, are unknown sequences which are to be determined by imposing suitable convergence conditions. For a1(n) = b1(m) = −1, a0(n) = b0(m) = 1 the operator (2.1) reduces to the ordinary bivariate operator (1.1). The total modulus of continuity ω(f ;δ1,δ2) of f ∈ C(D) is defined as ω(f ;δ1,δ2) := sup { |f (s,t)− f (x,y)| : s,t ∈ D, |s −x| ≤ δ1, |t −y| ≤ δ2 } , where δ1,δ2 > 0. The modulus ω(f ;δ1,δ2) satisfies the following properties. (i) ω(f ;δ1,δ2)→ 0, if δ1 → 0, δ2 → 0, (ii) |f (s,t)− f (x,y)| ≤ ω(f ;δ1,δ2) ( 1+ |s−x| δ1 )( 1+ |t−y| δ2 ) . The full modulus of continuity, ω(f ;δ) for f ∈ C(D) is defined as ω(f ;δ) := max√ (s−x)2+(t−y)2≤δ (x,y)∈D |f (s,t)− f (x,y)| . Other tools used to measure smoothness of a functions in C(D) are the partial modulus of continuity with respect to x and y and these are defined as ω(1)(f ;δ)= sup{|f (x1,y)− f (x2,y)| : y ∈ [0,1], |x1 −x2| ≤ δ} , ω(2)(f ;δ)= sup{|f (x,y1)− f (x,y2)| : x ∈ [0,1], |y1 −y2| ≤ δ} . 3. Preliminaries Let ei,j be the function ei,j(x,y) = x iy j. For a sequence Ln,m(f ;x,y) of bivariate linear positive operators we have the following extension of the Korovkin theorem in the square D. Theorem 3.1. [?] If {Ln,m} is a sequence of linear positive operators satisfying the conditions lim n,m→∞ ‖Ln,m(e00;x,y)−1‖C(D) =0, lim n,m→∞ ‖Ln,m(e10;x,y)−x‖C(D) =0, lim n,m→∞ ‖Ln,m(e01;x,y)−y‖C(D) =0, lim n,m→∞ ∥∥Ln,m(e201 +e210;x,y)− (x2 +y2)∥∥C(D) =0 then for any f ∈ C(D), lim n,m→∞ ‖Ln,m(f ;x,y)− f (x,y)‖C(D) =0. Lemma 3.1. The operator BM,1n,m verifies following identities: (1) BM,1n,m(e00,x,y)= (2a0(n)+a1(n))(2b0(m)+b1(m)), (2) BM,1n,m(e10,x,y)= [ x(2a0(n)+a1(n))+ (1−2x)(a0(n)+a1(n)) n ] (2b0(n)+b1(m)), Int. J. Anal. Appl. (2022), 20:60 5 (3) BM,1n,m(e01,x,y)= (2a0(n)+a1(n)) [ y(2b0(m)+b1(m))+ (1−2y)(b0(m)+b1(m)) m ] , (4) BM,1n,m(e11,x,y)= [ x(2a0(n)+a1(n))+ (1−2x)(a0(n)+a1(n)) n ] × [ y(2b0(m)+b1(m))+ (1−2y)(b0(m)+b1(m)) m ] , (5) BM,1n,m(e20,x,y)= [ x2(2a0(n)+a1(n))+ (4x−6x2)a0(n)+(3x−5x2)a1(n) n + (1−4x+4x2)(a1(n)+a0(n)) n2 ] (2b0(n)+b1), (6) BM,1n,m(e02,x,y)= (2a0(n)+a1(n)) [ y2(2b0(m)+b1(m))+ (4y−6y2)b0(m)+(3y−5y2)b1(m) m + (1−4y+4y2)(b1(m)+b0(m)) m2 ] , (7) BM,1n,m(e22,x,y)= [ x2(2a0(n)+a1(n))+ (4x−6x2)a0(n)+(3x−5x2)a1(n) n + (1−4x+4x2)(a1(n)+a0(n)) n2 ] [ y2(2b0(m)+b1(m))+ (4y−6y2)b0(m)+(3y−5y2)b1(m) m + (1−4y+4y2)(b1(m)+b0(m)) m2 ] , (8) BM,1n,m(e12,x,y)= [ x(2a0(n)+a1(n))+ (1−2x)(a0(n)+a1(n)) n ] [ y2(2b0(m)+b1(m))+ (4y−6y2)b0(m)+(3y−5y2)b1(m) m + (1−4y+4y2)(b1(m)+b0(m)) m2 ] , (9) BM,1n,m(e21,x,y)= [ x2(2a0(n)+a1(n))+ (4x−6x2)a0(n)+(3x−5x2)a1(n) n + (1−4x+4x2)(a1(n)+a0(n)) n2 ] × [ x(2a0(n)+a1(n))+ (1−2x)(a0(n)+a1(n)) n ] . Proof. The proof follows by straightforward calculations. � 4. Rate of Approximation In the next theorem we show that the linear positive operator sequence BM,1n,m defined by equation (2.1) converges uniformly to f with the help Theorem3.1 given by Volkov in [31]. Theorem 4.1. Let f ∈ B(D) and the sequences a0(n),a1(n),b0(m),b1(m) be convergent. If 2a0(n)+ a1(n)=2b0(m)+b1(m)=1 and the operator (2.1) is positive, then lim n,m→∞ BM,1n,m(f ;x,y)= f (x,y). Further, the convergence is uniform if f ∈ C(D). Proof. By lemma 3.1, for each i, j with 0≤ i ≤ 2,0≤ j ≤ 2 limn,m→∞BM,1n,m(eij,x,y)= eij, uniformly on D. Hence, lim n,m→∞ ∥∥BM,1n,m(e00;x,y)−1∥∥C(D) =0 lim n,m→∞ ∥∥BM,1n,m(e10;x,y)−x∥∥C(D) =0 lim n,m→∞ ∥∥BM,1n,m(e01;x,y)−y∥∥C(D) =0 lim n,m→∞ ∥∥Bn,m(e201 +e210;x,y)− (x2 +y2)∥∥C(D) =0 lim n,m→∞ ‖Bn,m(f ;x,y)− f (x,y)‖C(D) =0 Now the proof follows by the theorem 3.1. � 6 Int. J. Anal. Appl. (2022), 20:60 Theorem 4.2. (Voronovskaja Theorem)Let f ∈ C(D), fx, fy, fxx and fyy exist at a certain point (x,y) ∈ D and operator BM,1n,m(f ;x,y) be positive. If a1(n)+a0(n) = b1(m)+b0(m) = 0 and then we have, lim n,m→∞ 2n ( BM,1n,m(f ,x,y)− f (x,y) ) = ( (4x −6x2) lim n→∞ a0(n)+(3x −5x2) lim n→∞ a1(n)fxx ) + ( (4y −6y2) lim m→∞ b0(m)+(3y −5y2) lim m→∞ b1(m)fyy ) , where fxx and fyy are second order partial derivatives of f with respect to x and y respectively. Moreover, the above relation holds uniformly on D if fxx, fyy ∈ C(D). Proof. By Taylor’s formula, f ( k n , j m ) = f (x,y)+ (( k n −x ) ∂ ∂x + ( j m −y ) ∂ ∂y ) f (x,y) + 1 2! (( k n −x ) ∂ ∂x + ( j m −y ) ∂ ∂y )2 f (x,y) +ε ( k n , j m )(( k n −x ) ∂ ∂x + ( j m −y ) ∂ ∂y )2 f (x,y). The function ε(s,t) is bounded on D and lim(s,t)→(x,y))ε(s,t) = 0. Multiplying the above equation by PM,1 n,k (x)PM,1 n,j (y) and taking sum over i, j it follows that BM,1n,m(f ;x,y)= f (x,y)+B M,1 n,m((e10 −x);x,y)fx +B M,1 n,m((e01 −y);x,y)fy + 1 2 (BM,1n,m((e10 −x) 2;x,y)fxx +B M,1 n,m((e01 −y) 2;x,y)fyy +BM,1n,m(((e10 −x);x,y)((e01 −y);x,y)fxfy)) +BM,1n,m ( ε(e10,e01)(((e10,x,y)−x) fx +((e01,x,y)−y) fy)2 ;x,y ) . In view of the relation 2a0(n)+a1(n)=2b0(m)+b1(m)=1 and lemma 1, it follows that BM,1n,m(f ;x,y)= f (x,y)+ ( 1 2 (4a0(n)+3a1(n))x − (6a0(n)+5a1(n))x2 n ) fxx + ( 1 2 (4b0(m)+3b1(m))y − (6b0(m)+5b1(m))y2 m ) fyy +BM,1n,n ( ε(e10,e01)(((e10,x,y)−x) fx +((e01,x,y)−y) fy)2 ;x,y ) . (4.1) Since ε ( k n , j m ) is bounded on D and lim(s,t)→(x,y))ε(s,t)=0, proceeding in the standard manner, it follows that lim n,m→∞ BM,1n,m ( ε(e10,e01)(((e10,x,y)−x) fx +((e01,x,y)−y) fy)2 ;x,y ) =O ( 1 n , 1 m ) . (4.2) Now from equation (4.1) and (4.2) we obtain the desired result. � Int. J. Anal. Appl. (2022), 20:60 7 Theorem 4.3. Let BM,1n,m(f ;x,y) be the operator (2.1). If f ∈ C(D) and a0(n),a1(n),b0(m) and b1(m) are convergent sequences such that for each n and m, 2a0(n)+a1(n)=2b0(m)+b1(m)=1 then, ∣∣BM,1n,m(f ;x,y)− f (x,y)∣∣ ≤ 2 ( (1+3|a1(n)|)ω(1) ( f , n−1 2n √ n−2 ) +(1+3|b1(m)|)ω(2) ( f , m−1 2m √ m−2 )) . Proof. We have∣∣BM,1n,m(f ;x,y)− f (x,y)∣∣ ≤ ∣∣∣ n∑ k=0 m∑ j=0 ((a0(n)+a1(n))Pn−1,k(x)+(a0(n)+a1(n)(1−x))Pn−1,k−1(x)) ((b0(m)+b1(m))Pm−1,j(y)+(b0(m)+b1(m)(1−y))Pm−1,j−1(y))× × ( f ( k n , j m ) − f ( x, j m ))∣∣∣ + ∣∣∣ n∑ k=0 m∑ j=0 ((a0(n)+a1(n))Pn−1,k(x)+(a0(n)+a1(n)(1−x))Pn−1,k−1(x)) ((b0(m)+b1(m))Pm−1,j(y)+(b0(m)+b1(m)(1−y))Pm−1,j−1(y))× × ( f ( k n , j m ) − f ( k n ,y ))∣∣∣ ≤ ( |a0(n)+a1(n)| n∑ k=0 Pn−1,k(x)+ |(a0(n)+a1(n)(1−x))| n∑ k=0 Pn−1,k−1(x) ) × × ∣∣∣∣f ( k n , j m ) − f ( x, j m )∣∣∣∣ +  |b0(m)+b1(m)| m∑ j=0 Pm−1,j(y)+ |(b0(m)+b1(m)(1−y))| m∑ j=0 Pm−1,j−1(y)  × × ∣∣∣∣f ( k n , j m ) − f ( k n ,y )∣∣∣∣ ≤ |a0(n)+a1(n)| n∑ k=0 Pn−1,k(x)ω (1) (∣∣∣∣kn −x ∣∣∣∣ ) + |a0(n)+a1(n)(1−x)| n∑ k=0 Pn−1,k−1(x)ω (1) (∣∣∣∣kn −x ∣∣∣∣ ) + |b0(m)+b1(m)| m∑ j=0 Pm−1,j(y)ω (2) (∣∣∣∣ jm −y ∣∣∣∣ ) + |b0(m)+b1(m)(1−y)| m∑ j=0 Pm−1,j−1(y)ω (2) (∣∣∣∣ jm −y ∣∣∣∣ ) . Now using the property ω(f ,λδ)≤ (λ+1)ω(f ,δ), λ ≥ 0 8 Int. J. Anal. Appl. (2022), 20:60 of modulus of continuity, we get, ω(1) ( f ; ∣∣∣∣kn −x ∣∣∣∣ ) ≤ ( 1+ 2n √ n−2 n−1 ∣∣∣∣kn −x ∣∣∣∣ ) ω(1) ( f , n−1 2n √ n−2 ) . Similarly, ω(2) ( f ; ∣∣∣∣ jm −y ∣∣∣∣ ) ≤ ( 1+ 2m √ m−2 m−1 ∣∣∣∣ jm −y ∣∣∣∣ ) ω(2) ( f , m−1 2m √ m−2 ) . Therefore,∣∣BM,1n,m(f ;x,y)− f (x,y)∣∣ |a0(n)+a1(n)|ω(1) ( f , n−1 2n √ n−2 )( 1+ 2n √ n−2 n−1 n∑ k=0 Pn−1,k(x) ∣∣∣∣kn −x ∣∣∣∣ ) + |a0(n)+a1(n)(1−x)|ω(1) ( f , n−1 2n √ n−2 )( 1+ 2n √ n−2 n−1 n∑ k=0 Pn−1,k−1(x) ∣∣∣∣kn −x ∣∣∣∣ ) + |b0(m)+b1(m)|ω(2) ( f , m−1 2m √ m−2 )1+ 2m√m−2 m−1 m∑ j=0 Pm−1,j(y) ∣∣∣∣ jm −x ∣∣∣∣   + |b0(m)+b1(m)(1−y)|ω(2) ( f , m−1 2m √ m−2 )1+ 2m√m−2 m−1 m∑ j=0 Pm−1,j−1(y) ∣∣∣∣ jm −y ∣∣∣∣   . Now an application of Schwarz inequality and Lemma 1 yields, n∑ k=0 Pn−1,k(x) ∣∣∣∣kn −x ∣∣∣∣ ≤ ( n∑ k=0 Pn−1,k(x) ( k n −x )2)12 ( n∑ k=0 Pn−1,k(x) )1 2 = ( (2−n)x2 n2 + (n−1)x n2 )1 2 , and n∑ k=0 Pn−1,k−1(x) ∣∣∣∣kn −x ∣∣∣∣ ≤ ( n∑ k=0 Pn−1,k−1(x) ( k n −x )2)12 ( n∑ k=0 Pn−1,k−1(x) )1 2 = ( (2−n)x2 n2 + (n−3)x n2 + 1 n2 )1 2 . Since, max x∈[0,1],n>2 ( (2−n)x2 n2 + (n−1)x n2 ) = max x∈[0,1],n>2 ( (2−n)x2 n2 + (n−3)x n2 + 1 n2 ) = (n−1)2 4n2(n−2) therefore, n∑ k=0 Pn−1,k(x) ∣∣∣∣kn −x ∣∣∣∣ ≤ (n−1)2n√n−2, and similarly n∑ k=0 Pn−1,k−1(x) ∣∣∣∣kn −x ∣∣∣∣ ≤ (n−1)2n√n−2. Int. J. Anal. Appl. (2022), 20:60 9 In the same manner, m∑ j=0 Pm−1,j(y) ∣∣∣∣ jm −y ∣∣∣∣ ≤ (m−1)2m√m−2, m∑ j=0 Pm−1,j−1(y) ∣∣∣∣ jm −y ∣∣∣∣ ≤ (m−1)2m√m−2. Finally, using 2a0(n)+a1(n)=2b0(m)+b1(m)=1, we get∣∣BM,1n,m(f ;x,y)− f (x,y)∣∣ ≤ 4ω(1) ( f , n−1 2n √ n−2 )(∣∣∣∣a1(n)x + 1−a1(n)2 ∣∣∣∣+ ∣∣∣∣a1(n)(1−x)+ 1−a1(n)2 ∣∣∣∣ ) +4ω(2) ( f , m−1 2m √ m−2 )(∣∣∣∣b1(m)y + 1−b1(m)2 ∣∣∣∣+ ∣∣∣∣b1(m)(1−y)+ 1−b1(m)2 ∣∣∣∣ ) ≤ 2 [ (1+3|a1(n)|)ω(1) ( f , n−1 2n √ n−2 ) +(1+3|b1(m)|)ω(2) ( f , m−1 2m √ m−2 )] . The proof is now completed. � Remark 4.1. For a Lipschitz class we can established the degree of approximation by bivariate operator BM,1n,m(f ;x,y). For 0 < θ1 ≤ 1 and 0 < θ2 ≤ 1 and f ∈ B(D) the Lipschitz class LipM(θ1,θ2) is defined as follows LipM(θ1,θ2) := { f : |f (x1,x2)− f (y1,y2)| ≤ M |x1 −y1|θ1 |x2 −y2|θ2 } , where M > 0. So, if f (x,y) belongs to the class LipM(θ1,θ2) then, there holds the estimate∣∣BM,1n,m(f ;x,y)− f (x,y)∣∣ ≤ M(n,m)n θ12 m θ22 in view of the upper bounds ω(1) ( f , n−1 2n √ n−2 ) ≤ Cnθ1/2 and ω(2) ( f , m−1 2m √ m−2 ) ≤ Cmθ2/2. Here, the constant M(n,m)=2max{(1+3 |a1(n)|) ,(1+3 |b1(n)|)} . 5. The Second order operator BM,2n,m(f ;x,y) We extend our results to achieve 2nd order approximation on two variable. The proposed operator BM,2n (f ;x,y) is defined as BM,2n,m(f ;x,y)= n∑ k=0 m∑ j=0 PM,2 n,k (x)PM,2 m,j (y)f ( k n , j m ) , (5.1) where PM,2 n,k (x)= c(x,n)Pn−2,k(x)+d(x,n)Pn−2,k−1(x)+c(1−x,n)Pn−2,k−2(x), PM,2 m,j (y)= e(y,m)Pm−2(j)(y)+ f (y,m)Pm−2,j−1(y)+e(1−y,m)Pm−2(j −2)(y) and c(x,n)= c2(n)x 2 +c1(n)x +c0(n), d(x,n)= d0(n)x(1−x), e(y,m)= e2(m)y 2 +e1(m)y +e0, f (x,n)= f0(m)y(1−y). 10 Int. J. Anal. Appl. (2022), 20:60 As in the first order case, ci(n),ei(m), i =0,1,2, d0(n), f0(m) are the unknown sequences. Remark 5.1. If we put c2(n)= e2(m)=1,c1(n)= e1(m)=−2,c0(n)= e0(m)=1 and d0 = f0 =2 then the operator (5.1) becomes ordinary bivariate Bernstein operator (1.1). By straightforward calculation we obtain values of unknown sequences as follows. c2 = n 2 , c1 = −2−n 2 , c0 =1, d0 = n (5.2) e2 = m 2 , e1 = −2−m 2 , e0 =1, f0 = m. (5.3) Using these sequences we prove the following lemma. Lemma 5.1. For the operator BM,2n,m(f ,x,y) with the sequences chosen in (5.2)-(5.3), there holds the estimates (1) BM,2n,m(e00,x,y)=1, (2) BM,2n,m(e10,x,y)= x, (3) BM,2n,m(e01,x,y)= y, (4) BM,2n,m(e20,x,y)= x 2 + 2x(1−x) n2 , (5) BM,2n,m(e02,x,y)= y 2 + 2y(1−y) m2 , (6) BM,2n,m(e30,x,y)= x 3 + 2x(1−x)(5x−1) n2) − 6x(1−x)(2x−1) n3 , (7) BM,2n,m(e03,x,y)= y 3 + 2y(1−y)(5y−1) m2 − 6y(1−y)(2y−1) m3 . Proof. The proof follows by straightforward calculations. � Remark 5.2. By the lemma 5.1 we get improved rate of approximation for the functions eij, i, j =0(2) i.e. BM,2n,m((e01 − y),x,y) = 0,BM,2n,m((e10 − x),x,y) = 0,BM,2n,m((e20 − x2),x,y) = 2x(1−x) n2 and BM,2n,m((e02 −y2),x,y)= 2y(1−y) m2 . For twice differentiable function we obtain the following theorem. Theorem 5.1. Let BM,2n,m be the operator, f ∈ C3(D),n = m and let the sequences c2,c1,c0, e2,e1,e0,d0, f0 be bounded. Then, lim n→∞ n2(BM,2n,n (f ,x,y)− f (x,y))= x(1−x) ( fxx + (2x −1) 3 fxxx ) +y(1−y) ( fyy + (2y −1) 3 fyyy ) . Proof. We prove theorem for the cases when the operator BM,2n,n (f ,x,y) is positive. Since f ∈ C3(D), we can write Int. J. Anal. Appl. (2022), 20:60 11 f ( k n , j m ) = f (x,y)+ (( k n −x ) ∂ ∂x + ( j m −y ) ∂ ∂y ) f (x,y) + 1 2! (( k n −x ) ∂ ∂x + ( j m −y ) ∂ ∂y )2 f (x,y) + 1 3! (( k n −x ) ∂ ∂x + ( j m −y ) ∂ ∂y )3 f (x,y) + � ( k n , j m )(( k n −x ) ∂ ∂x + ( j m −y ) ∂ ∂y )3 f (x,y). (5.4) Multiplying equation (5.4) by PM,1 n,k (x)PM,1 n,j (y) and taking the sum over i, j we obtain BM,2n,n (f ;x,y)= f (x,y)+B M,2 n,n ((e10 −x);x,y)fx +B M,2 n,n ((e01 −y);x,y)fy + 1 2 (BM,2n,n ((e10 −x) 2;x,y)fxx +B M,2 n,n ((e01 −y) 2;x,y)fyy +BM,2n,n (((e10 −x);x,y)((e01 −y);x,y)fxfy) + 1 6 ( BM,2n,n ((e10 −x);x,y) 3fxxx +B M,2 n,n ((e01 −y) 3;x,y)fyyy ) +3BM,2n,n ((e10 −x);x,y)(e01 −y) 2;x,y)fxfyy +3BM,2n,n ((e10 −x) 2;x,y)(e01 −y);x,y)fxxfy) +BM,2n,n ( ε(e10,e01)(((e10,x,y)−x) fx +((e01,x,y)−y) fy)3 ;x,y ) here, BM,2n,n ( ε(e10,e01)(((e10,x,y)−x) fx +((e01,x,y)−y) fy)3 ;x,y ) ≤ � � is the error term which we have proved earlier for BM,1n,m, Now by lemma 5.1 we have, BM,2n,n (f ;x,y)− f (x,y) = 1 2 ( 2x(1−x) n2 fxx + 2y(1−y) n2 fyy ) + 1 6 ( 2x(1−x)(2x −1) n2 fxxx + 2y(1−y)(2y −1) n2 fyyy ) + � Finally, by standard computations, we have lim n→∞ n2(BM,2n,n (f ,x,y)− f (x,y))= x(1−x) ( fxx + (2x −1) 3 fxxx ) +y(1−y) ( fyy + (2x −1) 3 fyyy ) . � Remark 5.3. By theorem 5.1, for thrice differentiable functions, the degree of approximation is proved to be O(n−2) which is higher than the degree O(n−1) obtained by the ordinary operator Bn,n(f ;x,y). 12 Int. J. Anal. Appl. (2022), 20:60 6. The Third Order operator BM,3n,m(f ;x,y) The third order operator BM,3n,m(f ;x,y) on D is defined by BM,3n,m(f ;x,y)= n∑ k=0 m∑ j=0 f ( k n , j m ) PM,3 n,k (x)PM,3 m,j (y) (6.1) where, PM,3 n,k (x) = g(x,n)Pn−4,k(x)+ r(x,n)Pn−4,k−1(x)+ s(x,n)Pn−4,k−2(x) +r(1−x,n)Pn−4,k−3(x)+g(1−x,n)Pn−4,k−4(x), PM,3 m,j (y) = h(y,m)Pm−4,j(y)+ t(y,m)Pm−4,j−1(y)+w(y,m)Pm−4,j−2(y) +t(1−y,m)Pm−4,j−3(y)+g(1−y,n)Pm−4,j−4(y). and g(x,n)= g4(n)x 4 +g3(n)x 3 +g2(n)x 2 +g1(n)x +g0(n), r(x,n)= r4(n)x 4 + r3(n)x 3 + r2(n)x 2 + r1(n)x + r0(n), s(x,n)= s0(n)(x(1−x)2)2, h(y,m)= h4(m)y 4 +h3(m)y 3 +h2(m)y 2 +h1(m)y +h0(m), t(y,m)= t4(m)y4 + t3(m)y 3 + t2(m)y 2 + t1(m)y + t0(m), w(y,m)= w0(m)(y(1−y)2)2. Note 1. If we put g4(n) = h4(m) = 1,g3(n) = h3(m) = −4,g2(n) = h2(m) = 6,g1(n) = h1(m) = −4,g0(n) = h0(m) = 0, r4(n) = t4(m) = −4, r3(n) = t3(m) = 12, r2(n) = t2(m) = −12, r0(n) = t0(m) = 0 and s0(n) = w0(m) = 6,then the operator become the Bernstein operator (1.1) of two variable. By straightforward calculations we have the following values of the unknown sequences gi(n),hi(m), ri(n),ti(m), i =0,1,2,3,4. s0(n) and w0(m) for which the operator BM,3n,m(f ;x,y) attains the order O ( 1 n2 , 1 m2 ) . g4(n)=1+ 23 12 n+ 1 8 n2, g3(n)=−4− 143 n− 1 4 n2, g2(n)=6+ 10 3 n+ 1 8 n2, g1(n)=−4− 712n, g0(n)=1 r4(n)=−4− 233 n− 1 2 n2, r3(n)=12+17n+n 2, r2(n)=−12− 313 n− 1 2 n2 r1(n)=4+n, r0(n)=0, s0(n)=6+ 23 2 n+ 3 4 n2 h4(m)=1+ 23 12 m+ 1 8 m2, h3(m)=−4− 143 m− 1 4 m2, h2(m)=6+ 10 3 m+ 1 8 m2, h1(m)=−4− 712m, h0(m)=1 t4(m)=−4− 233 m− 1 2 m2, t3(m)=12+17m+m 2, t2(m)=−12− 313 m− 1 2 m2 t1(m)=4+m, t0(m)=0, w0(m)=6+ 23 2 m+ 3 4 m2. Int. J. Anal. Appl. (2022), 20:60 13 Lemma 6.1. The operator BM,3n,m(f ,x,y) together with the above values gi(n),hi(m), ri(n),ti(m), s0(n) and w0(m), i =0(4) verifies (1) BM,3n,m(e00,x,y)=1 (2) BM,3n,m(e10,x,y)= x (3) BM,3n,m(e01,x,y)= y (4) BM,3n,m(e20,x,y)= x 2 (5) BM,3n,m(e02,x,y)= y 2 (6) BM,3n,m(e30,x,y)= x 3 (7) BM,3n,m(e03,x,y)= y 3 (8) BM,3n,m(e 2 10 +e 2 01;x,y)= x 2 +y2. Proof. The proof follows by straightforward computations. � Proceeding along the lines of theorem 5.1, we obtain following theorem. Theorem 6.1. Let f ∈ C4(D) and together the values of gi,hi, ri,ti,si,wi, we have lim n,m→∞ (BM,2n,n (f ,x,y)− f (x,y))=O ( 1 n3 , 1 m3 ) . 7. Numerical Verification To verify the order of approximation by operator BM,3n,m(f ;x,y) we choose the function sinπx cosπy and for its images under the operator BM,3n,m for n =5,10 and m =5,10. Figure 1. Comparison of the function sinπx cosπy with BM,3n,m(f ,x,y) for n = m = 5,n = m =10. 14 Int. J. Anal. Appl. (2022), 20:60 Table 1 provides values of the function sinπx cosπy at equidistant nodes (x,y) with step size 0.2 and Tables 2, 3 provide values of the functions BM,35,5 (f ;x,y) and B M,3 10,10(f ;x,y) at the corresponding nodes. Table 1. f (x,y)= sinπx cosπy, (starred ∗ rows are of order 10−17) 0 0.2 0.4 0.6 0.8 1 0. 0. 0. 0. 0. 0. 0. 0.2 0.587785 0.475528 0.181636 -0.181636 -0.475528 -0.587785 0.4 0.951057 0.769421 0.293893 -0.293893 -0.769421 -0.951057 0.6 0.951057 0.769421 0.293893 -0.293893 -0.769421 -0.951057 0.8 0.587785 0.475528 0.181636 -0.181636 -0.475528 -0.587785 ∗1 12.2465 9.9076 3.78437 -3.78437 -9.9076 -12.2465 Table 2. BM,35,5 (f ;x,y) 0 0.2 0.4 0.6 0.8 1 0. 0. 0. 0. 0. 0. 0. 0.2 0.643288 0.527466 0.200931 -0.200931 -0.527466 -0.643288 0.4 1.03431 0.848086 0.323067 -0.323067 -0.848086 -1.03431 0.6 1.03431 0.848086 0.323067 -0.323067 -0.848086 -1.03431 0.8 0.643288 0.527466 0.200931 -0.200931 -0.527466 -0.643288 1. 0. 0. 0. 0. 0. 0. Table 3. BM,310,10(f ;x,y), (starred ∗ rows are of order 10−12) 0 0.2 0.4 0.6 0.8 1 0. 0. 0. 0. 0. 0. 0. 0.2 0.594148 0.482834 0.184553 -0.184553 -0.482834 -0.594148 0.4 0.961538 0.781392 0.298671 -0.298671 -0.781392 -0.961538 0.6 0.961538 0.781392 0.298671 -0.298671 -0.781392 -0.961538 0.8 0.594148 0.482834 0.184553 -0.184553 -0.482834 -0.594148 ∗1 5.38325 4.37469 1.67213 -1.67213 -4.37469 -5.38325 It follows that the maximum absolute errors by BM,3n,n (f ;x,y) for n =5,10 are found to be 0.0943346 at (0.5,0.919599) and 0.0134438 at (0.5,0.879425) respectively. Hence, in the figure 2 we compare f (x,y) with BM,3n,n (f ;x,y) for n =5 and n =10 in the neighbourhood of the point (0.5,0.8). Int. J. Anal. Appl. (2022), 20:60 15 Figure 2. Comparison of the function sinπx cosπy with BM,3n,m(f ,x,y) for n = m = 5,n = m =10. The next example compares BM,2n,n (f ;x,y) with B M,3 n,n (f ;x,y) for n =5 (See Fig. 3.) corresponding to the sufficiently smooth function e−x−y. It is observed that the accuracy increases as we move from BM,2n,n (f ;x,y) to B M,3 n,n (f ;x,y) for a fixed number of nodes. Figure 3. Comparison of the function e−x−y with BM,2n,m(f ,x,y), B M,3 n,m(f ,x,y) for n = 5,m =5. 16 Int. J. Anal. Appl. (2022), 20:60 Figure 4. Comparison of the function e−x−y with BM,2n,m(f ,x,y), B M,3 n,m(f ,x,y) for n = 5,m =5 in small rectangle The following tables show the values of f (x,y) and corresponding values under the operator BM,3n,m for n =5,10,15,20m =5,10,15,20 at different nodes with steps size 0.1. Table 4. f (x,y)= e−(x+y) 0 0.2 0.4 0.6 0.8 1 0 1. 0.818731 0.67032 0.548812 0.449329 0.367879 0.2 0.818731 0.67032 0.548812 0.449329 0.367879 0.301194 0.4 0.67032 0.548812 0.449329 0.367879 0.301194 0.246597 0.6 0.548812 0.449329 0.367879 0.301194 0.246597 0.201897 0.8 0.449329 0.367879 0.301194 0.246597 0.201897 0.165299 1 0.367879 0.301194 0.246597 0.201897 0.165299 0.135335 Table 5. BM,35,5 (f ;x,y) 0 0.2 0.4 0.6 0.8 1 0 1. 0.819147 0.670916 0.549393 0.449698 0.367879 0.2 0.819147 0.671002 0.549579 0.450034 0.368369 0.301347 0.4 0.670916 0.549579 0.450129 0.368597 0.30171 0.246816 0.6 0.549393 0.450034 0.368597 0.301833 0.247061 0.202111 0.8 0.449698 0.368369 0.30171 0.247061 0.202228 0.165435 1 0.367879 0.301347 0.246816 0.202111 0.165435 0.135335 Int. J. Anal. Appl. (2022), 20:60 17 Table 6. BM,310,10(f ;x,y) 0 0.2 0.4 0.6 0.8 1 0 1. 0.818788 0.670397 0.548882 0.44937 0.367879 0.2 0.818788 0.670414 0.548913 0.449418 0.367939 0.301215 0.4 0.670397 0.548913 0.449432 0.367969 0.301257 0.246625 0.6 0.548882 0.449418 0.367969 0.301271 0.246651 0.201922 0.8 0.44937 0.367939 0.301257 0.246651 0.201934 0.165314 1 0.367879 0.301215 0.246625 0.201922 0.165314 0.135335 Table 7. BM,320,20(f ;x,y) 0 0.2 0.4 0.6 0.8 1 0 1. 0.818738 0.67033 0.54882 0.449334 0.367879 0.2 0.818738 0.670332 0.548825 0.44934 0.367887 0.301197 0.4 0.67033 0.548825 0.449342 0.367891 0.301202 0.246601 0.6 0.54882 0.44934 0.367891 0.301204 0.246604 0.2019 0.8 0.449334 0.367887 0.301202 0.246604 0.201901 0.165301 1 0.367879 0.301197 0.246601 0.2019 0.165301 0.135335 The global errors of approximation of the function by BM,25,5 (f ,x,y)and B M,3 5,5 (f ,x,y)are0.00975603 and 0.000803997 which are obtained at (0.278459,0.278459) and (0.365884,0.365884) respectively. Here, again we choose small intervals for comparison in Fig 4 to show the three surfaces sufficiently separated. 8. Conclusion and Future Scope The surface plotting of a smooth function by a bivariate linear operators is constructed using the m,n nodes k/n, j/m k = 0,n, j = 0,m. The operator BM,jn,m(f ;x,y), j = 1,2,3 approximate such surfaces with the degree of convergence O(n−j), j = 1,2,3 which is significantly high in comparison of approximation by regular bivariate operator Bn,m(f ;x,y). We observe that the global absolute error decreases with the order O(n−1) from BM,2n,n (f ;x,y) to B M,3 n,n (f ;x,y) while keeping the degree of the polynomials same. On the other hand choosing large number of nodes can also improve the degree of approximation however, for a given fixed number of values it is better to apply the methods B M,j n,m(f ;x,y) for j = 2,3, ... It is of interest to determine exact class of the functions for which the methods BM,jn,m(f ;x,y) provide optimal degree of approximation. This can be considered as an open problem. Authors’ Contributions: All authors contributed equally and significantly in writing this paper. All authors read and approved the final manuscript. Conflicts of Interest: The authors declare that there are no conflicts of interest regarding the publi- cation of this paper. 18 Int. J. Anal. Appl. (2022), 20:60 References [1] J. Reinkenhof, Differentiation and Integration Using Bernstein’s Polynomials, Int. J. Numer. Meth. Eng. 11 (1977), 1627–1630. https://doi.org/10.1002/nme.1620111012. [2] J. Bustamante, Bernstein Operators and Their Properties, Springer, Cham, 2017. https://doi.org/10.1007/ 978-3-319-55402-0. [3] E. Voronovskaja, Détermination de la forme asymptotique d’approximation des fonctions par les polynômes de M. Bernstein, C.R. Acad. Sci. URSS, 4 (1932), 79–85. [4] G.G. Lorentz, Approximation of Functions, Athena Series, Holt, Rinehart and Winston, New York, (1966). [5] H. Khosravian-Arab, M. Dehghan, M.R. Eslahchi, A New Approach to Improve the Order of Approximation of the Bernstein Operators: Theory and Applications, Numer. Algorithms, 77 (2017), 111–150. https://doi.org/10. 1007/s11075-017-0307-z. [6] A.M. Acu, V. Gupta, G. Tachev, Better Numerical Approximation by Durrmeyer Type Operators, Results Math. 74 (2019), 90. https://doi.org/10.1007/s00025-019-1019-6. [7] A.R. Gairola, K.K. Singh, L.N. Mishra, Degree of Approximation by Certain Durrmeyer Type Operators, Discontin. Nonlinear. Complex. 11 (2022), 253–273. https://doi.org/10.5890/dnc.2022.06.006. [8] E.H. Kingsley, Bernstein Polynomials for Functions of Two Variables of Class C(k), Proc. Amer. Math. Soc. 2 (1951), 64–71. https://doi.org/10.2307/2032622. [9] T. Acar, A. Kajla, Degree of Approximation for Bivariate Generalized Bernstein Type Operators, Results Math. 73 (2018), 79. https://doi.org/10.1007/s00025-018-0838-1. [10] C. Ding, F. Cao, K-Functionals and Multivariate Bernstein Polynomials, J. Approx. Theory. 155 (2008), 125–135. https://doi.org/10.1016/j.jat.2008.03.011. [11] A. Fellhauer, Approximation of Smooth Functions Using Bernstein Polynomials in Multiple Variables, arXiv:1609.01940, (2016). https://doi.org/10.48550/arXiv.1609.01940. [12] C. Munoz, A. Narkawicz, Formalization of Bernstein Polynomials and Applications to Global Optimization, J. Autom. Reason. 51 (2012), 151–196. https://doi.org/10.1007/s10817-012-9256-3. [13] A. Pallini, Bernstein-Type Approximations of Smooth Functions, Statistica. 65 (2005), 169–191. https://doi. org/10.6092/ISSN.1973-2201/84. [14] T. Sauer, Multivariate Bernstein Polynomials and Convexity, Computer Aided Geom. Design. 8 (1991), 465–478. https://doi.org/10.1016/0167-8396(91)90031-6. [15] A.Yu. Veretennikov, E.V. Veretennikova, On Partial Derivatives of Multivariate Bernstein Polynomials, Sib. Adv. Math. 26 (2016), 294–305. https://doi.org/10.3103/s1055134416040039. [16] J. Wang, Z. Peng, S. Duan, J. Jing, Derivatives of Multivariate Bernstein Operators and Smoothness with Jacobi Weights, J. Appl. Math. 2012 (2012), 346132. https://doi.org/10.1155/2012/346132. [17] V.N. Mishra, K. Khatri, L.N. Mishra, Deepmala, Inverse Result in Simultaneous Approximation by Baskakov-Durrmeyer-Stancu Operators, J. Inequal. Appl. 2013 (2013), 586. https://doi.org/10.1186/ 1029-242x-2013-586. [18] R.B. Gandhi, Deepmala, V.N. Mishra, Local and Global Results for Modified Szász–Mirakjan Operators, Math. Methods Appl. Sci. 40 (2017), 2491–2504. https://doi.org/10.1002/mma.4171. [19] V.N. Mishra, K. Khatri, L.N. Mishra, Statistical Approximation by Kantorovich Type Discrete q−Beta Operators, Adv. Differ. Equ. 2013 (2013), 345. https://doi.org/10.1186/1687-1847-2013-345. [20] V.N. Mishra, K. Khatri, L.N. Mishra, On Simultaneous Approximation for Baskakov-Durrmeyer-Stancu Type Op- erators, J. Ultra Sci. Phys. Sci. 24 (2012), 567–577. [21] A.R. Gairola, Deepmala, L.N. Mishra, Rate of Approximation by Finite Iterates of q-Durrmeyer Operators, Proc. Natl. Acad. Sci., India, Sect. A Phys. Sci. 86 (2016), 229–234. https://doi.org/10.1007/s40010-016-0267-z. https://doi.org/10.1002/nme.1620111012 https://doi.org/10.1007/978-3-319-55402-0 https://doi.org/10.1007/978-3-319-55402-0 https://doi.org/10.1007/s11075-017-0307-z https://doi.org/10.1007/s11075-017-0307-z https://doi.org/10.1007/s00025-019-1019-6 https://doi.org/10.5890/dnc.2022.06.006 https://doi.org/10.2307/2032622 https://doi.org/10.1007/s00025-018-0838-1 https://doi.org/10.1016/j.jat.2008.03.011 https://doi.org/10.48550/arXiv.1609.01940 https://doi.org/10.1007/s10817-012-9256-3 https://doi.org/10.6092/ISSN.1973-2201/84 https://doi.org/10.6092/ISSN.1973-2201/84 https://doi.org/10.1016/0167-8396(91)90031-6 https://doi.org/10.3103/s1055134416040039 https://doi.org/10.1155/2012/346132 https://doi.org/10.1186/1029-242x-2013-586 https://doi.org/10.1186/1029-242x-2013-586 https://doi.org/10.1002/mma.4171 https://doi.org/10.1186/1687-1847-2013-345 https://doi.org/10.1007/s40010-016-0267-z Int. J. Anal. Appl. (2022), 20:60 19 [22] A.R. Gairola, Deepmala, L.N. Mishra, On the q-Derivatives of a Certain Linear Positive Operators, Iran. J. Sci. Technol. Trans. Sci. 42 (2017), 1409–1417. https://doi.org/10.1007/s40995-017-0227-8. [23] P.L. Butzer, On Two-Dimensional Bernstein Polynomials, Can. J. Math. 5 (1953), 107–113. https://doi.org/ 10.4153/cjm-1953-014-2. [24] O. T. Pop, About the Generalization of Voronovskaja’s Theorem for Bernstein Polynomials of Two Variables, Int. J. Pure Appl. Math. 38 (2007), 297-308. [25] D.D. Stancu, on Some Polynomials of Bernstein Type, Stud. Cerc. St. Mat. Iaşi, 11 (1960), 221-233. [26] X. Zhou, Approximation by Multivariate Bernstein Operators, Results. Math. 25 (1994), 166–191. https://doi. org/10.1007/bf03323150. [27] A. Bayad, T. Kim, S.H. Rim, Bernstein polynomials on Simplex, arXiv:1106.2482, (2011). https://doi.org/10. 48550/ARXIV.1106.2482. [28] N. Deo, N. Bhardwaj, Some Approximation Theorems for Multivariate Bernstein Operators, Southeast Asian Bull. Math. 34 (2010), 1023–1034. [29] T.D. Do, S. Waldron, Multivariate Bernstein Operators and Redundant Systems, J. Approx. Theory. 192 (2015), 215–233. https://doi.org/10.1016/j.jat.2014.12.001. [30] R.T. Farouki, T.N.T. Goodman, T. Sauer, Construction of Orthogonal Bases for Polynomials in Bernstein Form on Triangular and Simplex Domains, Computer Aided Geom. Design. 20 (2003), 209–230. https://doi.org/10. 1016/s0167-8396(03)00025-6. [31] V. I. Volkov, On the Convergence of Sequences of Linear Positive Operators in the Space of Continuous Functions of Two Variables, Dokl. Akad. Nauk SSSR, 115 (1957), 17–19. https://doi.org/10.1007/s40995-017-0227-8 https://doi.org/10.4153/cjm-1953-014-2 https://doi.org/10.4153/cjm-1953-014-2 https://doi.org/10.1007/bf03323150 https://doi.org/10.1007/bf03323150 https://doi.org/10.48550/ARXIV.1106.2482 https://doi.org/10.48550/ARXIV.1106.2482 https://doi.org/10.1016/j.jat.2014.12.001 https://doi.org/10.1016/s0167-8396(03)00025-6 https://doi.org/10.1016/s0167-8396(03)00025-6 1. Introduction 2. The first order operator BM,1n,m(f;x,y) 3. Preliminaries 4. Rate of Approximation 5. The Second order operator Bn,mM,2(f;x,y) 6. The Third Order operator Bn,mM,3(f;x,y) 7. Numerical Verification 8. Conclusion and Future Scope References